Revert "Remove gmp from AOS"
This reverts commit f37c97684f0910a3f241394549392f00145ab0f7.
We need gmp for SymEngine for symbolicmanipultion in C++
Change-Id: Ia13216d1715cf96944f7b4f422b7a799f921d4a4
Signed-off-by: Austin Schuh <austin.linux@gmail.com>
diff --git a/third_party/gmp/demos/pexpr.c b/third_party/gmp/demos/pexpr.c
new file mode 100644
index 0000000..d5009e9
--- /dev/null
+++ b/third_party/gmp/demos/pexpr.c
@@ -0,0 +1,1380 @@
+/* Program for computing integer expressions using the GNU Multiple Precision
+ Arithmetic Library.
+
+Copyright 1997, 1999-2002, 2005, 2008, 2012, 2015 Free Software Foundation, Inc.
+
+This program is free software; you can redistribute it and/or modify it under
+the terms of the GNU General Public License as published by the Free Software
+Foundation; either version 3 of the License, or (at your option) any later
+version.
+
+This program is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
+PARTICULAR PURPOSE. See the GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along with
+this program. If not, see https://www.gnu.org/licenses/. */
+
+
+/* This expressions evaluator works by building an expression tree (using a
+ recursive descent parser) which is then evaluated. The expression tree is
+ useful since we want to optimize certain expressions (like a^b % c).
+
+ Usage: pexpr [options] expr ...
+ (Assuming you called the executable `pexpr' of course.)
+
+ Command line options:
+
+ -b print output in binary
+ -o print output in octal
+ -d print output in decimal (the default)
+ -x print output in hexadecimal
+ -b<NUM> print output in base NUM
+ -t print timing information
+ -html output html
+ -wml output wml
+ -split split long lines each 80th digit
+*/
+
+/* Define LIMIT_RESOURCE_USAGE if you want to make sure the program doesn't
+ use up extensive resources (cpu, memory). Useful for the GMP demo on the
+ GMP web site, since we cannot load the server too much. */
+
+#include "pexpr-config.h"
+
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <setjmp.h>
+#include <signal.h>
+#include <ctype.h>
+
+#include <time.h>
+#include <sys/types.h>
+#include <sys/time.h>
+#if HAVE_SYS_RESOURCE_H
+#include <sys/resource.h>
+#endif
+
+#include "gmp.h"
+
+/* SunOS 4 and HPUX 9 don't define a canonical SIGSTKSZ, use a default. */
+#ifndef SIGSTKSZ
+#define SIGSTKSZ 4096
+#endif
+
+
+#define TIME(t,func) \
+ do { int __t0, __tmp; \
+ __t0 = cputime (); \
+ {func;} \
+ __tmp = cputime () - __t0; \
+ (t) = __tmp; \
+ } while (0)
+
+/* GMP version 1.x compatibility. */
+#if ! (__GNU_MP_VERSION >= 2)
+typedef MP_INT __mpz_struct;
+typedef __mpz_struct mpz_t[1];
+typedef __mpz_struct *mpz_ptr;
+#define mpz_fdiv_q mpz_div
+#define mpz_fdiv_r mpz_mod
+#define mpz_tdiv_q_2exp mpz_div_2exp
+#define mpz_sgn(Z) ((Z)->size < 0 ? -1 : (Z)->size > 0)
+#endif
+
+/* GMP version 2.0 compatibility. */
+#if ! (__GNU_MP_VERSION > 2 || __GNU_MP_VERSION_MINOR >= 1)
+#define mpz_swap(a,b) \
+ do { __mpz_struct __t; __t = *a; *a = *b; *b = __t;} while (0)
+#endif
+
+jmp_buf errjmpbuf;
+
+enum op_t {NOP, LIT, NEG, NOT, PLUS, MINUS, MULT, DIV, MOD, REM, INVMOD, POW,
+ AND, IOR, XOR, SLL, SRA, POPCNT, HAMDIST, GCD, LCM, SQRT, ROOT, FAC,
+ LOG, LOG2, FERMAT, MERSENNE, FIBONACCI, RANDOM, NEXTPRIME, BINOM,
+ TIMING};
+
+/* Type for the expression tree. */
+struct expr
+{
+ enum op_t op;
+ union
+ {
+ struct {struct expr *lhs, *rhs;} ops;
+ mpz_t val;
+ } operands;
+};
+
+typedef struct expr *expr_t;
+
+void cleanup_and_exit (int);
+
+char *skipspace (char *);
+void makeexp (expr_t *, enum op_t, expr_t, expr_t);
+void free_expr (expr_t);
+char *expr (char *, expr_t *);
+char *term (char *, expr_t *);
+char *power (char *, expr_t *);
+char *factor (char *, expr_t *);
+int match (char *, char *);
+int matchp (char *, char *);
+int cputime (void);
+
+void mpz_eval_expr (mpz_ptr, expr_t);
+void mpz_eval_mod_expr (mpz_ptr, expr_t, mpz_ptr);
+
+char *error;
+int flag_print = 1;
+int print_timing = 0;
+int flag_html = 0;
+int flag_wml = 0;
+int flag_splitup_output = 0;
+char *newline = "";
+gmp_randstate_t rstate;
+
+
+
+/* cputime() returns user CPU time measured in milliseconds. */
+#if ! HAVE_CPUTIME
+#if HAVE_GETRUSAGE
+int
+cputime (void)
+{
+ struct rusage rus;
+
+ getrusage (0, &rus);
+ return rus.ru_utime.tv_sec * 1000 + rus.ru_utime.tv_usec / 1000;
+}
+#else
+#if HAVE_CLOCK
+int
+cputime (void)
+{
+ if (CLOCKS_PER_SEC < 100000)
+ return clock () * 1000 / CLOCKS_PER_SEC;
+ return clock () / (CLOCKS_PER_SEC / 1000);
+}
+#else
+int
+cputime (void)
+{
+ return 0;
+}
+#endif
+#endif
+#endif
+
+
+int
+stack_downwards_helper (char *xp)
+{
+ char y;
+ return &y < xp;
+}
+int
+stack_downwards_p (void)
+{
+ char x;
+ return stack_downwards_helper (&x);
+}
+
+
+void
+setup_error_handler (void)
+{
+#if HAVE_SIGACTION
+ struct sigaction act;
+ act.sa_handler = cleanup_and_exit;
+ sigemptyset (&(act.sa_mask));
+#define SIGNAL(sig) sigaction (sig, &act, NULL)
+#else
+ struct { int sa_flags; } act;
+#define SIGNAL(sig) signal (sig, cleanup_and_exit)
+#endif
+ act.sa_flags = 0;
+
+ /* Set up a stack for signal handling. A typical cause of error is stack
+ overflow, and in such situation a signal can not be delivered on the
+ overflown stack. */
+#if HAVE_SIGALTSTACK
+ {
+ /* AIX uses stack_t, MacOS uses struct sigaltstack, various other
+ systems have both. */
+#if HAVE_STACK_T
+ stack_t s;
+#else
+ struct sigaltstack s;
+#endif
+ s.ss_sp = malloc (SIGSTKSZ);
+ s.ss_size = SIGSTKSZ;
+ s.ss_flags = 0;
+ if (sigaltstack (&s, NULL) != 0)
+ perror("sigaltstack");
+ act.sa_flags = SA_ONSTACK;
+ }
+#else
+#if HAVE_SIGSTACK
+ {
+ struct sigstack s;
+ s.ss_sp = malloc (SIGSTKSZ);
+ if (stack_downwards_p ())
+ s.ss_sp += SIGSTKSZ;
+ s.ss_onstack = 0;
+ if (sigstack (&s, NULL) != 0)
+ perror("sigstack");
+ act.sa_flags = SA_ONSTACK;
+ }
+#else
+#endif
+#endif
+
+#ifdef LIMIT_RESOURCE_USAGE
+ {
+ struct rlimit limit;
+
+ limit.rlim_cur = limit.rlim_max = 0;
+ setrlimit (RLIMIT_CORE, &limit);
+
+ limit.rlim_cur = 3;
+ limit.rlim_max = 4;
+ setrlimit (RLIMIT_CPU, &limit);
+
+ limit.rlim_cur = limit.rlim_max = 16 * 1024 * 1024;
+ setrlimit (RLIMIT_DATA, &limit);
+
+ getrlimit (RLIMIT_STACK, &limit);
+ limit.rlim_cur = 4 * 1024 * 1024;
+ setrlimit (RLIMIT_STACK, &limit);
+
+ SIGNAL (SIGXCPU);
+ }
+#endif /* LIMIT_RESOURCE_USAGE */
+
+ SIGNAL (SIGILL);
+ SIGNAL (SIGSEGV);
+#ifdef SIGBUS /* not in mingw */
+ SIGNAL (SIGBUS);
+#endif
+ SIGNAL (SIGFPE);
+ SIGNAL (SIGABRT);
+}
+
+int
+main (int argc, char **argv)
+{
+ struct expr *e;
+ int i;
+ mpz_t r;
+ int errcode = 0;
+ char *str;
+ int base = 10;
+
+ setup_error_handler ();
+
+ gmp_randinit (rstate, GMP_RAND_ALG_LC, 128);
+
+ {
+#if HAVE_GETTIMEOFDAY
+ struct timeval tv;
+ gettimeofday (&tv, NULL);
+ gmp_randseed_ui (rstate, tv.tv_sec + tv.tv_usec);
+#else
+ time_t t;
+ time (&t);
+ gmp_randseed_ui (rstate, t);
+#endif
+ }
+
+ mpz_init (r);
+
+ while (argc > 1 && argv[1][0] == '-')
+ {
+ char *arg = argv[1];
+
+ if (arg[1] >= '0' && arg[1] <= '9')
+ break;
+
+ if (arg[1] == 't')
+ print_timing = 1;
+ else if (arg[1] == 'b' && arg[2] >= '0' && arg[2] <= '9')
+ {
+ base = atoi (arg + 2);
+ if (base < 2 || base > 62)
+ {
+ fprintf (stderr, "error: invalid output base\n");
+ exit (-1);
+ }
+ }
+ else if (arg[1] == 'b' && arg[2] == 0)
+ base = 2;
+ else if (arg[1] == 'x' && arg[2] == 0)
+ base = 16;
+ else if (arg[1] == 'X' && arg[2] == 0)
+ base = -16;
+ else if (arg[1] == 'o' && arg[2] == 0)
+ base = 8;
+ else if (arg[1] == 'd' && arg[2] == 0)
+ base = 10;
+ else if (arg[1] == 'v' && arg[2] == 0)
+ {
+ printf ("pexpr linked to gmp %s\n", __gmp_version);
+ }
+ else if (strcmp (arg, "-html") == 0)
+ {
+ flag_html = 1;
+ newline = "<br>";
+ }
+ else if (strcmp (arg, "-wml") == 0)
+ {
+ flag_wml = 1;
+ newline = "<br/>";
+ }
+ else if (strcmp (arg, "-split") == 0)
+ {
+ flag_splitup_output = 1;
+ }
+ else if (strcmp (arg, "-noprint") == 0)
+ {
+ flag_print = 0;
+ }
+ else
+ {
+ fprintf (stderr, "error: unknown option `%s'\n", arg);
+ exit (-1);
+ }
+ argv++;
+ argc--;
+ }
+
+ for (i = 1; i < argc; i++)
+ {
+ int s;
+ int jmpval;
+
+ /* Set up error handler for parsing expression. */
+ jmpval = setjmp (errjmpbuf);
+ if (jmpval != 0)
+ {
+ fprintf (stderr, "error: %s%s\n", error, newline);
+ fprintf (stderr, " %s%s\n", argv[i], newline);
+ if (! flag_html)
+ {
+ /* ??? Dunno how to align expression position with arrow in
+ HTML ??? */
+ fprintf (stderr, " ");
+ for (s = jmpval - (long) argv[i]; --s >= 0; )
+ putc (' ', stderr);
+ fprintf (stderr, "^\n");
+ }
+
+ errcode |= 1;
+ continue;
+ }
+
+ str = expr (argv[i], &e);
+
+ if (str[0] != 0)
+ {
+ fprintf (stderr,
+ "error: garbage where end of expression expected%s\n",
+ newline);
+ fprintf (stderr, " %s%s\n", argv[i], newline);
+ if (! flag_html)
+ {
+ /* ??? Dunno how to align expression position with arrow in
+ HTML ??? */
+ fprintf (stderr, " ");
+ for (s = str - argv[i]; --s; )
+ putc (' ', stderr);
+ fprintf (stderr, "^\n");
+ }
+
+ errcode |= 1;
+ free_expr (e);
+ continue;
+ }
+
+ /* Set up error handler for evaluating expression. */
+ if (setjmp (errjmpbuf))
+ {
+ fprintf (stderr, "error: %s%s\n", error, newline);
+ fprintf (stderr, " %s%s\n", argv[i], newline);
+ if (! flag_html)
+ {
+ /* ??? Dunno how to align expression position with arrow in
+ HTML ??? */
+ fprintf (stderr, " ");
+ for (s = str - argv[i]; --s >= 0; )
+ putc (' ', stderr);
+ fprintf (stderr, "^\n");
+ }
+
+ errcode |= 2;
+ continue;
+ }
+
+ if (print_timing)
+ {
+ int t;
+ TIME (t, mpz_eval_expr (r, e));
+ printf ("computation took %d ms%s\n", t, newline);
+ }
+ else
+ mpz_eval_expr (r, e);
+
+ if (flag_print)
+ {
+ size_t out_len;
+ char *tmp, *s;
+
+ out_len = mpz_sizeinbase (r, base >= 0 ? base : -base) + 2;
+#ifdef LIMIT_RESOURCE_USAGE
+ if (out_len > 100000)
+ {
+ printf ("result is about %ld digits, not printing it%s\n",
+ (long) out_len - 3, newline);
+ exit (-2);
+ }
+#endif
+ tmp = malloc (out_len);
+
+ if (print_timing)
+ {
+ int t;
+ printf ("output conversion ");
+ TIME (t, mpz_get_str (tmp, base, r));
+ printf ("took %d ms%s\n", t, newline);
+ }
+ else
+ mpz_get_str (tmp, base, r);
+
+ out_len = strlen (tmp);
+ if (flag_splitup_output)
+ {
+ for (s = tmp; out_len > 80; s += 80)
+ {
+ fwrite (s, 1, 80, stdout);
+ printf ("%s\n", newline);
+ out_len -= 80;
+ }
+
+ fwrite (s, 1, out_len, stdout);
+ }
+ else
+ {
+ fwrite (tmp, 1, out_len, stdout);
+ }
+
+ free (tmp);
+ printf ("%s\n", newline);
+ }
+ else
+ {
+ printf ("result is approximately %ld digits%s\n",
+ (long) mpz_sizeinbase (r, base >= 0 ? base : -base),
+ newline);
+ }
+
+ free_expr (e);
+ }
+
+ mpz_clear (r);
+
+ exit (errcode);
+}
+
+char *
+expr (char *str, expr_t *e)
+{
+ expr_t e2;
+
+ str = skipspace (str);
+ if (str[0] == '+')
+ {
+ str = term (str + 1, e);
+ }
+ else if (str[0] == '-')
+ {
+ str = term (str + 1, e);
+ makeexp (e, NEG, *e, NULL);
+ }
+ else if (str[0] == '~')
+ {
+ str = term (str + 1, e);
+ makeexp (e, NOT, *e, NULL);
+ }
+ else
+ {
+ str = term (str, e);
+ }
+
+ for (;;)
+ {
+ str = skipspace (str);
+ switch (str[0])
+ {
+ case 'p':
+ if (match ("plus", str))
+ {
+ str = term (str + 4, &e2);
+ makeexp (e, PLUS, *e, e2);
+ }
+ else
+ return str;
+ break;
+ case 'm':
+ if (match ("minus", str))
+ {
+ str = term (str + 5, &e2);
+ makeexp (e, MINUS, *e, e2);
+ }
+ else
+ return str;
+ break;
+ case '+':
+ str = term (str + 1, &e2);
+ makeexp (e, PLUS, *e, e2);
+ break;
+ case '-':
+ str = term (str + 1, &e2);
+ makeexp (e, MINUS, *e, e2);
+ break;
+ default:
+ return str;
+ }
+ }
+}
+
+char *
+term (char *str, expr_t *e)
+{
+ expr_t e2;
+
+ str = power (str, e);
+ for (;;)
+ {
+ str = skipspace (str);
+ switch (str[0])
+ {
+ case 'm':
+ if (match ("mul", str))
+ {
+ str = power (str + 3, &e2);
+ makeexp (e, MULT, *e, e2);
+ break;
+ }
+ if (match ("mod", str))
+ {
+ str = power (str + 3, &e2);
+ makeexp (e, MOD, *e, e2);
+ break;
+ }
+ return str;
+ case 'd':
+ if (match ("div", str))
+ {
+ str = power (str + 3, &e2);
+ makeexp (e, DIV, *e, e2);
+ break;
+ }
+ return str;
+ case 'r':
+ if (match ("rem", str))
+ {
+ str = power (str + 3, &e2);
+ makeexp (e, REM, *e, e2);
+ break;
+ }
+ return str;
+ case 'i':
+ if (match ("invmod", str))
+ {
+ str = power (str + 6, &e2);
+ makeexp (e, REM, *e, e2);
+ break;
+ }
+ return str;
+ case 't':
+ if (match ("times", str))
+ {
+ str = power (str + 5, &e2);
+ makeexp (e, MULT, *e, e2);
+ break;
+ }
+ if (match ("thru", str))
+ {
+ str = power (str + 4, &e2);
+ makeexp (e, DIV, *e, e2);
+ break;
+ }
+ if (match ("through", str))
+ {
+ str = power (str + 7, &e2);
+ makeexp (e, DIV, *e, e2);
+ break;
+ }
+ return str;
+ case '*':
+ str = power (str + 1, &e2);
+ makeexp (e, MULT, *e, e2);
+ break;
+ case '/':
+ str = power (str + 1, &e2);
+ makeexp (e, DIV, *e, e2);
+ break;
+ case '%':
+ str = power (str + 1, &e2);
+ makeexp (e, MOD, *e, e2);
+ break;
+ default:
+ return str;
+ }
+ }
+}
+
+char *
+power (char *str, expr_t *e)
+{
+ expr_t e2;
+
+ str = factor (str, e);
+ while (str[0] == '!')
+ {
+ str++;
+ makeexp (e, FAC, *e, NULL);
+ }
+ str = skipspace (str);
+ if (str[0] == '^')
+ {
+ str = power (str + 1, &e2);
+ makeexp (e, POW, *e, e2);
+ }
+ return str;
+}
+
+int
+match (char *s, char *str)
+{
+ char *ostr = str;
+ int i;
+
+ for (i = 0; s[i] != 0; i++)
+ {
+ if (str[i] != s[i])
+ return 0;
+ }
+ str = skipspace (str + i);
+ return str - ostr;
+}
+
+int
+matchp (char *s, char *str)
+{
+ char *ostr = str;
+ int i;
+
+ for (i = 0; s[i] != 0; i++)
+ {
+ if (str[i] != s[i])
+ return 0;
+ }
+ str = skipspace (str + i);
+ if (str[0] == '(')
+ return str - ostr + 1;
+ return 0;
+}
+
+struct functions
+{
+ char *spelling;
+ enum op_t op;
+ int arity; /* 1 or 2 means real arity; 0 means arbitrary. */
+};
+
+struct functions fns[] =
+{
+ {"sqrt", SQRT, 1},
+#if __GNU_MP_VERSION >= 2
+ {"root", ROOT, 2},
+ {"popc", POPCNT, 1},
+ {"hamdist", HAMDIST, 2},
+#endif
+ {"gcd", GCD, 0},
+#if __GNU_MP_VERSION > 2 || __GNU_MP_VERSION_MINOR >= 1
+ {"lcm", LCM, 0},
+#endif
+ {"and", AND, 0},
+ {"ior", IOR, 0},
+#if __GNU_MP_VERSION > 2 || __GNU_MP_VERSION_MINOR >= 1
+ {"xor", XOR, 0},
+#endif
+ {"plus", PLUS, 0},
+ {"pow", POW, 2},
+ {"minus", MINUS, 2},
+ {"mul", MULT, 0},
+ {"div", DIV, 2},
+ {"mod", MOD, 2},
+ {"rem", REM, 2},
+#if __GNU_MP_VERSION >= 2
+ {"invmod", INVMOD, 2},
+#endif
+ {"log", LOG, 2},
+ {"log2", LOG2, 1},
+ {"F", FERMAT, 1},
+ {"M", MERSENNE, 1},
+ {"fib", FIBONACCI, 1},
+ {"Fib", FIBONACCI, 1},
+ {"random", RANDOM, 1},
+ {"nextprime", NEXTPRIME, 1},
+ {"binom", BINOM, 2},
+ {"binomial", BINOM, 2},
+ {"fac", FAC, 1},
+ {"fact", FAC, 1},
+ {"factorial", FAC, 1},
+ {"time", TIMING, 1},
+ {"", NOP, 0}
+};
+
+char *
+factor (char *str, expr_t *e)
+{
+ expr_t e1, e2;
+
+ str = skipspace (str);
+
+ if (isalpha (str[0]))
+ {
+ int i;
+ int cnt;
+
+ for (i = 0; fns[i].op != NOP; i++)
+ {
+ if (fns[i].arity == 1)
+ {
+ cnt = matchp (fns[i].spelling, str);
+ if (cnt != 0)
+ {
+ str = expr (str + cnt, &e1);
+ str = skipspace (str);
+ if (str[0] != ')')
+ {
+ error = "expected `)'";
+ longjmp (errjmpbuf, (int) (long) str);
+ }
+ makeexp (e, fns[i].op, e1, NULL);
+ return str + 1;
+ }
+ }
+ }
+
+ for (i = 0; fns[i].op != NOP; i++)
+ {
+ if (fns[i].arity != 1)
+ {
+ cnt = matchp (fns[i].spelling, str);
+ if (cnt != 0)
+ {
+ str = expr (str + cnt, &e1);
+ str = skipspace (str);
+
+ if (str[0] != ',')
+ {
+ error = "expected `,' and another operand";
+ longjmp (errjmpbuf, (int) (long) str);
+ }
+
+ str = skipspace (str + 1);
+ str = expr (str, &e2);
+ str = skipspace (str);
+
+ if (fns[i].arity == 0)
+ {
+ while (str[0] == ',')
+ {
+ makeexp (&e1, fns[i].op, e1, e2);
+ str = skipspace (str + 1);
+ str = expr (str, &e2);
+ str = skipspace (str);
+ }
+ }
+
+ if (str[0] != ')')
+ {
+ error = "expected `)'";
+ longjmp (errjmpbuf, (int) (long) str);
+ }
+
+ makeexp (e, fns[i].op, e1, e2);
+ return str + 1;
+ }
+ }
+ }
+ }
+
+ if (str[0] == '(')
+ {
+ str = expr (str + 1, e);
+ str = skipspace (str);
+ if (str[0] != ')')
+ {
+ error = "expected `)'";
+ longjmp (errjmpbuf, (int) (long) str);
+ }
+ str++;
+ }
+ else if (str[0] >= '0' && str[0] <= '9')
+ {
+ expr_t res;
+ char *s, *sc;
+
+ res = malloc (sizeof (struct expr));
+ res -> op = LIT;
+ mpz_init (res->operands.val);
+
+ s = str;
+ while (isalnum (str[0]))
+ str++;
+ sc = malloc (str - s + 1);
+ memcpy (sc, s, str - s);
+ sc[str - s] = 0;
+
+ mpz_set_str (res->operands.val, sc, 0);
+ *e = res;
+ free (sc);
+ }
+ else
+ {
+ error = "operand expected";
+ longjmp (errjmpbuf, (int) (long) str);
+ }
+ return str;
+}
+
+char *
+skipspace (char *str)
+{
+ while (str[0] == ' ')
+ str++;
+ return str;
+}
+
+/* Make a new expression with operation OP and right hand side
+ RHS and left hand side lhs. Put the result in R. */
+void
+makeexp (expr_t *r, enum op_t op, expr_t lhs, expr_t rhs)
+{
+ expr_t res;
+ res = malloc (sizeof (struct expr));
+ res -> op = op;
+ res -> operands.ops.lhs = lhs;
+ res -> operands.ops.rhs = rhs;
+ *r = res;
+ return;
+}
+
+/* Free the memory used by expression E. */
+void
+free_expr (expr_t e)
+{
+ if (e->op != LIT)
+ {
+ free_expr (e->operands.ops.lhs);
+ if (e->operands.ops.rhs != NULL)
+ free_expr (e->operands.ops.rhs);
+ }
+ else
+ {
+ mpz_clear (e->operands.val);
+ }
+}
+
+/* Evaluate the expression E and put the result in R. */
+void
+mpz_eval_expr (mpz_ptr r, expr_t e)
+{
+ mpz_t lhs, rhs;
+
+ switch (e->op)
+ {
+ case LIT:
+ mpz_set (r, e->operands.val);
+ return;
+ case PLUS:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ mpz_add (r, lhs, rhs);
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+ case MINUS:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ mpz_sub (r, lhs, rhs);
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+ case MULT:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ mpz_mul (r, lhs, rhs);
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+ case DIV:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ mpz_fdiv_q (r, lhs, rhs);
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+ case MOD:
+ mpz_init (rhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ mpz_abs (rhs, rhs);
+ mpz_eval_mod_expr (r, e->operands.ops.lhs, rhs);
+ mpz_clear (rhs);
+ return;
+ case REM:
+ /* Check if lhs operand is POW expression and optimize for that case. */
+ if (e->operands.ops.lhs->op == POW)
+ {
+ mpz_t powlhs, powrhs;
+ mpz_init (powlhs);
+ mpz_init (powrhs);
+ mpz_init (rhs);
+ mpz_eval_expr (powlhs, e->operands.ops.lhs->operands.ops.lhs);
+ mpz_eval_expr (powrhs, e->operands.ops.lhs->operands.ops.rhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ mpz_powm (r, powlhs, powrhs, rhs);
+ if (mpz_cmp_si (rhs, 0L) < 0)
+ mpz_neg (r, r);
+ mpz_clear (powlhs);
+ mpz_clear (powrhs);
+ mpz_clear (rhs);
+ return;
+ }
+
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ mpz_fdiv_r (r, lhs, rhs);
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+#if __GNU_MP_VERSION >= 2
+ case INVMOD:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ mpz_invert (r, lhs, rhs);
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+#endif
+ case POW:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ if (mpz_cmpabs_ui (lhs, 1) <= 0)
+ {
+ /* For 0^rhs and 1^rhs, we just need to verify that
+ rhs is well-defined. For (-1)^rhs we need to
+ determine (rhs mod 2). For simplicity, compute
+ (rhs mod 2) for all three cases. */
+ expr_t two, et;
+ two = malloc (sizeof (struct expr));
+ two -> op = LIT;
+ mpz_init_set_ui (two->operands.val, 2L);
+ makeexp (&et, MOD, e->operands.ops.rhs, two);
+ e->operands.ops.rhs = et;
+ }
+
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ if (mpz_cmp_si (rhs, 0L) == 0)
+ /* x^0 is 1 */
+ mpz_set_ui (r, 1L);
+ else if (mpz_cmp_si (lhs, 0L) == 0)
+ /* 0^y (where y != 0) is 0 */
+ mpz_set_ui (r, 0L);
+ else if (mpz_cmp_ui (lhs, 1L) == 0)
+ /* 1^y is 1 */
+ mpz_set_ui (r, 1L);
+ else if (mpz_cmp_si (lhs, -1L) == 0)
+ /* (-1)^y just depends on whether y is even or odd */
+ mpz_set_si (r, (mpz_get_ui (rhs) & 1) ? -1L : 1L);
+ else if (mpz_cmp_si (rhs, 0L) < 0)
+ /* x^(-n) is 0 */
+ mpz_set_ui (r, 0L);
+ else
+ {
+ unsigned long int cnt;
+ unsigned long int y;
+ /* error if exponent does not fit into an unsigned long int. */
+ if (mpz_cmp_ui (rhs, ~(unsigned long int) 0) > 0)
+ goto pow_err;
+
+ y = mpz_get_ui (rhs);
+ /* x^y == (x/(2^c))^y * 2^(c*y) */
+#if __GNU_MP_VERSION >= 2
+ cnt = mpz_scan1 (lhs, 0);
+#else
+ cnt = 0;
+#endif
+ if (cnt != 0)
+ {
+ if (y * cnt / cnt != y)
+ goto pow_err;
+ mpz_tdiv_q_2exp (lhs, lhs, cnt);
+ mpz_pow_ui (r, lhs, y);
+ mpz_mul_2exp (r, r, y * cnt);
+ }
+ else
+ mpz_pow_ui (r, lhs, y);
+ }
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+ pow_err:
+ error = "result of `pow' operator too large";
+ mpz_clear (lhs); mpz_clear (rhs);
+ longjmp (errjmpbuf, 1);
+ case GCD:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ mpz_gcd (r, lhs, rhs);
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+#if __GNU_MP_VERSION > 2 || __GNU_MP_VERSION_MINOR >= 1
+ case LCM:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ mpz_lcm (r, lhs, rhs);
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+#endif
+ case AND:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ mpz_and (r, lhs, rhs);
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+ case IOR:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ mpz_ior (r, lhs, rhs);
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+#if __GNU_MP_VERSION > 2 || __GNU_MP_VERSION_MINOR >= 1
+ case XOR:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ mpz_xor (r, lhs, rhs);
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+#endif
+ case NEG:
+ mpz_eval_expr (r, e->operands.ops.lhs);
+ mpz_neg (r, r);
+ return;
+ case NOT:
+ mpz_eval_expr (r, e->operands.ops.lhs);
+ mpz_com (r, r);
+ return;
+ case SQRT:
+ mpz_init (lhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ if (mpz_sgn (lhs) < 0)
+ {
+ error = "cannot take square root of negative numbers";
+ mpz_clear (lhs);
+ longjmp (errjmpbuf, 1);
+ }
+ mpz_sqrt (r, lhs);
+ return;
+#if __GNU_MP_VERSION > 2 || __GNU_MP_VERSION_MINOR >= 1
+ case ROOT:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ if (mpz_sgn (rhs) <= 0)
+ {
+ error = "cannot take non-positive root orders";
+ mpz_clear (lhs); mpz_clear (rhs);
+ longjmp (errjmpbuf, 1);
+ }
+ if (mpz_sgn (lhs) < 0 && (mpz_get_ui (rhs) & 1) == 0)
+ {
+ error = "cannot take even root orders of negative numbers";
+ mpz_clear (lhs); mpz_clear (rhs);
+ longjmp (errjmpbuf, 1);
+ }
+
+ {
+ unsigned long int nth = mpz_get_ui (rhs);
+ if (mpz_cmp_ui (rhs, ~(unsigned long int) 0) > 0)
+ {
+ /* If we are asked to take an awfully large root order, cheat and
+ ask for the largest order we can pass to mpz_root. This saves
+ some error prone special cases. */
+ nth = ~(unsigned long int) 0;
+ }
+ mpz_root (r, lhs, nth);
+ }
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+#endif
+ case FAC:
+ mpz_eval_expr (r, e->operands.ops.lhs);
+ if (mpz_size (r) > 1)
+ {
+ error = "result of `!' operator too large";
+ longjmp (errjmpbuf, 1);
+ }
+ mpz_fac_ui (r, mpz_get_ui (r));
+ return;
+#if __GNU_MP_VERSION >= 2
+ case POPCNT:
+ mpz_eval_expr (r, e->operands.ops.lhs);
+ { long int cnt;
+ cnt = mpz_popcount (r);
+ mpz_set_si (r, cnt);
+ }
+ return;
+ case HAMDIST:
+ { long int cnt;
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ cnt = mpz_hamdist (lhs, rhs);
+ mpz_clear (lhs); mpz_clear (rhs);
+ mpz_set_si (r, cnt);
+ }
+ return;
+#endif
+ case LOG2:
+ mpz_eval_expr (r, e->operands.ops.lhs);
+ { unsigned long int cnt;
+ if (mpz_sgn (r) <= 0)
+ {
+ error = "logarithm of non-positive number";
+ longjmp (errjmpbuf, 1);
+ }
+ cnt = mpz_sizeinbase (r, 2);
+ mpz_set_ui (r, cnt - 1);
+ }
+ return;
+ case LOG:
+ { unsigned long int cnt;
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ if (mpz_sgn (lhs) <= 0)
+ {
+ error = "logarithm of non-positive number";
+ mpz_clear (lhs); mpz_clear (rhs);
+ longjmp (errjmpbuf, 1);
+ }
+ if (mpz_cmp_ui (rhs, 256) >= 0)
+ {
+ error = "logarithm base too large";
+ mpz_clear (lhs); mpz_clear (rhs);
+ longjmp (errjmpbuf, 1);
+ }
+ cnt = mpz_sizeinbase (lhs, mpz_get_ui (rhs));
+ mpz_set_ui (r, cnt - 1);
+ mpz_clear (lhs); mpz_clear (rhs);
+ }
+ return;
+ case FERMAT:
+ {
+ unsigned long int t;
+ mpz_init (lhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ t = (unsigned long int) 1 << mpz_get_ui (lhs);
+ if (mpz_cmp_ui (lhs, ~(unsigned long int) 0) > 0 || t == 0)
+ {
+ error = "too large Mersenne number index";
+ mpz_clear (lhs);
+ longjmp (errjmpbuf, 1);
+ }
+ mpz_set_ui (r, 1);
+ mpz_mul_2exp (r, r, t);
+ mpz_add_ui (r, r, 1);
+ mpz_clear (lhs);
+ }
+ return;
+ case MERSENNE:
+ mpz_init (lhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ if (mpz_cmp_ui (lhs, ~(unsigned long int) 0) > 0)
+ {
+ error = "too large Mersenne number index";
+ mpz_clear (lhs);
+ longjmp (errjmpbuf, 1);
+ }
+ mpz_set_ui (r, 1);
+ mpz_mul_2exp (r, r, mpz_get_ui (lhs));
+ mpz_sub_ui (r, r, 1);
+ mpz_clear (lhs);
+ return;
+ case FIBONACCI:
+ { mpz_t t;
+ unsigned long int n, i;
+ mpz_init (lhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ if (mpz_sgn (lhs) <= 0 || mpz_cmp_si (lhs, 1000000000) > 0)
+ {
+ error = "Fibonacci index out of range";
+ mpz_clear (lhs);
+ longjmp (errjmpbuf, 1);
+ }
+ n = mpz_get_ui (lhs);
+ mpz_clear (lhs);
+
+#if __GNU_MP_VERSION > 2 || __GNU_MP_VERSION_MINOR >= 1
+ mpz_fib_ui (r, n);
+#else
+ mpz_init_set_ui (t, 1);
+ mpz_set_ui (r, 1);
+
+ if (n <= 2)
+ mpz_set_ui (r, 1);
+ else
+ {
+ for (i = 3; i <= n; i++)
+ {
+ mpz_add (t, t, r);
+ mpz_swap (t, r);
+ }
+ }
+ mpz_clear (t);
+#endif
+ }
+ return;
+ case RANDOM:
+ {
+ unsigned long int n;
+ mpz_init (lhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ if (mpz_sgn (lhs) <= 0 || mpz_cmp_si (lhs, 1000000000) > 0)
+ {
+ error = "random number size out of range";
+ mpz_clear (lhs);
+ longjmp (errjmpbuf, 1);
+ }
+ n = mpz_get_ui (lhs);
+ mpz_clear (lhs);
+ mpz_urandomb (r, rstate, n);
+ }
+ return;
+ case NEXTPRIME:
+ {
+ mpz_eval_expr (r, e->operands.ops.lhs);
+ mpz_nextprime (r, r);
+ }
+ return;
+ case BINOM:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_expr (lhs, e->operands.ops.lhs);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ {
+ unsigned long int k;
+ if (mpz_cmp_ui (rhs, ~(unsigned long int) 0) > 0)
+ {
+ error = "k too large in (n over k) expression";
+ mpz_clear (lhs); mpz_clear (rhs);
+ longjmp (errjmpbuf, 1);
+ }
+ k = mpz_get_ui (rhs);
+ mpz_bin_ui (r, lhs, k);
+ }
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+ case TIMING:
+ {
+ int t0;
+ t0 = cputime ();
+ mpz_eval_expr (r, e->operands.ops.lhs);
+ printf ("time: %d\n", cputime () - t0);
+ }
+ return;
+ default:
+ abort ();
+ }
+}
+
+/* Evaluate the expression E modulo MOD and put the result in R. */
+void
+mpz_eval_mod_expr (mpz_ptr r, expr_t e, mpz_ptr mod)
+{
+ mpz_t lhs, rhs;
+
+ switch (e->op)
+ {
+ case POW:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_mod_expr (lhs, e->operands.ops.lhs, mod);
+ mpz_eval_expr (rhs, e->operands.ops.rhs);
+ mpz_powm (r, lhs, rhs, mod);
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+ case PLUS:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_mod_expr (lhs, e->operands.ops.lhs, mod);
+ mpz_eval_mod_expr (rhs, e->operands.ops.rhs, mod);
+ mpz_add (r, lhs, rhs);
+ if (mpz_cmp_si (r, 0L) < 0)
+ mpz_add (r, r, mod);
+ else if (mpz_cmp (r, mod) >= 0)
+ mpz_sub (r, r, mod);
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+ case MINUS:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_mod_expr (lhs, e->operands.ops.lhs, mod);
+ mpz_eval_mod_expr (rhs, e->operands.ops.rhs, mod);
+ mpz_sub (r, lhs, rhs);
+ if (mpz_cmp_si (r, 0L) < 0)
+ mpz_add (r, r, mod);
+ else if (mpz_cmp (r, mod) >= 0)
+ mpz_sub (r, r, mod);
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+ case MULT:
+ mpz_init (lhs); mpz_init (rhs);
+ mpz_eval_mod_expr (lhs, e->operands.ops.lhs, mod);
+ mpz_eval_mod_expr (rhs, e->operands.ops.rhs, mod);
+ mpz_mul (r, lhs, rhs);
+ mpz_mod (r, r, mod);
+ mpz_clear (lhs); mpz_clear (rhs);
+ return;
+ default:
+ mpz_init (lhs);
+ mpz_eval_expr (lhs, e);
+ mpz_mod (r, lhs, mod);
+ mpz_clear (lhs);
+ return;
+ }
+}
+
+void
+cleanup_and_exit (int sig)
+{
+ switch (sig) {
+#ifdef LIMIT_RESOURCE_USAGE
+ case SIGXCPU:
+ printf ("expression took too long to evaluate%s\n", newline);
+ break;
+#endif
+ case SIGFPE:
+ printf ("divide by zero%s\n", newline);
+ break;
+ default:
+ printf ("expression required too much memory to evaluate%s\n", newline);
+ break;
+ }
+ exit (-2);
+}