Annotation of OpenXM_contrib/gmp/demos/pexpr.c, Revision 1.1
1.1 ! maekawa 1: /* Program for computing integer expressions using the GNU Multiple Precision
! 2: Arithmetic Library.
! 3:
! 4: Copyright (C) 1997, 1999, 2000 Free Software Foundation, Inc.
! 5:
! 6: This program is free software; you can redistribute it and/or modify it under
! 7: the terms of the GNU General Public License as published by the Free Software
! 8: Foundation; either version 2 of the License, or (at your option) any later
! 9: version.
! 10:
! 11: This program is distributed in the hope that it will be useful, but WITHOUT ANY
! 12: WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
! 13: PARTICULAR PURPOSE. See the GNU General Public License for more details.
! 14:
! 15: You should have received a copy of the GNU General Public License along with
! 16: this program; if not, write to the Free Software Foundation, Inc., 59 Temple
! 17: Place - Suite 330, Boston, MA 02111-1307, USA. */
! 18:
! 19:
! 20: /* This expressions evaluator works by building an expression tree (using a
! 21: recursive descent parser) which is then evaluated. The expression tree is
! 22: useful since we want to optimize certain expressions (like a^b % c).
! 23:
! 24: Usage: pexpr [options] expr ...
! 25: (Assuming you called the executable `pexpr' of course.)
! 26:
! 27: Command line options:
! 28:
! 29: -b print output in binary
! 30: -o print output in octal
! 31: -d print output in decimal (the default)
! 32: -x print output in hexadecimal
! 33: -<NUM> print output in base NUM
! 34: -t print timing information
! 35: -html output html
! 36: -nosplit do not split long lines each 60th digit
! 37: */
! 38:
! 39: /* Define LIMIT_RESOURCE_USAGE if you want to make sure the program doesn't
! 40: use up extensive resources (cpu, memory). Useful for the GMP demo on the
! 41: GMP web site, since we cannot load the server too much. */
! 42:
! 43: #ifdef LIMIT_RESOURCE_USAGE
! 44: #include <sys/types.h>
! 45: #include <sys/time.h>
! 46: #include <sys/resource.h>
! 47: #endif
! 48:
! 49: #include <string.h>
! 50: #include <stdio.h>
! 51: #include <stdlib.h>
! 52: #include <setjmp.h>
! 53: #include <signal.h>
! 54: #include <ctype.h>
! 55:
! 56: #include "gmp.h"
! 57:
! 58: /* GMP version 1.x compatibility. */
! 59: #if ! (__GNU_MP_VERSION >= 2)
! 60: typedef MP_INT __mpz_struct;
! 61: typedef __mpz_struct mpz_t[1];
! 62: typedef __mpz_struct *mpz_ptr;
! 63: #define mpz_fdiv_q mpz_div
! 64: #define mpz_fdiv_r mpz_mod
! 65: #define mpz_tdiv_q_2exp mpz_div_2exp
! 66: #define mpz_sgn(Z) ((Z)->size < 0 ? -1 : (Z)->size > 0)
! 67: #endif
! 68:
! 69: /* GMP version 2.0 compatibility. */
! 70: #if ! (__GNU_MP_VERSION > 2 || __GNU_MP_VERSION_MINOR >= 1)
! 71: #define mpz_swap(a,b) \
! 72: do { __mpz_struct __t; __t = *a; *a = *b; *b = __t;} while (0)
! 73: #endif
! 74:
! 75: jmp_buf errjmpbuf;
! 76:
! 77: enum op_t {NOP, LIT, NEG, NOT, PLUS, MINUS, MULT, DIV, MOD, REM, INVMOD, POW,
! 78: AND, IOR, XOR, SLL, SRA, POPCNT, HAMDIST, GCD, LCM, SQRT, ROOT, FAC,
! 79: LOG, LOG2, FERMAT, MERSENNE, FIBONACCI};
! 80:
! 81: /* Type for the expression tree. */
! 82: struct expr
! 83: {
! 84: enum op_t op;
! 85: union
! 86: {
! 87: struct {struct expr *lhs, *rhs;} ops;
! 88: mpz_t val;
! 89: } operands;
! 90: };
! 91:
! 92: typedef struct expr *expr_t;
! 93:
! 94: void cleanup_and_exit (int);
! 95:
! 96: char *skipspace (char *);
! 97: void makeexp (expr_t *, enum op_t, expr_t, expr_t);
! 98: void free_expr (expr_t);
! 99: char *expr (char *, expr_t *);
! 100: char *term (char *, expr_t *);
! 101: char *power (char *, expr_t *);
! 102: char *factor (char *, expr_t *);
! 103: int match (char *, char *);
! 104: int matchp (char *, char *);
! 105: int cputime (void);
! 106:
! 107: void mpz_eval_expr (mpz_ptr, expr_t);
! 108: void mpz_eval_mod_expr (mpz_ptr, expr_t, mpz_ptr);
! 109:
! 110: char *error;
! 111: int flag_print = 1;
! 112: int print_timing = 0;
! 113: int flag_html = 0;
! 114: int flag_splitup_output = 0;
! 115: char *newline = "";
! 116:
! 117: #ifdef _AIX
! 118: #define sigaltstack sigstack
! 119: #endif
! 120:
! 121: #if !defined(_WIN32) && !defined(__DJGPP__)
! 122: void
! 123: setup_error_handler ()
! 124: {
! 125: struct sigaction act;
! 126: struct sigaltstack sigstk;
! 127:
! 128: /* Set up a stack for signal handling. A typical cause of error is stack
! 129: overflow, and in such situation a signal can not be delivered on the
! 130: overflown stack. */
! 131: sigstk.ss_sp = malloc (SIGSTKSZ);
! 132: #ifndef _AIX
! 133: sigstk.ss_size = SIGSTKSZ;
! 134: sigstk.ss_flags = 0;
! 135: #endif /* ! _AIX */
! 136:
! 137: #ifndef _UNICOS
! 138: if (sigaltstack (&sigstk, 0) < 0)
! 139: perror("sigaltstack");
! 140: #endif
! 141:
! 142: /* Initialize structure for sigaction (called below). */
! 143: act.sa_handler = cleanup_and_exit;
! 144: sigemptyset (&(act.sa_mask));
! 145: act.sa_flags = SA_ONSTACK;
! 146:
! 147: #ifdef LIMIT_RESOURCE_USAGE
! 148: {
! 149: struct rlimit limit;
! 150:
! 151: limit.rlim_cur = limit.rlim_max = 0;
! 152: setrlimit (RLIMIT_CORE, &limit);
! 153:
! 154: limit.rlim_cur = 3;
! 155: limit.rlim_max = 4;
! 156: setrlimit (RLIMIT_CPU, &limit);
! 157:
! 158: limit.rlim_cur = limit.rlim_max = 4 * 1024 * 1024;
! 159: setrlimit (RLIMIT_DATA, &limit);
! 160:
! 161: getrlimit (RLIMIT_STACK, &limit);
! 162: limit.rlim_cur = 1 * 1024 * 1024;
! 163: setrlimit (RLIMIT_STACK, &limit);
! 164:
! 165: sigaction (SIGXCPU, &act, 0);
! 166: }
! 167: #endif /* LIMIT_RESOURCE_USAGE */
! 168:
! 169: sigaction (SIGILL, &act, 0);
! 170: sigaction (SIGSEGV, &act, 0);
! 171: sigaction (SIGBUS, &act, 0);
! 172: sigaction (SIGFPE, &act, 0);
! 173: }
! 174: #endif /* ! _WIN32 && ! __DJGPP__ */
! 175:
! 176: main (int argc, char **argv)
! 177: {
! 178: struct expr *e;
! 179: int i;
! 180: mpz_t r;
! 181: int errcode = 0;
! 182: char *str;
! 183: int base = 10;
! 184:
! 185: #if !defined(_WIN32) && !defined(__DJGPP__)
! 186: setup_error_handler ();
! 187: #endif
! 188:
! 189: mpz_init (r);
! 190:
! 191: while (argc > 1 && argv[1][0] == '-')
! 192: {
! 193: char *arg = argv[1];
! 194:
! 195: if (arg[1] >= '0' && arg[1] <= '9')
! 196: break;
! 197:
! 198: if (arg[1] == 't')
! 199: print_timing = 1;
! 200: else if (arg[1] == 'b' && arg[2] >= '0' && arg[2] <= '9')
! 201: {
! 202: base = atoi (arg + 2);
! 203: if (base < 2 || base > 36)
! 204: {
! 205: fprintf (stderr, "error: invalid output base\n");
! 206: exit (-1);
! 207: }
! 208: }
! 209: else if (arg[1] == 'b' && arg[2] == 0)
! 210: base = 2;
! 211: else if (arg[1] == 'x' && arg[2] == 0)
! 212: base = 16;
! 213: else if (arg[1] == 'o' && arg[2] == 0)
! 214: base = 8;
! 215: else if (arg[1] == 'd' && arg[2] == 0)
! 216: base = 10;
! 217: else if (strcmp (arg, "-html") == 0)
! 218: {
! 219: flag_html = 1;
! 220: newline = "<BR>";
! 221: }
! 222: else if (strcmp (arg, "-split") == 0)
! 223: {
! 224: flag_splitup_output = 1;
! 225: }
! 226: else if (strcmp (arg, "-noprint") == 0)
! 227: {
! 228: flag_print = 0;
! 229: }
! 230: else
! 231: {
! 232: fprintf (stderr, "error: unknown option `%s'\n", arg);
! 233: exit (-1);
! 234: }
! 235: argv++;
! 236: argc--;
! 237: }
! 238:
! 239: for (i = 1; i < argc; i++)
! 240: {
! 241: int s;
! 242: int jmpval;
! 243:
! 244: /* Set up error handler for parsing expression. */
! 245: jmpval = setjmp (errjmpbuf);
! 246: if (jmpval != 0)
! 247: {
! 248: fprintf (stderr, "error: %s%s\n", error, newline);
! 249: fprintf (stderr, " %s%s\n", argv[i], newline);
! 250: if (! flag_html)
! 251: {
! 252: /* ??? Dunno how to align expression position with arrow in
! 253: HTML ??? */
! 254: fprintf (stderr, " ");
! 255: for (s = jmpval - (long) argv[i]; --s >= 0; )
! 256: putc (' ', stderr);
! 257: fprintf (stderr, "^\n");
! 258: }
! 259:
! 260: errcode |= 1;
! 261: continue;
! 262: }
! 263:
! 264: str = expr (argv[i], &e);
! 265:
! 266: if (str[0] != 0)
! 267: {
! 268: fprintf (stderr,
! 269: "error: garbage where end of expression expected%s\n",
! 270: newline);
! 271: fprintf (stderr, " %s%s\n", argv[i], newline);
! 272: if (! flag_html)
! 273: {
! 274: /* ??? Dunno how to align expression position with arrow in
! 275: HTML ??? */
! 276: fprintf (stderr, " ");
! 277: for (s = str - argv[i]; --s; )
! 278: putc (' ', stderr);
! 279: fprintf (stderr, "^\n");
! 280: }
! 281:
! 282: errcode |= 1;
! 283: free_expr (e);
! 284: continue;
! 285: }
! 286:
! 287: /* Set up error handler for evaluating expression. */
! 288: if (setjmp (errjmpbuf))
! 289: {
! 290: fprintf (stderr, "error: %s%s\n", error, newline);
! 291: fprintf (stderr, " %s%s\n", argv[i], newline);
! 292: if (! flag_html)
! 293: {
! 294: /* ??? Dunno how to align expression position with arrow in
! 295: HTML ??? */
! 296: fprintf (stderr, " ");
! 297: for (s = str - argv[i]; --s >= 0; )
! 298: putc (' ', stderr);
! 299: fprintf (stderr, "^\n");
! 300: }
! 301:
! 302: errcode |= 2;
! 303: continue;
! 304: }
! 305:
! 306: {
! 307: int t0;
! 308:
! 309: if (print_timing)
! 310: t0 = cputime ();
! 311:
! 312: mpz_eval_expr (r, e);
! 313:
! 314: if (print_timing)
! 315: printf ("computation took %d ms%s\n", cputime () - t0, newline);
! 316: }
! 317:
! 318: if (flag_print)
! 319: {
! 320: size_t out_len;
! 321: char *tmp, *s;
! 322: int t0;
! 323:
! 324: out_len = mpz_sizeinbase (r, base) + 1;
! 325: tmp = malloc (out_len);
! 326:
! 327: if (print_timing)
! 328: t0 = cputime ();
! 329:
! 330: if (print_timing)
! 331: /* Print first half of message... */
! 332: printf ("output conversion ");
! 333:
! 334: mpz_get_str (tmp, -base, r);
! 335:
! 336: if (print_timing)
! 337: /* ...print 2nd half of message unless we caught a time limit
! 338: and therefore longjmp'ed */
! 339: printf ("took %d ms%s\n", cputime () - t0, newline);
! 340:
! 341: out_len = strlen (tmp);
! 342: if (flag_splitup_output)
! 343: {
! 344: for (s = tmp; out_len > 60; s += 60)
! 345: {
! 346: fwrite (s, 1, 60, stdout);
! 347: printf ("%s\n", newline);
! 348: out_len -= 60;
! 349: }
! 350:
! 351: fwrite (s, 1, out_len, stdout);
! 352: }
! 353: else
! 354: {
! 355: fwrite (tmp, 1, out_len, stdout);
! 356: }
! 357:
! 358: free (tmp);
! 359: printf ("%s\n", newline);
! 360: }
! 361: else
! 362: {
! 363: printf ("result is approximately %ld digits%s\n",
! 364: (long) mpz_sizeinbase (r, 10), newline);
! 365: }
! 366:
! 367: free_expr (e);
! 368: }
! 369:
! 370: exit (errcode);
! 371: }
! 372:
! 373: char *
! 374: expr (char *str, expr_t *e)
! 375: {
! 376: expr_t e2;
! 377:
! 378: str = skipspace (str);
! 379: if (str[0] == '+')
! 380: {
! 381: str = term (str + 1, e);
! 382: }
! 383: else if (str[0] == '-')
! 384: {
! 385: str = term (str + 1, e);
! 386: makeexp (e, NEG, *e, NULL);
! 387: }
! 388: else if (str[0] == '~')
! 389: {
! 390: str = term (str + 1, e);
! 391: makeexp (e, NOT, *e, NULL);
! 392: }
! 393: else
! 394: {
! 395: str = term (str, e);
! 396: }
! 397:
! 398: for (;;)
! 399: {
! 400: str = skipspace (str);
! 401: switch (str[0])
! 402: {
! 403: case 'p':
! 404: if (match ("plus", str))
! 405: {
! 406: str = term (str + 4, &e2);
! 407: makeexp (e, PLUS, *e, e2);
! 408: }
! 409: else
! 410: return str;
! 411: break;
! 412: case 'm':
! 413: if (match ("minus", str))
! 414: {
! 415: str = term (str + 5, &e2);
! 416: makeexp (e, MINUS, *e, e2);
! 417: }
! 418: else
! 419: return str;
! 420: break;
! 421: case '+':
! 422: str = term (str + 1, &e2);
! 423: makeexp (e, PLUS, *e, e2);
! 424: break;
! 425: case '-':
! 426: str = term (str + 1, &e2);
! 427: makeexp (e, MINUS, *e, e2);
! 428: break;
! 429: default:
! 430: return str;
! 431: }
! 432: }
! 433: }
! 434:
! 435: char *
! 436: term (char *str, expr_t *e)
! 437: {
! 438: expr_t e2;
! 439:
! 440: str = power (str, e);
! 441: for (;;)
! 442: {
! 443: str = skipspace (str);
! 444: switch (str[0])
! 445: {
! 446: case 'm':
! 447: if (match ("mul", str))
! 448: {
! 449: str = power (str + 3, &e2);
! 450: makeexp (e, MULT, *e, e2);
! 451: break;
! 452: }
! 453: if (match ("mod", str))
! 454: {
! 455: str = power (str + 3, &e2);
! 456: makeexp (e, MOD, *e, e2);
! 457: break;
! 458: }
! 459: return str;
! 460: case 'd':
! 461: if (match ("div", str))
! 462: {
! 463: str = power (str + 3, &e2);
! 464: makeexp (e, DIV, *e, e2);
! 465: break;
! 466: }
! 467: return str;
! 468: case 'r':
! 469: if (match ("rem", str))
! 470: {
! 471: str = power (str + 3, &e2);
! 472: makeexp (e, REM, *e, e2);
! 473: break;
! 474: }
! 475: return str;
! 476: case 'i':
! 477: if (match ("invmod", str))
! 478: {
! 479: str = power (str + 6, &e2);
! 480: makeexp (e, REM, *e, e2);
! 481: break;
! 482: }
! 483: return str;
! 484: case 't':
! 485: if (match ("times", str))
! 486: {
! 487: str = power (str + 5, &e2);
! 488: makeexp (e, MULT, *e, e2);
! 489: break;
! 490: }
! 491: if (match ("thru", str))
! 492: {
! 493: str = power (str + 4, &e2);
! 494: makeexp (e, DIV, *e, e2);
! 495: break;
! 496: }
! 497: if (match ("through", str))
! 498: {
! 499: str = power (str + 7, &e2);
! 500: makeexp (e, DIV, *e, e2);
! 501: break;
! 502: }
! 503: return str;
! 504: case '*':
! 505: str = power (str + 1, &e2);
! 506: makeexp (e, MULT, *e, e2);
! 507: break;
! 508: case '/':
! 509: str = power (str + 1, &e2);
! 510: makeexp (e, DIV, *e, e2);
! 511: break;
! 512: case '%':
! 513: str = power (str + 1, &e2);
! 514: makeexp (e, MOD, *e, e2);
! 515: break;
! 516: default:
! 517: return str;
! 518: }
! 519: }
! 520: }
! 521:
! 522: char *
! 523: power (char *str, expr_t *e)
! 524: {
! 525: expr_t e2;
! 526:
! 527: str = factor (str, e);
! 528: while (str[0] == '!')
! 529: {
! 530: str++;
! 531: makeexp (e, FAC, *e, NULL);
! 532: }
! 533: str = skipspace (str);
! 534: if (str[0] == '^')
! 535: {
! 536: str = power (str + 1, &e2);
! 537: makeexp (e, POW, *e, e2);
! 538: }
! 539: return str;
! 540: }
! 541:
! 542: int
! 543: match (char *s, char *str)
! 544: {
! 545: char *ostr = str;
! 546: int i;
! 547:
! 548: for (i = 0; s[i] != 0; i++)
! 549: {
! 550: if (str[i] != s[i])
! 551: return 0;
! 552: }
! 553: str = skipspace (str + i);
! 554: return str - ostr;
! 555: }
! 556:
! 557: int
! 558: matchp (char *s, char *str)
! 559: {
! 560: char *ostr = str;
! 561: int i;
! 562:
! 563: for (i = 0; s[i] != 0; i++)
! 564: {
! 565: if (str[i] != s[i])
! 566: return 0;
! 567: }
! 568: str = skipspace (str + i);
! 569: if (str[0] == '(')
! 570: return str - ostr + 1;
! 571: return 0;
! 572: }
! 573:
! 574: struct functions
! 575: {
! 576: char *spelling;
! 577: enum op_t op;
! 578: int arity; /* 1 or 2 means real arity; 0 means arbitrary. */
! 579: };
! 580:
! 581: struct functions fns[] =
! 582: {
! 583: {"sqrt", SQRT, 1},
! 584: #if __GNU_MP_VERSION >= 2
! 585: {"root", ROOT, 2},
! 586: {"popc", POPCNT, 1},
! 587: #endif
! 588: {"gcd", GCD, 0},
! 589: #if __GNU_MP_VERSION > 2 || __GNU_MP_VERSION_MINOR >= 1
! 590: {"lcm", LCM, 0},
! 591: #endif
! 592: {"and", AND, 0},
! 593: {"ior", IOR, 0},
! 594: #if __GNU_MP_VERSION > 2 || __GNU_MP_VERSION_MINOR >= 1
! 595: {"xor", XOR, 0},
! 596: #endif
! 597: {"plus", PLUS, 0},
! 598: {"minus", MINUS, 2},
! 599: {"mul", MULT, 0},
! 600: {"div", DIV, 2},
! 601: {"mod", MOD, 2},
! 602: {"rem", REM, 2},
! 603: #if __GNU_MP_VERSION >= 2
! 604: {"invmod", INVMOD, 2},
! 605: #endif
! 606: {"log", LOG, 2},
! 607: {"log2", LOG2, 1},
! 608: {"F", FERMAT, 1},
! 609: {"M", MERSENNE, 1},
! 610: {"fib", FIBONACCI, 1},
! 611: {"Fib", FIBONACCI, 1},
! 612: {"", NOP, 0}
! 613: };
! 614:
! 615: char *
! 616: factor (char *str, expr_t *e)
! 617: {
! 618: expr_t e1, e2;
! 619:
! 620: str = skipspace (str);
! 621:
! 622: if (isalpha (str[0]))
! 623: {
! 624: int i;
! 625: int cnt;
! 626:
! 627: for (i = 0; fns[i].op != NOP; i++)
! 628: {
! 629: if (fns[i].arity == 1)
! 630: {
! 631: cnt = matchp (fns[i].spelling, str);
! 632: if (cnt != 0)
! 633: {
! 634: str = expr (str + cnt, &e1);
! 635: str = skipspace (str);
! 636: if (str[0] != ')')
! 637: {
! 638: error = "expected `)'";
! 639: longjmp (errjmpbuf, (int) (long) str);
! 640: }
! 641: makeexp (e, fns[i].op, e1, NULL);
! 642: return str + 1;
! 643: }
! 644: }
! 645: }
! 646:
! 647: for (i = 0; fns[i].op != NOP; i++)
! 648: {
! 649: if (fns[i].arity != 1)
! 650: {
! 651: cnt = matchp (fns[i].spelling, str);
! 652: if (cnt != 0)
! 653: {
! 654: str = expr (str + cnt, &e1);
! 655: str = skipspace (str);
! 656:
! 657: if (str[0] != ',')
! 658: {
! 659: error = "expected `,' and another operand";
! 660: longjmp (errjmpbuf, (int) (long) str);
! 661: }
! 662:
! 663: str = skipspace (str + 1);
! 664: str = expr (str, &e2);
! 665: str = skipspace (str);
! 666:
! 667: if (fns[i].arity == 0)
! 668: {
! 669: while (str[0] == ',')
! 670: {
! 671: makeexp (&e1, fns[i].op, e1, e2);
! 672: str = skipspace (str + 1);
! 673: str = expr (str, &e2);
! 674: str = skipspace (str);
! 675: }
! 676: }
! 677:
! 678: if (str[0] != ')')
! 679: {
! 680: error = "expected `)'";
! 681: longjmp (errjmpbuf, (int) (long) str);
! 682: }
! 683:
! 684: makeexp (e, fns[i].op, e1, e2);
! 685: return str + 1;
! 686: }
! 687: }
! 688: }
! 689: }
! 690:
! 691: if (str[0] == '(')
! 692: {
! 693: str = expr (str + 1, e);
! 694: str = skipspace (str);
! 695: if (str[0] != ')')
! 696: {
! 697: error = "expected `)'";
! 698: longjmp (errjmpbuf, (int) (long) str);
! 699: }
! 700: str++;
! 701: }
! 702: else if (str[0] >= '0' && str[0] <= '9')
! 703: {
! 704: expr_t res;
! 705: char *s, *sc;
! 706:
! 707: res = malloc (sizeof (struct expr));
! 708: res -> op = LIT;
! 709: mpz_init (res->operands.val);
! 710:
! 711: s = str;
! 712: while (isalnum (str[0]))
! 713: str++;
! 714: sc = malloc (str - s + 1);
! 715: memcpy (sc, s, str - s);
! 716: sc[str - s] = 0;
! 717:
! 718: mpz_set_str (res->operands.val, sc, 0);
! 719: *e = res;
! 720: free (sc);
! 721: }
! 722: else
! 723: {
! 724: error = "operand expected";
! 725: longjmp (errjmpbuf, (int) (long) str);
! 726: }
! 727: return str;
! 728: }
! 729:
! 730: char *
! 731: skipspace (char *str)
! 732: {
! 733: while (str[0] == ' ')
! 734: str++;
! 735: return str;
! 736: }
! 737:
! 738: /* Make a new expression with operation OP and right hand side
! 739: RHS and left hand side lhs. Put the result in R. */
! 740: void
! 741: makeexp (expr_t *r, enum op_t op, expr_t lhs, expr_t rhs)
! 742: {
! 743: expr_t res;
! 744: res = malloc (sizeof (struct expr));
! 745: res -> op = op;
! 746: res -> operands.ops.lhs = lhs;
! 747: res -> operands.ops.rhs = rhs;
! 748: *r = res;
! 749: return;
! 750: }
! 751:
! 752: /* Free the memory used by expression E. */
! 753: void
! 754: free_expr (expr_t e)
! 755: {
! 756: if (e->op != LIT)
! 757: {
! 758: free_expr (e->operands.ops.lhs);
! 759: if (e->operands.ops.rhs != NULL)
! 760: free_expr (e->operands.ops.rhs);
! 761: }
! 762: else
! 763: {
! 764: mpz_clear (e->operands.val);
! 765: }
! 766: }
! 767:
! 768: /* Evaluate the expression E and put the result in R. */
! 769: void
! 770: mpz_eval_expr (mpz_ptr r, expr_t e)
! 771: {
! 772: mpz_t lhs, rhs;
! 773:
! 774: switch (e->op)
! 775: {
! 776: case LIT:
! 777: mpz_set (r, e->operands.val);
! 778: return;
! 779: case PLUS:
! 780: mpz_init (lhs); mpz_init (rhs);
! 781: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 782: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 783: mpz_add (r, lhs, rhs);
! 784: mpz_clear (lhs); mpz_clear (rhs);
! 785: return;
! 786: case MINUS:
! 787: mpz_init (lhs); mpz_init (rhs);
! 788: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 789: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 790: mpz_sub (r, lhs, rhs);
! 791: mpz_clear (lhs); mpz_clear (rhs);
! 792: return;
! 793: case MULT:
! 794: mpz_init (lhs); mpz_init (rhs);
! 795: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 796: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 797: mpz_mul (r, lhs, rhs);
! 798: mpz_clear (lhs); mpz_clear (rhs);
! 799: return;
! 800: case DIV:
! 801: mpz_init (lhs); mpz_init (rhs);
! 802: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 803: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 804: mpz_fdiv_q (r, lhs, rhs);
! 805: mpz_clear (lhs); mpz_clear (rhs);
! 806: return;
! 807: case MOD:
! 808: mpz_init (rhs);
! 809: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 810: mpz_abs (rhs, rhs);
! 811: mpz_eval_mod_expr (r, e->operands.ops.lhs, rhs);
! 812: mpz_clear (rhs);
! 813: return;
! 814: case REM:
! 815: /* Check if lhs operand is POW expression and optimize for that case. */
! 816: if (e->operands.ops.lhs->op == POW)
! 817: {
! 818: mpz_t powlhs, powrhs;
! 819: mpz_init (powlhs);
! 820: mpz_init (powrhs);
! 821: mpz_init (rhs);
! 822: mpz_eval_expr (powlhs, e->operands.ops.lhs->operands.ops.lhs);
! 823: mpz_eval_expr (powrhs, e->operands.ops.lhs->operands.ops.rhs);
! 824: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 825: mpz_powm (r, powlhs, powrhs, rhs);
! 826: if (mpz_cmp_si (rhs, 0L) < 0)
! 827: mpz_neg (r, r);
! 828: mpz_clear (powlhs);
! 829: mpz_clear (powrhs);
! 830: mpz_clear (rhs);
! 831: return;
! 832: }
! 833:
! 834: mpz_init (lhs); mpz_init (rhs);
! 835: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 836: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 837: mpz_fdiv_r (r, lhs, rhs);
! 838: mpz_clear (lhs); mpz_clear (rhs);
! 839: return;
! 840: #if __GNU_MP_VERSION >= 2
! 841: case INVMOD:
! 842: mpz_init (lhs); mpz_init (rhs);
! 843: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 844: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 845: mpz_invert (r, lhs, rhs);
! 846: mpz_clear (lhs); mpz_clear (rhs);
! 847: return;
! 848: #endif
! 849: case POW:
! 850: mpz_init (lhs); mpz_init (rhs);
! 851: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 852: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 853: if (mpz_cmp_si (rhs, 0L) == 0)
! 854: /* x^0 is 1 */
! 855: mpz_set_ui (r, 1L);
! 856: else if (mpz_cmp_si (lhs, 0L) == 0)
! 857: /* 0^y (where y != 0) is 0 */
! 858: mpz_set_ui (r, 0L);
! 859: else if (mpz_cmp_ui (lhs, 1L) == 0)
! 860: /* 1^y is 1 */
! 861: mpz_set_ui (r, 1L);
! 862: else if (mpz_cmp_si (lhs, -1L) == 0)
! 863: /* (-1)^y just depends on whether y is even or odd */
! 864: mpz_set_si (r, (mpz_get_ui (rhs) & 1) ? -1L : 1L);
! 865: else if (mpz_cmp_si (rhs, 0L) < 0)
! 866: /* x^(-n) is 0 */
! 867: mpz_set_ui (r, 0L);
! 868: else
! 869: {
! 870: unsigned long int cnt;
! 871: unsigned long int y;
! 872: /* error if exponent does not fit into an unsigned long int. */
! 873: if (mpz_cmp_ui (rhs, ~(unsigned long int) 0) > 0)
! 874: goto pow_err;
! 875:
! 876: y = mpz_get_ui (rhs);
! 877: /* x^y == (x/(2^c))^y * 2^(c*y) */
! 878: #if __GNU_MP_VERSION >= 2
! 879: cnt = mpz_scan1 (lhs, 0);
! 880: #else
! 881: cnt = 0;
! 882: #endif
! 883: if (cnt != 0)
! 884: {
! 885: if (y * cnt / cnt != y)
! 886: goto pow_err;
! 887: mpz_tdiv_q_2exp (lhs, lhs, cnt);
! 888: mpz_pow_ui (r, lhs, y);
! 889: mpz_mul_2exp (r, r, y * cnt);
! 890: }
! 891: else
! 892: mpz_pow_ui (r, lhs, y);
! 893: }
! 894: mpz_clear (lhs); mpz_clear (rhs);
! 895: return;
! 896: pow_err:
! 897: error = "result of `pow' operator too large";
! 898: mpz_clear (lhs); mpz_clear (rhs);
! 899: longjmp (errjmpbuf, 1);
! 900: case GCD:
! 901: mpz_init (lhs); mpz_init (rhs);
! 902: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 903: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 904: mpz_gcd (r, lhs, rhs);
! 905: mpz_clear (lhs); mpz_clear (rhs);
! 906: return;
! 907: #if __GNU_MP_VERSION > 2 || __GNU_MP_VERSION_MINOR >= 1
! 908: case LCM:
! 909: mpz_init (lhs); mpz_init (rhs);
! 910: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 911: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 912: mpz_lcm (r, lhs, rhs);
! 913: mpz_clear (lhs); mpz_clear (rhs);
! 914: return;
! 915: #endif
! 916: case AND:
! 917: mpz_init (lhs); mpz_init (rhs);
! 918: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 919: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 920: mpz_and (r, lhs, rhs);
! 921: mpz_clear (lhs); mpz_clear (rhs);
! 922: return;
! 923: case IOR:
! 924: mpz_init (lhs); mpz_init (rhs);
! 925: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 926: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 927: mpz_ior (r, lhs, rhs);
! 928: mpz_clear (lhs); mpz_clear (rhs);
! 929: return;
! 930: #if __GNU_MP_VERSION > 2 || __GNU_MP_VERSION_MINOR >= 1
! 931: case XOR:
! 932: mpz_init (lhs); mpz_init (rhs);
! 933: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 934: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 935: mpz_xor (r, lhs, rhs);
! 936: mpz_clear (lhs); mpz_clear (rhs);
! 937: return;
! 938: #endif
! 939: case NEG:
! 940: mpz_eval_expr (r, e->operands.ops.lhs);
! 941: mpz_neg (r, r);
! 942: return;
! 943: case NOT:
! 944: mpz_eval_expr (r, e->operands.ops.lhs);
! 945: mpz_com (r, r);
! 946: return;
! 947: case SQRT:
! 948: mpz_init (lhs);
! 949: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 950: if (mpz_sgn (lhs) < 0)
! 951: {
! 952: error = "cannot take square root of negative numbers";
! 953: mpz_clear (lhs);
! 954: longjmp (errjmpbuf, 1);
! 955: }
! 956: mpz_sqrt (r, lhs);
! 957: return;
! 958: #if __GNU_MP_VERSION > 2 || __GNU_MP_VERSION_MINOR >= 1
! 959: case ROOT:
! 960: mpz_init (lhs); mpz_init (rhs);
! 961: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 962: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 963: if (mpz_sgn (rhs) <= 0)
! 964: {
! 965: error = "cannot take non-positive root orders";
! 966: mpz_clear (lhs); mpz_clear (rhs);
! 967: longjmp (errjmpbuf, 1);
! 968: }
! 969: if (mpz_sgn (lhs) < 0 && (mpz_get_ui (rhs) & 1) == 0)
! 970: {
! 971: error = "cannot take even root orders of negative numbers";
! 972: mpz_clear (lhs); mpz_clear (rhs);
! 973: longjmp (errjmpbuf, 1);
! 974: }
! 975:
! 976: {
! 977: unsigned long int nth = mpz_get_ui (rhs);
! 978: if (mpz_cmp_ui (rhs, ~(unsigned long int) 0) > 0)
! 979: {
! 980: /* If we are asked to take an awfully large root order, cheat and
! 981: ask for the largest order we can pass to mpz_root. This saves
! 982: some error prone special cases. */
! 983: nth = ~(unsigned long int) 0;
! 984: }
! 985: mpz_root (r, lhs, nth);
! 986: }
! 987: mpz_clear (lhs); mpz_clear (rhs);
! 988: return;
! 989: #endif
! 990: case FAC:
! 991: mpz_eval_expr (r, e->operands.ops.lhs);
! 992: if (mpz_size (r) > 1)
! 993: {
! 994: error = "result of `!' operator too large";
! 995: longjmp (errjmpbuf, 1);
! 996: }
! 997: mpz_fac_ui (r, mpz_get_ui (r));
! 998: return;
! 999: #if __GNU_MP_VERSION >= 2
! 1000: case POPCNT:
! 1001: mpz_eval_expr (r, e->operands.ops.lhs);
! 1002: { unsigned long int cnt;
! 1003: cnt = mpz_popcount (r);
! 1004: mpz_set_ui (r, cnt);
! 1005: }
! 1006: return;
! 1007: #endif
! 1008: case LOG2:
! 1009: mpz_eval_expr (r, e->operands.ops.lhs);
! 1010: { unsigned long int cnt;
! 1011: if (mpz_sgn (r) <= 0)
! 1012: {
! 1013: error = "logarithm of non-positive number";
! 1014: longjmp (errjmpbuf, 1);
! 1015: }
! 1016: cnt = mpz_sizeinbase (r, 2);
! 1017: mpz_set_ui (r, cnt - 1);
! 1018: }
! 1019: return;
! 1020: case LOG:
! 1021: { unsigned long int cnt;
! 1022: mpz_init (lhs); mpz_init (rhs);
! 1023: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 1024: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 1025: if (mpz_sgn (lhs) <= 0)
! 1026: {
! 1027: error = "logarithm of non-positive number";
! 1028: mpz_clear (lhs); mpz_clear (rhs);
! 1029: longjmp (errjmpbuf, 1);
! 1030: }
! 1031: if (mpz_cmp_ui (rhs, 256) >= 0)
! 1032: {
! 1033: error = "logarithm base too large";
! 1034: mpz_clear (lhs); mpz_clear (rhs);
! 1035: longjmp (errjmpbuf, 1);
! 1036: }
! 1037: cnt = mpz_sizeinbase (lhs, mpz_get_ui (rhs));
! 1038: mpz_set_ui (r, cnt - 1);
! 1039: mpz_clear (lhs); mpz_clear (rhs);
! 1040: }
! 1041: return;
! 1042: case FERMAT:
! 1043: {
! 1044: unsigned long int t;
! 1045: mpz_init (lhs);
! 1046: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 1047: t = (unsigned long int) 1 << mpz_get_ui (lhs);
! 1048: if (mpz_cmp_ui (lhs, ~(unsigned long int) 0) > 0 || t == 0)
! 1049: {
! 1050: error = "too large Mersenne number index";
! 1051: mpz_clear (lhs);
! 1052: longjmp (errjmpbuf, 1);
! 1053: }
! 1054: mpz_set_ui (r, 1);
! 1055: mpz_mul_2exp (r, r, t);
! 1056: mpz_add_ui (r, r, 1);
! 1057: mpz_clear (lhs);
! 1058: }
! 1059: return;
! 1060: case MERSENNE:
! 1061: mpz_init (lhs);
! 1062: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 1063: if (mpz_cmp_ui (lhs, ~(unsigned long int) 0) > 0)
! 1064: {
! 1065: error = "too large Mersenne number index";
! 1066: mpz_clear (lhs);
! 1067: longjmp (errjmpbuf, 1);
! 1068: }
! 1069: mpz_set_ui (r, 1);
! 1070: mpz_mul_2exp (r, r, mpz_get_ui (lhs));
! 1071: mpz_sub_ui (r, r, 1);
! 1072: mpz_clear (lhs);
! 1073: return;
! 1074: case FIBONACCI:
! 1075: { mpz_t t;
! 1076: unsigned long int n, i;
! 1077: mpz_init (lhs);
! 1078: mpz_eval_expr (lhs, e->operands.ops.lhs);
! 1079: if (mpz_sgn (lhs) <= 0 || mpz_cmp_si (lhs, 1000000000) > 0)
! 1080: {
! 1081: error = "Fibonacci index out of range";
! 1082: mpz_clear (lhs);
! 1083: longjmp (errjmpbuf, 1);
! 1084: }
! 1085: n = mpz_get_ui (lhs);
! 1086: mpz_clear (lhs);
! 1087:
! 1088: #if __GNU_MP_VERSION > 2 || __GNU_MP_VERSION_MINOR >= 1
! 1089: mpz_fib_ui (r, n);
! 1090: #else
! 1091: mpz_init_set_ui (t, 1);
! 1092: mpz_set_ui (r, 1);
! 1093:
! 1094: if (n <= 2)
! 1095: mpz_set_ui (r, 1);
! 1096: else
! 1097: {
! 1098: for (i = 3; i <= n; i++)
! 1099: {
! 1100: mpz_add (t, t, r);
! 1101: mpz_swap (t, r);
! 1102: }
! 1103: }
! 1104: mpz_clear (t);
! 1105: #endif
! 1106: }
! 1107: return;
! 1108: default:
! 1109: abort ();
! 1110: }
! 1111: }
! 1112:
! 1113: /* Evaluate the expression E modulo MOD and put the result in R. */
! 1114: void
! 1115: mpz_eval_mod_expr (mpz_ptr r, expr_t e, mpz_ptr mod)
! 1116: {
! 1117: mpz_t lhs, rhs;
! 1118:
! 1119: switch (e->op)
! 1120: {
! 1121: case POW:
! 1122: mpz_init (lhs); mpz_init (rhs);
! 1123: mpz_eval_mod_expr (lhs, e->operands.ops.lhs, mod);
! 1124: mpz_eval_expr (rhs, e->operands.ops.rhs);
! 1125: mpz_powm (r, lhs, rhs, mod);
! 1126: mpz_clear (lhs); mpz_clear (rhs);
! 1127: return;
! 1128: case PLUS:
! 1129: mpz_init (lhs); mpz_init (rhs);
! 1130: mpz_eval_mod_expr (lhs, e->operands.ops.lhs, mod);
! 1131: mpz_eval_mod_expr (rhs, e->operands.ops.rhs, mod);
! 1132: mpz_add (r, lhs, rhs);
! 1133: if (mpz_cmp_si (r, 0L) < 0)
! 1134: mpz_add (r, r, mod);
! 1135: else if (mpz_cmp (r, mod) >= 0)
! 1136: mpz_sub (r, r, mod);
! 1137: mpz_clear (lhs); mpz_clear (rhs);
! 1138: return;
! 1139: case MINUS:
! 1140: mpz_init (lhs); mpz_init (rhs);
! 1141: mpz_eval_mod_expr (lhs, e->operands.ops.lhs, mod);
! 1142: mpz_eval_mod_expr (rhs, e->operands.ops.rhs, mod);
! 1143: mpz_sub (r, lhs, rhs);
! 1144: if (mpz_cmp_si (r, 0L) < 0)
! 1145: mpz_add (r, r, mod);
! 1146: else if (mpz_cmp (r, mod) >= 0)
! 1147: mpz_sub (r, r, mod);
! 1148: mpz_clear (lhs); mpz_clear (rhs);
! 1149: return;
! 1150: case MULT:
! 1151: mpz_init (lhs); mpz_init (rhs);
! 1152: mpz_eval_mod_expr (lhs, e->operands.ops.lhs, mod);
! 1153: mpz_eval_mod_expr (rhs, e->operands.ops.rhs, mod);
! 1154: mpz_mul (r, lhs, rhs);
! 1155: mpz_mod (r, r, mod);
! 1156: mpz_clear (lhs); mpz_clear (rhs);
! 1157: return;
! 1158: default:
! 1159: mpz_init (lhs);
! 1160: mpz_eval_expr (lhs, e);
! 1161: mpz_mod (r, lhs, mod);
! 1162: mpz_clear (lhs);
! 1163: return;
! 1164: }
! 1165: }
! 1166:
! 1167: void
! 1168: cleanup_and_exit (int sig)
! 1169: {
! 1170: #ifdef LIMIT_RESOURCE_USAGE
! 1171: if (sig == SIGXCPU)
! 1172: printf ("expression took too long time to evaluate%s\n", newline);
! 1173: else if (sig == SIGFPE)
! 1174: printf ("divide by zero%s\n", newline);
! 1175: else
! 1176: #endif
! 1177: printf ("expression required too much memory to evaluate%s\n", newline);
! 1178: exit (-2);
! 1179: }
! 1180:
! 1181: /* Return user CPU time measured in milliseconds. */
! 1182:
! 1183: #if defined (USG) || defined (__SVR4) || defined (_UNICOS) || defined (__hpux)
! 1184: #include <time.h>
! 1185:
! 1186: int
! 1187: cputime ()
! 1188: {
! 1189: if (CLOCKS_PER_SEC < 100000)
! 1190: return clock () * 1000 / CLOCKS_PER_SEC;
! 1191: return clock () / (CLOCKS_PER_SEC / 1000);
! 1192: }
! 1193: #else
! 1194: #include <sys/types.h>
! 1195: #include <sys/time.h>
! 1196: #include <sys/resource.h>
! 1197:
! 1198: int
! 1199: cputime ()
! 1200: {
! 1201: struct rusage rus;
! 1202:
! 1203: getrusage (0, &rus);
! 1204: return rus.ru_utime.tv_sec * 1000 + rus.ru_utime.tv_usec / 1000;
! 1205: }
! 1206: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to pexpr.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / demos