dropbear: libtommath/mtest/mpi.c comparison

comparison libtommath/mtest/mpi.c @ 1436:60fc6476e044

Update to libtommath v1.0

author	Matt Johnston <matt@ucc.asn.au>
date	Sat, 24 Jun 2017 22:37:14 +0800
parents	5ff8218bcee9
children

comparison

equal deleted inserted replaced

-:f849a5ca2efc
+:60fc6476e044
 by Michael J. Fromberger <[email protected]>
 Copyright (C) 1998 Michael J. Fromberger, All Rights Reserved
 Arbitrary precision integer arithmetic library
-$Id: mpi.c,v 1.2 2005/05/05 14:38:47 tom Exp $
+$Id$
 */
 #include "mpi.h"
 #include <stdlib.h>
 #include <string.h>
 #define DIAG(T,V) {fprintf(stderr,T);mp_print(V,stderr);fputc('\n',stderr);}
 #else
 #define DIAG(T,V)
 #endif
 /*
 If MP_LOGTAB is not defined, use the math library to compute the
 logarithms on the fly.  Otherwise, use the static table below.
 Pick which works best for your system.
 */
 #if MP_LOGTAB
 /* {{{ s_logv_2[] - log table for 2 in various bases */
 /*
 A table of the logs of 2 for various bases (the 0 and 1 entries of
 this table are meaningless and should not be referenced).
 This table is used to compute output lengths for the mp_toradix()
 function.  Since a number n in radix r takes up about log_r(n)
 digits, we estimate the output size by taking the least integer
 greater than log_r(n), where:
 log_r(n) = log_2(n) * log_r(2)
 This table, therefore, is a table of log_r(2) for 2 <= r <= 36,
 which are the output bases supported.
 */
 #include "logtab.h"
 /* }}} */
 };
 /* Value to digit maps for radix conversion   */
 /* s_dmap_1 - standard digits and letters */
 static const char *s_dmap_1 =
 "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";
 #if 0
 /* s_dmap_2 - base64 ordering for digits  */
 static const char *s_dmap_2 =
 /* }}} */
 /* {{{ Static function declarations */
 /*
 If MP_MACRO is false, these will be defined as actual functions;
 otherwise, suitable macro definitions will be used.  This works
 around the fact that ANSI C89 doesn't support an 'inline' keyword
 (although I hear C9x will ... about bloody time).  At present, the
 macro definitions are identical to the function bodies, but they'll
 }
 return MP_OKAY;
 CLEANUP:
 while(--pos >= 0)
 mp_clear(&mp[pos]);
 return res;
 } /* end mp_init_array() */
 usual
 */
 if(ALLOC(to) >= USED(from)) {
 s_mp_setz(DIGITS(to) + USED(from), ALLOC(to) - USED(from));
 s_mp_copy(DIGITS(from), DIGITS(to), USED(from));
 } else {
 if((tmp = s_mp_alloc(USED(from), sizeof(mp_digit))) == NULL)
 	return MP_MEM;
 s_mp_copy(DIGITS(from), tmp, USED(from));
 void   mp_clear_array(mp_int mp[], int count)
 {
 ARGCHK(mp != NULL && count > 0, MP_BADARG);
 while(--count >= 0)
 mp_clear(&mp[count]);
 } /* end mp_clear_array() */
 /* }}} */
 /* {{{ mp_zero(mp) */
 /*
 mp_zero(mp)
 Set mp to zero.  Does not change the allocated size of the structure,
 and therefore cannot fail (except on a bad argument, which we ignore)
 */
 void   mp_zero(mp_int *mp)
 for(ix = sizeof(long) - 1; ix >= 0; ix--) {
 if((res = s_mp_mul_2d(mp, CHAR_BIT)) != MP_OKAY)
 return res;
 res = s_mp_add_d(mp,
 		     (mp_digit)((v >> (ix * CHAR_BIT)) & UCHAR_MAX));
 if(res != MP_OKAY)
 return res;
 }
 ARGCHK(a != NULL && b != NULL, MP_BADARG);
 if((res = mp_copy(a, b)) != MP_OKAY)
 return res;
 if(s_mp_cmp_d(b, 0) == MP_EQ)
 SIGN(b) = MP_ZPOS;
 else
 SIGN(b) = (SIGN(b) == MP_NEG) ? MP_ZPOS : MP_NEG;
 return MP_OKAY;
 } /* end mp_neg() */
 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
 if(SIGN(a) == SIGN(b)) { /* same sign:  add values, keep sign */
 /* Commutativity of addition lets us do this in either order,
 so we avoid having to use a temporary even if the result
 is supposed to replace the output
 */
 if(c == b) {
 if((res = s_mp_add(c, a)) != MP_OKAY)
 	return res;
 } else {
 if(c != a && (res = mp_copy(a, c)) != MP_OKAY)
 	return res;
 if((res = s_mp_add(c, b)) != MP_OKAY)
 	return res;
 }
 } else if((cmp = s_mp_cmp(a, b)) > 0) {  /* different sign: a > b   */
 /* If the output is going to be clobbered, we will use a temporary
 variable; otherwise, we'll do it without touching the memory
 allocator at all, if possible
 */
 if(c == b) {
 mp_int  tmp;
 }
 s_mp_exch(&tmp, c);
 mp_clear(&tmp);
 } else {
 if(c != b && ((res = mp_copy(b, c)) != MP_OKAY))
 	return res;
 if((res = s_mp_sub(c, a)) != MP_OKAY)
 	return res;
 }
 return res;
 if((res = s_mp_mul(c, b)) != MP_OKAY)
 return res;
 }
 if(sgn == MP_ZPOS || s_mp_cmp_d(c, 0) == MP_EQ)
 SIGN(c) = MP_ZPOS;
 else
 SIGN(c) = sgn;
 return MP_OKAY;
 } /* end mp_mul() */
 /* }}} */
 if(r) {
 if((res = mp_copy(a, r)) != MP_OKAY)
 	return res;
 }
 if(q)
 mp_zero(q);
 return MP_OKAY;
 } else if(cmp == 0) {
 SIGN(&qtmp) = MP_ZPOS;
 if(s_mp_cmp_d(&rtmp, 0) == MP_EQ)
 SIGN(&rtmp) = MP_ZPOS;
 /* Copy output, if it is needed      */
 if(q)
 s_mp_exch(&qtmp, q);
 if(r)
 s_mp_exch(&rtmp, r);
 CLEANUP:
 mp_clear(&rtmp);
 mp_clear(&qtmp);
 mp_err mp_expt(mp_int *a, mp_int *b, mp_int *c)
 {
 mp_int   s, x;
 mp_err   res;
 mp_digit d;
-int      dig, bit;
+unsigned int bit, dig;
 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
 if(mp_cmp_z(b) < 0)
 return MP_RANGE;
 d = DIGIT(b, dig);
 /* Loop over bits of each non-maximal digit */
 for(bit = 0; bit < DIGIT_BIT; bit++) {
 if(d & 1) {
 	if((res = s_mp_mul(&s, &x)) != MP_OKAY)
 	  goto CLEANUP;
 }
 d >>= 1;
 if((res = s_mp_sqr(&x)) != MP_OKAY)
 	goto CLEANUP;
 }
 }
 d >>= 1;
 if((res = s_mp_sqr(&x)) != MP_OKAY)
 goto CLEANUP;
 }
 if(mp_iseven(b))
 SIGN(&s) = SIGN(a);
 res = mp_copy(&s, c);
 if(SIGN(m) == MP_NEG)
 return MP_RANGE;
 /*
 If |a| > m, we need to divide to get the remainder and take the
 absolute value.
 If |a| < m, we don't need to do any division, just copy and adjust
 the sign (if a is negative).
 If |a| == m, we can simply set the result to zero.
 that |a| != m, so we do those first.
 */
 if((mag = s_mp_cmp(a, m)) > 0) {
 if((res = mp_div(a, m, NULL, c)) != MP_OKAY)
 return res;
 if(SIGN(c) == MP_NEG) {
 if((res = mp_add(c, m, c)) != MP_OKAY)
 	return res;
 }
 if(mp_cmp_z(a) < 0) {
 if((res = mp_add(c, m, c)) != MP_OKAY)
 	return res;
 }
 } else {
 mp_zero(c);
 }
 /* Cannot take square root of a negative value */
 if(SIGN(a) == MP_NEG)
 return MP_RANGE;
 /* Special cases for zero and one, trivial     */
 if(mp_cmp_d(a, 0) == MP_EQ || mp_cmp_d(a, 1) == MP_EQ)
 return mp_copy(a, b);
 /* Initialize the temporaries we'll use below  */
 if((res = mp_init_size(&t, USED(a))) != MP_OKAY)
 return res;
 /* Compute an initial guess for the iteration as a itself */
 s_mp_exch(&x, b);
 CLEANUP:
 mp_clear(&x);
 X:
 mp_clear(&t);
 return res;
 } /* end mp_sqrt() */
 mp_exptmod(a, b, m, c)
 Compute c = (a ** b) mod m.  Uses a standard square-and-multiply
 method with modular reductions at each step. (This is basically the
 same code as mp_expt(), except for the addition of the reductions)
 The modular reductions are done using Barrett's algorithm (see
 s_mp_reduce() below for details)
 */
 mp_err mp_exptmod(mp_int *a, mp_int *b, mp_int *m, mp_int *c)
 {
 mp_int   s, x, mu;
 mp_err   res;
 mp_digit d, *db = DIGITS(b);
 mp_size  ub = USED(b);
-int      dig, bit;
+unsigned int bit, dig;
 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
 if(mp_cmp_z(b) < 0 || mp_cmp_z(m) <= 0)
 return MP_RANGE;
 goto MU;
 mp_set(&s, 1);
 /* mu = b^2k / m */
 s_mp_add_d(&mu, 1);
 s_mp_lshd(&mu, 2 * USED(m));
 if((res = mp_div(&mu, m, &mu, NULL)) != MP_OKAY)
 goto CLEANUP;
 /* Loop over digits of b in ascending order, except highest order */
 {
 mp_int  tmp;
 int     out;
 ARGCHK(a != NULL, MP_EQ);
 mp_init(&tmp); mp_set_int(&tmp, z);
 out = mp_cmp(a, &tmp);
 mp_clear(&tmp);
 return out;
 /* Initialize t */
 if(mp_isodd(&u)) {
 if((res = mp_copy(&v, &t)) != MP_OKAY)
 goto CLEANUP;
 /* t = -v */
 if(SIGN(&v) == MP_ZPOS)
 SIGN(&t) = MP_NEG;
 else
 SIGN(&t) = MP_ZPOS;
 } else {
 if((res = mp_copy(&u, &t)) != MP_OKAY)
 goto CLEANUP;
 }
 if(x)
 	if((res = mp_copy(&C, x)) != MP_OKAY) goto CLEANUP;
 if(y)
 	if((res = mp_copy(&D, y)) != MP_OKAY) goto CLEANUP;
 if(g)
 	if((res = mp_mul(&gx, &v, g)) != MP_OKAY) goto CLEANUP;
 break;
 }
 /*------------------------------------------------------------------------*/
 /* {{{ More I/O Functions */
 /* {{{ mp_read_signed_bin(mp, str, len) */
 /*
 mp_read_signed_bin(mp, str, len)
 Read in a raw value (base 256) into the given mp_int
 */
 return res;
 if((res = mp_add_d(mp, str[ix], mp)) != MP_OKAY)
 return res;
 }
 return MP_OKAY;
 } /* end mp_read_unsigned_bin() */
 /* }}} */
 /* {{{ mp_unsigned_bin_size(mp) */
 int     mp_unsigned_bin_size(mp_int *mp)
 {
 mp_digit   topdig;
 int        count;
 ARGCHK(mp != NULL, 0);
 return MP_OKAY;
 }
 /* Generate digits in reverse order */
 while(dp < end) {
-int      ix;
+unsigned int ix;
 d = *dp;
 for(ix = 0; ix < sizeof(mp_digit); ++ix) {
 *spos = d & UCHAR_MAX;
 d >>= CHAR_BIT;
 ++len;
 d >>= 1;
 }
 return len;
 } /* end mp_count_bits() */
 /* }}} */
 /* {{{ mp_read_radix(mp, str, radix) */
 {
 int     ix = 0, val = 0;
 mp_err  res;
 mp_sign sig = MP_ZPOS;
 ARGCHK(mp != NULL && str != NULL && radix >= 2 && radix <= MAX_RADIX,
 	 MP_BADARG);
 mp_zero(mp);
 /* Skip leading non-digit characters until a digit or '-' or '+' */
 while(str[ix] &&
 	(s_mp_tovalue(str[ix], radix) < 0) &&
 	str[ix] != '-' &&
 	str[ix] != '+') {
 ++ix;
 }
 /* {{{ mp_value_radix_size(num, qty, radix) */
 /* num = number of digits
 qty = number of bits per digit
 radix = target base
 Return the number of digits in the specified radix that would be
 needed to express 'num' digits of 'qty' bits each.
 */
 int    mp_value_radix_size(int num, int qty, int radix)
 {
 /* }}} */
 /* {{{ mp_toradix(mp, str, radix) */
-mp_err mp_toradix(mp_int *mp, unsigned char *str, int radix)
+mp_err mp_toradix(mp_int *mp, char *str, int radix)
 {
 int  ix, pos = 0;
 ARGCHK(mp != NULL && str != NULL, MP_BADARG);
 ARGCHK(radix > 1 && radix <= MAX_RADIX, MP_RANGE);
 str[pos--] = '\0';
 /* Reverse the digits and sign indicator     */
 ix = 0;
 while(ix < pos) {
-char tmp = str[ix];
+char _tmp = str[ix];
 str[ix] = str[pos];
-str[pos] = tmp;
+str[pos] = _tmp;
 ++ix;
 --pos;
 }
 mp_clear(&tmp);
 }
 return MP_OKAY;
 /* {{{ Arithmetic helpers */
 /* {{{ s_mp_lshd(mp, p) */
 /*
 Shift mp leftward by p digits, growing if needed, and zero-filling
 the in-shifted digits at the right end.  This is a convenient
 alternative to multiplication by powers of the radix
 */
 mp_err   s_mp_lshd(mp_int *mp, mp_size p)
 {
 mp_err   res;
 mp_size  pos;
 mp_digit *dp;
-int     ix;
+int ix;
 if(p == 0)
 return MP_OKAY;
 if((res = s_mp_pad(mp, USED(mp) + p)) != MP_OKAY)
 pos = USED(mp) - 1;
 dp = DIGITS(mp);
 /* Shift all the significant figures over as needed */
 for(ix = pos - p; ix >= 0; ix--)
 dp[ix + p] = dp[ix];
 /* Fill the bottom digits with zeroes */
-for(ix = 0; ix < p; ix++)
+for(ix = 0; (unsigned)ix < p; ix++)
 dp[ix] = 0;
 return MP_OKAY;
 } /* end s_mp_lshd() */
 /* }}} */
 /* {{{ s_mp_rshd(mp, p) */
 /*
 Shift mp rightward by p digits.  Maintains the invariant that
 digits above the precision are all zero.  Digits shifted off the
 end are lost.  Cannot fail.
 */
 /* {{{ s_mp_mul_2(mp) */
 mp_err s_mp_mul_2(mp_int *mp)
 {
-int      ix;
+unsigned int ix;
 mp_digit kin = 0, kout, *dp = DIGITS(mp);
 mp_err   res;
 /* Shift digits leftward by 1 bit */
 for(ix = 0; ix < USED(mp); ix++) {
 mp_err    s_mp_mul_2d(mp_int *mp, mp_digit d)
 {
 mp_err   res;
 mp_digit save, next, mask, *dp;
 mp_size  used;
-int      ix;
+unsigned int ix;
 if((res = s_mp_lshd(mp, d / DIGIT_BIT)) != MP_OKAY)
 return res;
 dp = DIGITS(mp); used = USED(mp);
 accomplish by multiplying a and b by a constant.  This constant is
 returned (so that it can be divided back out of the remainder at the
 end of the division process).
 We multiply by the smallest power of 2 that gives us a leading digit
 at least half the radix.  By choosing a power of 2, we simplify the
 multiplication and division steps to simple shifts.
 */
 mp_digit s_mp_norm(mp_int *a, mp_int *b)
 {
 mp_digit  t, d = 0;
 t = DIGIT(b, USED(b) - 1);
 while(t < (RADIX / 2)) {
 t <<= 1;
 ++d;
 }
 if(d != 0) {
 s_mp_mul_2d(a, d);
 s_mp_mul_2d(b, d);
 }
 /* If there is a precision increase, take care of it here; the above
 test guarantees we have enough storage to do this safely.
 */
 if(k) {
 dp[max] = k;
 USED(a) = max + 1;
 }
 s_mp_clamp(a);
 return MP_OKAY;
 } /* end s_mp_mul_d() */
 /* }}} */
 /* {{{ s_mp_div_d(mp, d, r) */
 *pa++ = ACCUM(w);
 w = CARRYOUT(w);
 }
 /* If we run out of 'b' digits before we're actually done, make
 sure the carries get propagated upward...
 */
 used = USED(a);
 while(w && ix < used) {
 w += *pa;
 *pa++ = ACCUM(w);
 }
 /* Clobber any leading zeroes we created    */
 s_mp_clamp(a);
 /*
 If there was a borrow out, then |b| > |a| in violation
 of our input invariant.  We've already done the work,
 but we'll at least complain about it...
 */
 if(w)
 #ifndef SHRT_MUL
 s_mp_mul(&q, m);
 s_mp_mod_2d(&q, (mp_digit)(DIGIT_BIT * (um + 1)));
 #else
 s_mp_mul_dig(&q, m, um + 1);
 #endif
 /* x = x - q */
 if((res = mp_sub(x, &q, x)) != MP_OKAY)
 goto CLEANUP;
 /* Outer loop:  Digits of b */
 pb = DIGITS(b);
 for(ix = 0; ix < ub; ++ix, ++pb) {
 if(*pb == 0)
 continue;
 /* Inner product:  Digits of a */
 pa = DIGITS(a);
 for(jx = 0; jx < ua; ++jx, ++pa) {
 mp_digit *pa, *pt;
 for(ix = 0; ix < len; ++ix, ++b) {
 if(*b == 0)
 continue;
 pa = a;
 for(jx = 0; jx < len; ++jx, ++pa) {
 pt = out + ix + jx;
 w = *b * *pa + k + *pt;
 *pt = ACCUM(w);
 overflow, we have to check explicitly for overflow conditions
 before they happen.
 */
 for(jx = ix + 1, pa2 = DIGITS(a) + jx; jx < used; ++jx, ++pa2) {
 mp_word  u = 0, v;
 /* Store this in a temporary to avoid indirections later */
 pt = pbt + ix + jx;
 /* Compute the multiplicative step */
 w = *pa1 * *pa2;
 /* Compute the additive step */
 v = *pt + k;
 /* If we do not already have an overflow carry, check to see
 	 if the addition will cause one, and set the carry out if so
 */
 u |= ((MP_WORD_MAX - v) < w);
 /* Add in the rest, again ignoring overflow */
 w += v;
 k = CARRYOUT(k);
 /* If we are carrying out, propagate the carry to the next digit
 in the output.  This may cascade, so we have to be somewhat
 circumspect -- but we will have enough precision in the output
 that we won't overflow
 */
 kx = 1;
 while(k) {
 k = pbt[ix + jx + kx] + 1;
 pbt[ix + jx + kx] = ACCUM(k);
 ix = USED(a) - 1;
 while(ix >= 0) {
 /* Find a partial substring of a which is at least b */
 while(s_mp_cmp(&rem, b) < 0 && ix >= 0) {
 if((res = s_mp_lshd(&rem, 1)) != MP_OKAY)
 	goto CLEANUP;
 if((res = s_mp_lshd(&quot, 1)) != MP_OKAY)
 	goto CLEANUP;
 s_mp_clamp(&rem);
 --ix;
 }
 /* If we didn't find one, we're finished dividing    */
 if(s_mp_cmp(&rem, b) < 0)
 break;
 /* Compute a guess for the next quotient digit       */
 q = DIGIT(&rem, USED(&rem) - 1);
 if(q <= DIGIT(b, USED(b) - 1) && USED(&rem) > 1)
 q = (q << DIGIT_BIT) | DIGIT(&rem, USED(&rem) - 2);
 /* See what that multiplies out to                   */
 mp_copy(b, &t);
 if((res = s_mp_mul_d(&t, q)) != MP_OKAY)
 goto CLEANUP;
 /*
 If it's too big, back it off.  We should not have to do this
 more than once, or, in rare cases, twice.  Knuth describes a
 method by which this could be reduced to a maximum of once, but
 I didn't implement that here.
 */
 */
 DIGIT(&quot, 0) = q;
 }
 /* Denormalize remainder                */
 if(d != 0)
 s_mp_div_2d(&rem, d);
 s_mp_clamp(&quot);
 s_mp_clamp(&rem);
 /* Copy quotient back to output         */
 s_mp_exch(&quot, a);
 /* Copy remainder back to output        */
 s_mp_exch(&rem, b);
 CLEANUP:
 mp_clear(&rem);
 bit = k % DIGIT_BIT;
 mp_zero(a);
 if((res = s_mp_pad(a, dig + 1)) != MP_OKAY)
 return res;
 DIGIT(a, dig) |= (1 << bit);
 return MP_OKAY;
 } /* end s_mp_2expt() */
 mp_digit *ap = DIGITS(a);
 if(ua > 1)
 return MP_GT;
 if(*ap < d)
 return MP_LT;
 else if(*ap > d)
 return MP_GT;
 else
 return MP_EQ;
 --dp; --ix;
 }
 return ((uv - 1) * DIGIT_BIT) + extra;
 }
 return -1;
 } /* end s_mp_ispow2() */
 expected to know what you're up to.
 */
 int      s_mp_tovalue(char ch, int r)
 {
 int    val, xch;
 if(r > 36)
 xch = ch;
 else
 xch = toupper(ch);
 val = xch - 'a' + 36;
 else if(xch == '+')
 val = 62;
 else if(xch == '/')
 val = 63;
 else
 return -1;
 if(val < 0 || val >= r)
 return -1;
 the value in the given radix.
 The results may be odd if you use a radix < 2 or > 64, you are
 expected to know what you're doing.
 */
 char     s_mp_todigit(int val, int r, int low)
 {
 char   ch;
 if(val < 0 || val >= r)
 /* }}} */
 /* {{{ s_mp_outlen(bits, radix) */
 /*
 Return an estimate for how long a string is needed to hold a radix
 r representation of a number with 'bits' significant bits.
 Does not include space for a sign or a NUL terminator.
 */
 /*------------------------------------------------------------------------*/
 /* HERE THERE BE DRAGONS                                                  */
 /* crc==4242132123, version==2, Sat Feb 02 06:43:52 2002 */
-/* $Source: /cvs/libtom/libtommath/mtest/mpi.c,v $ */
+/* $Source$ */
-/* $Revision: 1.2 $ */
+/* $Revision$ */
-/* $Date: 2005/05/05 14:38:47 $ */
+/* $Date$ */

Mercurial > dropbear

comparison libtommath/mtest/mpi.c @ 1436:60fc6476e044