commit 6223d322fd5c396f0c3f255b1fb332f7f6b7c74a
parent cbc841acb50d97dd48757a24afcf3bbe8af39c18
Author: Andres Navarro <canavarro82@gmail.com>
Date: Tue, 12 Apr 2011 15:13:54 -0300
Merged branch imath
Diffstat:
12 files changed, 3811 insertions(+), 331 deletions(-)
diff --git a/COPYRIGHT b/COPYRIGHT
@@ -3,12 +3,15 @@ klisp License
klisp is licensed under the terms of the MIT license reproduced below.
This means that klisp is free software and can be used for both academic
-and commercial purposes at absolutely no cost.
+and commercial purposes at absolutely no cost.
+The two projects whose code klisp uses, Lua & IMath, are also distributed
+under the MIT license.
===============================================================================
-Copyright (C) 2011 Andres Navarro
+klisp Parts: Copyright (C) 2011 Andres Navarro.
Lua Parts: Copyright (C) 1994-2010 Lua.org, PUC-Rio.
+IMath Parts: Copyright (C) 2002-2007 Michael J. Fromberger.
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
diff --git a/README b/README
@@ -7,7 +7,8 @@ README for klisp 0.1
the "Revised(-1) Report on the Kernel Programming Language", but
that probably won't happen for some time. It is written in C99 under
the MIT license. It draws heavily from the Lua interpreter source
- code & file structure.
+ code & file structure. It uses the IMath library for arbitrary sized
+ integers and (soon) rationals.
* What is the Kernel Programming Language?
----------------------------------------
diff --git a/src/Makefile b/src/Makefile
@@ -15,7 +15,7 @@ CORE_O= kobject.o ktoken.o kpair.o kstring.o ksymbol.o kread.o \
kgsymbols.o kgcontrol.o kgpairs_lists.o kgpair_mut.o kgenvironments.o \
kgenv_mut.o kgcombiners.o kgcontinuations.o kgencapsulations.o \
kgpromises.o kgkd_vars.o kgks_vars.o kgports.o kgchars.o kgnumbers.o \
- kgstrings.o
+ kgstrings.o imath.o
KRN_T= klisp
KRN_O= klisp.o
@@ -45,7 +45,7 @@ klisp.o: klisp.c klisp.h kobject.h kread.h kwrite.h klimits.h kstate.h kmem.h \
kobject.o: kobject.c kobject.h klisp.h
ktoken.o: ktoken.c ktoken.h kobject.h kstate.h kpair.h kstring.h ksymbol.h \
kerror.h klisp.h kinteger.h
-kinteger.o: kinteger.c kinteger.h kobject.h kstate.h kmem.h klisp.h
+kinteger.o: kinteger.c kinteger.h kobject.h kstate.h kmem.h klisp.h imath.h
kpair.o: kpair.c kpair.h kobject.h kstate.h kmem.h klisp.h
kstring.o: kstring.c kstring.h kobject.h kstate.h kmem.h klisp.h
# XXX: kpair.h because of use of list as symbol table
@@ -54,7 +54,7 @@ ksymbol.o: ksymbol.c ksymbol.h kobject.h kpair.h kstring.h kstate.h kmem.h \
kread.o: kread.c kread.h kobject.h ktoken.h kpair.h kstate.h kerror.h klisp.h \
kport.h
kwrite.o: kwrite.c kwrite.h kobject.h kpair.h kstring.h kstate.h kerror.h \
- klisp.h kport.h
+ klisp.h kport.h kinteger.h
kstate.o: kstate.c kstate.h klisp.h kobject.h kmem.h kstring.h klisp.h \
kground.h kenvironment.h kpair.h keval.h koperative.h kground.h \
krepl.h kcontinuation.h kapplicative.h kport.h ksymbol.h kport.h \
@@ -139,3 +139,4 @@ kgnumbers.o: kgnumbers.c kgnumbers.h kghelpers.h kstate.h klisp.h \
kgstrings.o: kgstrings.c kgstrings.h kghelpers.h kstate.h klisp.h \
kobject.h kerror.h kapplicative.h koperative.h kcontinuation.h \
ksymbol.h
+imath.o: kobject.h kstate.h kmem.h kerror.h
+\ No newline at end of file
diff --git a/src/imath.c b/src/imath.c
@@ -0,0 +1,3438 @@
+/*
+** imath.c
+** Arbitrary precision integer arithmetic routines.
+** See Copyright Notice in klisp.h
+*/
+
+/*
+** SOURCE NOTE: This is mostly from the IMath library, written by
+** M.J. Fromberger. It is adapted to klisp, mainly in the use of the
+** klisp allocator and fixing of digit size to 32 bits.
+** Imported from version (1.15) updated 01-Feb-2011 at 03:10 PM.
+*/
+
+#include "imath.h"
+
+#if DEBUG
+#include <stdio.h>
+#endif
+
+#include <stdlib.h>
+#include <string.h>
+#include <ctype.h>
+
+#include <assert.h>
+
+#if DEBUG
+#define STATIC /* public */
+#else
+#define STATIC static
+#endif
+
+/* Andres Navarro: klisp includes */
+#include "kobject.h"
+#include "kstate.h"
+#include "kmem.h"
+#include "kerror.h"
+
+/* {{{ Constants */
+
+const mp_result MP_OK = 0; /* no error, all is well */
+const mp_result MP_FALSE = 0; /* boolean false */
+const mp_result MP_TRUE = -1; /* boolean true */
+const mp_result MP_MEMORY = -2; /* out of memory */
+const mp_result MP_RANGE = -3; /* argument out of range */
+const mp_result MP_UNDEF = -4; /* result undefined */
+const mp_result MP_TRUNC = -5; /* output truncated */
+const mp_result MP_BADARG = -6; /* invalid null argument */
+const mp_result MP_MINERR = -6;
+
+const mp_sign MP_NEG = 1; /* value is strictly negative */
+const mp_sign MP_ZPOS = 0; /* value is non-negative */
+
+STATIC const char *s_unknown_err = "unknown result code";
+STATIC const char *s_error_msg[] = {
+ "error code 0",
+ "boolean true",
+ "out of memory",
+ "argument out of range",
+ "result undefined",
+ "output truncated",
+ "invalid argument",
+ NULL
+};
+
+/* }}} */
+
+/* Argument checking macros
+ Use CHECK() where a return value is required; NRCHECK() elsewhere */
+#define CHECK(TEST) assert(TEST)
+#define NRCHECK(TEST) assert(TEST)
+
+/* {{{ Logarithm table for computing output sizes */
+
+/* The ith entry of this table gives the value of log_i(2).
+
+ An integer value n requires ceil(log_i(n)) digits to be represented
+ in base i. Since it is easy to compute lg(n), by counting bits, we
+ can compute log_i(n) = lg(n) * log_i(2).
+
+ The use of this table eliminates a dependency upon linkage against
+ the standard math libraries.
+ */
+STATIC const double s_log2[] = {
+ 0.000000000, 0.000000000, 1.000000000, 0.630929754, /* 0 1 2 3 */
+ 0.500000000, 0.430676558, 0.386852807, 0.356207187, /* 4 5 6 7 */
+ 0.333333333, 0.315464877, 0.301029996, 0.289064826, /* 8 9 10 11 */
+ 0.278942946, 0.270238154, 0.262649535, 0.255958025, /* 12 13 14 15 */
+ 0.250000000, 0.244650542, 0.239812467, 0.235408913, /* 16 17 18 19 */
+ 0.231378213, 0.227670249, 0.224243824, 0.221064729, /* 20 21 22 23 */
+ 0.218104292, 0.215338279, 0.212746054, 0.210309918, /* 24 25 26 27 */
+ 0.208014598, 0.205846832, 0.203795047, 0.201849087, /* 28 29 30 31 */
+ 0.200000000, 0.198239863, 0.196561632, 0.194959022, /* 32 33 34 35 */
+ 0.193426404, /* 36 */
+};
+
+/* }}} */
+/* {{{ Various macros */
+
+/* Return the number of digits needed to represent a static value */
+#define MP_VALUE_DIGITS(V) \
+((sizeof(V)+(sizeof(mp_digit)-1))/sizeof(mp_digit))
+
+/* Round precision P to nearest word boundary */
+#define ROUND_PREC(P) ((mp_size)(2*(((P)+1)/2)))
+
+/* Set array P of S digits to zero */
+#define ZERO(P, S) \
+do{ \
+ mp_size i__ = (S) * sizeof(mp_digit); \
+ mp_digit *p__ = (P); \
+ memset(p__, 0, i__); \
+} while(0)
+
+/* Copy S digits from array P to array Q */
+#define COPY(P, Q, S) \
+do{ \
+ mp_size i__ = (S) * sizeof(mp_digit); \
+ mp_digit *p__ = (P), *q__ = (Q); \
+ memcpy(q__, p__, i__); \
+} while(0)
+
+/* Reverse N elements of type T in array A */
+#define REV(T, A, N) \
+do{ \
+ T *u_ = (A), *v_ = u_ + (N) - 1; \
+ while (u_ < v_) { \
+ T xch = *u_; \
+ *u_++ = *v_; \
+ *v_-- = xch; \
+ } \
+} while(0)
+
+#define CLAMP(Z) \
+do{ \
+ mp_int z_ = (Z); \
+ mp_size uz_ = MP_USED(z_); \
+ mp_digit *dz_ = MP_DIGITS(z_) + uz_ -1; \
+ while (uz_ > 1 && (*dz_-- == 0)) \
+ --uz_; \
+ MP_USED(z_) = uz_; \
+} while(0)
+
+/* Select min/max. Do not provide expressions for which multiple
+ evaluation would be problematic, e.g. x++ */
+#define MIN(A, B) ((B)<(A)?(B):(A))
+#define MAX(A, B) ((B)>(A)?(B):(A))
+
+/* Exchange lvalues A and B of type T, e.g.
+ SWAP(int, x, y) where x and y are variables of type int. */
+#define SWAP(T, A, B) \
+do{ \
+ T t_ = (A); \
+ A = (B); \
+ B = t_; \
+} while(0)
+
+/* Used to set up and access simple temp stacks within functions. */
+#define TEMP(K) (temp + (K))
+#define SETUP(E, C) \
+do{ \
+ if ((res = (E)) != MP_OK) \
+ goto CLEANUP; \
+ ++(C); \
+} while(0)
+
+/* Compare value to zero. */
+#define CMPZ(Z) \
+(((Z)->used==1&&(Z)->digits[0]==0)?0:((Z)->sign==MP_NEG)?-1:1)
+
+/* Multiply X by Y into Z, ignoring signs. Requires that Z have
+ enough storage preallocated to hold the result. */
+#define UMUL(K, X, Y, Z) \
+do{ \
+ mp_size ua_ = MP_USED(X), ub_ = MP_USED(Y); \
+ mp_size o_ = ua_ + ub_; \
+ ZERO(MP_DIGITS(Z), o_); \
+ (void) s_kmul((K), MP_DIGITS(X), MP_DIGITS(Y), MP_DIGITS(Z), ua_, ub_); \
+ MP_USED(Z) = o_; \
+ CLAMP(Z); \
+} while(0)
+
+/* Square X into Z. Requires that Z have enough storage to hold the
+ result. */
+#define USQR(K, X, Z) \
+do{ \
+ mp_size ua_ = MP_USED(X), o_ = ua_ + ua_; \
+ ZERO(MP_DIGITS(Z), o_); \
+ (void) s_ksqr((K), MP_DIGITS(X), MP_DIGITS(Z), ua_); \
+ MP_USED(Z) = o_; \
+ CLAMP(Z); \
+} while(0)
+
+#define UPPER_HALF(W) ((mp_word)((W) >> MP_DIGIT_BIT))
+#define LOWER_HALF(W) ((mp_digit)(W))
+#define HIGH_BIT_SET(W) ((W) >> (MP_WORD_BIT - 1))
+#define ADD_WILL_OVERFLOW(W, V) ((MP_WORD_MAX - (V)) < (W))
+
+/* }}} */
+/* {{{ Default configuration settings */
+
+/* Default number of digits allocated to a new mp_int */
+#if IMATH_TEST
+mp_size default_precision = MP_DEFAULT_PREC;
+#else
+STATIC const mp_size default_precision = MP_DEFAULT_PREC;
+#endif
+
+/* Minimum number of digits to invoke recursive multiply */
+#if IMATH_TEST
+mp_size multiply_threshold = MP_MULT_THRESH;
+#else
+STATIC const mp_size multiply_threshold = MP_MULT_THRESH;
+#endif
+
+/* }}} */
+
+/* Allocate a buffer of (at least) num digits, or return
+ NULL if that couldn't be done. */
+STATIC mp_digit *s_alloc(klisp_State *K, mp_size num);
+
+/* Release a buffer of digits allocated by s_alloc(). */
+STATIC void s_free(klisp_State *K, void *ptr, mp_size size);
+
+/* Insure that z has at least min digits allocated, resizing if
+ necessary. Returns true if successful, false if out of memory. */
+STATIC int s_pad(klisp_State *K, mp_int z, mp_size min);
+
+/* Fill in a "fake" mp_int on the stack with a given value */
+STATIC void s_fake(mp_int z, mp_small value, mp_digit vbuf[]);
+
+/* Compare two runs of digits of given length, returns <0, 0, >0 */
+STATIC int s_cdig(mp_digit *da, mp_digit *db, mp_size len);
+
+/* Pack the unsigned digits of v into array t */
+STATIC int s_vpack(mp_small v, mp_digit t[]);
+
+/* Compare magnitudes of a and b, returns <0, 0, >0 */
+STATIC int s_ucmp(mp_int a, mp_int b);
+
+/* Compare magnitudes of a and v, returns <0, 0, >0 */
+STATIC int s_vcmp(mp_int a, mp_small v);
+
+/* Unsigned magnitude addition; assumes dc is big enough.
+ Carry out is returned (no memory allocated). */
+STATIC mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc,
+ mp_size size_a, mp_size size_b);
+
+/* Unsigned magnitude subtraction. Assumes dc is big enough. */
+STATIC void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc,
+ mp_size size_a, mp_size size_b);
+
+/* Unsigned recursive multiplication. Assumes dc is big enough. */
+STATIC int s_kmul(klisp_State *K, mp_digit *da, mp_digit *db,
+ mp_digit *dc, mp_size size_a, mp_size size_b);
+
+/* Unsigned magnitude multiplication. Assumes dc is big enough. */
+STATIC void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc,
+ mp_size size_a, mp_size size_b);
+
+/* Unsigned recursive squaring. Assumes dc is big enough. */
+/* Andres Navarro: allocs temporaries */
+STATIC int s_ksqr(klisp_State *K, mp_digit *da, mp_digit *dc,
+ mp_size size_a);
+
+/* Unsigned magnitude squaring. Assumes dc is big enough. */
+STATIC void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a);
+
+/* Single digit addition. Assumes a is big enough. */
+STATIC void s_dadd(mp_int a, mp_digit b);
+
+/* Single digit multiplication. Assumes a is big enough. */
+STATIC void s_dmul(mp_int a, mp_digit b);
+
+/* Single digit multiplication on buffers; assumes dc is big enough. */
+STATIC void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc,
+ mp_size size_a);
+
+/* Single digit division. Replaces a with the quotient,
+ returns the remainder. */
+STATIC mp_digit s_ddiv(mp_int a, mp_digit b);
+
+/* Quick division by a power of 2, replaces z (no allocation) */
+STATIC void s_qdiv(mp_int z, mp_size p2);
+
+/* Quick remainder by a power of 2, replaces z (no allocation) */
+STATIC void s_qmod(mp_int z, mp_size p2);
+
+/* Quick multiplication by a power of 2, replaces z.
+ Allocates if necessary; returns false in case this fails. */
+STATIC int s_qmul(klisp_State *K, mp_int z, mp_size p2);
+
+/* Quick subtraction from a power of 2, replaces z.
+ Allocates if necessary; returns false in case this fails. */
+STATIC int s_qsub(klisp_State *K, mp_int z, mp_size p2);
+
+/* Return maximum k such that 2^k divides z. */
+STATIC int s_dp2k(mp_int z);
+
+/* Return k >= 0 such that z = 2^k, or -1 if there is no such k. */
+STATIC int s_isp2(mp_int z);
+
+/* Set z to 2^k. May allocate; returns false in case this fails. */
+STATIC int s_2expt(klisp_State *K, mp_int z, mp_small k);
+
+/* Normalize a and b for division, returns normalization constant */
+STATIC int s_norm(klisp_State *K, mp_int a, mp_int b);
+
+/* Compute constant mu for Barrett reduction, given modulus m, result
+ replaces z, m is untouched. */
+STATIC mp_result s_brmu(klisp_State *K, mp_int z, mp_int m);
+
+/* Reduce a modulo m, using Barrett's algorithm. */
+STATIC int s_reduce(klisp_State *K, mp_int x, mp_int m, mp_int mu,
+ mp_int q1, mp_int q2);
+
+/* Modular exponentiation, using Barrett reduction */
+STATIC mp_result s_embar(klisp_State *K, mp_int a, mp_int b, mp_int m,
+ mp_int mu, mp_int c);
+
+/* Unsigned magnitude division. Assumes |a| > |b|. Allocates
+ temporaries; overwrites a with quotient, b with remainder. */
+STATIC mp_result s_udiv(klisp_State *K, mp_int a, mp_int b);
+
+/* Compute the number of digits in radix r required to represent the
+ given value. Does not account for sign flags, terminators, etc. */
+STATIC int s_outlen(mp_int z, mp_size r);
+
+/* Guess how many digits of precision will be needed to represent a
+ radix r value of the specified number of digits. Returns a value
+ guaranteed to be no smaller than the actual number required. */
+STATIC mp_size s_inlen(int len, mp_size r);
+
+/* Convert a character to a digit value in radix r, or
+ -1 if out of range */
+STATIC int s_ch2val(char c, int r);
+
+/* Convert a digit value to a character */
+STATIC char s_val2ch(int v, int caps);
+
+/* Take 2's complement of a buffer in place */
+STATIC void s_2comp(unsigned char *buf, int len);
+
+/* Convert a value to binary, ignoring sign. On input, *limpos is the
+ bound on how many bytes should be written to buf; on output, *limpos
+ is set to the number of bytes actually written. */
+STATIC mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad);
+
+#if DEBUG
+/* Dump a representation of the mp_int to standard output */
+void s_print(char *tag, mp_int z);
+void s_print_buf(char *tag, mp_digit *buf, mp_size num);
+#endif
+
+/* {{{ mp_int_init(z) */
+
+mp_result mp_int_init(mp_int z)
+{
+ if(z == NULL)
+ return MP_BADARG;
+
+ z->single = 0;
+ z->digits = &(z->single);
+ z->alloc = 1;
+ z->used = 1;
+ z->sign = MP_ZPOS;
+
+ return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_alloc() */
+
+mp_int mp_int_alloc(klisp_State *K)
+{
+ mp_int out = klispM_new(K, mpz_t);
+
+ if(out != NULL)
+ mp_int_init(out);
+
+ return out;
+}
+
+/* }}} */
+
+/* {{{ mp_int_init_size(z, prec) */
+
+mp_result mp_int_init_size(klisp_State *K, mp_int z, mp_size prec)
+{
+ CHECK(z != NULL);
+
+ if(prec == 0)
+ prec = default_precision;
+ else if(prec == 1)
+ return mp_int_init(z);
+ else
+ prec = (mp_size) ROUND_PREC(prec);
+
+ if((MP_DIGITS(z) = s_alloc(K, prec)) == NULL)
+ return MP_MEMORY;
+
+ z->digits[0] = 0;
+ MP_USED(z) = 1;
+ MP_ALLOC(z) = prec;
+ MP_SIGN(z) = MP_ZPOS;
+
+ return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_init_copy(z, old) */
+
+mp_result mp_int_init_copy(klisp_State *K, mp_int z, mp_int old)
+{
+ mp_result res;
+ mp_size uold;
+
+ CHECK(z != NULL && old != NULL);
+
+ uold = MP_USED(old);
+ if(uold == 1) {
+ mp_int_init(z);
+ }
+ else {
+ mp_size target = MAX(uold, default_precision);
+
+ if((res = mp_int_init_size(K, z, target)) != MP_OK)
+ return res;
+ }
+
+ MP_USED(z) = uold;
+ MP_SIGN(z) = MP_SIGN(old);
+ COPY(MP_DIGITS(old), MP_DIGITS(z), uold);
+
+ return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_init_value(z, value) */
+
+mp_result mp_int_init_value(klisp_State *K, mp_int z, mp_small value)
+{
+ mpz_t vtmp;
+ mp_digit vbuf[MP_VALUE_DIGITS(value)];
+
+ s_fake(&vtmp, value, vbuf);
+ return mp_int_init_copy(K, z, &vtmp);
+}
+
+/* }}} */
+
+/* {{{ mp_int_set_value(z, value) */
+
+mp_result mp_int_set_value(klisp_State *K, mp_int z, mp_small value)
+{
+ mpz_t vtmp;
+ mp_digit vbuf[MP_VALUE_DIGITS(value)];
+
+ s_fake(&vtmp, value, vbuf);
+ return mp_int_copy(K, &vtmp, z);
+}
+
+/* }}} */
+
+/* {{{ mp_int_clear(z) */
+
+void mp_int_clear(klisp_State *K, mp_int z)
+{
+ if(z == NULL)
+ return;
+
+ if(MP_DIGITS(z) != NULL) {
+ if((void *) MP_DIGITS(z) != (void *) &MP_SINGLE(z))
+ s_free(K, MP_DIGITS(z), MP_ALLOC(z));
+ MP_DIGITS(z) = NULL;
+ }
+}
+
+/* }}} */
+
+/* {{{ mp_int_free(z) */
+
+void mp_int_free(klisp_State *K, mp_int z)
+{
+ NRCHECK(z != NULL);
+
+ mp_int_clear(K, z);
+ klispM_free(K, z); /* note: NOT s_free() */
+}
+
+/* }}} */
+
+/* {{{ mp_int_copy(a, c) */
+
+mp_result mp_int_copy(klisp_State *K, mp_int a, mp_int c)
+{
+ CHECK(a != NULL && c != NULL);
+
+ if(a != c) {
+ mp_size ua = MP_USED(a);
+ mp_digit *da, *dc;
+
+ if(!s_pad(K, c, ua))
+ return MP_MEMORY;
+
+ da = MP_DIGITS(a); dc = MP_DIGITS(c);
+ COPY(da, dc, ua);
+
+ MP_USED(c) = ua;
+ MP_SIGN(c) = MP_SIGN(a);
+ }
+
+ return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_swap(a, c) */
+
+void mp_int_swap(mp_int a, mp_int c)
+{
+ if(a != c) {
+ mpz_t tmp = *a;
+
+ *a = *c;
+ *c = tmp;
+ }
+}
+
+/* }}} */
+
+/* {{{ mp_int_zero(z) */
+
+void mp_int_zero(mp_int z)
+{
+ NRCHECK(z != NULL);
+
+ z->digits[0] = 0;
+ MP_USED(z) = 1;
+ MP_SIGN(z) = MP_ZPOS;
+}
+
+/* }}} */
+
+/* {{{ mp_int_abs(a, c) */
+
+mp_result mp_int_abs(klisp_State *K, mp_int a, mp_int c)
+{
+ mp_result res;
+
+ CHECK(a != NULL && c != NULL);
+
+ if((res = mp_int_copy(K, a, c)) != MP_OK)
+ return res;
+
+ MP_SIGN(c) = MP_ZPOS;
+ return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_neg(a, c) */
+
+mp_result mp_int_neg(klisp_State *K, mp_int a, mp_int c)
+{
+ mp_result res;
+
+ CHECK(a != NULL && c != NULL);
+
+ if((res = mp_int_copy(K, a, c)) != MP_OK)
+ return res;
+
+ if(CMPZ(c) != 0)
+ MP_SIGN(c) = 1 - MP_SIGN(a);
+
+ return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_add(a, b, c) */
+
+mp_result mp_int_add(klisp_State *K, mp_int a, mp_int b, mp_int c)
+{
+ mp_size ua, ub, uc, max;
+
+ CHECK(a != NULL && b != NULL && c != NULL);
+
+ ua = MP_USED(a); ub = MP_USED(b); uc = MP_USED(c);
+ max = MAX(ua, ub);
+
+ if(MP_SIGN(a) == MP_SIGN(b)) {
+ /* Same sign -- add magnitudes, preserve sign of addends */
+ mp_digit carry;
+
+ if(!s_pad(K, c, max))
+ return MP_MEMORY;
+
+ carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, ub);
+ uc = max;
+
+ if(carry) {
+ if(!s_pad(K, c, max + 1))
+ return MP_MEMORY;
+
+ c->digits[max] = carry;
+ ++uc;
+ }
+
+ MP_USED(c) = uc;
+ MP_SIGN(c) = MP_SIGN(a);
+
+ }
+ else {
+ /* Different signs -- subtract magnitudes, preserve sign of greater */
+ mp_int x, y;
+ int cmp = s_ucmp(a, b); /* magnitude comparision, sign ignored */
+
+ /* Set x to max(a, b), y to min(a, b) to simplify later code.
+ A special case yields zero for equal magnitudes.
+ */
+ if(cmp == 0) {
+ mp_int_zero(c);
+ return MP_OK;
+ }
+ else if(cmp < 0) {
+ x = b; y = a;
+ }
+ else {
+ x = a; y = b;
+ }
+
+ if(!s_pad(K, c, MP_USED(x)))
+ return MP_MEMORY;
+
+ /* Subtract smaller from larger */
+ s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y));
+ MP_USED(c) = MP_USED(x);
+ CLAMP(c);
+
+ /* Give result the sign of the larger */
+ MP_SIGN(c) = MP_SIGN(x);
+ }
+
+ return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_add_value(a, value, c) */
+
+mp_result mp_int_add_value(klisp_State *K, mp_int a, mp_small value, mp_int c)
+{
+ mpz_t vtmp;
+ mp_digit vbuf[MP_VALUE_DIGITS(value)];
+
+ s_fake(&vtmp, value, vbuf);
+
+ return mp_int_add(K, a, &vtmp, c);
+}
+
+/* }}} */
+
+/* {{{ mp_int_sub(a, b, c) */
+
+mp_result mp_int_sub(klisp_State *K, mp_int a, mp_int b, mp_int c)
+{
+ mp_size ua, ub, uc, max;
+
+ CHECK(a != NULL && b != NULL && c != NULL);
+
+ ua = MP_USED(a); ub = MP_USED(b); uc = MP_USED(c);
+ max = MAX(ua, ub);
+
+ if(MP_SIGN(a) != MP_SIGN(b)) {
+ /* Different signs -- add magnitudes and keep sign of a */
+ mp_digit carry;
+
+ if(!s_pad(K, c, max))
+ return MP_MEMORY;
+
+ carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, ub);
+ uc = max;
+
+ if(carry) {
+ if(!s_pad(K, c, max + 1))
+ return MP_MEMORY;
+
+ c->digits[max] = carry;
+ ++uc;
+ }
+
+ MP_USED(c) = uc;
+ MP_SIGN(c) = MP_SIGN(a);
+
+ }
+ else {
+ /* Same signs -- subtract magnitudes */
+ mp_int x, y;
+ mp_sign osign;
+ int cmp = s_ucmp(a, b);
+
+ if(!s_pad(K, c, max))
+ return MP_MEMORY;
+
+ if(cmp >= 0) {
+ x = a; y = b; osign = MP_ZPOS;
+ }
+ else {
+ x = b; y = a; osign = MP_NEG;
+ }
+
+ if(MP_SIGN(a) == MP_NEG && cmp != 0)
+ osign = 1 - osign;
+
+ s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y));
+ MP_USED(c) = MP_USED(x);
+ CLAMP(c);
+
+ MP_SIGN(c) = osign;
+ }
+
+ return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_sub_value(a, value, c) */
+
+mp_result mp_int_sub_value(klisp_State *K, mp_int a, mp_small value, mp_int c)
+{
+ mpz_t vtmp;
+ mp_digit vbuf[MP_VALUE_DIGITS(value)];
+
+ s_fake(&vtmp, value, vbuf);
+
+ return mp_int_sub(K, a, &vtmp, c);
+}
+
+/* }}} */
+
+/* {{{ mp_int_mul(a, b, c) */
+
+mp_result mp_int_mul(klisp_State *K, mp_int a, mp_int b, mp_int c)
+{
+ mp_digit *out;
+ mp_size osize, ua, ub, p = 0;
+ mp_sign osign;
+
+ CHECK(a != NULL && b != NULL && c != NULL);
+
+ /* If either input is zero, we can shortcut multiplication */
+ if(mp_int_compare_zero(a) == 0 || mp_int_compare_zero(b) == 0) {
+ mp_int_zero(c);
+ return MP_OK;
+ }
+
+ /* Output is positive if inputs have same sign, otherwise negative */
+ osign = (MP_SIGN(a) == MP_SIGN(b)) ? MP_ZPOS : MP_NEG;
+
+ /* If the output is not identical to any of the inputs, we'll write
+ the results directly; otherwise, allocate a temporary space. */
+ ua = MP_USED(a); ub = MP_USED(b);
+ osize = MAX(ua, ub);
+ osize = 4 * ((osize + 1) / 2);
+
+ if(c == a || c == b) {
+ p = ROUND_PREC(osize);
+ p = MAX(p, default_precision);
+
+ if((out = s_alloc(K, p)) == NULL)
+ return MP_MEMORY;
+ }
+ else {
+ if(!s_pad(K, c, osize))
+ return MP_MEMORY;
+
+ out = MP_DIGITS(c);
+ }
+ ZERO(out, osize);
+
+ if(!s_kmul(K, MP_DIGITS(a), MP_DIGITS(b), out, ua, ub))
+ return MP_MEMORY;
+
+ /* If we allocated a new buffer, get rid of whatever memory c was
+ already using, and fix up its fields to reflect that.
+ */
+ if(out != MP_DIGITS(c)) {
+ if((void *) MP_DIGITS(c) != (void *) &MP_SINGLE(c))
+ s_free(K, MP_DIGITS(c), MP_ALLOC(c));
+ MP_DIGITS(c) = out;
+ MP_ALLOC(c) = p;
+ }
+
+ MP_USED(c) = osize; /* might not be true, but we'll fix it ... */
+ CLAMP(c); /* ... right here */
+ MP_SIGN(c) = osign;
+
+ return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_mul_value(a, value, c) */
+
+mp_result mp_int_mul_value(klisp_State *K, mp_int a, mp_small value,
+ mp_int c)
+{
+ mpz_t vtmp;
+ mp_digit vbuf[MP_VALUE_DIGITS(value)];
+
+ s_fake(&vtmp, value, vbuf);
+
+ return mp_int_mul(K, a, &vtmp, c);
+}
+
+/* }}} */
+
+/* {{{ mp_int_mul_pow2(a, p2, c) */
+
+mp_result mp_int_mul_pow2(klisp_State *K, mp_int a, mp_small p2, mp_int c)
+{
+ mp_result res;
+ CHECK(a != NULL && c != NULL && p2 >= 0);
+
+ if((res = mp_int_copy(K, a, c)) != MP_OK)
+ return res;
+
+ if(s_qmul(K, c, (mp_size) p2))
+ return MP_OK;
+ else
+ return MP_MEMORY;
+}
+
+/* }}} */
+
+/* {{{ mp_int_sqr(a, c) */
+
+mp_result mp_int_sqr(klisp_State *K, mp_int a, mp_int c)
+{
+ mp_digit *out;
+ mp_size osize, p = 0;
+
+ CHECK(a != NULL && c != NULL);
+
+ /* Get a temporary buffer big enough to hold the result */
+ osize = (mp_size) 4 * ((MP_USED(a) + 1) / 2);
+ if(a == c) {
+ p = ROUND_PREC(osize);
+ p = MAX(p, default_precision);
+
+ if((out = s_alloc(K, p)) == NULL)
+ return MP_MEMORY;
+ }
+ else {
+ if(!s_pad(K, c, osize))
+ return MP_MEMORY;
+
+ out = MP_DIGITS(c);
+ }
+ ZERO(out, osize);
+
+ s_ksqr(K, MP_DIGITS(a), out, MP_USED(a));
+
+ /* Get rid of whatever memory c was already using, and fix up its
+ fields to reflect the new digit array it's using
+ */
+ if(out != MP_DIGITS(c)) {
+ if((void *) MP_DIGITS(c) != (void *) &MP_SINGLE(c))
+ s_free(K, MP_DIGITS(c), MP_ALLOC(c));
+ MP_DIGITS(c) = out;
+ MP_ALLOC(c) = p;
+ }
+
+ MP_USED(c) = osize; /* might not be true, but we'll fix it ... */
+ CLAMP(c); /* ... right here */
+ MP_SIGN(c) = MP_ZPOS;
+
+ return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_div(a, b, q, r) */
+
+mp_result mp_int_div(klisp_State *K, mp_int a, mp_int b, mp_int q, mp_int r)
+{
+ int cmp, last = 0, lg;
+ mp_result res = MP_OK;
+ mpz_t temp[2];
+ mp_int qout, rout;
+ mp_sign sa = MP_SIGN(a), sb = MP_SIGN(b);
+
+ CHECK(a != NULL && b != NULL && q != r);
+
+ if(CMPZ(b) == 0)
+ return MP_UNDEF;
+ else if((cmp = s_ucmp(a, b)) < 0) {
+ /* If |a| < |b|, no division is required:
+ q = 0, r = a
+ */
+ if(r && (res = mp_int_copy(K, a, r)) != MP_OK)
+ return res;
+
+ if(q)
+ mp_int_zero(q);
+
+ return MP_OK;
+ }
+ else if(cmp == 0) {
+ /* If |a| = |b|, no division is required:
+ q = 1 or -1, r = 0
+ */
+ if(r)
+ mp_int_zero(r);
+
+ if(q) {
+ mp_int_zero(q);
+ q->digits[0] = 1;
+
+ if(sa != sb)
+ MP_SIGN(q) = MP_NEG;
+ }
+
+ return MP_OK;
+ }
+
+ /* When |a| > |b|, real division is required. We need someplace to
+ store quotient and remainder, but q and r are allowed to be NULL
+ or to overlap with the inputs.
+ */
+ if((lg = s_isp2(b)) < 0) {
+ if(q && b != q) {
+ if((res = mp_int_copy(K, a, q)) != MP_OK)
+ goto CLEANUP;
+ else
+ qout = q;
+ }
+ else {
+ qout = TEMP(last);
+ SETUP(mp_int_init_copy(K, TEMP(last), a), last);
+ }
+
+ if(r && a != r) {
+ if((res = mp_int_copy(K, b, r)) != MP_OK)
+ goto CLEANUP;
+ else
+ rout = r;
+ }
+ else {
+ rout = TEMP(last);
+ SETUP(mp_int_init_copy(K, TEMP(last), b), last);
+ }
+
+ if((res = s_udiv(K, qout, rout)) != MP_OK) goto CLEANUP;
+ }
+ else {
+ if(q && (res = mp_int_copy(K, a, q)) != MP_OK) goto CLEANUP;
+ if(r && (res = mp_int_copy(K, a, r)) != MP_OK) goto CLEANUP;
+
+ if(q) s_qdiv(q, (mp_size) lg); qout = q;
+ if(r) s_qmod(r, (mp_size) lg); rout = r;
+ }
+
+ /* Recompute signs for output */
+ if(rout) {
+ MP_SIGN(rout) = sa;
+ if(CMPZ(rout) == 0)
+ MP_SIGN(rout) = MP_ZPOS;
+ }
+ if(qout) {
+ MP_SIGN(qout) = (sa == sb) ? MP_ZPOS : MP_NEG;
+ if(CMPZ(qout) == 0)
+ MP_SIGN(qout) = MP_ZPOS;
+ }
+
+ if(q && (res = mp_int_copy(K, qout, q)) != MP_OK) goto CLEANUP;
+ if(r && (res = mp_int_copy(K, rout, r)) != MP_OK) goto CLEANUP;
+
+ CLEANUP:
+ while(--last >= 0)
+ mp_int_clear(K, TEMP(last));
+
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_mod(a, m, c) */
+
+mp_result mp_int_mod(klisp_State *K, mp_int a, mp_int m, mp_int c)
+{
+ mp_result res;
+ mpz_t tmp;
+ mp_int out;
+
+ if(m == c) {
+ mp_int_init(&tmp);
+ out = &tmp;
+ }
+ else {
+ out = c;
+ }
+
+ if((res = mp_int_div(K, a, m, NULL, out)) != MP_OK)
+ goto CLEANUP;
+
+ if(CMPZ(out) < 0)
+ res = mp_int_add(K, out, m, c);
+ else
+ res = mp_int_copy(K, out, c);
+
+ CLEANUP:
+ if(out != c)
+ mp_int_clear(K, &tmp);
+
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_div_value(a, value, q, r) */
+
+mp_result mp_int_div_value(klisp_State *K, mp_int a, mp_small value,
+ mp_int q, mp_small *r)
+{
+ mpz_t vtmp, rtmp;
+ mp_digit vbuf[MP_VALUE_DIGITS(value)];
+ mp_result res;
+
+ mp_int_init(&rtmp);
+ s_fake(&vtmp, value, vbuf);
+
+ if((res = mp_int_div(K, a, &vtmp, q, &rtmp)) != MP_OK)
+ goto CLEANUP;
+
+ if(r)
+ (void) mp_int_to_int(&rtmp, r); /* can't fail */
+
+ CLEANUP:
+ mp_int_clear(K, &rtmp);
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_div_pow2(a, p2, q, r) */
+
+mp_result mp_int_div_pow2(klisp_State *K, mp_int a, mp_small p2, mp_int q,
+ mp_int r)
+{
+ mp_result res = MP_OK;
+
+ CHECK(a != NULL && p2 >= 0 && q != r);
+
+ if(q != NULL && (res = mp_int_copy(K, a, q)) == MP_OK)
+ s_qdiv(q, (mp_size) p2);
+
+ if(res == MP_OK && r != NULL && (res = mp_int_copy(K, a, r)) == MP_OK)
+ s_qmod(r, (mp_size) p2);
+
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_expt(a, b, c) */
+
+mp_result mp_int_expt(klisp_State *K, mp_int a, mp_small b, mp_int c)
+{
+ mpz_t t;
+ mp_result res;
+ unsigned int v = abs(b);
+
+ CHECK(c != NULL);
+
+ CHECK(b >= 0 && c != NULL);
+
+ if((res = mp_int_init_copy(K, &t, a)) != MP_OK)
+ return res;
+
+ (void) mp_int_set_value(K, c, 1);
+ while(v != 0) {
+ if(v & 1) {
+ if((res = mp_int_mul(K, c, &t, c)) != MP_OK)
+ goto CLEANUP;
+ }
+
+ v >>= 1;
+ if(v == 0) break;
+
+ if((res = mp_int_sqr(K, &t, &t)) != MP_OK)
+ goto CLEANUP;
+ }
+
+ CLEANUP:
+ mp_int_clear(K, &t);
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_expt_value(a, b, c) */
+
+mp_result mp_int_expt_value(klisp_State *K, mp_small a, mp_small b, mp_int c)
+{
+ mpz_t t;
+ mp_result res;
+ unsigned int v = abs(b);
+
+ CHECK(b >= 0 && c != NULL);
+
+ if((res = mp_int_init_value(K, &t, a)) != MP_OK)
+ return res;
+
+ (void) mp_int_set_value(K, c, 1);
+ while(v != 0) {
+ if(v & 1) {
+ if((res = mp_int_mul(K, c, &t, c)) != MP_OK)
+ goto CLEANUP;
+ }
+
+ v >>= 1;
+ if(v == 0) break;
+
+ if((res = mp_int_sqr(K, &t, &t)) != MP_OK)
+ goto CLEANUP;
+ }
+
+ CLEANUP:
+ mp_int_clear(K, &t);
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_expt_full(a, b, c) */
+
+mp_result mp_int_expt_full(klisp_State *K, mp_int a, mp_int b, mp_int c)
+{
+ mpz_t t;
+ mp_result res;
+ int ix, jx;
+
+ CHECK(a != NULL && b != NULL && c != NULL);
+
+ if ((res = mp_int_init_copy(K, &t, a)) != MP_OK)
+ return res;
+
+ (void) mp_int_set_value(K, c, 1);
+ for (ix = 0; ix < MP_USED(b); ++ix) {
+ mp_digit d = b->digits[ix];
+
+ for (jx = 0; jx < MP_DIGIT_BIT; ++jx) {
+ if (d & 1) {
+ if ((res = mp_int_mul(K, c, &t, c)) != MP_OK)
+ goto CLEANUP;
+ }
+
+ d >>= 1;
+ if (d == 0 && ix + 1 == MP_USED(b))
+ break;
+ if ((res = mp_int_sqr(K, &t, &t)) != MP_OK)
+ goto CLEANUP;
+ }
+ }
+
+ CLEANUP:
+ mp_int_clear(K, &t);
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_compare(a, b) */
+
+int mp_int_compare(mp_int a, mp_int b)
+{
+ mp_sign sa;
+
+ CHECK(a != NULL && b != NULL);
+
+ sa = MP_SIGN(a);
+ if(sa == MP_SIGN(b)) {
+ int cmp = s_ucmp(a, b);
+
+ /* If they're both zero or positive, the normal comparison
+ applies; if both negative, the sense is reversed. */
+ if(sa == MP_ZPOS)
+ return cmp;
+ else
+ return -cmp;
+
+ }
+ else {
+ if(sa == MP_ZPOS)
+ return 1;
+ else
+ return -1;
+ }
+}
+
+/* }}} */
+
+/* {{{ mp_int_compare_unsigned(a, b) */
+
+int mp_int_compare_unsigned(mp_int a, mp_int b)
+{
+ NRCHECK(a != NULL && b != NULL);
+
+ return s_ucmp(a, b);
+}
+
+/* }}} */
+
+/* {{{ mp_int_compare_zero(z) */
+
+int mp_int_compare_zero(mp_int z)
+{
+ NRCHECK(z != NULL);
+
+ if(MP_USED(z) == 1 && z->digits[0] == 0)
+ return 0;
+ else if(MP_SIGN(z) == MP_ZPOS)
+ return 1;
+ else
+ return -1;
+}
+
+/* }}} */
+
+/* {{{ mp_int_compare_value(z, value) */
+
+int mp_int_compare_value(mp_int z, mp_small value)
+{
+ mp_sign vsign = (value < 0) ? MP_NEG : MP_ZPOS;
+ int cmp;
+
+ CHECK(z != NULL);
+
+ if(vsign == MP_SIGN(z)) {
+ cmp = s_vcmp(z, value);
+
+ if(vsign == MP_ZPOS)
+ return cmp;
+ else
+ return -cmp;
+ }
+ else {
+ if(value < 0)
+ return 1;
+ else
+ return -1;
+ }
+}
+
+/* }}} */
+
+/* {{{ mp_int_exptmod(a, b, m, c) */
+
+mp_result mp_int_exptmod(klisp_State *K, mp_int a, mp_int b, mp_int m,
+ mp_int c)
+{
+ mp_result res;
+ mp_size um;
+ mpz_t temp[3];
+ mp_int s;
+ int last = 0;
+
+ CHECK(a != NULL && b != NULL && c != NULL && m != NULL);
+
+ /* Zero moduli and negative exponents are not considered. */
+ if(CMPZ(m) == 0)
+ return MP_UNDEF;
+ if(CMPZ(b) < 0)
+ return MP_RANGE;
+
+ um = MP_USED(m);
+ SETUP(mp_int_init_size(K, TEMP(0), 2 * um), last);
+ SETUP(mp_int_init_size(K, TEMP(1), 2 * um), last);
+
+ if(c == b || c == m) {
+ SETUP(mp_int_init_size(K, TEMP(2), 2 * um), last);
+ s = TEMP(2);
+ }
+ else {
+ s = c;
+ }
+
+ if((res = mp_int_mod(K, a, m, TEMP(0))) != MP_OK) goto CLEANUP;
+
+ if((res = s_brmu(K, TEMP(1), m)) != MP_OK) goto CLEANUP;
+
+ if((res = s_embar(K, TEMP(0), b, m, TEMP(1), s)) != MP_OK)
+ goto CLEANUP;
+
+ res = mp_int_copy(K, s, c);
+
+ CLEANUP:
+ while(--last >= 0)
+ mp_int_clear(K, TEMP(last));
+
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_exptmod_evalue(a, value, m, c) */
+
+mp_result mp_int_exptmod_evalue(klisp_State *K, mp_int a, mp_small value,
+ mp_int m, mp_int c)
+{
+ mpz_t vtmp;
+ mp_digit vbuf[MP_VALUE_DIGITS(value)];
+
+ s_fake(&vtmp, value, vbuf);
+
+ return mp_int_exptmod(K, a, &vtmp, m, c);
+}
+
+/* }}} */
+
+/* {{{ mp_int_exptmod_bvalue(v, b, m, c) */
+
+mp_result mp_int_exptmod_bvalue(klisp_State *K, mp_small value, mp_int b,
+ mp_int m, mp_int c)
+{
+ mpz_t vtmp;
+ mp_digit vbuf[MP_VALUE_DIGITS(value)];
+
+ s_fake(&vtmp, value, vbuf);
+
+ return mp_int_exptmod(K, &vtmp, b, m, c);
+}
+
+/* }}} */
+
+/* {{{ mp_int_exptmod_known(a, b, m, mu, c) */
+
+mp_result mp_int_exptmod_known(klisp_State *K, mp_int a, mp_int b, mp_int m,
+ mp_int mu, mp_int c)
+{
+ mp_result res;
+ mp_size um;
+ mpz_t temp[2];
+ mp_int s;
+ int last = 0;
+
+ CHECK(a && b && m && c);
+
+ /* Zero moduli and negative exponents are not considered. */
+ if(CMPZ(m) == 0)
+ return MP_UNDEF;
+ if(CMPZ(b) < 0)
+ return MP_RANGE;
+
+ um = MP_USED(m);
+ SETUP(mp_int_init_size(K, TEMP(0), 2 * um), last);
+
+ if(c == b || c == m) {
+ SETUP(mp_int_init_size(K, TEMP(1), 2 * um), last);
+ s = TEMP(1);
+ }
+ else {
+ s = c;
+ }
+
+ if((res = mp_int_mod(K, a, m, TEMP(0))) != MP_OK) goto CLEANUP;
+
+ if((res = s_embar(K, TEMP(0), b, m, mu, s)) != MP_OK)
+ goto CLEANUP;
+
+ res = mp_int_copy(K, s, c);
+
+ CLEANUP:
+ while(--last >= 0)
+ mp_int_clear(K, TEMP(last));
+
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_redux_const(m, c) */
+
+mp_result mp_int_redux_const(klisp_State *K, mp_int m, mp_int c)
+{
+ CHECK(m != NULL && c != NULL && m != c);
+
+ return s_brmu(K, c, m);
+}
+
+/* }}} */
+
+/* {{{ mp_int_invmod(a, m, c) */
+
+mp_result mp_int_invmod(klisp_State *K, mp_int a, mp_int m, mp_int c)
+{
+ mp_result res;
+ mp_sign sa;
+ int last = 0;
+ mpz_t temp[2];
+
+ CHECK(a != NULL && m != NULL && c != NULL);
+
+ if(CMPZ(a) == 0 || CMPZ(m) <= 0)
+ return MP_RANGE;
+
+ sa = MP_SIGN(a); /* need this for the result later */
+
+ for(last = 0; last < 2; ++last)
+ mp_int_init(TEMP(last));
+
+ if((res = mp_int_egcd(K, a, m, TEMP(0), TEMP(1), NULL)) != MP_OK)
+ goto CLEANUP;
+
+ if(mp_int_compare_value(TEMP(0), 1) != 0) {
+ res = MP_UNDEF;
+ goto CLEANUP;
+ }
+
+ /* It is first necessary to constrain the value to the proper range */
+ if((res = mp_int_mod(K, TEMP(1), m, TEMP(1))) != MP_OK)
+ goto CLEANUP;
+
+ /* Now, if 'a' was originally negative, the value we have is
+ actually the magnitude of the negative representative; to get the
+ positive value we have to subtract from the modulus. Otherwise,
+ the value is okay as it stands.
+ */
+ if(sa == MP_NEG)
+ res = mp_int_sub(K, m, TEMP(1), c);
+ else
+ res = mp_int_copy(K, TEMP(1), c);
+
+ CLEANUP:
+ while(--last >= 0)
+ mp_int_clear(K, TEMP(last));
+
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_gcd(a, b, c) */
+
+/* Binary GCD algorithm due to Josef Stein, 1961 */
+mp_result mp_int_gcd(klisp_State *K, mp_int a, mp_int b, mp_int c)
+{
+ int ca, cb, k = 0;
+ mpz_t u, v, t;
+ mp_result res;
+
+ CHECK(a != NULL && b != NULL && c != NULL);
+
+ ca = CMPZ(a);
+ cb = CMPZ(b);
+ if(ca == 0 && cb == 0)
+ return MP_UNDEF;
+ else if(ca == 0)
+ return mp_int_abs(K, b, c);
+ else if(cb == 0)
+ return mp_int_abs(K, a, c);
+
+ mp_int_init(&t);
+ if((res = mp_int_init_copy(K, &u, a)) != MP_OK)
+ goto U;
+ if((res = mp_int_init_copy(K, &v, b)) != MP_OK)
+ goto V;
+
+ MP_SIGN(&u) = MP_ZPOS; MP_SIGN(&v) = MP_ZPOS;
+
+ { /* Divide out common factors of 2 from u and v */
+ int div2_u = s_dp2k(&u), div2_v = s_dp2k(&v);
+
+ k = MIN(div2_u, div2_v);
+ s_qdiv(&u, (mp_size) k);
+ s_qdiv(&v, (mp_size) k);
+ }
+
+ if(mp_int_is_odd(&u)) {
+ if((res = mp_int_neg(K, &v, &t)) != MP_OK)
+ goto CLEANUP;
+ }
+ else {
+ if((res = mp_int_copy(K, &u, &t)) != MP_OK)
+ goto CLEANUP;
+ }
+
+ for(;;) {
+ s_qdiv(&t, s_dp2k(&t));
+
+ if(CMPZ(&t) > 0) {
+ if((res = mp_int_copy(K, &t, &u)) != MP_OK)
+ goto CLEANUP;
+ }
+ else {
+ if((res = mp_int_neg(K, &t, &v)) != MP_OK)
+ goto CLEANUP;
+ }
+
+ if((res = mp_int_sub(K, &u, &v, &t)) != MP_OK)
+ goto CLEANUP;
+
+ if(CMPZ(&t) == 0)
+ break;
+ }
+
+ if((res = mp_int_abs(K, &u, c)) != MP_OK)
+ goto CLEANUP;
+ if(!s_qmul(K, c, (mp_size) k))
+ res = MP_MEMORY;
+
+ CLEANUP:
+ mp_int_clear(K, &v);
+ V: mp_int_clear(K, &u);
+ U: mp_int_clear(K, &t);
+
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_egcd(a, b, c, x, y) */
+
+/* This is the binary GCD algorithm again, but this time we keep track
+ of the elementary matrix operations as we go, so we can get values
+ x and y satisfying c = ax + by.
+ */
+mp_result mp_int_egcd(klisp_State *K, mp_int a, mp_int b, mp_int c,
+ mp_int x, mp_int y)
+{
+ int k, last = 0, ca, cb;
+ mpz_t temp[8];
+ mp_result res;
+
+ CHECK(a != NULL && b != NULL && c != NULL &&
+ (x != NULL || y != NULL));
+
+ ca = CMPZ(a);
+ cb = CMPZ(b);
+ if(ca == 0 && cb == 0)
+ return MP_UNDEF;
+ else if(ca == 0) {
+ if((res = mp_int_abs(K, b, c)) != MP_OK) return res;
+ mp_int_zero(x); (void) mp_int_set_value(K, y, 1); return MP_OK;
+ }
+ else if(cb == 0) {
+ if((res = mp_int_abs(K, a, c)) != MP_OK) return res;
+ (void) mp_int_set_value(K, x, 1); mp_int_zero(y); return MP_OK;
+ }
+
+ /* Initialize temporaries:
+ A:0, B:1, C:2, D:3, u:4, v:5, ou:6, ov:7 */
+ for(last = 0; last < 4; ++last)
+ mp_int_init(TEMP(last));
+ TEMP(0)->digits[0] = 1;
+ TEMP(3)->digits[0] = 1;
+
+ SETUP(mp_int_init_copy(K, TEMP(4), a), last);
+ SETUP(mp_int_init_copy(K, TEMP(5), b), last);
+
+ /* We will work with absolute values here */
+ MP_SIGN(TEMP(4)) = MP_ZPOS;
+ MP_SIGN(TEMP(5)) = MP_ZPOS;
+
+ { /* Divide out common factors of 2 from u and v */
+ int div2_u = s_dp2k(TEMP(4)), div2_v = s_dp2k(TEMP(5));
+
+ k = MIN(div2_u, div2_v);
+ s_qdiv(TEMP(4), k);
+ s_qdiv(TEMP(5), k);
+ }
+
+ SETUP(mp_int_init_copy(K, TEMP(6), TEMP(4)), last);
+ SETUP(mp_int_init_copy(K, TEMP(7), TEMP(5)), last);
+
+ for(;;) {
+ while(mp_int_is_even(TEMP(4))) {
+ s_qdiv(TEMP(4), 1);
+
+ if(mp_int_is_odd(TEMP(0)) || mp_int_is_odd(TEMP(1))) {
+ if((res = mp_int_add(K, TEMP(0), TEMP(7), TEMP(0))) != MP_OK)
+ goto CLEANUP;
+ if((res = mp_int_sub(K, TEMP(1), TEMP(6), TEMP(1))) != MP_OK)
+ goto CLEANUP;
+ }
+
+ s_qdiv(TEMP(0), 1);
+ s_qdiv(TEMP(1), 1);
+ }
+
+ while(mp_int_is_even(TEMP(5))) {
+ s_qdiv(TEMP(5), 1);
+
+ if(mp_int_is_odd(TEMP(2)) || mp_int_is_odd(TEMP(3))) {
+ if((res = mp_int_add(K, TEMP(2), TEMP(7), TEMP(2))) != MP_OK)
+ goto CLEANUP;
+ if((res = mp_int_sub(K, TEMP(3), TEMP(6), TEMP(3))) != MP_OK)
+ goto CLEANUP;
+ }
+
+ s_qdiv(TEMP(2), 1);
+ s_qdiv(TEMP(3), 1);
+ }
+
+ if(mp_int_compare(TEMP(4), TEMP(5)) >= 0) {
+ if((res = mp_int_sub(K, TEMP(4), TEMP(5), TEMP(4))) != MP_OK) goto CLEANUP;
+ if((res = mp_int_sub(K, TEMP(0), TEMP(2), TEMP(0))) != MP_OK) goto CLEANUP;
+ if((res = mp_int_sub(K, TEMP(1), TEMP(3), TEMP(1))) != MP_OK) goto CLEANUP;
+ }
+ else {
+ if((res = mp_int_sub(K, TEMP(5), TEMP(4), TEMP(5))) != MP_OK) goto CLEANUP;
+ if((res = mp_int_sub(K, TEMP(2), TEMP(0), TEMP(2))) != MP_OK) goto CLEANUP;
+ if((res = mp_int_sub(K, TEMP(3), TEMP(1), TEMP(3))) != MP_OK) goto CLEANUP;
+ }
+
+ if(CMPZ(TEMP(4)) == 0) {
+ if(x && (res = mp_int_copy(K, TEMP(2), x)) != MP_OK) goto CLEANUP;
+ if(y && (res = mp_int_copy(K, TEMP(3), y)) != MP_OK) goto CLEANUP;
+ if(c) {
+ if(!s_qmul(K, TEMP(5), k)) {
+ res = MP_MEMORY;
+ goto CLEANUP;
+ }
+
+ res = mp_int_copy(K, TEMP(5), c);
+ }
+
+ break;
+ }
+ }
+
+ CLEANUP:
+ while(--last >= 0)
+ mp_int_clear(K, TEMP(last));
+
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_lcm(a, b, c) */
+
+mp_result mp_int_lcm(klisp_State *K, mp_int a, mp_int b, mp_int c)
+{
+ mpz_t lcm;
+ mp_result res;
+
+ CHECK(a != NULL && b != NULL && c != NULL);
+
+ /* Since a * b = gcd(a, b) * lcm(a, b), we can compute
+ lcm(a, b) = (a / gcd(a, b)) * b.
+
+ This formulation insures everything works even if the input
+ variables share space.
+ */
+ if((res = mp_int_init(&lcm)) != MP_OK)
+ return res;
+ if((res = mp_int_gcd(K, a, b, &lcm)) != MP_OK)
+ goto CLEANUP;
+ if((res = mp_int_div(K, a, &lcm, &lcm, NULL)) != MP_OK)
+ goto CLEANUP;
+ if((res = mp_int_mul(K, &lcm, b, &lcm)) != MP_OK)
+ goto CLEANUP;
+
+ res = mp_int_copy(K, &lcm, c);
+
+ CLEANUP:
+ mp_int_clear(K, &lcm);
+
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_divisible_value(a, v) */
+
+int mp_int_divisible_value(klisp_State *K, mp_int a, mp_small v)
+{
+ mp_small rem = 0;
+
+ if(mp_int_div_value(K, a, v, NULL, &rem) != MP_OK)
+ return 0;
+
+ return rem == 0;
+}
+
+/* }}} */
+
+/* {{{ mp_int_is_pow2(z) */
+
+int mp_int_is_pow2(mp_int z)
+{
+ CHECK(z != NULL);
+
+ return s_isp2(z);
+}
+
+/* }}} */
+
+/* {{{ mp_int_root(a, b, c) */
+
+/* Implementation of Newton's root finding method, based loosely on a
+ patch contributed by Hal Finkel <half@halssoftware.com>
+ modified by M. J. Fromberger.
+ */
+mp_result mp_int_root(klisp_State *K, mp_int a, mp_small b, mp_int c)
+{
+ mp_result res = MP_OK;
+ mpz_t temp[5];
+ int last = 0;
+ int flips = 0;
+
+ CHECK(a != NULL && c != NULL && b > 0);
+
+ if(b == 1) {
+ return mp_int_copy(K, a, c);
+ }
+ if(MP_SIGN(a) == MP_NEG) {
+ if(b % 2 == 0)
+ return MP_UNDEF; /* root does not exist for negative a with even b */
+ else
+ flips = 1;
+ }
+
+ SETUP(mp_int_init_copy(K, TEMP(last), a), last);
+ SETUP(mp_int_init_copy(K, TEMP(last), a), last);
+ SETUP(mp_int_init(TEMP(last)), last);
+ SETUP(mp_int_init(TEMP(last)), last);
+ SETUP(mp_int_init(TEMP(last)), last);
+
+ (void) mp_int_abs(K, TEMP(0), TEMP(0));
+ (void) mp_int_abs(K, TEMP(1), TEMP(1));
+
+ for(;;) {
+ if((res = mp_int_expt(K, TEMP(1), b, TEMP(2))) != MP_OK)
+ goto CLEANUP;
+
+ if(mp_int_compare_unsigned(TEMP(2), TEMP(0)) <= 0)
+ break;
+
+ if((res = mp_int_sub(K, TEMP(2), TEMP(0), TEMP(2))) != MP_OK)
+ goto CLEANUP;
+ if((res = mp_int_expt(K, TEMP(1), b - 1, TEMP(3))) != MP_OK)
+ goto CLEANUP;
+ if((res = mp_int_mul_value(K, TEMP(3), b, TEMP(3))) != MP_OK)
+ goto CLEANUP;
+ if((res = mp_int_div(K, TEMP(2), TEMP(3), TEMP(4), NULL)) != MP_OK)
+ goto CLEANUP;
+ if((res = mp_int_sub(K, TEMP(1), TEMP(4), TEMP(4))) != MP_OK)
+ goto CLEANUP;
+
+ if(mp_int_compare_unsigned(TEMP(1), TEMP(4)) == 0) {
+ if((res = mp_int_sub_value(K, TEMP(4), 1, TEMP(4))) != MP_OK)
+ goto CLEANUP;
+ }
+ if((res = mp_int_copy(K, TEMP(4), TEMP(1))) != MP_OK)
+ goto CLEANUP;
+ }
+
+ if((res = mp_int_copy(K, TEMP(1), c)) != MP_OK)
+ goto CLEANUP;
+
+ /* If the original value of a was negative, flip the output sign. */
+ if(flips)
+ (void) mp_int_neg(K, c, c); /* cannot fail */
+
+ CLEANUP:
+ while(--last >= 0)
+ mp_int_clear(K, TEMP(last));
+
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_to_int(z, out) */
+
+mp_result mp_int_to_int(mp_int z, mp_small *out)
+{
+ mp_usmall uv = 0;
+ mp_size uz;
+ mp_digit *dz;
+ mp_sign sz;
+
+ CHECK(z != NULL);
+
+ /* Make sure the value is representable as an int */
+ sz = MP_SIGN(z);
+ if((sz == MP_ZPOS && mp_int_compare_value(z, MP_SMALL_MAX) > 0) ||
+ mp_int_compare_value(z, MP_SMALL_MIN) < 0)
+ return MP_RANGE;
+
+ uz = MP_USED(z);
+ dz = MP_DIGITS(z) + uz - 1;
+
+ while(uz > 0) {
+ uv <<= MP_DIGIT_BIT/2;
+ uv = (uv << (MP_DIGIT_BIT/2)) | *dz--;
+ --uz;
+ }
+
+ if(out)
+ *out = (sz == MP_NEG) ? -(mp_small)uv : (mp_small)uv;
+
+ return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_to_uint(z, *out) */
+
+mp_result mp_int_to_uint(mp_int z, mp_usmall *out)
+{
+ mp_usmall uv = 0;
+ mp_size uz;
+ mp_digit *dz;
+ mp_sign sz;
+
+ CHECK(z != NULL);
+
+ /* Make sure the value is representable as an int */
+ sz = MP_SIGN(z);
+ if(!(sz == MP_ZPOS && mp_int_compare_value(z, UINT_MAX) <= 0))
+ return MP_RANGE;
+
+ uz = MP_USED(z);
+ dz = MP_DIGITS(z) + uz - 1;
+
+ while(uz > 0) {
+ uv <<= MP_DIGIT_BIT/2;
+ uv = (uv << (MP_DIGIT_BIT/2)) | *dz--;
+ --uz;
+ }
+
+ if(out)
+ *out = uv;
+
+ return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_to_string(z, radix, str, limit) */
+
+mp_result mp_int_to_string(klisp_State *K, mp_int z, mp_size radix,
+ char *str, int limit)
+{
+ mp_result res;
+ int cmp = 0;
+
+ CHECK(z != NULL && str != NULL && limit >= 2);
+
+ if(radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
+ return MP_RANGE;
+
+ if(CMPZ(z) == 0) {
+ *str++ = s_val2ch(0, 1);
+ }
+ else {
+ mpz_t tmp;
+ char *h, *t;
+
+ if((res = mp_int_init_copy(K, &tmp, z)) != MP_OK)
+ return res;
+
+ if(MP_SIGN(z) == MP_NEG) {
+ *str++ = '-';
+ --limit;
+ }
+ h = str;
+
+ /* Generate digits in reverse order until finished or limit reached */
+ for(/* */; limit > 0; --limit) {
+ mp_digit d;
+
+ if((cmp = CMPZ(&tmp)) == 0)
+ break;
+
+ d = s_ddiv(&tmp, (mp_digit)radix);
+ *str++ = s_val2ch(d, 1);
+ }
+ t = str - 1;
+
+ /* Put digits back in correct output order */
+ while(h < t) {
+ char tc = *h;
+ *h++ = *t;
+ *t-- = tc;
+ }
+
+ mp_int_clear(K, &tmp);
+ }
+
+ *str = '\0';
+ if(cmp == 0)
+ return MP_OK;
+ else
+ return MP_TRUNC;
+}
+
+/* }}} */
+
+/* {{{ mp_int_string_len(z, radix) */
+
+mp_result mp_int_string_len(mp_int z, mp_size radix)
+{
+ int len;
+
+ CHECK(z != NULL);
+
+ if(radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
+ return MP_RANGE;
+
+ len = s_outlen(z, radix) + 1; /* for terminator */
+
+ /* Allow for sign marker on negatives */
+ if(MP_SIGN(z) == MP_NEG)
+ len += 1;
+
+ return len;
+}
+
+/* }}} */
+
+/* {{{ mp_int_read_string(z, radix, *str) */
+
+/* Read zero-terminated string into z */
+mp_result mp_int_read_string(klisp_State *K, mp_int z, mp_size radix,
+ const char *str)
+{
+ return mp_int_read_cstring(K, z, radix, str, NULL);
+}
+
+/* }}} */
+
+/* {{{ mp_int_read_cstring(z, radix, *str, **end) */
+
+mp_result mp_int_read_cstring(klisp_State *K, mp_int z, mp_size radix,
+ const char *str, char **end)
+{
+ int ch;
+
+ CHECK(z != NULL && str != NULL);
+
+ if(radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
+ return MP_RANGE;
+
+ /* Skip leading whitespace */
+ while(isspace((int)*str))
+ ++str;
+
+ /* Handle leading sign tag (+/-, positive default) */
+ switch(*str) {
+ case '-':
+ MP_SIGN(z) = MP_NEG;
+ ++str;
+ break;
+ case '+':
+ ++str; /* fallthrough */
+ default:
+ MP_SIGN(z) = MP_ZPOS;
+ break;
+ }
+
+ /* Skip leading zeroes */
+ while((ch = s_ch2val(*str, radix)) == 0)
+ ++str;
+
+ /* Make sure there is enough space for the value */
+ if(!s_pad(K, z, s_inlen(strlen(str), radix)))
+ return MP_MEMORY;
+
+ MP_USED(z) = 1; z->digits[0] = 0;
+
+ while(*str != '\0' && ((ch = s_ch2val(*str, radix)) >= 0)) {
+ s_dmul(z, (mp_digit)radix);
+ s_dadd(z, (mp_digit)ch);
+ ++str;
+ }
+
+ CLAMP(z);
+
+ /* Override sign for zero, even if negative specified. */
+ if(CMPZ(z) == 0)
+ MP_SIGN(z) = MP_ZPOS;
+
+ if(end != NULL)
+ *end = (char *)str;
+
+ /* Return a truncation error if the string has unprocessed
+ characters remaining, so the caller can tell if the whole string
+ was done */
+ if(*str != '\0')
+ return MP_TRUNC;
+ else
+ return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_count_bits(z) */
+
+mp_result mp_int_count_bits(mp_int z)
+{
+ mp_size nbits = 0, uz;
+ mp_digit d;
+
+ CHECK(z != NULL);
+
+ uz = MP_USED(z);
+ if(uz == 1 && z->digits[0] == 0)
+ return 1;
+
+ --uz;
+ nbits = uz * MP_DIGIT_BIT;
+ d = z->digits[uz];
+
+ while(d != 0) {
+ d >>= 1;
+ ++nbits;
+ }
+
+ return nbits;
+}
+
+/* }}} */
+
+/* {{{ mp_int_to_binary(z, buf, limit) */
+
+mp_result mp_int_to_binary(klisp_State *K, mp_int z, unsigned char *buf,
+ int limit)
+{
+ static const int PAD_FOR_2C = 1;
+
+ mp_result res;
+ int limpos = limit;
+
+ CHECK(z != NULL && buf != NULL);
+
+ res = s_tobin(z, buf, &limpos, PAD_FOR_2C);
+
+ if(MP_SIGN(z) == MP_NEG)
+ s_2comp(buf, limpos);
+
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_int_read_binary(z, buf, len) */
+
+mp_result mp_int_read_binary(klisp_State *K, mp_int z, unsigned char *buf,
+ int len)
+{
+ mp_size need, i;
+ unsigned char *tmp;
+ mp_digit *dz;
+
+ CHECK(z != NULL && buf != NULL && len > 0);
+
+ /* Figure out how many digits are needed to represent this value */
+ need = ((len * CHAR_BIT) + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT;
+ if(!s_pad(K, z, need))
+ return MP_MEMORY;
+
+ mp_int_zero(z);
+
+ /* If the high-order bit is set, take the 2's complement before
+ reading the value (it will be restored afterward) */
+ if(buf[0] >> (CHAR_BIT - 1)) {
+ MP_SIGN(z) = MP_NEG;
+ s_2comp(buf, len);
+ }
+
+ dz = MP_DIGITS(z);
+ for(tmp = buf, i = len; i > 0; --i, ++tmp) {
+ s_qmul(K, z, (mp_size) CHAR_BIT);
+ *dz |= *tmp;
+ }
+
+ /* Restore 2's complement if we took it before */
+ if(MP_SIGN(z) == MP_NEG)
+ s_2comp(buf, len);
+
+ return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_binary_len(z) */
+
+mp_result mp_int_binary_len(mp_int z)
+{
+ mp_result res = mp_int_count_bits(z);
+ int bytes = mp_int_unsigned_len(z);
+
+ if(res <= 0)
+ return res;
+
+ bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT;
+
+ /* If the highest-order bit falls exactly on a byte boundary, we
+ need to pad with an extra byte so that the sign will be read
+ correctly when reading it back in. */
+ if(bytes * CHAR_BIT == res)
+ ++bytes;
+
+ return bytes;
+}
+
+/* }}} */
+
+/* {{{ mp_int_to_unsigned(z, buf, limit) */
+
+mp_result mp_int_to_unsigned(klisp_State *K, mp_int z, unsigned char *buf,
+ int limit)
+{
+ static const int NO_PADDING = 0;
+
+ CHECK(z != NULL && buf != NULL);
+
+ return s_tobin(z, buf, &limit, NO_PADDING);
+}
+
+/* }}} */
+
+/* {{{ mp_int_read_unsigned(z, buf, len) */
+
+mp_result mp_int_read_unsigned(klisp_State *K, mp_int z, unsigned char *buf, int len)
+{
+ mp_size need, i;
+ unsigned char *tmp;
+ mp_digit *dz;
+
+ CHECK(z != NULL && buf != NULL && len > 0);
+
+ /* Figure out how many digits are needed to represent this value */
+ need = ((len * CHAR_BIT) + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT;
+ if(!s_pad(K, z, need))
+ return MP_MEMORY;
+
+ mp_int_zero(z);
+
+ dz = MP_DIGITS(z);
+ for(tmp = buf, i = len; i > 0; --i, ++tmp) {
+ (void) s_qmul(K, z, CHAR_BIT);
+ *dz |= *tmp;
+ }
+
+ return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_int_unsigned_len(z) */
+
+mp_result mp_int_unsigned_len(mp_int z)
+{
+ mp_result res = mp_int_count_bits(z);
+ int bytes;
+
+ if(res <= 0)
+ return res;
+
+ bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT;
+
+ return bytes;
+}
+
+/* }}} */
+
+/* {{{ mp_error_string(res) */
+
+const char *mp_error_string(mp_result res)
+{
+ int ix;
+ if(res > 0)
+ return s_unknown_err;
+
+ res = -res;
+ for(ix = 0; ix < res && s_error_msg[ix] != NULL; ++ix)
+ ;
+
+ if(s_error_msg[ix] != NULL)
+ return s_error_msg[ix];
+ else
+ return s_unknown_err;
+}
+
+/* }}} */
+
+/*------------------------------------------------------------------------*/
+/* Private functions for internal use. These make assumptions. */
+
+/* {{{ s_alloc(num) */
+
+STATIC mp_digit *s_alloc(klisp_State *K, mp_size num)
+{
+ mp_digit *out = klispM_malloc(K, num * sizeof(mp_digit));
+
+ assert(out != NULL); /* for debugging */
+#if DEBUG > 1
+ {
+ mp_digit v = (mp_digit) 0xdeadbeef;
+ int ix;
+
+ for(ix = 0; ix < num; ++ix)
+ out[ix] = v;
+ }
+#endif
+
+ return out;
+}
+
+/* }}} */
+
+/* {{{ s_realloc(old, osize, nsize) */
+
+STATIC mp_digit *s_realloc(klisp_State *K, mp_digit *old, mp_size osize,
+ mp_size nsize)
+{
+#if DEBUG > 1
+ mp_digit *new = s_alloc(K, nsize);
+ int ix;
+
+ for(ix = 0; ix < nsize; ++ix)
+ new[ix] = (mp_digit) 0xdeadbeef;
+
+ memcpy(new, old, osize * sizeof(mp_digit));
+#else
+ mp_digit *new = klispM_realloc_(K, old, osize * sizeof(mp_digit),
+ nsize * sizeof(mp_digit));
+
+ assert(new != NULL); /* for debugging */
+#endif
+ return new;
+}
+
+/* }}} */
+
+/* {{{ s_free(ptr) */
+
+STATIC void s_free(klisp_State *K, void *ptr, mp_size size)
+{
+ klispM_freemem(K, ptr, size * sizeof(mp_digit));
+}
+
+/* }}} */
+
+/* {{{ s_pad(z, min) */
+
+STATIC int s_pad(klisp_State *K, mp_int z, mp_size min)
+{
+ if(MP_ALLOC(z) < min) {
+ mp_size nsize = ROUND_PREC(min);
+ mp_digit *tmp;
+
+ if((void *)z->digits == (void *)&(z->single)) {
+ if((tmp = s_alloc(K, nsize)) == NULL)
+ return 0;
+
+ COPY(MP_DIGITS(z), tmp, MP_USED(z));
+ } else if((tmp = s_realloc(K, MP_DIGITS(z), MP_ALLOC(z), nsize)) == NULL)
+ return 0;
+
+ MP_DIGITS(z) = tmp;
+ MP_ALLOC(z) = nsize;
+ }
+
+ return 1;
+}
+
+/* }}} */
+
+/* {{{ s_fake(z, value, vbuf) */
+
+STATIC void s_fake(mp_int z, mp_small value, mp_digit vbuf[])
+{
+ mp_size uv = (mp_size) s_vpack(value, vbuf);
+
+ z->used = uv;
+ z->alloc = MP_VALUE_DIGITS(value);
+ z->sign = (value < 0) ? MP_NEG : MP_ZPOS;
+ z->digits = vbuf;
+}
+
+/* }}} */
+
+/* {{{ s_cdig(da, db, len) */
+
+STATIC int s_cdig(mp_digit *da, mp_digit *db, mp_size len)
+{
+ mp_digit *dat = da + len - 1, *dbt = db + len - 1;
+
+ for(/* */; len != 0; --len, --dat, --dbt) {
+ if(*dat > *dbt)
+ return 1;
+ else if(*dat < *dbt)
+ return -1;
+ }
+
+ return 0;
+}
+
+/* }}} */
+
+/* {{{ s_vpack(v, t[]) */
+
+STATIC int s_vpack(mp_small v, mp_digit t[])
+{
+ mp_usmall uv = (mp_usmall) ((v < 0) ? -v : v);
+ int ndig = 0;
+
+ if(uv == 0)
+ t[ndig++] = 0;
+ else {
+ while(uv != 0) {
+ t[ndig++] = (mp_digit) uv;
+ uv >>= MP_DIGIT_BIT/2;
+ uv >>= MP_DIGIT_BIT/2;
+ }
+ }
+
+ return ndig;
+}
+
+/* }}} */
+
+/* {{{ s_ucmp(a, b) */
+
+STATIC int s_ucmp(mp_int a, mp_int b)
+{
+ mp_size ua = MP_USED(a), ub = MP_USED(b);
+
+ if(ua > ub)
+ return 1;
+ else if(ub > ua)
+ return -1;
+ else
+ return s_cdig(MP_DIGITS(a), MP_DIGITS(b), ua);
+}
+
+/* }}} */
+
+/* {{{ s_vcmp(a, v) */
+
+STATIC int s_vcmp(mp_int a, mp_small v)
+{
+ mp_digit vdig[MP_VALUE_DIGITS(v)];
+ int ndig = 0;
+ mp_size ua = MP_USED(a);
+
+ ndig = s_vpack(v, vdig);
+
+ if(ua > ndig)
+ return 1;
+ else if(ua < ndig)
+ return -1;
+ else
+ return s_cdig(MP_DIGITS(a), vdig, ndig);
+}
+
+/* }}} */
+
+/* {{{ s_uadd(da, db, dc, size_a, size_b) */
+
+STATIC mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc,
+ mp_size size_a, mp_size size_b)
+{
+ mp_size pos;
+ mp_word w = 0;
+
+ /* Insure that da is the longer of the two to simplify later code */
+ if(size_b > size_a) {
+ SWAP(mp_digit *, da, db);
+ SWAP(mp_size, size_a, size_b);
+ }
+
+ /* Add corresponding digits until the shorter number runs out */
+ for(pos = 0; pos < size_b; ++pos, ++da, ++db, ++dc) {
+ w = w + (mp_word) *da + (mp_word) *db;
+ *dc = LOWER_HALF(w);
+ w = UPPER_HALF(w);
+ }
+
+ /* Propagate carries as far as necessary */
+ for(/* */; pos < size_a; ++pos, ++da, ++dc) {
+ w = w + *da;
+
+ *dc = LOWER_HALF(w);
+ w = UPPER_HALF(w);
+ }
+
+ /* Return carry out */
+ return (mp_digit)w;
+}
+
+/* }}} */
+
+/* {{{ s_usub(da, db, dc, size_a, size_b) */
+
+STATIC void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc,
+ mp_size size_a, mp_size size_b)
+{
+ mp_size pos;
+ mp_word w = 0;
+
+ /* We assume that |a| >= |b| so this should definitely hold */
+ assert(size_a >= size_b);
+
+ /* Subtract corresponding digits and propagate borrow */
+ for(pos = 0; pos < size_b; ++pos, ++da, ++db, ++dc) {
+ w = ((mp_word)MP_DIGIT_MAX + 1 + /* MP_RADIX */
+ (mp_word)*da) - w - (mp_word)*db;
+
+ *dc = LOWER_HALF(w);
+ w = (UPPER_HALF(w) == 0);
+ }
+
+ /* Finish the subtraction for remaining upper digits of da */
+ for(/* */; pos < size_a; ++pos, ++da, ++dc) {
+ w = ((mp_word)MP_DIGIT_MAX + 1 + /* MP_RADIX */
+ (mp_word)*da) - w;
+
+ *dc = LOWER_HALF(w);
+ w = (UPPER_HALF(w) == 0);
+ }
+
+ /* If there is a borrow out at the end, it violates the precondition */
+ assert(w == 0);
+}
+
+/* }}} */
+
+/* {{{ s_kmul(da, db, dc, size_a, size_b) */
+
+STATIC int s_kmul(klisp_State *K, mp_digit *da, mp_digit *db,
+ mp_digit *dc, mp_size size_a, mp_size size_b)
+{
+ mp_size bot_size;
+
+ /* Make sure b is the smaller of the two input values */
+ if(size_b > size_a) {
+ SWAP(mp_digit *, da, db);
+ SWAP(mp_size, size_a, size_b);
+ }
+
+ /* Insure that the bottom is the larger half in an odd-length split;
+ the code below relies on this being true.
+ */
+ bot_size = (size_a + 1) / 2;
+
+ /* If the values are big enough to bother with recursion, use the
+ Karatsuba algorithm to compute the product; otherwise use the
+ normal multiplication algorithm
+ */
+ if(multiply_threshold &&
+ size_a >= multiply_threshold &&
+ size_b > bot_size) {
+
+ mp_digit *t1, *t2, *t3, carry;
+
+ mp_digit *a_top = da + bot_size;
+ mp_digit *b_top = db + bot_size;
+
+ mp_size at_size = size_a - bot_size;
+ mp_size bt_size = size_b - bot_size;
+ mp_size buf_size = 2 * bot_size;
+
+ /* Do a single allocation for all three temporary buffers needed;
+ each buffer must be big enough to hold the product of two
+ bottom halves, and one buffer needs space for the completed
+ product; twice the space is plenty.
+ */
+ if((t1 = s_alloc(K, 4 * buf_size)) == NULL) return 0;
+ t2 = t1 + buf_size;
+ t3 = t2 + buf_size;
+ ZERO(t1, 4 * buf_size);
+
+ /* t1 and t2 are initially used as temporaries to compute the inner product
+ (a1 + a0)(b1 + b0) = a1b1 + a1b0 + a0b1 + a0b0
+ */
+ carry = s_uadd(da, a_top, t1, bot_size, at_size); /* t1 = a1 + a0 */
+ t1[bot_size] = carry;
+
+ carry = s_uadd(db, b_top, t2, bot_size, bt_size); /* t2 = b1 + b0 */
+ t2[bot_size] = carry;
+
+ /* t3 = t1 * t2 */
+ (void) s_kmul(K, t1, t2, t3, bot_size + 1, bot_size + 1);
+
+ /* Now we'll get t1 = a0b0 and t2 = a1b1, and subtract them out so that
+ we're left with only the pieces we want: t3 = a1b0 + a0b1
+ */
+ ZERO(t1, buf_size);
+ ZERO(t2, buf_size);
+ (void) s_kmul(K, da, db, t1, bot_size, bot_size); /* t1 = a0 * b0 */
+ (void) s_kmul(K, a_top, b_top, t2, at_size, bt_size); /* t2 = a1 * b1 */
+
+ /* Subtract out t1 and t2 to get the inner product */
+ s_usub(t3, t1, t3, buf_size + 2, buf_size);
+ s_usub(t3, t2, t3, buf_size + 2, buf_size);
+
+ /* Assemble the output value */
+ COPY(t1, dc, buf_size);
+ carry = s_uadd(t3, dc + bot_size, dc + bot_size,
+ buf_size + 1, buf_size);
+ assert(carry == 0);
+
+ carry = s_uadd(t2, dc + 2*bot_size, dc + 2*bot_size,
+ buf_size, buf_size);
+ assert(carry == 0);
+
+ /* note t2 and t3 are just internal pointers to t1 */
+ s_free(K, t1, 4 * buf_size);
+ }
+ else {
+ s_umul(da, db, dc, size_a, size_b);
+ }
+
+ return 1;
+}
+
+/* }}} */
+
+/* {{{ s_umul(da, db, dc, size_a, size_b) */
+
+STATIC void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc,
+ mp_size size_a, mp_size size_b)
+{
+ mp_size a, b;
+ mp_word w;
+
+ for(a = 0; a < size_a; ++a, ++dc, ++da) {
+ mp_digit *dct = dc;
+ mp_digit *dbt = db;
+
+ if(*da == 0)
+ continue;
+
+ w = 0;
+ for(b = 0; b < size_b; ++b, ++dbt, ++dct) {
+ w = (mp_word)*da * (mp_word)*dbt + w + (mp_word)*dct;
+
+ *dct = LOWER_HALF(w);
+ w = UPPER_HALF(w);
+ }
+
+ *dct = (mp_digit)w;
+ }
+}
+
+/* }}} */
+
+/* {{{ s_ksqr(da, dc, size_a) */
+
+STATIC int s_ksqr(klisp_State *K, mp_digit *da, mp_digit *dc,
+ mp_size size_a)
+{
+ if(multiply_threshold && size_a > multiply_threshold) {
+ mp_size bot_size = (size_a + 1) / 2;
+ mp_digit *a_top = da + bot_size;
+ mp_digit *t1, *t2, *t3, carry;
+ mp_size at_size = size_a - bot_size;
+ mp_size buf_size = 2 * bot_size;
+
+ if((t1 = s_alloc(K, 4 * buf_size)) == NULL) return 0;
+ t2 = t1 + buf_size;
+ t3 = t2 + buf_size;
+ ZERO(t1, 4 * buf_size);
+
+ (void) s_ksqr(K, da, t1, bot_size); /* t1 = a0 ^ 2 */
+ (void) s_ksqr(K, a_top, t2, at_size); /* t2 = a1 ^ 2 */
+
+ (void) s_kmul(K, da, a_top, t3, bot_size, at_size); /* t3 = a0 * a1 */
+
+ /* Quick multiply t3 by 2, shifting left (can't overflow) */
+ {
+ int i, top = bot_size + at_size;
+ mp_word w, save = 0;
+
+ for(i = 0; i < top; ++i) {
+ w = t3[i];
+ w = (w << 1) | save;
+ t3[i] = LOWER_HALF(w);
+ save = UPPER_HALF(w);
+ }
+ t3[i] = LOWER_HALF(save);
+ }
+
+ /* Assemble the output value */
+ COPY(t1, dc, 2 * bot_size);
+ carry = s_uadd(t3, dc + bot_size, dc + bot_size,
+ buf_size + 1, buf_size);
+ assert(carry == 0);
+
+ carry = s_uadd(t2, dc + 2*bot_size, dc + 2*bot_size,
+ buf_size, buf_size);
+ assert(carry == 0);
+
+ /* note that t2 and t2 are internal pointers only */
+ s_free(K, t1, 4 * buf_size);
+ }
+ else {
+ s_usqr(da, dc, size_a);
+ }
+
+ return 1;
+}
+
+/* }}} */
+
+/* {{{ s_usqr(da, dc, size_a) */
+
+STATIC void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a)
+{
+ mp_size i, j;
+ mp_word w;
+
+ for(i = 0; i < size_a; ++i, dc += 2, ++da) {
+ mp_digit *dct = dc, *dat = da;
+
+ if(*da == 0)
+ continue;
+
+ /* Take care of the first digit, no rollover */
+ w = (mp_word)*dat * (mp_word)*dat + (mp_word)*dct;
+ *dct = LOWER_HALF(w);
+ w = UPPER_HALF(w);
+ ++dat; ++dct;
+
+ for(j = i + 1; j < size_a; ++j, ++dat, ++dct) {
+ mp_word t = (mp_word)*da * (mp_word)*dat;
+ mp_word u = w + (mp_word)*dct, ov = 0;
+
+ /* Check if doubling t will overflow a word */
+ if(HIGH_BIT_SET(t))
+ ov = 1;
+
+ w = t + t;
+
+ /* Check if adding u to w will overflow a word */
+ if(ADD_WILL_OVERFLOW(w, u))
+ ov = 1;
+
+ w += u;
+
+ *dct = LOWER_HALF(w);
+ w = UPPER_HALF(w);
+ if(ov) {
+ w += MP_DIGIT_MAX; /* MP_RADIX */
+ ++w;
+ }
+ }
+
+ w = w + *dct;
+ *dct = (mp_digit)w;
+ while((w = UPPER_HALF(w)) != 0) {
+ ++dct; w = w + *dct;
+ *dct = LOWER_HALF(w);
+ }
+
+ assert(w == 0);
+ }
+}
+
+/* }}} */
+
+/* {{{ s_dadd(a, b) */
+
+STATIC void s_dadd(mp_int a, mp_digit b)
+{
+ mp_word w = 0;
+ mp_digit *da = MP_DIGITS(a);
+ mp_size ua = MP_USED(a);
+
+ w = (mp_word)*da + b;
+ *da++ = LOWER_HALF(w);
+ w = UPPER_HALF(w);
+
+ for(ua -= 1; ua > 0; --ua, ++da) {
+ w = (mp_word)*da + w;
+
+ *da = LOWER_HALF(w);
+ w = UPPER_HALF(w);
+ }
+
+ if(w) {
+ *da = (mp_digit)w;
+ MP_USED(a) += 1;
+ }
+}
+
+/* }}} */
+
+/* {{{ s_dmul(a, b) */
+
+STATIC void s_dmul(mp_int a, mp_digit b)
+{
+ mp_word w = 0;
+ mp_digit *da = MP_DIGITS(a);
+ mp_size ua = MP_USED(a);
+
+ while(ua > 0) {
+ w = (mp_word)*da * b + w;
+ *da++ = LOWER_HALF(w);
+ w = UPPER_HALF(w);
+ --ua;
+ }
+
+ if(w) {
+ *da = (mp_digit)w;
+ MP_USED(a) += 1;
+ }
+}
+
+/* }}} */
+
+/* {{{ s_dbmul(da, b, dc, size_a) */
+
+STATIC void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, mp_size size_a)
+{
+ mp_word w = 0;
+
+ while(size_a > 0) {
+ w = (mp_word)*da++ * (mp_word)b + w;
+
+ *dc++ = LOWER_HALF(w);
+ w = UPPER_HALF(w);
+ --size_a;
+ }
+
+ if(w)
+ *dc = LOWER_HALF(w);
+}
+
+/* }}} */
+
+/* {{{ s_ddiv(da, d, dc, size_a) */
+
+STATIC mp_digit s_ddiv(mp_int a, mp_digit b)
+{
+ mp_word w = 0, qdigit;
+ mp_size ua = MP_USED(a);
+ mp_digit *da = MP_DIGITS(a) + ua - 1;
+
+ for(/* */; ua > 0; --ua, --da) {
+ w = (w << MP_DIGIT_BIT) | *da;
+
+ if(w >= b) {
+ qdigit = w / b;
+ w = w % b;
+ }
+ else {
+ qdigit = 0;
+ }
+
+ *da = (mp_digit)qdigit;
+ }
+
+ CLAMP(a);
+ return (mp_digit)w;
+}
+
+/* }}} */
+
+/* {{{ s_qdiv(z, p2) */
+
+STATIC void s_qdiv(mp_int z, mp_size p2)
+{
+ mp_size ndig = p2 / MP_DIGIT_BIT, nbits = p2 % MP_DIGIT_BIT;
+ mp_size uz = MP_USED(z);
+
+ if(ndig) {
+ mp_size mark;
+ mp_digit *to, *from;
+
+ if(ndig >= uz) {
+ mp_int_zero(z);
+ return;
+ }
+
+ to = MP_DIGITS(z); from = to + ndig;
+
+ for(mark = ndig; mark < uz; ++mark)
+ *to++ = *from++;
+
+ MP_USED(z) = uz - ndig;
+ }
+
+ if(nbits) {
+ mp_digit d = 0, *dz, save;
+ mp_size up = MP_DIGIT_BIT - nbits;
+
+ uz = MP_USED(z);
+ dz = MP_DIGITS(z) + uz - 1;
+
+ for(/* */; uz > 0; --uz, --dz) {
+ save = *dz;
+
+ *dz = (*dz >> nbits) | (d << up);
+ d = save;
+ }
+
+ CLAMP(z);
+ }
+
+ if(MP_USED(z) == 1 && z->digits[0] == 0)
+ MP_SIGN(z) = MP_ZPOS;
+}
+
+/* }}} */
+
+/* {{{ s_qmod(z, p2) */
+
+STATIC void s_qmod(mp_int z, mp_size p2)
+{
+ mp_size start = p2 / MP_DIGIT_BIT + 1, rest = p2 % MP_DIGIT_BIT;
+ mp_size uz = MP_USED(z);
+ mp_digit mask = (1 << rest) - 1;
+
+ if(start <= uz) {
+ MP_USED(z) = start;
+ z->digits[start - 1] &= mask;
+ CLAMP(z);
+ }
+}
+
+/* }}} */
+
+/* {{{ s_qmul(z, p2) */
+
+STATIC int s_qmul(klisp_State *K, mp_int z, mp_size p2)
+{
+ mp_size uz, need, rest, extra, i;
+ mp_digit *from, *to, d;
+
+ if(p2 == 0)
+ return 1;
+
+ uz = MP_USED(z);
+ need = p2 / MP_DIGIT_BIT; rest = p2 % MP_DIGIT_BIT;
+
+ /* Figure out if we need an extra digit at the top end; this occurs
+ if the topmost `rest' bits of the high-order digit of z are not
+ zero, meaning they will be shifted off the end if not preserved */
+ extra = 0;
+ if(rest != 0) {
+ mp_digit *dz = MP_DIGITS(z) + uz - 1;
+
+ if((*dz >> (MP_DIGIT_BIT - rest)) != 0)
+ extra = 1;
+ }
+
+ if(!s_pad(K, z, uz + need + extra))
+ return 0;
+
+ /* If we need to shift by whole digits, do that in one pass, then
+ to back and shift by partial digits.
+ */
+ if(need > 0) {
+ from = MP_DIGITS(z) + uz - 1;
+ to = from + need;
+
+ for(i = 0; i < uz; ++i)
+ *to-- = *from--;
+
+ ZERO(MP_DIGITS(z), need);
+ uz += need;
+ }
+
+ if(rest) {
+ d = 0;
+ for(i = need, from = MP_DIGITS(z) + need; i < uz; ++i, ++from) {
+ mp_digit save = *from;
+
+ *from = (*from << rest) | (d >> (MP_DIGIT_BIT - rest));
+ d = save;
+ }
+
+ d >>= (MP_DIGIT_BIT - rest);
+ if(d != 0) {
+ *from = d;
+ uz += extra;
+ }
+ }
+
+ MP_USED(z) = uz;
+ CLAMP(z);
+
+ return 1;
+}
+
+/* }}} */
+
+/* {{{ s_qsub(z, p2) */
+
+/* Compute z = 2^p2 - |z|; requires that 2^p2 >= |z|
+ The sign of the result is always zero/positive.
+ */
+STATIC int s_qsub(klisp_State *K, mp_int z, mp_size p2)
+{
+ mp_digit hi = (1 << (p2 % MP_DIGIT_BIT)), *zp;
+ mp_size tdig = (p2 / MP_DIGIT_BIT), pos;
+ mp_word w = 0;
+
+ if(!s_pad(K, z, tdig + 1))
+ return 0;
+
+ for(pos = 0, zp = MP_DIGITS(z); pos < tdig; ++pos, ++zp) {
+ w = ((mp_word) MP_DIGIT_MAX + 1) - w - (mp_word)*zp;
+
+ *zp = LOWER_HALF(w);
+ w = UPPER_HALF(w) ? 0 : 1;
+ }
+
+ w = ((mp_word) MP_DIGIT_MAX + 1 + hi) - w - (mp_word)*zp;
+ *zp = LOWER_HALF(w);
+
+ assert(UPPER_HALF(w) != 0); /* no borrow out should be possible */
+
+ MP_SIGN(z) = MP_ZPOS;
+ CLAMP(z);
+
+ return 1;
+}
+
+/* }}} */
+
+/* {{{ s_dp2k(z) */
+
+STATIC int s_dp2k(mp_int z)
+{
+ int k = 0;
+ mp_digit *dp = MP_DIGITS(z), d;
+
+ if(MP_USED(z) == 1 && *dp == 0)
+ return 1;
+
+ while(*dp == 0) {
+ k += MP_DIGIT_BIT;
+ ++dp;
+ }
+
+ d = *dp;
+ while((d & 1) == 0) {
+ d >>= 1;
+ ++k;
+ }
+
+ return k;
+}
+
+/* }}} */
+
+/* {{{ s_isp2(z) */
+
+STATIC int s_isp2(mp_int z)
+{
+ mp_size uz = MP_USED(z), k = 0;
+ mp_digit *dz = MP_DIGITS(z), d;
+
+ while(uz > 1) {
+ if(*dz++ != 0)
+ return -1;
+ k += MP_DIGIT_BIT;
+ --uz;
+ }
+
+ d = *dz;
+ while(d > 1) {
+ if(d & 1)
+ return -1;
+ ++k; d >>= 1;
+ }
+
+ return (int) k;
+}
+
+/* }}} */
+
+/* {{{ s_2expt(z, k) */
+
+STATIC int s_2expt(klisp_State *K, mp_int z, mp_small k)
+{
+ mp_size ndig, rest;
+ mp_digit *dz;
+
+ ndig = (k + MP_DIGIT_BIT) / MP_DIGIT_BIT;
+ rest = k % MP_DIGIT_BIT;
+
+ if(!s_pad(K, z, ndig))
+ return 0;
+
+ dz = MP_DIGITS(z);
+ ZERO(dz, ndig);
+ *(dz + ndig - 1) = (1 << rest);
+ MP_USED(z) = ndig;
+
+ return 1;
+}
+
+/* }}} */
+
+/* {{{ s_norm(a, b) */
+
+STATIC int s_norm(klisp_State *K, mp_int a, mp_int b)
+{
+ mp_digit d = b->digits[MP_USED(b) - 1];
+ int k = 0;
+
+ while(d < (mp_digit) (1 << (MP_DIGIT_BIT - 1))) { /* d < (MP_RADIX / 2) */
+ d <<= 1;
+ ++k;
+ }
+
+ /* These multiplications can't fail */
+ if(k != 0) {
+ (void) s_qmul(K, a, (mp_size) k);
+ (void) s_qmul(K, b, (mp_size) k);
+ }
+
+ return k;
+}
+
+/* }}} */
+
+/* {{{ s_brmu(z, m) */
+
+STATIC mp_result s_brmu(klisp_State *K, mp_int z, mp_int m)
+{
+ mp_size um = MP_USED(m) * 2;
+
+ if(!s_pad(K, z, um))
+ return MP_MEMORY;
+
+ s_2expt(K, z, MP_DIGIT_BIT * um);
+ return mp_int_div(K, z, m, z, NULL);
+}
+
+/* }}} */
+
+/* {{{ s_reduce(x, m, mu, q1, q2) */
+
+STATIC int s_reduce(klisp_State *K, mp_int x, mp_int m, mp_int mu,
+ mp_int q1, mp_int q2)
+{
+ mp_size um = MP_USED(m), umb_p1, umb_m1;
+
+ umb_p1 = (um + 1) * MP_DIGIT_BIT;
+ umb_m1 = (um - 1) * MP_DIGIT_BIT;
+
+ if(mp_int_copy(K, x, q1) != MP_OK)
+ return 0;
+
+ /* Compute q2 = floor((floor(x / b^(k-1)) * mu) / b^(k+1)) */
+ s_qdiv(q1, umb_m1);
+ UMUL(K, q1, mu, q2);
+ s_qdiv(q2, umb_p1);
+
+ /* Set x = x mod b^(k+1) */
+ s_qmod(x, umb_p1);
+
+ /* Now, q is a guess for the quotient a / m.
+ Compute x - q * m mod b^(k+1), replacing x. This may be off
+ by a factor of 2m, but no more than that.
+ */
+ UMUL(K, q2, m, q1);
+ s_qmod(q1, umb_p1);
+ (void) mp_int_sub(K, x, q1, x); /* can't fail */
+
+ /* The result may be < 0; if it is, add b^(k+1) to pin it in the
+ proper range. */
+ if((CMPZ(x) < 0) && !s_qsub(K, x, umb_p1))
+ return 0;
+
+ /* If x > m, we need to back it off until it is in range.
+ This will be required at most twice. */
+ if(mp_int_compare(x, m) >= 0) {
+ (void) mp_int_sub(K, x, m, x);
+ if(mp_int_compare(x, m) >= 0)
+ (void) mp_int_sub(K, x, m, x);
+ }
+
+ /* At this point, x has been properly reduced. */
+ return 1;
+}
+
+/* }}} */
+
+/* {{{ s_embar(a, b, m, mu, c) */
+
+/* Perform modular exponentiation using Barrett's method, where mu is
+ the reduction constant for m. Assumes a < m, b > 0. */
+STATIC mp_result s_embar(klisp_State *K, mp_int a, mp_int b, mp_int m,
+ mp_int mu, mp_int c)
+{
+ mp_digit *db, *dbt, umu, d;
+ mpz_t temp[3];
+ mp_result res;
+ int last = 0;
+
+ umu = MP_USED(mu); db = MP_DIGITS(b); dbt = db + MP_USED(b) - 1;
+
+ while(last < 3) {
+ SETUP(mp_int_init_size(K, TEMP(last), 4 * umu), last);
+ ZERO(MP_DIGITS(TEMP(last - 1)), MP_ALLOC(TEMP(last - 1)));
+ }
+
+ (void) mp_int_set_value(K, c, 1);
+
+ /* Take care of low-order digits */
+ while(db < dbt) {
+ int i;
+
+ for(d = *db, i = MP_DIGIT_BIT; i > 0; --i, d >>= 1) {
+ if(d & 1) {
+ /* The use of a second temporary avoids allocation */
+ UMUL(K, c, a, TEMP(0));
+ if(!s_reduce(K, TEMP(0), m, mu, TEMP(1), TEMP(2))) {
+ res = MP_MEMORY; goto CLEANUP;
+ }
+ mp_int_copy(K, TEMP(0), c);
+ }
+
+
+ USQR(K, a, TEMP(0));
+ assert(MP_SIGN(TEMP(0)) == MP_ZPOS);
+ if(!s_reduce(K, TEMP(0), m, mu, TEMP(1), TEMP(2))) {
+ res = MP_MEMORY; goto CLEANUP;
+ }
+ assert(MP_SIGN(TEMP(0)) == MP_ZPOS);
+ mp_int_copy(K, TEMP(0), a);
+
+
+ }
+
+ ++db;
+ }
+
+ /* Take care of highest-order digit */
+ d = *dbt;
+ for(;;) {
+ if(d & 1) {
+ UMUL(K, c, a, TEMP(0));
+ if(!s_reduce(K, TEMP(0), m, mu, TEMP(1), TEMP(2))) {
+ res = MP_MEMORY; goto CLEANUP;
+ }
+ mp_int_copy(K, TEMP(0), c);
+ }
+
+ d >>= 1;
+ if(!d) break;
+
+ USQR(K, a, TEMP(0));
+ if(!s_reduce(K, TEMP(0), m, mu, TEMP(1), TEMP(2))) {
+ res = MP_MEMORY; goto CLEANUP;
+ }
+ (void) mp_int_copy(K, TEMP(0), a);
+ }
+
+ CLEANUP:
+ while(--last >= 0)
+ mp_int_clear(K, TEMP(last));
+
+ return res;
+}
+
+/* }}} */
+
+/* {{{ s_udiv(a, b) */
+
+/* Precondition: a >= b and b > 0
+ Postcondition: a' = a / b, b' = a % b
+ */
+STATIC mp_result s_udiv(klisp_State *K, mp_int a, mp_int b)
+{
+ mpz_t q, r, t;
+ mp_size ua, ub, qpos = 0;
+ mp_digit *da, btop;
+ mp_result res = MP_OK;
+ int k, skip = 0;
+
+ /* Force signs to positive */
+ MP_SIGN(a) = MP_ZPOS;
+ MP_SIGN(b) = MP_ZPOS;
+
+ /* Normalize, per Knuth */
+ k = s_norm(K, a, b);
+
+ ua = MP_USED(a); ub = MP_USED(b); btop = b->digits[ub - 1];
+ if((res = mp_int_init_size(K, &q, ua)) != MP_OK) return res;
+ if((res = mp_int_init_size(K, &t, ua + 1)) != MP_OK) goto CLEANUP;
+
+ da = MP_DIGITS(a);
+ r.digits = da + ua - 1; /* The contents of r are shared with a */
+ r.used = 1;
+ r.sign = MP_ZPOS;
+ r.alloc = MP_ALLOC(a);
+ ZERO(t.digits, t.alloc);
+
+ /* Solve for quotient digits, store in q.digits in reverse order */
+ while(r.digits >= da) {
+ assert(qpos <= q.alloc);
+
+ if(s_ucmp(b, &r) > 0) {
+ r.digits -= 1;
+ r.used += 1;
+
+ if(++skip > 1 && qpos > 0)
+ q.digits[qpos++] = 0;
+
+ CLAMP(&r);
+ }
+ else {
+ mp_word pfx = r.digits[r.used - 1];
+ mp_word qdigit;
+
+ if(r.used > 1 && pfx <= btop) {
+ pfx <<= MP_DIGIT_BIT / 2;
+ pfx <<= MP_DIGIT_BIT / 2;
+ pfx |= r.digits[r.used - 2];
+ }
+
+ qdigit = pfx / btop;
+ if(qdigit > MP_DIGIT_MAX) {
+ qdigit = MP_DIGIT_MAX;
+ }
+
+ s_dbmul(MP_DIGITS(b), (mp_digit) qdigit, t.digits, ub);
+ t.used = ub + 1; CLAMP(&t);
+ while(s_ucmp(&t, &r) > 0) {
+ --qdigit;
+ (void) mp_int_sub(K, &t, b, &t); /* cannot fail */
+ }
+
+ s_usub(r.digits, t.digits, r.digits, r.used, t.used);
+ CLAMP(&r);
+
+ q.digits[qpos++] = (mp_digit) qdigit;
+ ZERO(t.digits, t.used);
+ skip = 0;
+ }
+ }
+
+ /* Put quotient digits in the correct order, and discard extra zeroes */
+ q.used = qpos;
+ REV(mp_digit, q.digits, qpos);
+ CLAMP(&q);
+
+ /* Denormalize the remainder */
+ CLAMP(a);
+ if(k != 0)
+ s_qdiv(a, k);
+
+ mp_int_copy(K, a, b); /* ok: 0 <= r < b */
+ mp_int_copy(K, &q, a); /* ok: q <= a */
+
+ mp_int_clear(K, &t);
+ CLEANUP:
+ mp_int_clear(K, &q);
+ return res;
+}
+
+/* }}} */
+
+/* {{{ s_outlen(z, r) */
+
+STATIC int s_outlen(mp_int z, mp_size r)
+{
+ mp_result bits;
+ double raw;
+
+ assert(r >= MP_MIN_RADIX && r <= MP_MAX_RADIX);
+
+ bits = mp_int_count_bits(z);
+ raw = (double)bits * s_log2[r];
+
+ return (int)(raw + 0.999999);
+}
+
+/* }}} */
+
+/* {{{ s_inlen(len, r) */
+
+STATIC mp_size s_inlen(int len, mp_size r)
+{
+ double raw = (double)len / s_log2[r];
+ mp_size bits = (mp_size)(raw + 0.5);
+
+ return (mp_size)((bits + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT);
+}
+
+/* }}} */
+
+/* {{{ s_ch2val(c, r) */
+
+STATIC int s_ch2val(char c, int r)
+{
+ int out;
+
+ if(isdigit((unsigned char) c))
+ out = c - '0';
+ else if(r > 10 && isalpha((unsigned char) c))
+ out = toupper(c) - 'A' + 10;
+ else
+ return -1;
+
+ return (out >= r) ? -1 : out;
+}
+
+/* }}} */
+
+/* {{{ s_val2ch(v, caps) */
+
+STATIC char s_val2ch(int v, int caps)
+{
+ assert(v >= 0);
+
+ if(v < 10)
+ return v + '0';
+ else {
+ char out = (v - 10) + 'a';
+
+ if(caps)
+ return toupper(out);
+ else
+ return out;
+ }
+}
+
+/* }}} */
+
+/* {{{ s_2comp(buf, len) */
+
+STATIC void s_2comp(unsigned char *buf, int len)
+{
+ int i;
+ unsigned short s = 1;
+
+ for(i = len - 1; i >= 0; --i) {
+ unsigned char c = ~buf[i];
+
+ s = c + s;
+ c = s & UCHAR_MAX;
+ s >>= CHAR_BIT;
+
+ buf[i] = c;
+ }
+
+ /* last carry out is ignored */
+}
+
+/* }}} */
+
+/* {{{ s_tobin(z, buf, *limpos) */
+
+STATIC mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad)
+{
+ mp_size uz;
+ mp_digit *dz;
+ int pos = 0, limit = *limpos;
+
+ uz = MP_USED(z); dz = MP_DIGITS(z);
+ while(uz > 0 && pos < limit) {
+ mp_digit d = *dz++;
+ int i;
+
+ for(i = sizeof(mp_digit); i > 0 && pos < limit; --i) {
+ buf[pos++] = (unsigned char)d;
+ d >>= CHAR_BIT;
+
+ /* Don't write leading zeroes */
+ if(d == 0 && uz == 1)
+ i = 0; /* exit loop without signaling truncation */
+ }
+
+ /* Detect truncation (loop exited with pos >= limit) */
+ if(i > 0) break;
+
+ --uz;
+ }
+
+ if(pad != 0 && (buf[pos - 1] >> (CHAR_BIT - 1))) {
+ if(pos < limit)
+ buf[pos++] = 0;
+ else
+ uz = 1;
+ }
+
+ /* Digits are in reverse order, fix that */
+ REV(unsigned char, buf, pos);
+
+ /* Return the number of bytes actually written */
+ *limpos = pos;
+
+ return (uz == 0) ? MP_OK : MP_TRUNC;
+}
+
+/* }}} */
+
+/* {{{ s_print(tag, z) */
+
+#if DEBUG
+void s_print(char *tag, mp_int z)
+{
+ int i;
+
+ fprintf(stderr, "%s: %c ", tag,
+ (MP_SIGN(z) == MP_NEG) ? '-' : '+');
+
+ for(i = MP_USED(z) - 1; i >= 0; --i)
+ fprintf(stderr, "%0*X", (int)(MP_DIGIT_BIT / 4), z->digits[i]);
+
+ fputc('\n', stderr);
+
+}
+
+void s_print_buf(char *tag, mp_digit *buf, mp_size num)
+{
+ int i;
+
+ fprintf(stderr, "%s: ", tag);
+
+ for(i = num - 1; i >= 0; --i)
+ fprintf(stderr, "%0*X", (int)(MP_DIGIT_BIT / 4), buf[i]);
+
+ fputc('\n', stderr);
+}
+#endif
+
+/* }}} */
+
+/* HERE THERE BE DRAGONS */
diff --git a/src/imath.h b/src/imath.h
@@ -0,0 +1,298 @@
+/*
+** imath.h
+** Arbitrary precision integer arithmetic routines.
+** See Copyright Notice in klisp.h
+*/
+
+/*
+** SOURCE NOTE: This is mostly from the IMath library, written by
+** M.J. Fromberger. It is adapted to klisp, mainly in the use of the
+** klisp allocator and fixing of digit size to 32 bits.
+** Imported from version (1.15) updated 01-Feb-2011 at 03:10 PM.
+*/
+
+#ifndef IMATH_H_
+#define IMATH_H_
+
+#include <limits.h>
+
+/* Andres Navarro: use c99 constants */
+#ifndef USE_C99
+#define USE_C99 1
+#endif
+
+/* Andres Navarro: klisp includes */
+#include "kobject.h"
+#include "kstate.h"
+
+#ifdef USE_C99
+#include <stdint.h>
+#endif
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#if USE_C99
+typedef unsigned char mp_sign;
+typedef uint32_t mp_size;
+typedef int mp_result;
+typedef int32_t mp_small; /* must be a signed type */
+typedef uint32_t mp_usmall; /* must be an unsigned type */
+typedef uint32_t mp_digit;
+typedef uint64_t mp_word;
+#else /* USE_C99 */
+typedef unsigned char mp_sign;
+typedef unsigned int mp_size;
+typedef int mp_result;
+typedef long mp_small; /* must be a signed type */
+typedef unsigned long mp_usmall; /* must be an unsigned type */
+#ifdef USE_LONG_LONG
+typedef unsigned int mp_digit;
+typedef unsigned long long mp_word;
+#else /* USE_LONG_LONG */
+typedef unsigned short mp_digit;
+typedef unsigned int mp_word;
+#endif /* USE_LONG_LONG */
+#endif /* USE_C99 */
+
+/* Andres Navarro: Use kobject type instead */
+/*
+typedef struct mpz {
+ mp_digit single;
+ mp_digit *digits;
+ mp_size alloc;
+ mp_size used;
+ mp_sign sign;
+} mpz_t, *mp_int;
+*/
+typedef Bigint mpz_t, *mp_int;
+
+#define MP_SINGLE(Z) ((Z)->single) /* added to correct check in mp_int_clear */
+#define MP_DIGITS(Z) ((Z)->digits)
+#define MP_ALLOC(Z) ((Z)->alloc)
+#define MP_USED(Z) ((Z)->used)
+#define MP_SIGN(Z) ((Z)->sign)
+
+extern const mp_result MP_OK;
+extern const mp_result MP_FALSE;
+extern const mp_result MP_TRUE;
+extern const mp_result MP_MEMORY;
+extern const mp_result MP_RANGE;
+extern const mp_result MP_UNDEF;
+extern const mp_result MP_TRUNC;
+extern const mp_result MP_BADARG;
+extern const mp_result MP_MINERR;
+
+#define MP_DIGIT_BIT (sizeof(mp_digit) * CHAR_BIT)
+#define MP_WORD_BIT (sizeof(mp_word) * CHAR_BIT)
+/* Andres Navarro: USE_C99 */
+#ifdef USE_C99
+#define MP_SMALL_MIN INT32_MIN
+#define MP_SMALL_MAX INT32_MAX
+#define MP_USMALL_MIN UINT32_MIN
+#define MP_USMALL_MAX UINT32_MAX
+#define MP_DIGIT_MAX ((uint64_t) UINT32_MAX)
+#define MP_WORD_MAX UINT64_MAX
+#else /* USE_C99 */
+#define MP_SMALL_MIN LONG_MIN
+#define MP_SMALL_MAX LONG_MAX
+#define MP_USMALL_MIN ULONG_MIN
+#define MP_USMALL_MAX ULONG_MAX
+#ifdef USE_LONG_LONG
+# ifndef ULONG_LONG_MAX
+# ifdef ULLONG_MAX
+# define ULONG_LONG_MAX ULLONG_MAX
+# else
+# error "Maximum value of unsigned long long not defined!"
+# endif
+# endif
+# define MP_DIGIT_MAX (UINT_MAX * 1ULL)
+# define MP_WORD_MAX ULONG_LONG_MAX
+#else /* USE_LONG_LONG */
+# define MP_DIGIT_MAX (USHRT_MAX * 1UL)
+# define MP_WORD_MAX (UINT_MAX * 1UL)
+#endif /* USE_LONG_LONG */
+#endif /* USE_C99 */
+
+#define MP_MIN_RADIX 2
+#define MP_MAX_RADIX 36
+
+/* Values with fewer than this many significant digits use the
+ standard multiplication algorithm; otherwise, a recursive algorithm
+ is used. Choose a value to suit your platform.
+ */
+#define MP_MULT_THRESH 22
+
+#define MP_DEFAULT_PREC 8 /* default memory allocation, in digits */
+
+extern const mp_sign MP_NEG;
+extern const mp_sign MP_ZPOS;
+
+#define mp_int_is_odd(Z) ((Z)->digits[0] & 1)
+#define mp_int_is_even(Z) !((Z)->digits[0] & 1)
+
+/* NOTE: this doesn't use the allocator */
+mp_result mp_int_init(mp_int z);
+mp_int mp_int_alloc(klisp_State *K);
+mp_result mp_int_init_size(klisp_State *K, mp_int z, mp_size prec);
+mp_result mp_int_init_copy(klisp_State *K, mp_int z, mp_int old);
+mp_result mp_int_init_value(klisp_State *K, mp_int z, mp_small value);
+mp_result mp_int_set_value(klisp_State *K, mp_int z, mp_small value);
+void mp_int_clear(klisp_State *K, mp_int z);
+void mp_int_free(klisp_State *K, mp_int z);
+
+mp_result mp_int_copy(klisp_State *K, mp_int a, mp_int c); /* c = a */
+/* NOTE: this doesn't use the allocator */
+void mp_int_swap(mp_int a, mp_int c); /* swap a, c */
+/* NOTE: this doesn't use the allocator */
+void mp_int_zero(mp_int z); /* z = 0 */
+mp_result mp_int_abs(klisp_State *K, mp_int a, mp_int c); /* c = |a| */
+mp_result mp_int_neg(klisp_State *K, mp_int a, mp_int c); /* c = -a */
+/* c = a + b */
+mp_result mp_int_add(klisp_State *K, mp_int a, mp_int b, mp_int c);
+mp_result mp_int_add_value(klisp_State *K, mp_int a, mp_small value,
+ mp_int c);
+/* c = a - b */
+mp_result mp_int_sub(klisp_State *K, mp_int a, mp_int b, mp_int c);
+mp_result mp_int_sub_value(klisp_State *K, mp_int a, mp_small value,
+ mp_int c);
+/* c = a * b */
+mp_result mp_int_mul(klisp_State *K, mp_int a, mp_int b, mp_int c);
+mp_result mp_int_mul_value(klisp_State *K, mp_int a, mp_small value,
+ mp_int c);
+mp_result mp_int_mul_pow2(klisp_State *K, mp_int a, mp_small p2, mp_int c);
+mp_result mp_int_sqr(klisp_State *K, mp_int a, mp_int c); /* c = a * a */
+/* q = a / b */
+/* r = a % b */
+mp_result mp_int_div(klisp_State *K, mp_int a, mp_int b, mp_int q,
+ mp_int r);
+/* q = a / value */
+/* r = a % value */
+mp_result mp_int_div_value(klisp_State *K, mp_int a, mp_small value,
+ mp_int q, mp_small *r);
+/* q = a / 2^p2 */
+/* r = q % 2^p2 */
+mp_result mp_int_div_pow2(klisp_State *K, mp_int a, mp_small p2,
+ mp_int q, mp_int r);
+/* c = a % m */
+mp_result mp_int_mod(klisp_State *K, mp_int a, mp_int m, mp_int c);
+#define mp_int_mod_value(K, A, V, R) \
+ mp_int_div_value((K), (A), (V), 0, (R))
+/* c = a^b */
+mp_result mp_int_expt(klisp_State *K, mp_int a, mp_small b, mp_int c);
+/* c = a^b */
+mp_result mp_int_expt_value(klisp_State *K, mp_small a, mp_small b, mp_int c);
+/* c = a^b */
+mp_result mp_int_expt_full(klisp_State *K, mp_int a, mp_int b, mp_int c);
+
+/* NOTE: this doesn't use the allocator */
+int mp_int_compare(mp_int a, mp_int b); /* a <=> b */
+/* NOTE: this doesn't use the allocator */
+int mp_int_compare_unsigned(mp_int a, mp_int b); /* |a| <=> |b| */
+/* NOTE: this doesn't use the allocator */
+int mp_int_compare_zero(mp_int z); /* a <=> 0 */
+/* NOTE: this doesn't use the allocator */
+int mp_int_compare_value(mp_int z, mp_small value); /* a <=> v */
+
+/* Returns true if v|a, false otherwise (including errors) */
+int mp_int_divisible_value(klisp_State *K, mp_int a, mp_small v);
+
+/* NOTE: this doesn't use the allocator */
+/* Returns k >= 0 such that z = 2^k, if one exists; otherwise < 0 */
+int mp_int_is_pow2(mp_int z);
+
+mp_result mp_int_exptmod(klisp_State *K, mp_int a, mp_int b, mp_int m,
+ mp_int c); /* c = a^b (mod m) */
+mp_result mp_int_exptmod_evalue(klisp_State *K, mp_int a, mp_small value,
+ mp_int m, mp_int c); /* c = a^v (mod m) */
+mp_result mp_int_exptmod_bvalue(klisp_State *K, mp_small value, mp_int b,
+ mp_int m, mp_int c); /* c = v^b (mod m) */
+mp_result mp_int_exptmod_known(klisp_State *K, mp_int a, mp_int b,
+ mp_int m, mp_int mu,
+ mp_int c); /* c = a^b (mod m) */
+mp_result mp_int_redux_const(klisp_State *K, mp_int m, mp_int c);
+
+/* c = 1/a (mod m) */
+mp_result mp_int_invmod(klisp_State *K, mp_int a, mp_int m, mp_int c);
+
+/* c = gcd(a, b) */
+mp_result mp_int_gcd(klisp_State *K, mp_int a, mp_int b, mp_int c);
+
+/* c = gcd(a, b) */
+/* c = ax + by */
+mp_result mp_int_egcd(klisp_State *K, mp_int a, mp_int b, mp_int c,
+ mp_int x, mp_int y);
+
+/* c = lcm(a, b) */
+mp_result mp_int_lcm(klisp_State *K, mp_int a, mp_int b, mp_int c);
+
+/* c = floor(a^{1/b}) */
+mp_result mp_int_root(klisp_State *K, mp_int a, mp_small b, mp_int c);
+/* c = floor(sqrt(a)) */
+#define mp_int_sqrt(K, a, c) mp_int_root((K), a, 2, c)
+
+/* Convert to a small int, if representable; else MP_RANGE */
+/* NOTE: this doesn't use the allocator */
+mp_result mp_int_to_int(mp_int z, mp_small *out);
+/* NOTE: this doesn't use the allocator */
+mp_result mp_int_to_uint(mp_int z, mp_usmall *out);
+
+/* Convert to nul-terminated string with the specified radix, writing at
+ most limit characters including the nul terminator */
+mp_result mp_int_to_string(klisp_State *K, mp_int z, mp_size radix,
+ char *str, int limit);
+
+/* Return the number of characters required to represent
+ z in the given radix. May over-estimate. */
+/* NOTE: this doesn't use the allocator */
+mp_result mp_int_string_len(mp_int z, mp_size radix);
+
+/* Read zero-terminated string into z */
+mp_result mp_int_read_string(klisp_State *K, mp_int z, mp_size radix,
+ const char *str);
+mp_result mp_int_read_cstring(klisp_State *K, mp_int z, mp_size radix,
+ const char *str, char **end);
+
+/* Return the number of significant bits in z */
+/* NOTE: this doesn't use the allocator */
+mp_result mp_int_count_bits(mp_int z);
+
+/* Convert z to two's complement binary, writing at most limit bytes */
+mp_result mp_int_to_binary(klisp_State *K, mp_int z, unsigned char *buf,
+ int limit);
+
+/* Read a two's complement binary value into z from the given buffer */
+mp_result mp_int_read_binary(klisp_State *K, mp_int z, unsigned char *buf,
+ int len);
+
+/* Return the number of bytes required to represent z in binary. */
+/* NOTE: this doesn't use the allocator */
+mp_result mp_int_binary_len(mp_int z);
+
+/* Convert z to unsigned binary, writing at most limit bytes */
+mp_result mp_int_to_unsigned(klisp_State *K, mp_int z, unsigned char *buf,
+ int limit);
+
+/* Read an unsigned binary value into z from the given buffer */
+mp_result mp_int_read_unsigned(klisp_State *K, mp_int z, unsigned char *buf,
+ int len);
+
+/* Return the number of bytes required to represent z as unsigned output */
+/* NOTE: this doesn't use the allocator */
+mp_result mp_int_unsigned_len(mp_int z);
+
+/* Return a statically allocated string describing error code res */
+/* NOTE: this doesn't use the allocator */
+const char *mp_error_string(mp_result res);
+
+#if DEBUG
+void s_print(char *tag, mp_int z);
+void s_print_buf(char *tag, mp_digit *buf, mp_size num);
+#endif
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* end IMATH_H_ */
diff --git a/src/kinteger.c b/src/kinteger.c
@@ -14,21 +14,6 @@
#include "kstate.h"
#include "kmem.h"
-#define bind_iter kbind_bigint_iter
-#define iter_has_next kbigint_iter_has_more
-#define iter_next kbigint_iter_next
-#define iter_update_last kbigint_iter_update_last
-
-#define LOG_BASE(base) (log(2.0) / log(base))
-
-Bigint_Node *make_new_node(klisp_State *K, uint32_t digit)
-{
- Bigint_Node *node = klispM_new(K, Bigint_Node);
- node->digit = digit;
- node->next_xor_prev = (uintptr_t) 0; /* NULL ^ NULL: 0 */
- return node;
-}
-
/* for now used only for reading */
/* NOTE: is uint to allow INT32_MIN as positive argument in read */
TValue kbigint_new(klisp_State *K, bool sign, uint32_t digit)
@@ -43,16 +28,14 @@ TValue kbigint_new(klisp_State *K, bool sign, uint32_t digit)
new_bigint->flags = 0;
/* bigint specific fields */
-
- /* GC: root bigint */
- /* put dummy value to work if garbage collections happens while allocating
- node */
- new_bigint->sign_size = 0;
- new_bigint->first = new_bigint->last = NULL;
-
- Bigint_Node *node = make_new_node(K, digit);
- new_bigint->first = new_bigint->last = node;
- new_bigint->sign_size = sign? -1 : 1;
+ /* If later changed to alloc obj:
+ GC: root bigint & put dummy value to work if garbage collections
+ happens while allocating array */
+ new_bigint->single = digit;
+ new_bigint->digits = &(new_bigint->single);
+ new_bigint->alloc = 1;
+ new_bigint->used = 1;
+ new_bigint->sign = sign? MP_NEG : MP_ZPOS;
return gc2bigint(new_bigint);
}
@@ -60,19 +43,9 @@ TValue kbigint_new(klisp_State *K, bool sign, uint32_t digit)
/* used in write to destructively get the digits */
TValue kbigint_copy(klisp_State *K, TValue src)
{
- Bigint *srcb = tv2bigint(src);
- /* iterate in little endian mode */
- bind_iter(iter, srcb, false);
- uint32_t digit = iter_next(iter);
- /* GC: root copy */
- TValue copy = kbigint_new(K, kbigint_sign(srcb), digit);
- Bigint *copyb = tv2bigint(copy);
-
- while(iter_has_next(iter)) {
- uint32_t digit = iter_next(iter);
- kbigint_add_node(copyb, make_new_node(K, digit));
- }
-
+ TValue copy = kbigint_new(K, false, 0);
+ /* arguments are in reverse order with respect to mp_int_copy */
+ UNUSED(mp_int_init_copy(K, tv2bigint(copy), tv2bigint(src)));
return copy;
}
@@ -83,227 +56,100 @@ void kbigint_add_digit(klisp_State *K, TValue tv_bigint, int32_t base,
int32_t digit)
{
Bigint *bigint = tv2bigint(tv_bigint);
- /* iterate in little endian mode */
- bind_iter(iter, bigint, false);
-
- uint64_t carry = digit;
-
- while(iter_has_next(iter)) {
- uint32_t cur = iter_next(iter);
- /* I hope the compiler understand that this should be a
- 32bits x 32bits = 64bits multiplication instruction */
- uint64_t res_carry = (uint64_t) cur * base + carry;
- carry = res_carry >> 32;
- uint32_t new_cur = (uint32_t) res_carry;
- iter_update_last(iter, new_cur);
- }
-
- if (carry != 0) {
- /* must add one node to the bigint */
- kbigint_add_node(bigint, make_new_node(K, (uint32_t) carry));
- }
+ UNUSED(mp_int_mul_value(K, bigint, base, bigint));
+ UNUSED(mp_int_add_value(K, bigint, digit, bigint));
}
/* This is used by the writer to get the digits of a number
tv_bigint must be positive */
int32_t kbigint_remove_digit(klisp_State *K, TValue tv_bigint, int32_t base)
{
- assert(kbigint_has_digits(K, tv_bigint));
-
+ UNUSED(K);
Bigint *bigint = tv2bigint(tv_bigint);
-
- /* iterate in big endian mode */
- bind_iter(iter, bigint, true);
-
- uint32_t result = 0;
- uint32_t rest = 0;
- uint32_t divisor = base;
-
- while(iter_has_next(iter)) {
- uint64_t dividend = (((uint64_t) rest) << 32) |
- (uint64_t) iter_next(iter);
-
- if (dividend >= divisor) { /* avoid division if possible */
- /* I hope the compiler calculates this only once and
- is able to get that this is a 64bits by 32bits division
- instruction */
- result = (uint32_t) (dividend / divisor);
- rest = (uint32_t) (dividend % divisor);
- } else {
- result = 0;
- rest = (uint32_t) dividend;
- }
- iter_update_last(iter, result);
- }
-
- /* rest now contains the last digit & result the value of the top node */
-
- /* adjust the node list, at most the bigint should lose a node */
- if (bigint->first->digit == 0) {
- Bigint_Node *node = kbigint_remove_node(bigint);
- klispM_free(K, node);
- }
-
- return (int32_t) rest;
+ int32_t r;
+ UNUSED(mp_int_div_value(K, bigint, base, bigint, &r));
+ return r;
}
/* This is used by write to test if there is any digit left to print */
bool kbigint_has_digits(klisp_State *K, TValue tv_bigint)
{
UNUSED(K);
- return kbigint_size(tv2bigint(tv_bigint)) != 0;
+ return (mp_int_compare_zero(tv2bigint(tv_bigint)) != 0);
}
/* Mutate the bigint to have the opposite sign, used in read,
write and abs */
-void kbigint_invert_sign(TValue tv_bigint)
+void kbigint_invert_sign(klisp_State *K, TValue tv_bigint)
{
Bigint *bigint = tv2bigint(tv_bigint);
- bigint->sign_size = -bigint->sign_size;
+ UNUSED(mp_int_neg(K, bigint, bigint));
}
/* this is used by write to estimate the number of chars necessary to
print the number */
int32_t kbigint_print_size(TValue tv_bigint, int32_t base)
{
- /* count the number of bits and estimate using the log of
- the base */
- Bigint *bigint = tv2bigint(tv_bigint);
-
- int32_t num_bits = 0;
- uint32_t first_digit = bigint->first->digit;
- while(first_digit != 0) {
- ++num_bits;
- first_digit >>= 1;
- }
- num_bits += 32 * (kbigint_size(bigint)) - 2 ;
- /* add 1.5 for safety */
- return (int32_t)(LOG_BASE(base) * num_bits + 1.0);
+ return mp_int_string_len(tv2bigint(tv_bigint), base);
}
bool kbigint_eqp(TValue tv_bigint1, TValue tv_bigint2)
{
- Bigint *bigint1 = tv2bigint(tv_bigint1);
- Bigint *bigint2 = tv2bigint(tv_bigint2);
-
- if (bigint1->sign_size != bigint2->sign_size)
- return false;
-
- /* iterate in big endian mode */
- bind_iter(iter1, bigint1, true);
- bind_iter(iter2, bigint2, true);
-
- while(iter_has_next(iter1)) {
- uint32_t digit1 = iter_next(iter1);
- uint32_t digit2 = iter_next(iter2);
- if (digit1 != digit2)
- return false;
- }
- return true;
+ return (mp_int_compare(tv2bigint(tv_bigint1),
+ tv2bigint(tv_bigint2)) == 0);
}
bool kbigint_ltp(TValue tv_bigint1, TValue tv_bigint2)
{
- Bigint *bigint1 = tv2bigint(tv_bigint1);
- Bigint *bigint2 = tv2bigint(tv_bigint2);
-
- /* first take care of the easy sign cases */
- if (kbigint_negp(bigint1)) {
- if (kbigint_posp(bigint2)) {
- return true;
- } else {
- /* if both are negative reverse the order to compare
- as positive */
- Bigint *temp = bigint1;
- bigint1 = bigint2;
- bigint2 = temp;
- /* swap the tvalues just in case */
- TValue tv_temp = tv_bigint1;
- tv_bigint1 = tv_bigint2;
- tv_bigint2 = tv_temp;
- }
- } else if (kbigint_negp(bigint2)) {
- return false;
- }
-
- /* the the easy size cases */
- int32_t size1 = kbigint_size(bigint1);
- int32_t size2 = kbigint_size(bigint2);
-
- if (size1 < size2)
- return true;
- else if (size1 > size2)
- return false;
-
- /* size and sign equal, iterate in big endian mode */
- bind_iter(iter1, bigint1, true);
- bind_iter(iter2, bigint2, true);
-
- while(iter_has_next(iter1) && iter_has_next(iter2)) {
- uint32_t digit1 = iter_next(iter1);
- uint32_t digit2 = iter_next(iter2);
- if (digit1 < digit2)
- return true;
- else if (digit1 > digit2)
- return false;
- /* if equal we keep comparing */
- }
-
- return false;
+ return (mp_int_compare(tv2bigint(tv_bigint1),
+ tv2bigint(tv_bigint2)) < 0);
}
bool kbigint_lep(TValue tv_bigint1, TValue tv_bigint2)
{
- /* a <= b == !(a > b) == !(b < a) */
- return !kbigint_ltp(tv_bigint2, tv_bigint1);
+ return (mp_int_compare(tv2bigint(tv_bigint1),
+ tv2bigint(tv_bigint2)) <= 0);
}
bool kbigint_gtp(TValue tv_bigint1, TValue tv_bigint2)
{
- /* a > b == (b < a) */
- return kbigint_ltp(tv_bigint2, tv_bigint1);
+ return (mp_int_compare(tv2bigint(tv_bigint1),
+ tv2bigint(tv_bigint2)) > 0);
}
bool kbigint_gep(TValue tv_bigint1, TValue tv_bigint2)
{
- /* a >= b == !(a < b) */
- return !kbigint_ltp(tv_bigint1, tv_bigint2);
+ return (mp_int_compare(tv2bigint(tv_bigint1),
+ tv2bigint(tv_bigint2)) >= 0);
}
bool kbigint_negativep(TValue tv_bigint)
{
- return kbigint_negp(tv2bigint(tv_bigint));
+ return (mp_int_compare_zero(tv2bigint(tv_bigint)) < 0);
}
-/* unlike the positive? applicative this would return true on zero,
- but zero is never represented as a bigint so there is no problem */
-/* Bigints constructed from fixints could be, but we don't care about
- zero returning positive in other place than in positive? */
bool kbigint_positivep(TValue tv_bigint)
{
- return kbigint_posp(tv2bigint(tv_bigint));
+ return (mp_int_compare_zero(tv2bigint(tv_bigint)) > 0);
}
bool kbigint_oddp(TValue tv_bigint)
{
- Bigint *bigint = tv2bigint(tv_bigint);
- return ((bigint->last->digit) & 1) != 0;
+ return mp_int_is_odd(tv2bigint(tv_bigint));
}
bool kbigint_evenp(TValue tv_bigint)
{
- Bigint *bigint = tv2bigint(tv_bigint);
- return ((bigint->last->digit) & 1) == 0;
+ return mp_int_is_even(tv2bigint(tv_bigint));
}
TValue kbigint_abs(klisp_State *K, TValue tv_bigint)
{
- if (kbigint_positivep(tv_bigint)) {
- return tv_bigint;
+ if (kbigint_negativep(tv_bigint)) {
+ TValue copy = kbigint_new(K, false, 0);
+ UNUSED(mp_int_abs(K, tv2bigint(tv_bigint), tv2bigint(copy)));
+ return copy;
} else {
- TValue res = kbigint_copy(K, tv_bigint);
- kbigint_invert_sign(res);
- return res;
+ return tv_bigint;
}
}
-
diff --git a/src/kinteger.h b/src/kinteger.h
@@ -13,6 +13,7 @@
#include "kobject.h"
#include "kstate.h"
+#include "imath.h"
/* for now used only for reading */
/* NOTE: is uint and has flag to allow INT32_MIN as positive argument */
@@ -21,22 +22,23 @@ TValue kbigint_new(klisp_State *K, bool sign, uint32_t digit);
/* used in write to destructively get the digits */
TValue kbigint_copy(klisp_State *K, TValue src);
+/* XXX/TODO: rewrite this to use IMath */
+
/* Create a stack allocated bigints from a fixint,
useful for mixed operations, relatively light weight compared
to creating it in the heap and burdening the gc */
#define kbind_bigint(name, fixint) \
int32_t (KUNIQUE_NAME(i)) = ivalue(fixint); \
- Bigint_Node KUNIQUE_NAME(node); \
- (KUNIQUE_NAME(node)).digit = ({ \
+ Bigint KUNIQUE_NAME(bigint); \
+ (KUNIQUE_NAME(bigint)).single = ({ \
int64_t temp = (KUNIQUE_NAME(i)); \
(uint32_t) ((temp < 0)? -temp : temp); \
}); \
- /* NULL ^ NULL: 0 */ \
- (KUNIQUE_NAME(node)).next_xor_prev = (uintptr_t) 0; \
- Bigint KUNIQUE_NAME(bigint); \
- (KUNIQUE_NAME(bigint)).first = &(KUNIQUE_NAME(node)); \
- (KUNIQUE_NAME(bigint)).last = &(KUNIQUE_NAME(node)); \
- (KUNIQUE_NAME(bigint)).sign_size = (KUNIQUE_NAME(i)) < 0? -1 : 1; \
+ (KUNIQUE_NAME(bigint)).digits = &((KUNIQUE_NAME(bigint)).single); \
+ (KUNIQUE_NAME(bigint)).alloc = 1; \
+ (KUNIQUE_NAME(bigint)).used = 1; \
+ (KUNIQUE_NAME(bigint)).sign = (KUNIQUE_NAME(i)) < 0? \
+ MP_NEG : MP_ZPOS; \
Bigint *name = &(KUNIQUE_NAME(bigint));
/* This can be used prior to calling a bigint functions
@@ -77,7 +79,7 @@ bool kbigint_evenp(TValue tv_bigint);
TValue kbigint_abs(klisp_State *K, TValue tv_bigint);
/* Mutate the bigint to have the opposite sign, used in read & write */
-void kbigint_invert_sign(TValue tv_bigint);
+void kbigint_invert_sign(klisp_State *K, TValue tv_bigint);
/* this is used by write to estimate the number of chars necessary to
print the number */
diff --git a/src/klisp.h b/src/klisp.h
@@ -39,6 +39,7 @@ void klisp_close (klisp_State *K);
/******************************************************************************
* Copyright (C) 2011 Andres Navarro. All rights reserved.
* Lua parts: Copyright (C) 1994-2010 Lua.org, PUC-Rio. All rights reserved.
+* IMath Parts: Copyright (C) 2002-2007 Michael J. Fromberger.
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
diff --git a/src/kobject.h b/src/kobject.h
@@ -255,50 +255,16 @@ typedef __attribute__((aligned (8))) union {
** Individual heap-allocated values
*/
-/*
-** Node Structure for the list of 'digits' in bignums.
-** The node list is a xor encoded doubly linked list of big endian 'digits'
-** in base 2^32.
-** It is a linked list instead of say, an array because it is immutable and
-** there is no need to index the list.
-** It is doubly linked because some operations (like + and *) work from least
-** significant to most significant and some others (like <? and >?) work the
-** other way around.
-** It is xor encoded because, well... I never had a chance to actually use
-** this and this is a perfect example of where it is actually useful: Lists
-** are to be traversed in both directions, the node is small so the space
-** gain is real (2/3 the space in 32 bits, 3/5 in 64 bits) and all the code
-** using the pointers is encapsulated in a (relatively) small module.
-*/
-
-typedef struct __attribute__ ((__packed__)) {
- uint32_t digit;
- uintptr_t next_xor_prev;
-} Bigint_Node;
-
typedef struct __attribute__ ((__packed__)) {
CommonHeader;
- Bigint_Node *first;
- Bigint_Node *last;
- int32_t sign_size;
+/* These are all from IMath (XXX: find a way to use mp_types directly) */
+ uint32_t single;
+ uint32_t *digits;
+ uint32_t alloc;
+ uint32_t used;
+ unsigned char sign;
} Bigint;
-/* macros to access size/sign */
-#define kbigint_sign(b_) ((b_)->sign_size < 0)
-#define kbigint_negp(b_) (kbigint_sign(b_))
-#define kbigint_posp(b_) (!kbigint_sign(b_))
-#define kbigint_size(b_) ({ int32_t ss = (b_)->sign_size; \
- ss < 0? -ss : ss;})
-
-/* Java Style Iterator for the xor list */
-typedef struct {
- Bigint_Node *cur;
- Bigint_Node *next;
-} Bigint_Iter;
-
-#define kxor_next(cur, prev) \
- ((Bigint_Node *) ((cur)->next_xor_prev ^ (uintptr_t) (prev)))
-
/* REFACTOR: move these macros somewhere else */
/* NOTE: The use of the intermediate KCONCAT is needed to allow
expansion of the __LINE__ macro. */
@@ -306,78 +272,6 @@ typedef struct {
#define KCONCAT(a, b) KCONCAT_(a, b)
#define KUNIQUE_NAME(prefix) KCONCAT(prefix, __LINE__ )
-#define kbind_bigint_iter(name, bigint, be) \
- Bigint_Iter KUNIQUE_NAME(iter); \
- Bigint_Iter *(name) = &(KUNIQUE_NAME(iter)); \
- kbigint_iter_init((name), (bigint), (be));
-
-inline void kbigint_iter_init(Bigint_Iter *i, Bigint *bigint, bool big_endian)
-{
- i->cur = (Bigint_Node *) NULL;
- i->next = big_endian? bigint->first : bigint->last;
-}
-
-inline bool kbigint_iter_has_more(Bigint_Iter *i)
-{
- return i->next != NULL;
-}
-
-inline uint32_t kbigint_iter_next(Bigint_Iter *i)
-{
- assert(kbigint_iter_has_more(i));
-
- Bigint_Node *next_next = kxor_next(i->next, i->cur);
- i->cur = i->next;
- i->next = next_next;
- return i->cur->digit;
-}
-
-/* this is needed for add_digit */
-inline void kbigint_iter_update_last(Bigint_Iter *i, uint32_t digit)
-{
- assert(i->cur != NULL);
- i->cur->digit = digit;
-}
-
-/* Helper for adding nodes to the head of the list.
- Neither can be NULL.
- The second argument should have NULL as previous for this to work */
-inline void kbigint_node_cons(Bigint_Node *head, Bigint_Node *tail)
-{
- head->next_xor_prev = (uintptr_t) tail ^ (uintptr_t) NULL;
- tail->next_xor_prev = (tail->next_xor_prev ^ (uintptr_t) NULL) ^
- (uintptr_t) head;
-}
-
-/* used in add_digit, bigint has at least one node */
-inline void kbigint_add_node(Bigint *bigint, Bigint_Node *node)
-{
- kbigint_node_cons(node, bigint->first);
- bigint->first = node;
- bigint->sign_size += bigint->sign_size < 0? -1 : 1;
-}
-
-/* Helper for removing a node from the head of the list.
- The argument should have NULL as previous for this to work */
-inline Bigint_Node *kbigint_remove_node(Bigint *bigint)
-{
- Bigint_Node *head = bigint->first;
- assert(head != NULL);
- Bigint_Node *tail = (Bigint_Node *) (head->next_xor_prev ^
- (uintptr_t) NULL);
- if (tail == NULL) { /* last node removed */
- bigint->sign_size = 0;
- bigint->first = bigint->last = (Bigint_Node *) NULL;
- } else {
- tail->next_xor_prev = (tail->next_xor_prev ^ (uintptr_t) head) ^
- (uintptr_t) NULL;
- bigint->first = tail;
- bigint->sign_size -= bigint->sign_size < 0? -1 : 1;
- }
- head->next_xor_prev = 0; /* NULL ^ NULL: 0 */
- return head;
-}
-
typedef struct __attribute__ ((__packed__)) {
CommonHeader;
TValue mark; /* for cycle/sharing aware algorithms */
diff --git a/src/kstate.c b/src/kstate.c
@@ -35,6 +35,8 @@
#include "kgpairs_lists.h" /* for creating list_app */
+#include "imath.h" /* for memory freeing */
+
/*
** State creation and destruction
*/
@@ -456,14 +458,7 @@ void klisp_close (klisp_State *K)
switch(type) {
case K_TBIGINT: {
- Bigint *bigint = (Bigint *)obj;
- int size = kbigint_size(bigint);
- Bigint_Node *node;
- while(size--) {
- node = kbigint_remove_node(bigint);
- klispM_free(K, node);
- }
- klispM_free(K, bigint);
+ mp_int_free(K, (Bigint *)obj);
break;
}
case K_TPAIR:
diff --git a/src/ktoken.c b/src/ktoken.c
@@ -399,7 +399,7 @@ TValue ktok_read_number(klisp_State *K, bool is_pos)
return i2tv(fixint);
} else {
if (!is_pos)
- kbigint_invert_sign(bigint_res);
+ kbigint_invert_sign(K, bigint_res);
return bigint_res;
}
}
diff --git a/src/kwrite.c b/src/kwrite.c
@@ -50,7 +50,7 @@ void kw_print_bigint(klisp_State *K, TValue bigint)
TValue copy = kbigint_copy(K, bigint);
/* must work with positive bigint to get the digits */
if (kbigint_negativep(bigint))
- kbigint_invert_sign(copy);
+ kbigint_invert_sign(K, copy);
while(kbigint_has_digits(K, copy)) {
int32_t digit = kbigint_remove_digit(K, copy, 10);
