diff options
Diffstat (limited to 'lisp/mp')
| -rw-r--r-- | lisp/mp/mp.c | 822 | ||||
| -rw-r--r-- | lisp/mp/mp.h | 435 | ||||
| -rw-r--r-- | lisp/mp/mpi.c | 1656 | ||||
| -rw-r--r-- | lisp/mp/mpr.c | 436 |
4 files changed, 3349 insertions, 0 deletions
diff --git a/lisp/mp/mp.c b/lisp/mp/mp.c new file mode 100644 index 0000000..78b7a0e --- /dev/null +++ b/lisp/mp/mp.c @@ -0,0 +1,822 @@ +/* + * Copyright (c) 2002 by The XFree86 Project, Inc. + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE XFREE86 PROJECT BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, + * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF + * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * Except as contained in this notice, the name of the XFree86 Project shall + * not be used in advertising or otherwise to promote the sale, use or other + * dealings in this Software without prior written authorization from the + * XFree86 Project. + * + * Author: Paulo César Pereira de Andrade + */ + +/* $XFree86: xc/programs/xedit/lisp/mp/mp.c,v 1.3 2002/11/20 07:44:43 paulo Exp $ */ + +#include "mp.h" + +/* + * TODO: + * o Optimize squaring + * o Write better division code and move from mpi.c to here + * o Make multiplication code don't required memory to be zeroed + * + The first step is easy, just multiply the low word, + * then the high word, that may overlap with the result + * of the first multiply (in case of carry), and then + * just make sure carry is properly propagated in the + * subsequent multiplications. + * + Some code needs also to be rewritten because some + * intermediate addition code in mp_mul, mp_karatsuba_mul, + * and mp_toom_mul is assuming the memory is zeroed. + */ + +/* + * Prototypes + */ + /* out of memory handler */ +static void mp_outmem(void); + + /* memory allocation fallback functions */ +static void *_mp_malloc(size_t); +static void *_mp_calloc(size_t, size_t); +static void *_mp_realloc(void*, size_t); +static void _mp_free(void*); + +/* + * Initialization + */ +static mp_malloc_fun __mp_malloc = _mp_malloc; +static mp_calloc_fun __mp_calloc = _mp_calloc; +static mp_realloc_fun __mp_realloc = _mp_realloc; +static mp_free_fun __mp_free = _mp_free; + +/* + * Implementation + */ +static void +mp_outmem(void) +{ + fprintf(stderr, "out of memory in MP library.\n"); + exit(1); +} + +static void * +_mp_malloc(size_t size) +{ + return (malloc(size)); +} + +void * +mp_malloc(size_t size) +{ + void *pointer = (*__mp_malloc)(size); + + if (pointer == NULL) + mp_outmem(); + + return (pointer); +} + +mp_malloc_fun +mp_set_malloc(mp_malloc_fun fun) +{ + mp_malloc_fun old = __mp_malloc; + + __mp_malloc = fun; + + return (old); +} + +static void * +_mp_calloc(size_t nmemb, size_t size) +{ + return (calloc(nmemb, size)); +} + +void * +mp_calloc(size_t nmemb, size_t size) +{ + void *pointer = (*__mp_calloc)(nmemb, size); + + if (pointer == NULL) + mp_outmem(); + + return (pointer); +} + +mp_calloc_fun +mp_set_calloc(mp_calloc_fun fun) +{ + mp_calloc_fun old = __mp_calloc; + + __mp_calloc = fun; + + return (old); +} + +static void * +_mp_realloc(void *old, size_t size) +{ + return (realloc(old, size)); +} + +void * +mp_realloc(void *old, size_t size) +{ + void *pointer = (*__mp_realloc)(old, size); + + if (pointer == NULL) + mp_outmem(); + + return (pointer); +} + +mp_realloc_fun +mp_set_realloc(mp_realloc_fun fun) +{ + mp_realloc_fun old = __mp_realloc; + + __mp_realloc = fun; + + return (old); +} + +static void +_mp_free(void *pointer) +{ + free(pointer); +} + +void +mp_free(void *pointer) +{ + (*__mp_free)(pointer); +} + +mp_free_fun +mp_set_free(mp_free_fun fun) +{ + mp_free_fun old = __mp_free; + + __mp_free = fun; + + return (old); +} + +long +mp_add(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2) +{ + BNI value; /* intermediate result */ + BNS carry; /* carry flag */ + long size; /* result size */ + + if (len1 < len2) + MP_SWAP(op1, op2, len1, len2); + + /* unroll start of loop */ + value = op1[0] + op2[0]; + rop[0] = value; + carry = value >> BNSBITS; + + /* add op1 and op2 */ + for (size = 1; size < len2; size++) { + value = op1[size] + op2[size] + carry; + rop[size] = value; + carry = value >> BNSBITS; + } + if (rop != op1) { + for (; size < len1; size++) { + value = op1[size] + carry; + rop[size] = value; + carry = value >> BNSBITS; + } + } + else { + /* if rop == op1, than just adjust carry */ + for (; carry && size < len1; size++) { + value = op1[size] + carry; + rop[size] = value; + carry = value >> BNSBITS; + } + size = len1; + } + if (carry) + rop[size++] = carry; + + return (size); +} + +long +mp_sub(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2) +{ + long svalue; /* intermediate result */ + BNS carry; /* carry flag */ + long size; /* result size */ + + /* special case */ + if (op1 == op2) { + rop[0] = 0; + + return (1); + } + + /* unroll start of loop */ + svalue = op1[0] - op2[0]; + rop[0] = svalue; + carry = svalue < 0; + + /* subtracts op2 from op1 */ + for (size = 1; size < len2; size++) { + svalue = (long)(op1[size]) - op2[size] - carry; + rop[size] = svalue; + carry = svalue < 0; + } + if (rop != op1) { + for (; size < len1; size++) { + svalue = op1[size] - carry; + rop[size] = svalue; + carry = svalue < 0; + } + } + else { + /* if rop == op1, than just adjust carry */ + for (; carry && size < len1; size++) { + svalue = op1[size] - carry; + rop[size] = svalue; + carry = svalue < 0; + } + size = len1; + } + + /* calculate result size */ + while (size > 1 && rop[size - 1] == 0) + --size; + + return (size); +} + +long +mp_lshift(BNS *rop, BNS *op, BNI len, long shift) +{ + long i, size; + BNI words, bits; /* how many word and bit shifts */ + + words = shift / BNSBITS; + bits = shift % BNSBITS; + size = len + words; + + if (bits) { + BNS hi, lo; + BNI carry; + int adj; + + for (i = 1, carry = CARRY >> 1; carry; i++, carry >>= 1) + if (op[len - 1] & carry) + break; + adj = (bits + (BNSBITS - i)) / BNSBITS; + size += adj; + + lo = hi = op[0]; + rop[words] = lo << bits; + for (i = 1; i < len; i++) { + hi = op[i]; + rop[words + i] = hi << bits | (lo >> (BNSBITS - bits)); + lo = hi; + } + if (adj) + rop[size - 1] = hi >> (BNSBITS - bits); + } + else + memmove(rop + size - len, op, sizeof(BNS) * len); + + if (words) + memset(rop, '\0', sizeof(BNS) * words); + + return (size); +} + +long +mp_rshift(BNS *rop, BNS *op, BNI len, long shift) +{ + int adj = 0; + long i, size; + BNI words, bits; /* how many word and bit shifts */ + + words = shift / BNSBITS; + bits = shift % BNSBITS; + size = len - words; + + if (bits) { + BNS hi, lo; + BNI carry; + + for (i = 0, carry = CARRY >> 1; carry; i++, carry >>= 1) + if (op[len - 1] & carry) + break; + adj = (bits + i) / BNSBITS; + if (size - adj == 0) { + rop[0] = 0; + + return (1); + } + + hi = lo = op[words + size - 1]; + rop[size - 1] = hi >> bits; + for (i = size - 2; i >= 0; i--) { + lo = op[words + i]; + rop[i] = (lo >> bits) | (hi << (BNSBITS - bits)); + hi = lo; + } + if (adj) + rop[0] |= lo << (BNSBITS - bits); + } + else + memmove(rop, op + len - size, size * sizeof(BNS)); + + return (size - adj); +} + + /* rop must be a pointer to len1 + len2 elements + * rop cannot be either op1 or op2 + * rop must be all zeros */ +long +mp_base_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2) +{ + long i, j; /* counters */ + BNI value; /* intermediate result */ + BNS carry; /* carry value */ + long size = len1 + len2; + + /* simple optimization: first pass does not need to deference rop[i+j] */ + if (op1[0]) { + value = (BNI)(op1[0]) * op2[0]; + rop[0] = value; + carry = (BNS)(value >> BNSBITS); + for (j = 1; j < len2; j++) { + value = (BNI)(op1[0]) * op2[j] + carry; + rop[j] = value; + carry = (BNS)(value >> BNSBITS); + } + rop[j] = carry; + } + + /* do the multiplication */ + for (i = 1; i < len1; i++) { + if (op1[i]) { + /* unrool loop initialization */ + value = (BNI)(op1[i]) * op2[0] + rop[i]; + rop[i] = value; + carry = (BNS)(value >> BNSBITS); + /* multiply */ + for (j = 1; j < len2; j++) { + value = (BNI)(op1[i]) * op2[j] + rop[i + j] + carry; + rop[i + j] = value; + carry = (BNS)(value >> BNSBITS); + } + rop[i + j] = carry; + } + } + + if (size > 1 && rop[size - 1] == 0) + --size; + + return (size); +} + + /* Karatsuba method + * t + ((a0 + a1) (b0 + b1) - t - u) x + ux² + * where t = a0b0 and u = a1b1 + * + * Karatsuba method reduces the number of multiplications. Example: + * Square a 40 length number + * instead of a plain 40*40 = 1600 multiplies/adds, it does: + * 20*20+20*20+20*20 = 1200 + * but since it is recursive, every 20*20=400 is reduced to + * 10*10+10*10+10*10=300 + * and so on. + * The multiplication by x and x² is a just a shift, as it is a + * power of two, and is implemented below by just writting at the + * correct offset */ +long +mp_karatsuba_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2) +{ + BNI x; /* shift count */ + BNI la0, la1, lb0, lb1; /* length of a0, a1, b0, and b1 */ + BNS *t; /* temporary memory for t product */ + BNS *u; /* temporary memory for u product */ + BNS *r; /* pointer to rop */ + long xlen, tlen, ulen; + + /* calculate value of x, that is 2^(BNSBITS*x) */ + if (len1 >= len2) + x = (len1 + 1) >> 1; + else + x = (len2 + 1) >> 1; + + /* calculate length of operands */ + la0 = x; + la1 = len1 - x; + lb0 = x; + lb1 = len2 - x; + + /* allocate buffer for t and (a0 + a1) */ + tlen = la0 + lb0; + t = mp_malloc(sizeof(BNS) * tlen); + + /* allocate buffer for u and (b0 + b1) */ + if (la1 + lb1 < lb0 + lb1 + 1) + ulen = lb0 + lb1 + 1; + else + ulen = la1 + lb1; + u = mp_malloc(sizeof(BNS) * ulen); + + /* calculate a0 + a1, store result in t */ + tlen = mp_add(t, op1, op1 + x, la0, la1); + + /* calculate b0 + b1, store result in u */ + ulen = mp_add(u, op2, op2 + x, lb0, lb1); + + /* store (a0 + a1) * (b0 + b1) in rop */ + + r = rop + x; /* multiplied by 2^(BNSBITS*x) */ + xlen = mp_mul(r, t, u, tlen, ulen); + + /* must zero t and u memory, this is required for mp_mul */ + + /* calculate t = a0 * b0 */ + tlen = la0 + lb0; + memset(t, '\0', sizeof(BNS) * tlen); + tlen = mp_mul(t, op1, op2, la0, lb0); + + /* calculate u = a1 * b1 */ + ulen = la1 + lb1; + memset(u, '\0', sizeof(BNS) * ulen); + ulen = mp_mul(u, op1 + x, op2 + x, la1, lb1); + + /* subtract t from partial result */ + xlen = mp_sub(r, r, t, xlen, tlen); + + /* subtract u form partial result */ + xlen = mp_sub(r, r, u, xlen, ulen); + + /* add ux^2 to partial result */ + + r = rop + (x << 1); /* multiplied by x^2 = 2^(BNSBITS*x*2) */ + xlen = len1 + len2; + xlen = mp_add(r, r, u, xlen, ulen); + + /* now add t to final result */ + xlen = mp_add(rop, rop, t, xlen, tlen); + + mp_free(t); + mp_free(u); + + if (xlen > 1 && rop[xlen - 1] == 0) + --xlen; + + return (xlen); +} + + /* Toom method (partially based on GMP documentation) + * Evaluation at k = [ 0 1/2 1 2 oo ] + * U(x) = (U2k + U1)k + U0 + * V(x) = (V2k + V1)k + V0 + * W(x) = U(x)V(x) + * + * Sample: + * 123 * 456 + * + * EVALUATION: + * U(0) = (1*0+2)*0+3 => 3 + * U(1) = 1+(2+3*2)*2 => 17 + * U(2) = 1+2+3 => 6 + * U(3) = (1*2+2)*2+3 => 11 + * U(4) = 1+(2+3*0)*0 => 1 + * + * V(0) = (4*0+5)*0+6 => 6 + * V(1) = 4+(5+6*2)*2 => 38 + * V(2) = 4+5+6 => 15 + * V(3) = (4*2+5)*2+6 => 32 + * V(4) = 4+(5+6*0)*0 => 4 + * + * U = [ 3 17 6 11 1 ] + * V = [ 6 38 15 32 4 ] + * W = [ 18 646 90 352 4 ] + * + * After that, we have: + * a = 18 (w0 already known) + * b = 16w0 + 8w1 + 4w2 + 2w3 + w4 + * c = w0 + w1 + w2 + w3 + w4 + * d = w0 + 2w1 + 4w2 + 8w3 + 16w4 + * e = 4 (w4 already known) + * + * INTERPOLATION: + * b = b -16a - e (354) + * c = c - a - e (68) + * d = d - a - 16e (270) + * + * w = (b + d) - 8c = (10w1+8w2+10w3) - (8w1+8w2+8w3) = 2w1+2w3 + * w = 2c - w (56) + * b = b/2 = 4w1+w+w3 + * b = b-c = 4w1+w+w3 - w1+w2+w3 = 3w1+w2 + * c = w/2 (w2 = 28) + * b = b-c = 3w1+c - c = 3w1 + * b = b/3 (w1 = 27) + * d = d/2 + * d = d-b-w = b+w+4w3 - b-w = 4w3 + * d = d/4 (w3 = 13) + * + * RESULT: + * w4*10^4 + w3*10³ + w2*10² + w1*10 + w0 + * 40000 + 13000 + 2800 + 270 + 18 + * 10 is the base where the calculation was done + * + * This sample uses small numbers, so it does not show the + * advantage of the method. But for example (in base 10), when squaring + * 123456789012345678901234567890 + * The normal method would do 30*30=900 multiplications + * Karatsuba method would do 15*15*3=675 multiplications + * Toom method would do 10*10*5=500 multiplications + * Toom method has a larger overhead if compared with Karatsuba method, + * due to evaluation and interpolation, so it should be used for larger + * numbers, so that the computation time of evaluation/interpolation + * would be smaller than the time spent using other methods. + * + * Note that Karatsuba method can be seen as a special case of + * Toom method, i.e: + * U1U0 * V1V0 + * with k = [ 0 1 oo ] + * U = [ U0 U1+U0 U1 ] + * V = [ V0 V1+V0 V1 ] + * W = [ U0*V0 (U1+U0)*(V1+V0) (U1+V1) ] + * + * w0 = U0*V0 + * w = (U1+U0)*(V1+V0) + * w2 = (U1*V1) + * + * w1 = w - w0 - w2 + * w2x² + w1x + w0 + * + * See Knuth's Seminumerical Algorithms for a sample implemention + * using 4 stacks and k = [ 0 1 2 3 ... ], based on the size of the + * input. + */ +long +mp_toom_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2) +{ + long size, xsize, i; + BNI value; /* used in division */ + BNS carry; + BNI x; /* shift count */ + BNI l1, l2; + BNI al, bl, cl, dl, el, Ul[3], Vl[3]; + BNS *a, *b, *c, *d, *e, *U[3], *V[3]; + + /* x is the base i.e. 2^(BNSBITS*x) */ + x = (len1 + len2 + 4) / 6; + l1 = len1 - (x << 1); /* length of remaining piece of op1 */ + l2 = len2 - (x << 1); /* length of remaining piece of op2 */ + + /* allocate memory for storing U and V */ + U[0] = mp_malloc(sizeof(BNS) * (x + 2)); + V[0] = mp_malloc(sizeof(BNS) * (x + 2)); + U[1] = mp_malloc(sizeof(BNS) * (x + 1)); + V[1] = mp_malloc(sizeof(BNS) * (x + 1)); + U[2] = mp_malloc(sizeof(BNS) * (x + 2)); + V[2] = mp_malloc(sizeof(BNS) * (x + 2)); + + /* EVALUATE U AND V */ + + /* Numbers are in the format U2x²+U1x+U0 and V2x²+V1x+V0 */ + + /* U[0] = U2+U1*2+U0*4 */ + + /* store U1*2 in U[1], this value is used twice */ + Ul[1] = mp_lshift(U[1], op1 + x, x, 1); + + /* store U0*4 in U[0] */ + Ul[0] = mp_lshift(U[0], op1, x, 2); + /* add U1*2 to U[0] */ + Ul[0] = mp_add(U[0], U[0], U[1], Ul[0], Ul[1]); + /* add U2 to U[0] */ + Ul[0] = mp_add(U[0], U[0], op1 + x + x, Ul[0], l1); + + /* U[2] = U2*4+U1*2+U0 */ + + /* store U2*4 in U[2] */ + Ul[2] = mp_lshift(U[2], op1 + x + x, l1, 2); + /* add U1*2 to U[2] */ + Ul[2] = mp_add(U[2], U[2], U[1], Ul[2], Ul[1]); + /* add U0 to U[2] */ + Ul[2] = mp_add(U[2], U[2], op1, Ul[2], x); + + /* U[1] = U2+U1+U0 */ + + Ul[1] = mp_add(U[1], op1, op1 + x, x, x); + Ul[1] = mp_add(U[1], U[1], op1 + x + x, Ul[1], l1); + + + /* Evaluate V[x], same code as U[x] */ + Vl[1] = mp_lshift(V[1], op2 + x, x, 1); + Vl[0] = mp_lshift(V[0], op2, x, 2); + Vl[0] = mp_add(V[0], V[0], V[1], Vl[0], Vl[1]); + Vl[0] = mp_add(V[0], V[0], op2 + x + x, Vl[0], l2); + Vl[2] = mp_lshift(V[2], op2 + x + x, l2, 2); + Vl[2] = mp_add(V[2], V[2], V[1], Vl[2], Vl[1]); + Vl[2] = mp_add(V[2], V[2], op2, Vl[2], x); + Vl[1] = mp_add(V[1], op2, op2 + x, x, x); + Vl[1] = mp_add(V[1], V[1], op2 + x + x, Vl[1], l2); + + + /* MULTIPLY U[] AND V[] */ + + /* calculate (U2+U1*2+U0*4) * (V2+V1*2+V0*4) */ + b = mp_calloc(1, sizeof(BNS) * (Ul[0] * Vl[0])); + bl = mp_mul(b, U[0], V[0], Ul[0], Vl[0]); + mp_free(U[0]); + mp_free(V[0]); + + /* calculate (U2+U1+U0) * (V2+V1+V0) */ + c = mp_calloc(1, sizeof(BNS) * (Ul[1] * Vl[1])); + cl = mp_mul(c, U[1], V[1], Ul[1], Vl[1]); + mp_free(U[1]); + mp_free(V[1]); + + /* calculate (U2*4+U1*2+U0) * (V2*4+V1*2+V0) */ + d = mp_calloc(1, sizeof(BNS) * (Ul[2] * Vl[2])); + dl = mp_mul(d, U[2], V[2], Ul[2], Vl[2]); + mp_free(U[2]); + mp_free(V[2]); + + /* calculate U0 * V0 */ + a = mp_calloc(1, sizeof(BNS) * (x + x)); + al = mp_mul(a, op1, op2, x, x); + + /* calculate U2 * V2 */ + e = mp_calloc(1, sizeof(BNS) * (l1 + l2)); + el = mp_mul(e, op1 + x + x, op2 + x + x, l1, l2); + + + /* INTERPOLATE COEFFICIENTS */ + + /* b = b - 16a - e */ + size = mp_lshift(rop, a, al, 4); + bl = mp_sub(b, b, rop, bl, size); + bl = mp_sub(b, b, e, bl, el); + + /* c = c - a - e*/ + cl = mp_sub(c, c, a, cl, al); + cl = mp_sub(c, c, e, cl, el); + + /* d = d - a - 16e */ + dl = mp_sub(d, d, a, dl, al); + size = mp_lshift(rop, e, el, 4); + dl = mp_sub(d, d, rop, dl, size); + + /* w = (b + d) - 8c */ + size = mp_add(rop, b, d, bl, dl); + xsize = mp_lshift(rop + size, c, cl, 3); /* rop has enough storage */ + size = mp_sub(rop, rop, rop + size, size, xsize); + + /* w = 2c - w*/ + xsize = mp_lshift(rop + size, c, cl, 1); + size = mp_sub(rop, rop + size, rop, xsize, size); + + /* b = b/2 */ + bl = mp_rshift(b, b, bl, 1); + + /* b = b - c */ + bl = mp_sub(b, b, c, bl, cl); + + /* c = w / 2 */ + cl = mp_rshift(c, rop, size, 1); + + /* b = b - c */ + bl = mp_sub(b, b, c, bl, cl); + + /* b = b/3 */ + /* maybe the most expensive calculation */ + i = bl - 1; + value = b[i]; + b[i] = value / 3; + for (--i; i >= 0; i--) { + carry = value % 3; + value = ((BNI)carry << BNSBITS) + b[i]; + b[i] = (BNS)(value / 3); + } + + /* d = d/2 */ + dl = mp_rshift(d, d, dl, 1); + + /* d = d - b - w */ + dl = mp_sub(d, d, b, dl, bl); + dl = mp_sub(d, d, rop, dl, size); + + /* d = d/4 */ + dl = mp_rshift(d, d, dl, 2); + + + /* STORE RESULT IN ROP */ + /* first clear memory used as temporary variable w and 8c */ + memset(rop, '\0', sizeof(BNS) * (len1 + len2)); + + i = x * 4; + xsize = (len1 + len2) - i; + size = mp_add(rop + i, rop + i, e, xsize, el) + i; + i = x * 3; + xsize = size - i; + size = mp_add(rop + i, rop + i, d, xsize, dl) + i; + i = x * 2; + xsize = size - i; + size = mp_add(rop + i, rop + i, c, xsize, cl) + i; + i = x; + xsize = size - i; + size = mp_add(rop + i, rop + i, b, xsize, bl) + i; + size = mp_add(rop, rop, a, size, al); + + mp_free(e); + mp_free(d); + mp_free(c); + mp_free(b); + mp_free(a); + + if (size > 1 && rop[size - 1] == 0) + --size; + + return (size); +} + +long +mp_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2) +{ + if (len1 < len2) + MP_SWAP(op1, op2, len1, len2); + + if (len1 < KARATSUBA || len2 < KARATSUBA) + return (mp_base_mul(rop, op1, op2, len1, len2)); + else if (len1 < TOOM && len2 < TOOM && len2 > ((len1 + 1) >> 1)) + return (mp_karatsuba_mul(rop, op1, op2, len1, len2)); + else if (len1 >= TOOM && len2 >= TOOM && (len2 + 2) / 3 == (len1 + 2) / 3) + return (mp_toom_mul(rop, op1, op2, len1, len2)); + else { + long xsize, psize, isize; + BNS *ptr; + + /* adjust index pointer and estimated size of result */ + isize = 0; + xsize = len1 + len2; + mp_mul(rop, op1, op2, len2, len2); + /* adjust pointers */ + len1 -= len2; + op1 += len2; + + /* allocate buffer for intermediate multiplications */ + if (len1 > len2) + ptr = mp_calloc(1, sizeof(BNS) * (len2 + len2)); + else + ptr = mp_calloc(1, sizeof(BNS) * (len1 + len2)); + + /* loop multiplying len2 size operands at a time */ + while (len1 >= len2) { + isize += len2; + psize = mp_mul(ptr, op1, op2, len2, len2); + mp_add(rop + isize, rop + isize, ptr, xsize - isize, psize); + len1 -= len2; + op1 += len2; + + /* multiplication routines require zeroed memory */ + memset(ptr, '\0', sizeof(BNS) * (MIN(len1, len2) + len2)); + } + + /* len1 was not a multiple of len2 */ + if (len1) { + isize += len2; + psize = mp_mul(ptr, op2, op1, len2, len1); + mp_add(rop + isize, rop + isize, ptr, xsize, psize); + } + + /* adjust result size */ + if (rop[xsize - 1] == 0) + --xsize; + + mp_free(ptr); + + return (xsize); + } +} diff --git a/lisp/mp/mp.h b/lisp/mp/mp.h new file mode 100644 index 0000000..88f1b24 --- /dev/null +++ b/lisp/mp/mp.h @@ -0,0 +1,435 @@ +/* + * Copyright (c) 2002 by The XFree86 Project, Inc. + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE XFREE86 PROJECT BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, + * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF + * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * Except as contained in this notice, the name of the XFree86 Project shall + * not be used in advertising or otherwise to promote the sale, use or other + * dealings in this Software without prior written authorization from the + * XFree86 Project. + * + * Author: Paulo César Pereira de Andrade + */ + +/* $XFree86: xc/programs/xedit/lisp/mp/mp.h,v 1.6 2003/01/12 03:55:51 tsi Exp $ */ + +#include <stdio.h> +#include <math.h> +#ifdef sun +#include <ieeefp.h> +#endif +#include <float.h> +#include <stdlib.h> +#include <limits.h> +#include <ctype.h> +#include <string.h> + +#ifndef __mp_h_ +#define __mp_h_ + +#ifdef __GNUC__ +#define INLINE __inline__ +#else +#define INLINE /**/ +#endif + +/* this normally is better for multiplication and also + * simplify addition loops putting the larger value first */ +#define MP_SWAP(op1, op2, len1, len2) { \ + BNS *top = op1; \ + BNI tlen = len1; \ + \ + op1 = op2; \ + len1 = len2; \ + op2 = top; \ + len2 = tlen; \ +} + +/* + * At least this length to use Karatsuba multiplication method + */ +#define KARATSUBA 32 + +/* + * At least this length to use Toom multiplication method + */ +#define TOOM 128 + +#if ULONG_MAX > 4294967295UL + /* sizeof(long) == 8 and sizeof(int) == 4 */ +# define BNI unsigned long +# define BNS unsigned int +# define MINSLONG 0x8000000000000000UL +# define CARRY 0x100000000 +# define LMASK 0xffffffff00000000UL +# define SMASK 0x00000000ffffffffUL +# define BNIBITS 64 +# define BNSBITS 32 +# ifndef LONG64 +# define LONG64 +# endif +#else + /* sizeof(long) == 4 and sizeof(short) == 2 */ +# define BNI unsigned long +# define BNS unsigned short +# define MINSLONG 0x80000000UL +# define CARRY 0x10000 +# define LMASK 0xffff0000UL +# define SMASK 0x0000ffffUL +# define BNIBITS 32 +# define BNSBITS 16 +#endif + +#ifdef MAX +#undef MAX +#endif +#define MAX(a, b) ((a) > (b) ? (a) : (b)) + +#ifdef MIN +#undef MIN +#endif +#define MIN(a, b) ((a) < (b) ? (a) : (b)) + +/* + * Types + */ +typedef struct _mpi { + unsigned int size : 31; + unsigned int sign : 1; + BNI alloc; + BNS *digs; /* LSF format */ +} mpi; + +typedef struct _mpr { + mpi num; + mpi den; +} mpr; + +typedef void *(*mp_malloc_fun)(size_t); +typedef void *(*mp_calloc_fun)(size_t, size_t); +typedef void *(*mp_realloc_fun)(void*, size_t); +typedef void (*mp_free_fun)(void*); + +/* + * Prototypes + */ +/* GENERIC FUNCTIONS */ + /* memory allocation wrappers */ +void *mp_malloc(size_t size); +void *mp_calloc(size_t nmemb, size_t size); +void *mp_realloc(void *pointer, size_t size); +void mp_free(void *pointer); +mp_malloc_fun mp_set_malloc(mp_malloc_fun); +mp_calloc_fun mp_set_calloc(mp_calloc_fun); +mp_realloc_fun mp_set_realloc(mp_realloc_fun); +mp_free_fun mp_set_free(mp_free_fun); + + /* adds op1 and op2, stores result in rop + * rop must pointer to at least len1 + len2 + 1 elements + * rop can be either op1 or op2 */ +long mp_add(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2); + + /* subtracts op2 from op1, stores result in rop + * rop must pointer to at least len1 + len2 elements + * op1 must be >= op2 + * rop can be either op1 or op2 */ +long mp_sub(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2); + + /* shift op to the left shift bits + * rop must have enough storage for result + * rop can be op */ +long mp_lshift(BNS *rop, BNS *op, BNI len, long shift); + + /* shift op to the right shift bits + * shift must be positive + * rop can be op */ +long mp_rshift(BNS *rop, BNS *op, BNI len, long shift); + + /* use simple generic multiplication method + * rop cannot be the same as op1 or op2 + * rop must be zeroed + * op1 can be op2 */ +long mp_base_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2); + + /* use Karatsuba method + * MIN(len1, len2) must be larger than (MAX(len1, len2) + 1) >> 1 + * MAX(len1, len2) should be at least 2 + * rop cannot be the same as op1 or op2 + * rop must be zeroed + * op1 can be op2 */ +long mp_karatsuba_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2); + + /* use Toom method + * len1 / 3 should be equal to len2 / 3 + * len1 / 3 should be at least 1 + * rop cannot be the same as op1 or op2 + * rop must be zeroed + * op1 can be op2 */ +long mp_toom_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2); + + /* chooses the available multiplication methods based on it's input + * rop must be a pointer to len1 + len2 elements + * rop cannot be the same as op1 or op2 + * rop must be zeroed + * op1 can be op2 */ +long mp_mul(BNS *rop, BNS *op1, BNS *op2, BNI len1, BNI len2); + +/* INTEGER FUNCTIONS */ + /* initialize op and set it to 0 */ +void mpi_init(mpi *op); + + /* clear memory associated to op */ +void mpi_clear(mpi *op); + + /* set rop to the value of op */ +void mpi_set(mpi *rop, mpi *op); + + /* set rop to the value of si */ +void mpi_seti(mpi *rop, long si); + + /* set rop to the floor(fabs(d)) */ +void mpi_setd(mpi *rop, double d); + + /* initialize rop to number representation in str in the given base. + * leading zeros are skipped. + * if sign present, it is processed. + * base must be in the range 2 to 36. */ +void mpi_setstr(mpi *rop, char *str, int base); + + /* adds two mp integers */ +void mpi_add(mpi *rop, mpi *op1, mpi *op2); + + /* adds op1 and op2 */ +void mpi_addi(mpi *rop, mpi *op1, long op2); + + /* subtracts two mp integers */ +void mpi_sub(mpi *rop, mpi *op1, mpi *op2); + + /* subtracts op2 from op1 */ +void mpi_subi(mpi *rop, mpi *op1, long op2); + + /* multiply two mp integers */ +void mpi_mul(mpi *rop, mpi *op1, mpi *op2); + + /* multiply op1 by op2 */ +void mpi_muli(mpi *rop, mpi *op1, long op2); + + /* divides num by den and sets rop to result */ +void mpi_div(mpi *rop, mpi *num, mpi *den); + + /* divides num by den and sets rop to the remainder */ +void mpi_rem(mpi *rop, mpi *num, mpi *den); + + /* divides num by den, sets quotient to qrop and remainder to rrop + * qrop is truncated towards zero. + * qrop and rrop are optional + * qrop and rrop cannot be the same variable */ +void mpi_divqr(mpi *qrop, mpi *rrop, mpi *num, mpi *den); + + /* divides num by then and stores result in rop */ +void mpi_divi(mpi *rop, mpi *num, long den); + + /* divides num by den and returns remainder */ +long mpi_remi(mpi *num, long den); + + /* divides num by den + * stores quotient in qrop and returns remainder */ +long mpi_divqri(mpi *qrop, mpi *num, long den); + + /* sets rop to num modulo den */ +void mpi_mod(mpi *rop, mpi *num, mpi *den); + + /* returns num modulo den */ +long mpi_modi(mpi *num, long den); + + /* sets rop to the greatest common divisor of num and den + * result is always positive */ +void mpi_gcd(mpi *rop, mpi *num, mpi *den); + + /* sets rop to the least common multiple of num and den + * result is always positive */ +void mpi_lcm(mpi *rop, mpi *num, mpi *den); + + /* sets rop to op raised to exp */ +void mpi_pow(mpi *rop, mpi *op, unsigned long exp); + + /* sets rop to the integer part of the nth root of op. + * returns 1 if result is exact, 0 otherwise */ +int mpi_root(mpi *rop, mpi *op, unsigned long nth); + + /* sets rop to the integer part of the square root of op. + * returns 1 if result is exact, 0 otherwise */ +int mpi_sqrt(mpi *rop, mpi *op); + + /* bit shift, left if shift positive, right if negative + * a fast way to multiply and divide by powers of two */ +void mpi_ash(mpi *rop, mpi *op, long shift); + + /* sets rop to op1 logand op2 */ +void mpi_and(mpi *rop, mpi *op1, mpi *op2); + + /* sets rop to op1 logior op2 */ +void mpi_ior(mpi *rop, mpi *op1, mpi *op2); + + /* sets rop to op1 logxor op2 */ +void mpi_xor(mpi *rop, mpi *op1, mpi *op2); + + /* sets rop to one's complement of op */ +void mpi_com(mpi *rop, mpi *op); + + /* sets rop to -op */ +void mpi_neg(mpi *rop, mpi *op); + + /* sets rop to the absolute value of op */ +void mpi_abs(mpi *rop, mpi *op); + + /* compares op1 and op2 + * returns >0 if op1 > op2, 0 if op1 = op2, and <0 if op1 < op2 */ +int mpi_cmp(mpi *op1, mpi *op2); + + /* mpi_cmp with a long integer operand */ +int mpi_cmpi(mpi *op1, long op2); + + /* compares absolute value of op1 and op2 + * returns >0 if abs(op1) > abs(op2), 0 if abs(op1) = abs(op2), + * and <0 if abs(op1) < abs(op2) */ +int mpi_cmpabs(mpi *op1, mpi *op2); + + /* mpi_cmpabs with a long integer operand */ +int mpi_cmpabsi(mpi *op1, long op2); + + /* returns 1 if op1 > 0, 0 if op1 = 0, and -1 if op1 < 0 */ +int mpi_sgn(mpi *op); + + /* fastly swaps contents of op1 and op2 */ +void mpi_swap(mpi *op1, mpi *op2); + + /* returns 1 if op fits in a signed long int, 0 otherwise */ +int mpi_fiti(mpi *op); + + /* converts mp integer to long int + * to know if the value will fit, call mpi_fiti */ +long mpi_geti(mpi *op); + + /* convert mp integer to double */ +double mpi_getd(mpi *op); + + /* returns exact number of characters to represent mp integer + * in given base, excluding sign and ending null character. + * base must be in the range 2 to 36 */ +unsigned long mpi_getsize(mpi *op, int base); + + /* returns pointer to string with representation of mp integer + * if str is not NULL, it must have enough space to store integer + * representation, if NULL a newly allocated string is returned. + * base must be in the range 2 to 36 */ +char *mpi_getstr(char *str, mpi *op, int base); + +/* RATIO FUNCTIONS */ +#define mpr_num(op) (&((op)->num)) +#define mpr_den(op) (&((op)->den)) + + /* initialize op and set it to 0/1 */ +void mpr_init(mpr *op); + + /* clear memory associated to op */ +void mpr_clear(mpr *op); + + /* set rop to the value of op */ +void mpr_set(mpr *rop, mpr *op); + + /* set rop to num/den */ +void mpr_seti(mpr *rop, long num, long den); + + /* set rop to the value of d */ +void mpr_setd(mpr *rop, double d); + + /* initialize rop to number representation in str in the given base. + * leading zeros are skipped. + * if sign present, it is processed. + * base must be in the range 2 to 36. */ +void mpr_setstr(mpr *rop, char *str, int base); + + /* remove common factors of op */ +void mpr_canonicalize(mpr *op); + + /* adds two mp rationals */ +void mpr_add(mpr *rop, mpr *op1, mpr *op2); + + /* adds op1 and op2 */ +void mpr_addi(mpr *rop, mpr *op1, long op2); + + /* subtracts two mp rationals */ +void mpr_sub(mpr *rop, mpr *op1, mpr *op2); + + /* subtracts op2 from op1 */ +void mpr_subi(mpr *rop, mpr *op1, long op2); + + /* multiply two mp rationals */ +void mpr_mul(mpr *rop, mpr *op1, mpr *op2); + + /* multiply op1 by op2 */ +void mpr_muli(mpr *rop, mpr *op1, long op2); + + /* divide two mp rationals */ +void mpr_div(mpr *rop, mpr *op1, mpr *op2); + + /* divides op1 by op2 */ +void mpr_divi(mpr *rop, mpr *op1, long op2); + + /* sets rop to 1/op */ +void mpr_inv(mpr *rop, mpr *op); + + /* sets rop to -op */ +void mpr_neg(mpr *rop, mpr *op); + + /* sets rop to the absolute value of op */ +void mpr_abs(mpr *rop, mpr *op); + + /* compares op1 and op2 + * returns >0 if op1 > op2, 0 if op1 = op2, and <0 if op1 < op2 */ +int mpr_cmp(mpr *op1, mpr *op2); + + /* mpr_cmp with a long integer operand */ +int mpr_cmpi(mpr *op1, long op2); + + /* compares absolute value of op1 and op2 + * returns >0 if abs(op1) > abs(op2), 0 if abs(op1) = abs(op2), + * and <0 if abs(op1) < abs(op2) */ +int mpr_cmpabs(mpr *op1, mpr *op2); + + /* mpr_cmpabs with a long integer operand */ +int mpr_cmpabsi(mpr *op1, long op2); + + /* fastly swaps contents of op1 and op2 */ +void mpr_swap(mpr *op1, mpr *op2); + + /* returns 1 if op fits in a signed long int, 0 otherwise */ +int mpr_fiti(mpr *op); + + /* convert mp rational to double */ +double mpr_getd(mpr *op); + + /* returns pointer to string with representation of mp rational + * if str is not NULL, it must have enough space to store rational + * representation, if NULL a newly allocated string is returned. + * base must be in the range 2 to 36 */ +char *mpr_getstr(char *str, mpr *op, int base); + +#endif /* __mp_h_ */ diff --git a/lisp/mp/mpi.c b/lisp/mp/mpi.c new file mode 100644 index 0000000..506dc7e --- /dev/null +++ b/lisp/mp/mpi.c @@ -0,0 +1,1656 @@ +/* + * Copyright (c) 2002 by The XFree86 Project, Inc. + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE XFREE86 PROJECT BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, + * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF + * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * Except as contained in this notice, the name of the XFree86 Project shall + * not be used in advertising or otherwise to promote the sale, use or other + * dealings in this Software without prior written authorization from the + * XFree86 Project. + * + * Author: Paulo César Pereira de Andrade + */ + +/* $XFree86: xc/programs/xedit/lisp/mp/mpi.c,v 1.12 2002/11/20 07:44:43 paulo Exp $ */ + +#include "mp.h" + +/* + * Prototypes + */ + /* do the hard work of mpi_add and mpi_sub */ +static void mpi_addsub(mpi *rop, mpi *op1, mpi *op2, int sub); + + /* logical functions implementation */ +static INLINE BNS mpi_logic(BNS op1, BNS op2, BNS op); +static void mpi_log(mpi *rop, mpi *op1, mpi *op2, BNS op); + + /* internal mpi_seti, whithout memory allocation */ +static void _mpi_seti(mpi *rop, long si); + +/* + * Initialization + */ +static BNS onedig[1] = { 1 }; +static mpi mpone = { 1, 1, 0, (BNS*)&onedig }; + +/* + * Implementation + */ +void +mpi_init(mpi *op) +{ + op->sign = 0; + op->size = op->alloc = 1; + op->digs = mp_malloc(sizeof(BNS)); + op->digs[0] = 0; +} + +void +mpi_clear(mpi *op) +{ + op->sign = 0; + op->size = op->alloc = 0; + mp_free(op->digs); +} + +void +mpi_set(mpi *rop, mpi *op) +{ + if (rop != op) { + if (rop->alloc < op->size) { + rop->digs = mp_realloc(rop->digs, sizeof(BNS) * op->size); + rop->alloc = op->size; + } + rop->size = op->size; + memcpy(rop->digs, op->digs, sizeof(BNS) * op->size); + rop->sign = op->sign; + } +} + +void +mpi_seti(mpi *rop, long si) +{ + unsigned long ui; + int sign = si < 0; + int size; + + if (si == MINSLONG) { + ui = MINSLONG; + size = 2; + } + else { + if (sign) + ui = -si; + else + ui = si; + if (ui < CARRY) + size = 1; + else + size = 2; + } + + if (rop->alloc < size) { + rop->digs = mp_realloc(rop->digs, sizeof(BNS) * size); + rop->alloc = size; + } + rop->size = size; + + /* store data in small mp integer */ + rop->digs[0] = (BNS)ui; + if (size > 1) + rop->digs[1] = (BNS)(ui >> BNSBITS); + rop->size = size; + + /* adjust result sign */ + rop->sign = sign; +} + +static void +_mpi_seti(mpi *rop, long si) +{ + unsigned long ui; + int sign = si < 0; + int size; + + if (si == MINSLONG) { + ui = MINSLONG; + size = 2; + } + else { + if (sign) + ui = -si; + else + ui = si; + if (ui < CARRY) + size = 1; + else + size = 2; + } + + rop->digs[0] = (BNS)ui; + if (size > 1) + rop->digs[1] = (BNS)(ui >> BNSBITS); + rop->size = size; + + rop->sign = sign; +} + +void +mpi_setd(mpi *rop, double d) +{ + long i; + double mantissa; + int shift, exponent; + BNI size; + + if (isnan(d)) + d = 0.0; + else if (!finite(d)) + d = copysign(1.0, d) * DBL_MAX; + + /* check if number is larger than 1 */ + if (fabs(d) < 1.0) { + rop->digs[0] = 0; + rop->size = 1; + rop->sign = d < 0.0; + + return; + } + + mantissa = frexp(d, &exponent); + if (mantissa < 0) + mantissa = -mantissa; + + size = (exponent + (BNSBITS - 1)) / BNSBITS; + shift = BNSBITS - (exponent & (BNSBITS - 1)); + + /* adjust amount of memory */ + if (rop->alloc < size) { + rop->digs = mp_realloc(rop->digs, sizeof(BNS) * size); + rop->alloc = size; + } + rop->size = size; + + /* adjust the exponent */ + if (shift < BNSBITS) + mantissa = ldexp(mantissa, -shift); + + /* convert double */ + for (i = size - 1; i >= 0 && mantissa != 0.0; i--) { + mantissa = ldexp(mantissa, BNSBITS); + rop->digs[i] = (BNS)mantissa; + mantissa -= rop->digs[i]; + } + for (; i >= 0; i--) + rop->digs[i] = 0; + + /* normalize */ + if (size > 1 && rop->digs[size - 1] == 0) + --rop->size; + + rop->sign = d < 0.0; +} + +/* how many BNS in the given base, log(base) / log(CARRY) */ +#ifdef LONG64 +static double str_bases[37] = { + 0.0000000000000000, 0.0000000000000000, 0.0312500000000000, + 0.0495300781475362, 0.0625000000000000, 0.0725602529652301, + 0.0807800781475362, 0.0877298413143002, 0.0937500000000000, + 0.0990601562950723, 0.1038102529652301, 0.1081072380824156, + 0.1120300781475362, 0.1156387411919092, 0.1189798413143002, + 0.1220903311127662, 0.1250000000000000, 0.1277332137890731, + 0.1303101562950723, 0.1327477347951120, 0.1350602529652300, + 0.1372599194618363, 0.1393572380824156, 0.1413613111267817, + 0.1432800781475362, 0.1451205059304602, 0.1468887411919092, + 0.1485902344426084, 0.1502298413143002, 0.1518119060977367, + 0.1533403311127662, 0.1548186346995899, 0.1562500000000000, + 0.1576373162299517, 0.1589832137890731, 0.1602900942795302, + 0.1615601562950723, +}; +#else +static double str_bases[37] = { + 0.0000000000000000, 0.0000000000000000, 0.0625000000000000, + 0.0990601562950723, 0.1250000000000000, 0.1451205059304602, + 0.1615601562950723, 0.1754596826286003, 0.1875000000000000, + 0.1981203125901446, 0.2076205059304602, 0.2162144761648311, + 0.2240601562950723, 0.2312774823838183, 0.2379596826286003, + 0.2441806622255325, 0.2500000000000000, 0.2554664275781462, + 0.2606203125901445, 0.2654954695902241, 0.2701205059304602, + 0.2745198389236725, 0.2787144761648311, 0.2827226222535633, + 0.2865601562950723, 0.2902410118609203, 0.2937774823838183, + 0.2971804688852168, 0.3004596826286003, 0.3036238121954733, + 0.3066806622255324, 0.3096372693991797, 0.3125000000000000, + 0.3152746324599034, 0.3179664275781462, 0.3205801885590604, + 0.3231203125901446, +}; +#endif + +void +mpi_setstr(mpi *rop, char *str, int base) +{ + long i; /* counter */ + int sign; /* result sign */ + BNI carry; /* carry value */ + BNI value; /* temporary value */ + BNI size; /* size of result */ + char *ptr; /* end of valid input */ + + /* initialization */ + sign = 0; + carry = 0; + + /* skip leading spaces */ + while (isspace(*str)) + ++str; + + /* check if sign supplied */ + if (*str == '-') { + sign = 1; + ++str; + } + else if (*str == '+') + ++str; + + /* skip leading zeros */ + while (*str == '0') + ++str; + + ptr = str; + while (*ptr) { + if (*ptr >= '0' && *ptr <= '9') { + if (*ptr - '0' >= base) + break; + } + else if (*ptr >= 'A' && *ptr <= 'Z') { + if (*ptr - 'A' + 10 >= base) + break; + } + else if (*ptr >= 'a' && *ptr <= 'z') { + if (*ptr - 'a' + 10 >= base) + break; + } + else + break; + ++ptr; + } + + /* resulting size */ + size = (ptr - str) * str_bases[base] + 1; + + /* make sure rop has enough storage */ + if (rop->alloc < size) { + rop->digs = mp_realloc(rop->digs, size * sizeof(BNS)); + rop->alloc = size; + } + rop->size = size; + + /* initialize rop to zero */ + memset(rop->digs, '\0', size * sizeof(BNS)); + + /* set result sign */ + rop->sign = sign; + + /* convert string */ + for (; str < ptr; str++) { + value = *str; + if (islower(value)) + value = toupper(value); + value = value > '9' ? value - 'A' + 10 : value - '0'; + value += rop->digs[0] * base; + carry = value >> BNSBITS; + rop->digs[0] = (BNS)value; + for (i = 1; i < size; i++) { + value = (BNI)rop->digs[i] * base + carry; + carry = value >> BNSBITS; + rop->digs[i] = (BNS)value; + } + } + + /* normalize */ + if (rop->size > 1 && rop->digs[rop->size - 1] == 0) + --rop->size; +} + +void +mpi_add(mpi *rop, mpi *op1, mpi *op2) +{ + mpi_addsub(rop, op1, op2, 0); +} + +void +mpi_addi(mpi *rop, mpi *op1, long op2) +{ + BNS digs[2]; + mpi op; + + op.digs = (BNS*)digs; + _mpi_seti(&op, op2); + + mpi_addsub(rop, op1, &op, 0); +} + +void +mpi_sub(mpi *rop, mpi *op1, mpi *op2) +{ + mpi_addsub(rop, op1, op2, 1); +} + +void +mpi_subi(mpi *rop, mpi *op1, long op2) +{ + BNS digs[2]; + mpi op; + + op.digs = (BNS*)digs; + _mpi_seti(&op, op2); + + mpi_addsub(rop, op1, &op, 1); +} + +static void +mpi_addsub(mpi *rop, mpi *op1, mpi *op2, int sub) +{ + long xlen; /* maximum result size */ + + if (sub ^ (op1->sign == op2->sign)) { + /* plus one for possible carry */ + xlen = MAX(op1->size, op2->size) + 1; + if (rop->alloc < xlen) { + rop->digs = mp_realloc(rop->digs, sizeof(BNS) * xlen); + rop->alloc = xlen; + } + rop->size = mp_add(rop->digs, op1->digs, op2->digs, + op1->size, op2->size); + rop->sign = op1->sign; + } + else { + long cmp; /* check for larger operator */ + + cmp = mpi_cmpabs(op1, op2); + if (cmp == 0) { + rop->digs[0] = 0; + rop->size = 1; + rop->sign = 0; + } + else { + xlen = MAX(op1->size, op2->size); + if (rop->alloc < xlen) { + rop->digs = mp_realloc(rop->digs, sizeof(BNS) * xlen); + rop->alloc = xlen; + } + if (cmp > 0) { + rop->size = mp_sub(rop->digs, op1->digs, op2->digs, + op1->size, op2->size); + rop->sign = op1->sign; + } + else { + rop->size = mp_sub(rop->digs, op2->digs, op1->digs, + op2->size, op1->size); + rop->sign = sub ^ op2->sign; + } + } + } +} + +void +mpi_mul(mpi *rop, mpi *op1, mpi *op2) +{ + int sign; /* sign flag */ + BNS *digs; /* result data */ + long xsize; /* result size */ + + /* get result sign */ + sign = op1->sign ^ op2->sign; + + /* check for special cases */ + if (op1->size == 1) { + if (*op1->digs == 0) { + /* multiply by 0 */ + mpi_seti(rop, 0); + return; + } + else if (*op1->digs == 1) { + /* multiply by +-1 */ + if (rop->alloc < op2->size) { + rop->digs = mp_realloc(rop->digs, sizeof(BNS) * op2->size); + rop->alloc = op2->size; + } + rop->size = op2->size; + memmove(rop->digs, op2->digs, sizeof(BNS) * op2->size); + rop->sign = op2->size > 1 || *op2->digs ? sign : 0; + + return; + } + } + else if (op2->size == 1) { + if (*op2->digs == 0) { + /* multiply by 0 */ + mpi_seti(rop, 0); + return; + } + else if (*op2->digs == 1) { + /* multiply by +-1 */ + if (rop->alloc < op1->size) { + rop->digs = mp_realloc(rop->digs, sizeof(BNS) * op1->size); + rop->alloc = op1->size; + } + rop->size = op1->size; + memmove(rop->digs, op1->digs, sizeof(BNS) * op1->size); + rop->sign = op1->size > 1 || *op1->digs ? sign : 0; + + return; + } + } + + /* allocate result data and set it to zero */ + xsize = op1->size + op2->size; + if (rop->digs == op1->digs || rop->digs == op2->digs) + /* rop is also an operand */ + digs = mp_calloc(1, sizeof(BNS) * xsize); + else { + if (rop->alloc < xsize) { + rop->digs = mp_realloc(rop->digs, sizeof(BNS) * xsize); + rop->alloc = xsize; + } + digs = rop->digs; + memset(digs, '\0', sizeof(BNS) * xsize); + } + + /* multiply operands */ + xsize = mp_mul(digs, op1->digs, op2->digs, op1->size, op2->size); + + /* store result in rop */ + if (digs != rop->digs) { + /* if rop was an operand, free old data */ + mp_free(rop->digs); + rop->digs = digs; + } + rop->size = xsize; + + /* set result sign */ + rop->sign = sign; +} + +void +mpi_muli(mpi *rop, mpi *op1, long op2) +{ + BNS digs[2]; + mpi op; + + op.digs = (BNS*)digs; + _mpi_seti(&op, op2); + + mpi_mul(rop, op1, &op); +} + +void +mpi_div(mpi *rop, mpi *num, mpi *den) +{ + mpi_divqr(rop, NULL, num, den); +} + +void +mpi_rem(mpi *rop, mpi *num, mpi *den) +{ + mpi_divqr(NULL, rop, num, den); +} + +/* + * Could/should be changed to not allocate qdigs if qrop is NULL + * Performance wouldn't suffer too much with a test on every loop iteration. + */ +void +mpi_divqr(mpi *qrop, mpi *rrop, mpi *num, mpi *den) +{ + long i, j; /* counters */ + int qsign; /* sign of quotient */ + int rsign; /* sign of remainder */ + BNI qsize; /* size of quotient */ + BNI rsize; /* size of remainder */ + BNS qest; /* estimative of quotient value */ + BNS *qdigs, *rdigs; /* work copy or result */ + BNS *ndigs, *ddigs; /* work copy or divisor and dividend */ + BNI value; /* temporary result */ + long svalue; /* signed temporary result (2's complement) */ + BNS carry, scarry, denorm; /* carry and normalization */ + BNI dpos, npos; /* offsets in data */ + + /* get signs */ + rsign = num->sign; + qsign = rsign ^ den->sign; + + /* check for special case */ + if (num->size < den->size) { + /* quotient is zero and remainder is numerator */ + if (rrop && rrop->digs != num->digs) { + if (rrop->alloc < num->size) { + rrop->digs = mp_realloc(rrop->digs, sizeof(BNS) * num->size); + rrop->alloc = num->size; + } + rrop->size = num->size; + memcpy(rrop->digs, num->digs, sizeof(BNS) * num->size); + rrop->sign = rsign; + } + if (qrop) + mpi_seti(qrop, 0); + + return; + } + + /* estimate result sizes */ + rsize = den->size; + qsize = num->size - den->size + 1; + + /* offsets */ + npos = num->size - 1; + dpos = den->size - 1; + + /* allocate space for quotient and remainder */ + if (qrop == NULL || qrop->digs == num->digs || qrop->digs == den->digs) + qdigs = mp_calloc(1, sizeof(BNS) * qsize); + else { + if (qrop->alloc < qsize) { + qrop->digs = mp_realloc(qrop->digs, sizeof(BNS) * qsize); + qrop->alloc = qsize; + } + memset(qrop->digs, '\0', sizeof(BNS) * qsize); + qdigs = qrop->digs; + } + if (rrop) { + if (rrop->digs == num->digs || rrop->digs == den->digs) + rdigs = mp_calloc(1, sizeof(BNS) * rsize); + else { + if (rrop->alloc < rsize) { + rrop->digs = mp_realloc(rrop->digs, sizeof(BNS) * rsize); + rrop->alloc = rsize; + } + memset(rrop->digs, '\0', sizeof(BNS) * rsize); + rdigs = rrop->digs; + } + } + else + rdigs = NULL; /* fix gcc warning */ + + /* special case, only one word in divisor */ + if (dpos == 0) { + for (carry = 0, i = npos; i >= 0; i--) { + value = ((BNI)carry << BNSBITS) + num->digs[i]; + qdigs[i] = (BNS)(value / den->digs[0]); + carry = (BNS)(value % den->digs[0]); + } + if (rrop) + rdigs[0] = carry; + + goto mpi_divqr_done; + } + + /* make work copy of numerator */ + ndigs = mp_malloc(sizeof(BNS) * (num->size + 1)); + /* allocate one extra word an update offset */ + memcpy(ndigs, num->digs, sizeof(BNS) * num->size); + ndigs[num->size] = 0; + ++npos; + + /* normalize */ + denorm = (BNS)((BNI)CARRY / ((BNI)(den->digs[dpos]) + 1)); + + if (denorm > 1) { + /* i <= num->size because ndigs has an extra word */ + for (carry = 0, i = 0; i <= num->size; i++) { + value = ndigs[i] * (BNI)denorm + carry; + ndigs[i] = (BNS)value; + carry = (BNS)(value >> BNSBITS); + } + /* make work copy of denominator */ + ddigs = mp_malloc(sizeof(BNS) * den->size); + memcpy(ddigs, den->digs, sizeof(BNS) * den->size); + for (carry = 0, i = 0; i < den->size; i++) { + value = ddigs[i] * (BNI)denorm + carry; + ddigs[i] = (BNS)value; + carry = (BNS)(value >> BNSBITS); + } + } + else + /* only allocate copy of denominator if going to change it */ + ddigs = den->digs; + + /* divide mp integers */ + for (j = qsize - 1; j >= 0; j--, npos--) { + /* estimate quotient */ + if (ndigs[npos] == ddigs[dpos]) + qest = (BNS)SMASK; + else + qest = (BNS)((((BNI)(ndigs[npos]) << 16) + ndigs[npos - 1]) / + ddigs[dpos]); + + while ((value = ((BNI)(ndigs[npos]) << 16) + ndigs[npos - 1] - + qest * (BNI)(ddigs[dpos])) < CARRY && + ddigs[dpos - 1] * (BNI)qest > + (value << BNSBITS) + ndigs[npos - 2]) + --qest; + + /* multiply and subtract */ + carry = scarry = 0; + for (i = 0; i < den->size; i++) { + value = qest * (BNI)ddigs[i] + carry; + carry = (BNS)(value >> BNSBITS); + svalue = (long)ndigs[npos - dpos + i - 1] - (long)(value & SMASK) - + (long)scarry; + ndigs[npos - dpos + i - 1] = (BNS)svalue; + scarry = svalue < 0; + } + + svalue = (long)ndigs[npos] - (long)(carry & SMASK) - (long)scarry; + ndigs[npos] = (BNS)svalue; + + if (svalue & LMASK) { + /* quotient too big */ + --qest; + carry = 0; + for (i = 0; i < den->size; i++) { + value = ndigs[npos - dpos + i - 1] + (BNI)carry + (BNI)ddigs[i]; + ndigs[npos - dpos + i - 1] = (BNS)value; + carry = (BNS)(value >> BNSBITS); + } + ndigs[npos] += carry; + } + + qdigs[j] = qest; + } + + /* calculate remainder */ + if (rrop) { + for (carry = 0, j = dpos; j >= 0; j--) { + value = ((BNI)carry << BNSBITS) + ndigs[j]; + rdigs[j] = (BNS)(value / denorm); + carry = (BNS)(value % denorm); + } + } + + mp_free(ndigs); + if (ddigs != den->digs) + mp_free(ddigs); + +mpi_divqr_done: + if (rrop) { + if (rrop->digs != rdigs) + mp_free(rrop->digs); + /* normalize remainder */ + for (i = rsize - 1; i >= 0; i--) + if (rdigs[i] != 0) + break; + if (i != rsize - 1) { + if (i < 0) { + rsign = 0; + rsize = 1; + } + else + rsize = i + 1; + } + rrop->digs = rdigs; + rrop->sign = rsign; + rrop->size = rsize; + } + + /* normalize quotient */ + if (qrop) { + if (qrop->digs != qdigs) + mp_free(qrop->digs); + for (i = qsize - 1; i >= 0; i--) + if (qdigs[i] != 0) + break; + if (i != qsize - 1) { + if (i < 0) { + qsign = 0; + qsize = 1; + } + else + qsize = i + 1; + } + qrop->digs = qdigs; + qrop->sign = qsign; + qrop->size = qsize; + } + else + mp_free(qdigs); +} + +long +mpi_divqri(mpi *qrop, mpi *num, long den) +{ + BNS ddigs[2]; + mpi dop, rrop; + long remainder; + + dop.digs = (BNS*)ddigs; + _mpi_seti(&dop, den); + + memset(&rrop, '\0', sizeof(mpi)); + mpi_init(&rrop); + mpi_divqr(qrop, &rrop, num, &dop); + remainder = rrop.digs[0]; + if (rrop.size > 1) + remainder |= (BNI)(rrop.digs[1]) << BNSBITS; + if (rrop.sign) + remainder = -remainder; + mpi_clear(&rrop); + + return (remainder); +} + +void +mpi_divi(mpi *rop, mpi *num, long den) +{ + BNS ddigs[2]; + mpi dop; + + dop.digs = (BNS*)ddigs; + _mpi_seti(&dop, den); + + mpi_divqr(rop, NULL, num, &dop); +} + +long +mpi_remi(mpi *num, long den) +{ + return (mpi_divqri(NULL, num, den)); +} + +void +mpi_mod(mpi *rop, mpi *num, mpi *den) +{ + mpi_rem(rop, num, den); + if (num->sign ^ den->sign) + mpi_add(rop, rop, den); +} + +long +mpi_modi(mpi *num, long den) +{ + long remainder; + + remainder = mpi_remi(num, den); + if (num->sign ^ (den < 0)) + remainder += den; + + return (remainder); +} + +void +mpi_gcd(mpi *rop, mpi *num, mpi *den) +{ + long cmp; + mpi rem; + + /* check if result already given */ + cmp = mpi_cmpabs(num, den); + + /* check if num is equal to den or if num is zero */ + if (cmp == 0 || (num->size == 1 && num->digs[0] == 0)) { + mpi_set(rop, den); + rop->sign = 0; + return; + } + /* check if den is not zero */ + if (den->size == 1 && den->digs[0] == 0) { + mpi_set(rop, num); + rop->sign = 0; + return; + } + + /* don't call mpi_init, relies on realloc(0, size) == malloc(size) */ + memset(&rem, '\0', sizeof(mpi)); + + /* if num larger than den */ + if (cmp > 0) { + mpi_rem(&rem, num, den); + if (rem.size == 1 && rem.digs[0] == 0) { + /* exact division */ + mpi_set(rop, den); + rop->sign = 0; + mpi_clear(&rem); + return; + } + mpi_set(rop, den); + } + else { + mpi_rem(&rem, den, num); + if (rem.size == 1 && rem.digs[0] == 0) { + /* exact division */ + mpi_set(rop, num); + rop->sign = 0; + mpi_clear(&rem); + return; + } + mpi_set(rop, num); + } + + /* loop using positive values */ + rop->sign = rem.sign = 0; + + /* cannot optimize this inverting rem/rop assignment earlier + * because rop mais be an operand */ + mpi_swap(rop, &rem); + + /* Euclides algorithm */ + for (;;) { + mpi_rem(&rem, &rem, rop); + if (rem.size == 1 && rem.digs[0] == 0) + break; + mpi_swap(rop, &rem); + } + mpi_clear(&rem); +} + +void +mpi_lcm(mpi *rop, mpi *num, mpi *den) +{ + mpi gcd; + + /* check for zero operand */ + if ((num->size == 1 && num->digs[0] == 0) || + (den->size == 1 && den->digs[0] == 0)) { + rop->digs[0] = 0; + rop->sign = 0; + return; + } + + /* don't call mpi_init, relies on realloc(0, size) == malloc(size) */ + memset(&gcd, '\0', sizeof(mpi)); + + mpi_gcd(&gcd, num, den); + mpi_div(&gcd, den, &gcd); + mpi_mul(rop, &gcd, num); + rop->sign = 0; + + mpi_clear(&gcd); +} + +void +mpi_pow(mpi *rop, mpi *op, unsigned long exp) +{ + mpi zop, top; + + if (exp == 2) { + mpi_mul(rop, op, op); + return; + } + /* check for op**0 */ + else if (exp == 0) { + rop->digs[0] = 1; + rop->size = 1; + rop->sign = 0; + return; + } + else if (exp == 1) { + mpi_set(rop, op); + return; + } + else if (op->size == 1) { + if (op->digs[0] == 0) { + mpi_seti(rop, 0); + return; + } + else if (op->digs[0] == 1) { + mpi_seti(rop, op->sign && (exp & 1) ? -1 : 1); + return; + } + } + + memset(&zop, '\0', sizeof(mpi)); + memset(&top, '\0', sizeof(mpi)); + mpi_set(&zop, op); + mpi_set(&top, op); + for (--exp; exp; exp >>= 1) { + if (exp & 1) + mpi_mul(&zop, &top, &zop); + mpi_mul(&top, &top, &top); + } + + mpi_clear(&top); + rop->sign = zop.sign; + mp_free(rop->digs); + rop->digs = zop.digs; + rop->size = zop.size; +} + +/* Find integer root of given number using the iteration + * x{n+1} = ((K-1) * x{n} + N / x{n}^(K-1)) / K + */ +int +mpi_root(mpi *rop, mpi *op, unsigned long nth) +{ + long bits, cmp; + int exact; + int sign; + mpi *r, t, temp, quot, old, rem; + + sign = op->sign; + + /* divide by zero op**1/0 */ + if (nth == 0) { + int one = 1, zero = 0; + one = one / zero; + } + /* result is complex */ + if (sign && !(nth & 1)) { + int one = 1, zero = 0; + one = one / zero; + } + + /* special case op**1/1 = op */ + if (nth == 1) { + mpi_set(rop, op); + return (1); + } + + bits = mpi_getsize(op, 2) - 2; + + if (bits < 0 || bits / nth == 0) { + /* integral root is surely less than 2 */ + exact = op->size == 1 && (op->digs[0] == 1 || op->digs[0] == 0); + mpi_seti(rop, sign ? -1 : op->digs[0] == 0 ? 0 : 1); + + return (exact == 1); + } + + /* initialize */ + if (rop != op) + r = rop; + else { + r = &t; + memset(r, '\0', sizeof(mpi)); + } + memset(&temp, '\0', sizeof(mpi)); + memset(", '\0', sizeof(mpi)); + memset(&old, '\0', sizeof(mpi)); + memset(&rem, '\0', sizeof(mpi)); + + if (sign) + r->sign = 0; + + /* root aproximation */ + mpi_ash(r, op, -(bits - (bits / nth))); + + for (;;) { + mpi_pow(&temp, r, nth - 1); + mpi_divqr(", &rem, op, &temp); + cmp = mpi_cmpabs(r, "); + if (cmp == 0) { + exact = mpi_cmpi(&rem, 0) == 0; + break; + } + else if (cmp < 0) { + if (mpi_cmpabs(r, &old) == 0) { + exact = 0; + break; + } + mpi_set(&old, r); + } + mpi_muli(&temp, r, nth - 1); + mpi_add(", ", &temp); + mpi_divi(r, ", nth); + } + + mpi_clear(&temp); + mpi_clear("); + mpi_clear(&old); + mpi_clear(&rem); + if (r != rop) { + mpi_set(rop, r); + mpi_clear(r); + } + rop->sign = sign; + + return (exact); +} + +/* + * Find square root using the iteration: + * x{n+1} = (x{n}+N/x{n})/2 + */ +int +mpi_sqrt(mpi *rop, mpi *op) +{ + long bits, cmp; + int exact; + mpi *r, t, quot, rem, old; + + /* result is complex */ + if (op->sign) { + int one = 1, zero = 0; + one = one / zero; + } + + bits = mpi_getsize(op, 2) - 2; + + if (bits < 2) { + /* integral root is surely less than 2 */ + exact = op->size == 1 && (op->digs[0] == 1 || op->digs[0] == 0); + mpi_seti(rop, op->digs[0] == 0 ? 0 : 1); + + return (exact == 1); + } + + /* initialize */ + if (rop != op) + r = rop; + else { + r = &t; + memset(r, '\0', sizeof(mpi)); + } + memset(", '\0', sizeof(mpi)); + memset(&rem, '\0', sizeof(mpi)); + memset(&old, '\0', sizeof(mpi)); + + /* root aproximation */ + mpi_ash(r, op, -(bits - (bits / 2))); + + for (;;) { + if (mpi_cmpabs(r, &old) == 0) { + exact = 0; + break; + } + mpi_divqr(", &rem, op, r); + cmp = mpi_cmpabs(", r); + if (cmp == 0) { + exact = mpi_cmpi(&rem, 0) == 0; + break; + } + else if (cmp > 0 && rem.size == 1 && rem.digs[0] == 0) + mpi_subi(", ", 1); + mpi_set(&old, r); + mpi_add(r, r, "); + mpi_ash(r, r, -1); + } + mpi_clear("); + mpi_clear(&rem); + mpi_clear(&old); + if (r != rop) { + mpi_set(rop, r); + mpi_clear(r); + } + + return (exact); +} + +void +mpi_ash(mpi *rop, mpi *op, long shift) +{ + long i; /* counter */ + long xsize; /* maximum result size */ + BNS *digs; + + /* check for 0 shift, multiply/divide by 1 */ + if (shift == 0) { + if (rop != op) { + if (rop->alloc < op->size) { + rop->digs = mp_realloc(rop->digs, sizeof(BNS) * op->size); + rop->alloc = op->size; + } + rop->size = op->size; + memcpy(rop->digs, op->digs, sizeof(BNS) * op->size); + } + + return; + } + else if (op->size == 1 && op->digs[0] == 0) { + rop->sign = 0; + rop->size = 1; + rop->digs[0] = 0; + + return; + } + + /* check shift and initialize */ + if (shift > 0) + xsize = op->size + (shift / BNSBITS) + 1; + else { + xsize = op->size - ((-shift) / BNSBITS); + if (xsize <= 0) { + rop->size = 1; + rop->sign = op->sign; + rop->digs[0] = op->sign ? 1 : 0; + + return; + } + } + + /* allocate/adjust memory for result */ + if (rop == op) + digs = mp_malloc(sizeof(BNS) * xsize); + else { + if (rop->alloc < xsize) { + rop->digs = mp_realloc(rop->digs, sizeof(BNS) * xsize); + rop->alloc = xsize; + } + digs = rop->digs; + } + + /* left shift, multiply by power of two */ + if (shift > 0) + rop->size = mp_lshift(digs, op->digs, op->size, shift); + /* right shift, divide by power of two */ + else { + long carry = 0; + + if (op->sign) { + BNI words, bits; + + words = -shift / BNSBITS; + bits = -shift % BNSBITS; + for (i = 0; i < words; i++) + carry |= op->digs[xsize + i]; + if (!carry) { + for (i = 0; i < bits; i++) + if (op->digs[op->size - xsize] & (1 << i)) { + carry = 1; + break; + } + } + } + rop->size = mp_rshift(digs, op->digs, op->size, -shift); + + if (carry) + /* emulates two's complement subtracting 1 from the result */ + rop->size = mp_add(digs, digs, mpone.digs, rop->size, 1); + } + + if (rop->digs != digs) { + mp_free(rop->digs); + rop->alloc = rop->size; + rop->digs = digs; + } + rop->sign = op->sign; +} + +static INLINE BNS +mpi_logic(BNS op1, BNS op2, BNS op) +{ + switch (op) { + case '&': + return (op1 & op2); + case '|': + return (op1 | op2); + case '^': + return (op1 ^ op2); + } + + return (SMASK); +} + +static void +mpi_log(mpi *rop, mpi *op1, mpi *op2, BNS op) +{ + long i; /* counter */ + long c, c1, c2; /* carry */ + BNS *digs, *digs1, *digs2; /* pointers to mp data */ + BNI size, size1, size2; + BNS sign, sign1, sign2; + BNS n, n1, n2; /* logical operands */ + BNI sum; + + /* initialize */ + size1 = op1->size; + size2 = op2->size; + + sign1 = op1->sign ? SMASK : 0; + sign2 = op2->sign ? SMASK : 0; + + sign = mpi_logic(sign1, sign2, op); + + size = MAX(size1, size2); + if (sign) + ++size; + if (rop->alloc < size) { + rop->digs = mp_realloc(rop->digs, sizeof(BNS) * size); + rop->alloc = size; + } + + digs = rop->digs; + digs1 = op1->digs; + digs2 = op2->digs; + + c = c1 = c2 = 1; + + /* apply logical operation */ + for (i = 0; i < size; i++) { + if (i >= size1) + n1 = sign1; + else if (sign1) { + sum = (BNI)(BNS)(~digs1[i]) + c1; + c1 = (long)(sum >> BNSBITS); + n1 = (BNS)sum; + } + else + n1 = digs1[i]; + + if (i >= size2) + n2 = sign2; + else if (sign2) { + sum = (BNI)(BNS)(~digs2[i]) + c2; + c2 = (long)(sum >> BNSBITS); + n2 = (BNS)sum; + } + else + n2 = digs2[i]; + + n = mpi_logic(n1, n2, op); + if (sign) { + sum = (BNI)(BNS)(~n) + c; + c = (long)(sum >> BNSBITS); + digs[i] = (BNS)sum; + } + else + digs[i] = n; + } + + /* normalize */ + for (i = size - 1; i >= 0; i--) + if (digs[i] != 0) + break; + if (i != size - 1) { + if (i < 0) { + sign = 0; + size = 1; + } + else + size = i + 1; + } + + rop->sign = sign != 0; + rop->size = size; +} + +void +mpi_and(mpi *rop, mpi *op1, mpi *op2) +{ + mpi_log(rop, op1, op2, '&'); +} + +void +mpi_ior(mpi *rop, mpi *op1, mpi *op2) +{ + mpi_log(rop, op1, op2, '|'); +} + +void +mpi_xor(mpi *rop, mpi *op1, mpi *op2) +{ + mpi_log(rop, op1, op2, '^'); +} + +void +mpi_com(mpi *rop, mpi *op) +{ + static BNS digs[1] = { 1 }; + static mpi one = { 0, 1, 1, (BNS*)&digs }; + + mpi_log(rop, rop, &one, '^'); +} + +void +mpi_neg(mpi *rop, mpi *op) +{ + int sign = op->sign ^ 1; + + if (rop->digs != op->digs) { + if (rop->alloc < op->size) { + rop->digs = mp_realloc(rop->digs, sizeof(BNS) * rop->size); + rop->alloc = op->size; + } + rop->size = op->size; + memcpy(rop->digs, op->digs, sizeof(BNS) * rop->size); + } + + rop->sign = sign; +} + +void +mpi_abs(mpi *rop, mpi *op) +{ + if (rop->digs != op->digs) { + if (rop->alloc < op->size) { + rop->digs = mp_realloc(rop->digs, sizeof(BNS) * rop->size); + rop->alloc = op->size; + } + rop->size = op->size; + memcpy(rop->digs, op->digs, sizeof(BNS) * rop->size); + } + + rop->sign = 0; +} + +int +mpi_cmp(mpi *op1, mpi *op2) +{ + if (op1->sign ^ op2->sign) + return (op1->sign ? -1 : 1); + + if (op1->size == op2->size) { + long i, cmp = 0; + + for (i = op1->size - 1; i >= 0; i--) + if ((cmp = (long)op1->digs[i] - (long)op2->digs[i]) != 0) + break; + + return (cmp == 0 ? 0 : (cmp < 0) ^ op1->sign ? -1 : 1); + } + + return ((op1->size < op2->size) ^ op1->sign ? -1 : 1); +} + +int +mpi_cmpi(mpi *op1, long op2) +{ + long cmp; + + if (op1->size > 2) + return (op1->sign ? -1 : 1); + + cmp = op1->digs[0]; + if (op1->size == 2) { + cmp |= (long)op1->digs[1] << BNSBITS; + if (cmp == MINSLONG) + return (op2 == cmp && op1->sign ? 0 : op1->sign ? -1 : 1); + } + if (op1->sign) + cmp = -cmp; + + return (cmp - op2); +} + +int +mpi_cmpabs(mpi *op1, mpi *op2) +{ + if (op1->size == op2->size) { + long i, cmp = 0; + + for (i = op1->size - 1; i >= 0; i--) + if ((cmp = (long)op1->digs[i] - (long)op2->digs[i]) != 0) + break; + + return (cmp); + } + + return ((op1->size < op2->size) ? -1 : 1); +} + +int +mpi_cmpabsi(mpi *op1, long op2) +{ + unsigned long cmp; + + if (op1->size > 2) + return (1); + + cmp = op1->digs[0]; + if (op1->size == 2) + cmp |= (unsigned long)op1->digs[1] << BNSBITS; + + return (cmp > op2 ? 1 : cmp == op2 ? 0 : -1); +} + +int +mpi_sgn(mpi *op) +{ + return (op->sign ? -1 : op->size > 1 || op->digs[0] ? 1 : 0); +} + +void +mpi_swap(mpi *op1, mpi *op2) +{ + if (op1 != op2) { + mpi swap; + + memcpy(&swap, op1, sizeof(mpi)); + memcpy(op1, op2, sizeof(mpi)); + memcpy(op2, &swap, sizeof(mpi)); + } +} + +int +mpi_fiti(mpi *op) +{ + if (op->size == 1) + return (1); + else if (op->size == 2) { + unsigned long value = ((BNI)(op->digs[1]) << BNSBITS) | op->digs[0]; + + if (value & MINSLONG) + return (op->sign && value == MINSLONG) ? 1 : 0; + + return (1); + } + + return (0); +} + +long +mpi_geti(mpi *op) +{ + long value; + + value = op->digs[0]; + if (op->size > 1) + value |= (BNI)(op->digs[1]) << BNSBITS; + + return (op->sign && value != MINSLONG ? -value : value); +} + +double +mpi_getd(mpi *op) +{ + long i, len; + double d = 0.0; + int exponent; + +#define FLOATDIGS sizeof(double) / sizeof(BNS) + + switch (op->size) { + case 2: + d = (BNI)(op->digs[1]) << BNSBITS; + case 1: + d += op->digs[0]; + return (op->sign ? -d : d); + default: + break; + } + + for (i = 0, len = op->size; len > 0 && i < FLOATDIGS; i++) + d = ldexp(d, BNSBITS) + op->digs[--len]; + d = frexp(d, &exponent); + if (len > 0) + exponent += len * BNSBITS; + + if (d == 0.0) + return (d); + + d = ldexp(d, exponent); + + return (op->sign ? -d : d); +} + +/* how many digits in a given base, floor(log(CARRY) / log(base)) */ +#ifdef LONG64 +static char dig_bases[37] = { + 0, 0, 32, 20, 16, 13, 12, 11, 10, 10, 9, 9, 8, 8, 8, 8, + 8, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6, 6, + 6, 6, 6, 6, 6, +}; +#else +static char dig_bases[37] = { + 0, 0, 16, 10, 8, 6, 6, 5, 5, 5, 4, 4, 4, 4, 4, 4, + 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, +}; +#endif + +/* how many digits per bit in a given base, log(2) / log(base) */ +static double bit_bases[37] = { + 0.0000000000000000, 0.0000000000000000, 1.0000000000000000, + 0.6309297535714575, 0.5000000000000000, 0.4306765580733931, + 0.3868528072345416, 0.3562071871080222, 0.3333333333333334, + 0.3154648767857287, 0.3010299956639811, 0.2890648263178878, + 0.2789429456511298, 0.2702381544273197, 0.2626495350371936, + 0.2559580248098155, 0.2500000000000000, 0.2446505421182260, + 0.2398124665681315, 0.2354089133666382, 0.2313782131597592, + 0.2276702486969530, 0.2242438242175754, 0.2210647294575037, + 0.2181042919855316, 0.2153382790366965, 0.2127460535533632, + 0.2103099178571525, 0.2080145976765095, 0.2058468324604344, + 0.2037950470905062, 0.2018490865820999, 0.2000000000000000, + 0.1982398631705605, 0.1965616322328226, 0.1949590218937863, + 0.1934264036172708, +}; + +/* normalization base for string conversion, pow(base, dig_bases[base]) & ~CARRY */ +#ifdef LONG64 +static BNS big_bases[37] = { + 0x00000001, 0x00000001, 0x00000000, 0xCFD41B91, 0x00000000, 0x48C27395, + 0x81BF1000, 0x75DB9C97, 0x40000000, 0xCFD41B91, 0x3B9ACA00, 0x8C8B6D2B, + 0x19A10000, 0x309F1021, 0x57F6C100, 0x98C29B81, 0x00000000, 0x18754571, + 0x247DBC80, 0x3547667B, 0x4C4B4000, 0x6B5A6E1D, 0x94ACE180, 0xCAF18367, + 0x0B640000, 0x0E8D4A51, 0x1269AE40, 0x17179149, 0x1CB91000, 0x23744899, + 0x2B73A840, 0x34E63B41, 0x40000000, 0x4CFA3CC1, 0x5C13D840, 0x6D91B519, + 0x81BF1000, +}; +#else +static BNS big_bases[37] = { + 0x0001, 0x0001, 0x0000, 0xE6A9, 0x0000, 0x3D09, 0xB640, 0x41A7, 0x8000, + 0xE6A9, 0x2710, 0x3931, 0x5100, 0x6F91, 0x9610, 0xC5C1, 0x0000, 0x1331, + 0x16C8, 0x1ACB, 0x1F40, 0x242D, 0x2998, 0x2F87, 0x3600, 0x3D09, 0x44A8, + 0x4CE3, 0x55C0, 0x5F45, 0x6978, 0x745F, 0x8000, 0x8C61, 0x9988, 0xA77B, + 0xb640, +}; +#endif + +unsigned long +mpi_getsize(mpi *op, int base) +{ + unsigned long value, bits; + + value = op->digs[op->size - 1]; + + /* count leading bits */ + if (value) { + unsigned long count, carry; + + for (count = 0, carry = CARRY >> 1; carry; count++, carry >>= 1) + if (value & carry) + break; + + bits = BNSBITS - count; + } + else + bits = 0; + + return ((bits + (op->size - 1) * BNSBITS) * bit_bases[base] + 1); +} + +char * +mpi_getstr(char *str, mpi *op, int base) +{ + long i; /* counter */ + BNS *digs, *xdigs; /* copy of op data */ + BNI size; /* size of op */ + BNI digits; /* digits per word in given base */ + BNI bigbase; /* big base of given base */ + BNI strsize; /* size of resulting string */ + char *cp; /* pointer in str for conversion */ + + /* initialize */ + size = op->size; + strsize = mpi_getsize(op, base) + op->sign + 1; + + if (str == NULL) + str = mp_malloc(strsize); + + /* check for zero */ + if (size == 1 && op->digs[0] == 0) { + str[0] = '0'; + str[1] = '\0'; + + return (str); + } + + digits = dig_bases[base]; + bigbase = big_bases[base]; + + cp = str + strsize; + *--cp = '\0'; + + /* make copy of op data and adjust digs */ + xdigs = mp_malloc(size * sizeof(BNS)); + memcpy(xdigs, op->digs, size * sizeof(unsigned short)); + digs = xdigs + size - 1; + + /* convert to string */ + for (;;) { + long count = -1; + BNI value; + BNS quotient, remainder = 0; + + /* if power of two base */ + if ((base & (base - 1)) == 0) { + for (i = 0; i < size; i++) { + quotient = remainder; + remainder = digs[-i]; + digs[-i] = quotient; + if (count < 0 && quotient) + count = i; + } + } + else { + for (i = 0; i < size; i++) { + value = digs[-i] + ((BNI)remainder << BNSBITS); + quotient = (BNS)(value / bigbase); + remainder = (BNS)(value % bigbase); + digs[-i] = quotient; + if (count < 0 && quotient) + count = i; + } + } + quotient = remainder; + for (i = 0; i < digits; i++) { + if (quotient == 0 && count < 0) + break; + remainder = quotient % base; + quotient /= base; + *--cp = remainder < 10 ? remainder + '0' : remainder - 10 + 'A'; + } + if (count < 0) + break; + digs -= count; + size -= count; + } + + /* adjust sign */ + if (op->sign) + *--cp = '-'; + + /* remove any extra characters */ + if (cp > str) + strcpy(str, cp); + + mp_free(xdigs); + + return (str); +} diff --git a/lisp/mp/mpr.c b/lisp/mp/mpr.c new file mode 100644 index 0000000..8b26fe0 --- /dev/null +++ b/lisp/mp/mpr.c @@ -0,0 +1,436 @@ +/* + * Copyright (c) 2002 by The XFree86 Project, Inc. + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE XFREE86 PROJECT BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, + * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF + * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + * + * Except as contained in this notice, the name of the XFree86 Project shall + * not be used in advertising or otherwise to promote the sale, use or other + * dealings in this Software without prior written authorization from the + * XFree86 Project. + * + * Author: Paulo César Pereira de Andrade + */ + +/* $XFree86: xc/programs/xedit/lisp/mp/mpr.c,v 1.2 2002/11/08 08:01:00 paulo Exp $ */ + +#include "mp.h" + +/* + * TODO: + * Implement a fast gcd and divexact for integers, so that this code + * could be changed to do intermediary calculations faster using smaller + * numbers. + */ + +/* + * Prototypes + */ + /* do the hard work of mpr_add and mpr_sub */ +static void mpr_addsub(mpr *rop, mpr *op1, mpr *op2, int sub); + + /* do the hard work of mpr_cmp and mpr_cmpabs */ +static int mpr_docmp(mpr *op1, mpr *op2, int sign); + +/* + * Implementation + */ +void +mpr_init(mpr *op) +{ + op->num.digs = mp_malloc(sizeof(BNS)); + op->num.sign = 0; + op->num.size = op->num.alloc = 1; + op->num.digs[0] = 0; + + op->den.digs = mp_malloc(sizeof(BNS)); + op->den.sign = 0; + op->den.size = op->den.alloc = 1; + op->den.digs[0] = 1; +} + +void +mpr_clear(mpr *op) +{ + op->num.sign = 0; + op->num.size = op->num.alloc = 0; + mp_free(op->num.digs); + + op->den.sign = 0; + op->den.size = op->den.alloc = 0; + mp_free(op->den.digs); +} + +void +mpr_set(mpr *rop, mpr *op) +{ + if (rop != op) { + mpi_set(mpr_num(rop), mpr_num(op)); + mpi_set(mpr_den(rop), mpr_den(op)); + } +} + +void +mpr_seti(mpr *rop, long num, long den) +{ + mpi_seti(mpr_num(rop), num); + mpi_seti(mpr_den(rop), den); +} + +void +mpr_setd(mpr *rop, double d) +{ + double val, num; + int e, sign; + + /* initialize */ + if (d < 0) { + sign = 1; + val = -d; + } + else { + sign = 0; + val = d; + } + + val = frexp(val, &e); + while (modf(val, &num) != 0.0 && val <= DBL_MAX / 2.0) { + --e; + val *= 2.0; + } + if (e >= 0) { + mpi_setd(mpr_num(rop), d); + mpi_seti(mpr_den(rop), 1); + } + else { + mpi_setd(mpr_num(rop), sign ? -num : num); + mpi_setd(mpr_den(rop), ldexp(1.0, -e)); + } +} + +void +mpr_setstr(mpr *rop, char *str, int base) +{ + char *slash = strchr(str, '/'); + + mpi_setstr(mpr_num(rop), str, base); + if (slash != NULL) + mpi_setstr(mpr_den(rop), slash + 1, base); + else + mpi_seti(mpr_den(rop), 1); +} + +void +mpr_canonicalize(mpr *op) +{ + mpi gcd; + + memset(&gcd, '\0', sizeof(mpi)); + + mpi_gcd(&gcd, mpr_num(op), mpr_den(op)); + if (mpi_cmpabsi(&gcd, 1)) { + mpi_div(mpr_num(op), mpr_num(op), &gcd); + mpi_div(mpr_den(op), mpr_den(op), &gcd); + } + + if (op->den.sign) { + op->num.sign = !op->num.sign; + op->den.sign = 0; + } + + mpi_clear(&gcd); +} + +void +mpr_add(mpr *rop, mpr *op1, mpr *op2) +{ + mpr_addsub(rop, op1, op2, 0); +} + +void +mpr_addi(mpr *rop, mpr *op1, long op2) +{ + mpi prod; + + memset(&prod, '\0', sizeof(mpi)); + + mpi_muli(&prod, mpr_den(op1), op2); + mpi_add(mpr_num(rop), mpr_num(op1), &prod); + mpi_clear(&prod); +} + +void +mpr_sub(mpr *rop, mpr *op1, mpr *op2) +{ + mpr_addsub(rop, op1, op2, 1); +} + +void +mpr_subi(mpr *rop, mpr *op1, long op2) +{ + mpi prod; + + memset(&prod, '\0', sizeof(mpi)); + + mpi_muli(&prod, mpr_den(op1), op2); + mpi_sub(mpr_num(rop), mpr_num(op1), &prod); + mpi_clear(&prod); +} + +static void +mpr_addsub(mpr *rop, mpr *op1, mpr *op2, int sub) +{ + mpi prod1, prod2; + + memset(&prod1, '\0', sizeof(mpi)); + memset(&prod2, '\0', sizeof(mpi)); + + mpi_mul(&prod1, mpr_num(op1), mpr_den(op2)); + mpi_mul(&prod2, mpr_num(op2), mpr_den(op1)); + + if (sub) + mpi_sub(mpr_num(rop), &prod1, &prod2); + else + mpi_add(mpr_num(rop), &prod1, &prod2); + + mpi_clear(&prod1); + mpi_clear(&prod2); + + mpi_mul(mpr_den(rop), mpr_den(op1), mpr_den(op2)); +} + +void +mpr_mul(mpr *rop, mpr *op1, mpr *op2) +{ + /* check if temporary storage is required */ + if (op1 == op2 && rop == op1) { + mpi prod; + + memset(&prod, '\0', sizeof(mpi)); + + mpi_mul(&prod, mpr_num(op1), mpr_num(op2)); + mpi_mul(mpr_den(rop), mpr_den(op1), mpr_den(op2)); + mpi_set(mpr_num(rop), &prod); + + mpi_clear(&prod); + } + else { + mpi_mul(mpr_num(rop), mpr_num(op1), mpr_num(op2)); + mpi_mul(mpr_den(rop), mpr_den(op1), mpr_den(op2)); + } +} + +void +mpr_muli(mpr *rop, mpr *op1, long op2) +{ + mpi_muli(mpr_num(rop), mpr_num(op1), op2); +} + +void +mpr_div(mpr *rop, mpr *op1, mpr *op2) +{ + /* check if temporary storage is required */ + if (op1 == op2 && rop == op1) { + mpi prod; + + memset(&prod, '\0', sizeof(mpi)); + + mpi_mul(&prod, mpr_num(op1), mpr_den(op2)); + mpi_mul(mpr_den(rop), mpr_num(op2), mpr_den(op1)); + mpi_set(mpr_num(rop), &prod); + + mpi_clear(&prod); + } + else { + mpi_mul(mpr_num(rop), mpr_num(op1), mpr_den(op2)); + mpi_mul(mpr_den(rop), mpr_num(op2), mpr_den(op1)); + } +} + +void +mpr_divi(mpr *rop, mpr *op1, long op2) +{ + mpi_muli(mpr_den(rop), mpr_den(op1), op2); +} + +void +mpr_inv(mpr *rop, mpr *op) +{ + if (rop == op) + mpi_swap(mpr_num(op), mpr_den(op)); + else { + mpi_set(mpr_num(rop), mpr_den(op)); + mpi_set(mpr_den(rop), mpr_num(op)); + } +} + +void +mpr_neg(mpr *rop, mpr *op) +{ + mpi_neg(mpr_num(rop), mpr_num(op)); + mpi_set(mpr_den(rop), mpr_den(op)); +} + +void +mpr_abs(mpr *rop, mpr *op) +{ + if (mpr_num(op)->sign) + mpi_neg(mpr_num(rop), mpr_num(op)); + else + mpi_set(mpr_num(rop), mpr_num(op)); + + /* op may not be canonicalized */ + if (mpr_den(op)->sign) + mpi_neg(mpr_den(rop), mpr_den(op)); + else + mpi_set(mpr_den(rop), mpr_den(op)); +} + +int +mpr_cmp(mpr *op1, mpr *op2) +{ + return (mpr_docmp(op1, op2, 1)); +} + +int +mpr_cmpi(mpr *op1, long op2) +{ + int cmp; + mpr rat; + + mpr_init(&rat); + mpi_seti(mpr_num(&rat), op2); + cmp = mpr_docmp(op1, &rat, 1); + mpr_clear(&rat); + + return (cmp); +} + +int +mpr_cmpabs(mpr *op1, mpr *op2) +{ + return (mpr_docmp(op1, op2, 0)); +} + +int +mpr_cmpabsi(mpr *op1, long op2) +{ + int cmp; + mpr rat; + + mpr_init(&rat); + mpi_seti(mpr_num(&rat), op2); + cmp = mpr_docmp(op1, &rat, 0); + mpr_clear(&rat); + + return (cmp); +} + +static int +mpr_docmp(mpr *op1, mpr *op2, int sign) +{ + int cmp, neg; + mpi prod1, prod2; + + neg = 0; + if (sign) { + /* if op1 is negative */ + if (mpr_num(op1)->sign ^ mpr_den(op1)->sign) { + /* if op2 is positive */ + if (!(mpr_num(op2)->sign ^ mpr_den(op2)->sign)) + return (-1); + else + neg = 1; + } + /* if op2 is negative */ + else if (mpr_num(op2)->sign ^ mpr_den(op2)->sign) + return (1); + /* else same sign */ + } + + /* if denominators are equal, compare numerators */ + if (mpi_cmpabs(mpr_den(op1), mpr_den(op2)) == 0) { + cmp = mpi_cmpabs(mpr_num(op1), mpr_num(op2)); + if (cmp == 0) + return (0); + if (sign && neg) + return (cmp < 0 ? 1 : -1); + return (cmp); + } + + memset(&prod1, '\0', sizeof(mpi)); + memset(&prod2, '\0', sizeof(mpi)); + + /* "divide" op1 by op2 + * if result is smaller than 1, op1 is smaller than op2 */ + mpi_mul(&prod1, mpr_num(op1), mpr_den(op2)); + mpi_mul(&prod2, mpr_num(op2), mpr_den(op1)); + + cmp = mpi_cmpabs(&prod1, &prod2); + + mpi_clear(&prod1); + mpi_clear(&prod2); + + if (sign && neg) + return (cmp < 0 ? 1 : -1); + return (cmp); +} + +void +mpr_swap(mpr *op1, mpr *op2) +{ + if (op1 != op2) { + mpr swap; + + memcpy(&swap, op1, sizeof(mpr)); + memcpy(op1, op2, sizeof(mpr)); + memcpy(op2, &swap, sizeof(mpr)); + } +} + +int +mpr_fiti(mpr *op) +{ + return (mpi_fiti(mpr_num(op)) && mpi_fiti(mpr_den(op))); +} + +double +mpr_getd(mpr *op) +{ + return (mpi_getd(mpr_num(op)) / mpi_getd(mpr_den(op))); +} + +char * +mpr_getstr(char *str, mpr *op, int base) +{ + int len; + + if (str == NULL) { + len = mpi_getsize(mpr_num(op), base) + mpr_num(op)->sign + 1 + + mpi_getsize(mpr_den(op), base) + mpr_den(op)->sign + 1; + + str = mp_malloc(len); + } + + (void)mpi_getstr(str, mpr_num(op), base); + len = strlen(str); + str[len] = '/'; + (void)mpi_getstr(str + len + 1, mpr_den(op), base); + + return (str); +} |