![[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|
Return to pexpr.c CVS log | Up to [local] / OpenXM_contrib / gmp / demos |
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