diff options
Diffstat (limited to 'lisp/mp/mpi.c')
| -rw-r--r-- | lisp/mp/mpi.c | 1656 |
1 files changed, 1656 insertions, 0 deletions
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); +} |