/* * 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" #ifdef __UNIXOS2__ # define finite(x) isfinite(x) #endif /* * 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); }