/* Implementation of bignums by Henrik.Johansson@Nexus.Comm.SE in 1991 * version 1.2 */ /* Modifications in this file (Nov 1999, Rainer Wichmann): * - a few compiler warnings fixed: * gcc insists on 'if ( ( b_eq_qr = ((b == q) || (b == r)) ) )' * instead of 'if ( b_eq_qr = ((b == q) || (b == r)) )' * - code that is not used in samhain is #defined out by: * #if (defined (SH_WITH_CLIENT) || defined (SH_WITH_SERVER)) * or * #ifdef UNUSED_CODE * - the error message (error_string[]) is enclosed in N_(),_() macros * - the following four #includes have been added */ #include "config_xor.h" #include #ifdef HAVE_MEMORY_H #include #endif #include "samhain.h" #include "bignum.h" #if (!defined(HAVE_LIBGMP) || !defined(HAVE_GMP_H)) && (defined (SH_WITH_CLIENT) || defined (SH_WITH_SERVER)) #define DIGIT BIGNUM_DIGIT #define DIGIT2 BIGNUM_TWO_DIGITS #ifndef TRUE #define TRUE (1 == 1) #endif #ifndef FALSE #define FALSE (1 != 1) #endif #define MIN_ALLOC ((sizeof(long) / sizeof(DIGIT)) << 1) /* Don't expect the "sign" field to have the right value (BIG_SIGN_0) at * all times, so compare real content of bignum with zerop. * Maybe it would be better to be sure at all times that the sign is right? */ #define zerop(BIG) ((*(BIG)->dp == 0) && ((BIG)->dgs_used == 1)) #define uonep(BIG) ((*(BIG)->dp == 1) && ((BIG)->dgs_used == 1)) #define NEGATE_SIGN(SGN) (-(SGN)) #define DIGIT_BITS BIGNUM_DIGIT_BITS #define DIGIT2_BITS BIGNUM_TWO_DIGITS_BITS #define DIGIT_PART(N) ((DIGIT)((N) & ((((DIGIT2)1) << DIGIT_BITS) - 1))) #define RESULT_MINUSP(RES) (((RES) & ((DIGIT2)1 << (DIGIT2_BITS - 1))) != 0) #define SHIFT_DIGIT_DOWN(N) (DIGIT_PART((N) >> DIGIT_BITS)) #define SHIFT_DIGIT_UP(N) ((DIGIT2)((N) << DIGIT_BITS)) #define DIGIT_BASE (((DIGIT2)1) << DIGIT_BITS) #define MAX_DIGIT ((((DIGIT2)1) << DIGIT_BITS) - 1) #define uint unsigned int #define ulong unsigned long char *last_string = NULL; char *big_end_string = NULL; ulong length_last_string; static char big_base_digits1[] = N_("00112233445566778899AaBbCcDdEeFfGgHhIiJjKkLlMmNnOoPpQqRrSsTtUuVvWwXxYyZz"); static char error_string[] = N_("(big_errno != 0, something has gone wrong)"); static char big_base_digits[73] = "\0"; bignum big_one; bignum tmp_string, tmp_add, tmp_mul, tmp_q, tmp_r, tmp_rand, tmp_round; DIGIT *tmp_digits; ulong tmp_ulong; /* #define COUNT_ALLOC */ #ifdef COUNT_ALLOC ulong dgs_alloc = 0, digits_freed = 0; #endif int big_errno = BIG_OK; double log2tbl[] = { -1.0, -1.0, 1.000000000000000, 1.584961891174317, 2.000000000000002, 2.321928024291994, 2.584960937500003, 2.807353973388675, 3.000000000000004, 3.169923782348637, 3.321928024291997, 3.459430694580083, 3.584960937500005, 3.700439453125006, 3.807353973388678, 3.906888961791999, 4.000000000000014, 4.087459564208999, 4.169921875000016, 4.247924804687517, 4.321926116943377, 4.392314910888691, 4.459430694580098, 4.523559570312520, 4.584960937500021, 4.643856048584007, 4.700439453125023, 4.754886627197290, 4.807353973388697, 4.857978820800807, 4.906887054443386, 4.954193115234403, 5.000000000000028, 5.044391632080107, 5.087459564209015, 5.129280090332062, 5.169921875000032 }; struct digit_blocks { int digCnt; DIGIT dig_base; } big_block[37]; /* The code uses index 2 to 36 */ /* ---------------------------------------------------------------------- * All static utility functions are placed first in this file. */ #ifndef HAVE_STRCHR static char * strchr (char * str_ptr, char ch) { do { if (*str_ptr == ch) { return str_ptr; } } while (*str_ptr++ != '\0'); return NULL; } #endif #if 0 extern int printf(); static void print_digits(prefix, a) char *prefix; bignum *a; { unsigned long i; ulong tmp; printf("%s", prefix); if (a->sign == BIG_SIGN_MINUS) { printf("- "); } else { if (a->sign == BIG_SIGN_PLUS) { printf("+ "); } else { printf("+/- "); } } for (i = a->dgs_used - 1; i > 0; i--) { tmp = a->dp[i]; printf("%lu, ", tmp); } tmp = *a->dp; printf("%lu\n", tmp); } #else #define print_digits(prefix, big) /* */ #endif static void init_digit_blocks(void) { uint digcnt, base; DIGIT maxdigit, curdigit, tmp; for (base = 2; base <= 36; base++) { tmp = ((1U << (DIGIT_BITS - 1)) / base); maxdigit = tmp * 2; curdigit = 1; digcnt = 0; while (curdigit < maxdigit) { digcnt++; curdigit *= base; } big_block[base].digCnt = digcnt; big_block[base].dig_base = curdigit; } } #ifdef COUNT_ALLOC static void free_digits(DIGIT * dp, ulong count) { /* Mar 4 2000 R. Wichmann: check for dp == NULL */ if (!dp) return; digits_freed += count; free((char *)dp); dp = NULL; } #else /* June 5 2000 R. Wichmann: check for dp == NULL */ /* #define free_digits(DP,CNT) free((char *)DP) */ static void free_digits(DIGIT * dp, ulong count) { /* Mar 4 2000 R. Wichmann: check for dp == NULL */ if (!dp) return; (void) count; free((char *)dp); dp = NULL; } #endif static DIGIT * allocate_digits(ulong alloclen) { DIGIT *digit_ptr; if (big_errno != BIG_OK) { return NULL; } #ifdef COUNT_ALLOC dgs_alloc += alloclen; #endif digit_ptr = calloc(alloclen, sizeof(DIGIT)); if (digit_ptr == NULL) { big_errno = BIG_MEMERR; return NULL; } return digit_ptr; } static bigerr_t newsize(DIGIT **dp_p, ulong *cursize_p, ulong minsz, ulong newsz) { if (*cursize_p >= minsz) { return big_errno; } free_digits(*dp_p, *cursize_p); if (newsz < MIN_ALLOC) { newsz = MIN_ALLOC; } if ((*dp_p = allocate_digits(newsz)) == NULL) { return big_errno; } *cursize_p = newsz; return big_errno; } /* `memcpy' uses an `int' as counter. If an int can't hold a big enough * counter, then `digits_cpy' becomes a function (and a little slower, * probably. Unfortunately `memcpy' takes a signed counter, but _that_ * size of numbers are almost too big! Or isn't it? */ #ifdef MEMCPY_LONG_COUNTER /* extern char *memcpy(); */ #define digits_cpy(dst, src, count) memcpy((char *)(dst), (char *)(src), \ (count) * sizeof(DIGIT)) #else static void digits_cpy(DIGIT *dst, DIGIT *src, ulong count) { while (count-- > 0) { *dst++ = *src++; } } #endif static DIGIT * copy_digits(DIGIT *src, ulong cp_count, ulong alloclen) { DIGIT *dst; if (big_errno != BIG_OK) { return NULL; } if ((dst = allocate_digits(alloclen)) == NULL) { return NULL; } digits_cpy(dst, src, cp_count); return dst; } static void add_digit(bignum * a, DIGIT d) { DIGIT2 res = d; DIGIT *last_digit, *digit_ptr = a->dp, *digvect; last_digit = digit_ptr + a->dgs_used; while (digit_ptr < last_digit) { res = *digit_ptr + res; *digit_ptr = DIGIT_PART(res); res = SHIFT_DIGIT_DOWN(res); digit_ptr++; if (res == 0) { break; } } if (res != 0) { if (a->dgs_used < a->dgs_alloc) { *digit_ptr = DIGIT_PART(res); } else { if ((digvect = copy_digits(a->dp, a->dgs_used, a->dgs_used + 4)) == NULL) { return; /* big_errno will be set to last error */ } digvect[a->dgs_used] = DIGIT_PART(res); free_digits(a->dp, a->dgs_alloc); a->dgs_alloc = a->dgs_used + 4; a->dp = digvect; } a->dgs_used++; } } static DIGIT vect_div_digit(DIGIT *digvect, ulong *len, DIGIT d) { DIGIT *digit_ptr = digvect + *len - 1, newval; DIGIT2 rest = 0; if (d == 0) { big_errno = BIG_DIV_ZERO; return -1; } if (d == 1) { return 0; } while (digit_ptr >= digvect) { rest = (DIGIT2)SHIFT_DIGIT_UP(rest); newval = DIGIT_PART((*digit_ptr + rest) / d); rest = (*digit_ptr + rest) % d; *digit_ptr = newval; digit_ptr--; } if (*len > 1) { if (digvect[*len - 1] == 0) { *len -= 1; } } return DIGIT_PART(rest); } static DIGIT udiv_digit(bignum *a, DIGIT d) { DIGIT res; res = vect_div_digit(a->dp, &a->dgs_used, d); if (zerop(a)) { a->sign = BIG_SIGN_0; } return res; } static DIGIT vect_mul_digit(DIGIT *digvect, ulong len, DIGIT x) { DIGIT *digit_ptr = digvect + len; DIGIT2 res = 0; while (digvect < digit_ptr) { res += *digvect * (DIGIT2)x; *digvect = DIGIT_PART(res); res = SHIFT_DIGIT_DOWN(res); digvect++; } return DIGIT_PART(res); } static void umul_digit(bignum *a, DIGIT x) { DIGIT ovfl, *digvect; ovfl = vect_mul_digit(a->dp, a->dgs_used, x); if (ovfl != 0) { if (a->dgs_used < a->dgs_alloc) { a->dp[a->dgs_used] = ovfl; } else { digvect = copy_digits(a->dp, a->dgs_used, a->dgs_used + 4); digvect[a->dgs_used] = ovfl; free_digits(a->dp, a->dgs_alloc); a->dgs_alloc = a->dgs_used + 4; a->dp = digvect; } a->dgs_used++; } } static int ucompare_digits(bignum *a, bignum *b) { DIGIT *a_ptr, *b_ptr; int retval = 0; if (a->dgs_used == b->dgs_used) { a_ptr = a->dp + a->dgs_used - 1; b_ptr = b->dp + b->dgs_used - 1; while ((*a_ptr == *b_ptr) && (a_ptr >= a->dp)) { a_ptr--; b_ptr--; } if (a_ptr < a->dp) { return retval; } else { if (retval == 0) retval = (*a_ptr > *b_ptr) ? 1 : -1; } return retval; } return (a->dgs_used > b->dgs_used) ? 1 : -1; } static void uadd_digits(bignum *a, bignum *b, bignum *c) { DIGIT *dp_x, *dp_y, *dp_z, *res_dp, *end_x, *end_y; ulong len_x, len_y; DIGIT2 res = 0; if (a->dgs_used > b->dgs_used) { dp_x = a->dp; len_x = a->dgs_used; dp_y = b->dp; len_y = b->dgs_used; } else { dp_x = b->dp; len_x = b->dgs_used; dp_y = a->dp; len_y = a->dgs_used; } end_x = dp_x + len_x; end_y = dp_y + len_y; if (c->dgs_alloc >= len_x) { dp_z = c->dp; } else { if (newsize(&tmp_add.dp, &tmp_add.dgs_alloc, len_x, len_x + 4) != BIG_OK) { return; } dp_z = tmp_add.dp; } res_dp = dp_z; while (dp_y < end_y) { res += ((DIGIT2)*dp_x++) + *dp_y++; *res_dp++ = DIGIT_PART(res); res = SHIFT_DIGIT_DOWN(res); } while (dp_x < end_x) { res += *dp_x++; *res_dp++ = DIGIT_PART(res); res = SHIFT_DIGIT_DOWN(res); } if (res != 0) { *res_dp++ = DIGIT_PART(res); } if (dp_z != c->dp) { tmp_digits = c->dp; c->dp = tmp_add.dp; tmp_add.dp = tmp_digits; tmp_ulong = c->dgs_alloc; c->dgs_alloc = tmp_add.dgs_alloc; tmp_add.dgs_alloc = tmp_ulong; } c->dgs_used = res_dp - c->dp; } static void usub_digits(bignum *a, bignum *b, bignum *c) { DIGIT *dp_x, *dp_y, *dp_z, *res_dp, *end_x, *end_y; ulong len_x, len_y; DIGIT2 res = 0; dp_x = a->dp; len_x = a->dgs_used; dp_y = b->dp; len_y = b->dgs_used; end_x = dp_x + len_x; end_y = dp_y + len_y; if (c->dgs_alloc >= len_x) { dp_z = c->dp; } else { if (newsize(&tmp_add.dp, &tmp_add.dgs_alloc, len_x, len_x + 2) != BIG_OK) { return; } dp_z = tmp_add.dp; } res_dp = dp_z; while (dp_y < end_y) { res = ((DIGIT2)*dp_x++) - *dp_y++ - RESULT_MINUSP(res); *res_dp++ = DIGIT_PART(res); } while (dp_x < end_x) { res = *dp_x++ - RESULT_MINUSP(res); *res_dp++ = DIGIT_PART(res); } #ifdef BIG_CHECK_LIMITS if (RESULT_MINUSP(res) != 0) { big_errno = BIG_ALGERR; return; } #endif while ((*--res_dp == 0) && (res_dp > dp_z)) { /* Move pointer down until we reach a non-zero */ } if (dp_z != c->dp) { tmp_digits = c->dp; c->dp = tmp_add.dp; tmp_add.dp = tmp_digits; tmp_ulong = tmp_add.dgs_alloc; tmp_add.dgs_alloc = c->dgs_alloc; c->dgs_alloc = tmp_ulong; } c->dgs_used = res_dp - dp_z + 1; } /* * This (pseudo) random number generator is not very good. It has a long * period [ 2 ^ (DIGIT_BITS * 2) (?)] before it starts with the same sequence * again, but the lower bits are "less random" than they should be. I've * solved it by using word of double length, and returning the upper bits. * Please mail me if you know how to make it into a "more random" generator. * One important thing though: it will have to reasonably fast. * The good thing with this one is that it is very portable, but that doesn't * help much when you want _good_ random numbers. * Henrik.Johansson@Nexus.Comm.SE */ #ifdef UNUSED_CODE static DIGIT rand(void) { static DIGIT2 x1 = 33, x2 = 45; /* Just give them two starting values */ if (x2 = BIG_RAND_A2 * x2 + BIG_RAND_C2, x2 == 45) { x1 = BIG_RAND_A1 * x1 + BIG_RAND_C1; } return SHIFT_DIGIT_DOWN(x1 + x2); /* Skip least significant bits */ } /* UNUSED_CODE */ #endif /* ---------------------------------------------------------------------- * All external functions comes here. */ bigerr_t big_init_pkg() { strcpy (big_base_digits, _(big_base_digits1)); /* known to fit */ init_digit_blocks(); big_create(&tmp_string); big_create(&tmp_add); big_create(&tmp_mul); big_create(&tmp_q); big_create(&tmp_r); big_create(&tmp_rand); big_create(&big_one); big_set_long((long)1, &big_one); length_last_string = 10; if ((last_string = calloc(1,length_last_string)) == NULL) { big_errno = BIG_MEMERR; } return big_errno; } void big_release_pkg() { big_destroy(&tmp_string); big_destroy(&tmp_add); big_destroy(&tmp_mul); big_destroy(&tmp_q); big_destroy(&tmp_r); big_destroy(&tmp_round); big_destroy(&big_one); if (last_string != NULL) free(last_string); #ifdef COUNT_ALLOC printf("Allocated digits: %lu\n", dgs_alloc); printf("Freed digits: %lu\n", digits_freed); #endif } bigerr_t big_create(a) bignum *a; { /* Make sure that big_create alway returns a NULL bignum * that can safely be destroyed. */ if (a != NULL) { a->dp = NULL; } if (big_errno != BIG_OK) { return big_errno; } if (a == NULL) { big_errno = BIG_ARGERR; return big_errno; } a->sign = BIG_SIGN_0; a->dgs_used = 1; if ((a->dp = allocate_digits((long)sizeof(long))) == NULL) { return big_errno; } a->dgs_alloc = sizeof(long); *a->dp = 0; return big_errno; } void big_destroy(a) bignum *a; { /* Mar 4, 2000 R. Wichmann: check for a == NULL */ /* free_digits() checks for a->dp == NULL */ if (a != NULL) free_digits(a->dp, a->dgs_alloc); } ulong big_bitcount(a) bignum *a; { int bits = 0; DIGIT high_digit; if (big_errno != BIG_OK) { return 0; } high_digit = a->dp[a->dgs_used - 1]; while (high_digit != 0) { bits++; high_digit >>= 1; } return DIGIT_BITS * (a->dgs_used - 1) + bits; } bigerr_t big_set_big(a, b) bignum *a; bignum *b; { if ((big_errno != BIG_OK) || (a == b)) { return big_errno; } if (newsize(&b->dp, &b->dgs_alloc, a->dgs_used, a->dgs_used) != BIG_OK) { return big_errno; } b->dgs_used = a->dgs_used; b->sign = a->sign; digits_cpy(b->dp, a->dp, a->dgs_used); return big_errno; } void big_set_ulong(n, a) ulong n; bignum *a; { unsigned int i; if (big_errno != BIG_OK) { return; } if (n == 0) { *a->dp = 0; a->dgs_used = 1; a->sign = BIG_SIGN_0; } else { a->dgs_used = 0; for (i = 0; i < sizeof(long) / sizeof(DIGIT); i++) { if (n == 0) { break; } a->dgs_used++; a->dp[i] = DIGIT_PART(n); n = SHIFT_DIGIT_DOWN(n); } a->sign = BIG_SIGN_PLUS; } print_digits("set_(u)long: a = ", a); } void big_set_long(n, a) long n; bignum *a; { ulong m; m = (n < 0) ? - n : n; big_set_ulong(m, a); if (!(n >= 0)) { a->sign = BIG_SIGN_MINUS; } } bigerr_t big_set_string(numstr, base, a) char *numstr; int base; bignum *a; { char *src, *maxdigit, *chrptr; DIGIT dig_base, dig_sum, last_base; int cnt, maxcnt; if (big_errno != BIG_OK) { return big_errno; } big_end_string = numstr; if ((base < 2) || (base > 36) || (numstr == NULL) || (a == NULL)) { big_errno = BIG_ARGERR; return big_errno; } src = numstr; maxdigit = big_base_digits + (base << 1); maxcnt = big_block[base].digCnt; dig_base = big_block[base].dig_base; a->dgs_used = 1; *a->dp = 0; a->sign = BIG_SIGN_PLUS; /* Assume it will be positive */ while (strchr(" \t\n\r", *src) != NULL) /* Skip whitespaces */ { src++; } if ((*src == '-') || (*src == '+')) /* First non-white is a sign? */ { a->sign = (*src == '-') ? BIG_SIGN_MINUS : BIG_SIGN_PLUS; src++; } chrptr = strchr(big_base_digits, *src); if ((chrptr == NULL) || (chrptr >= maxdigit)) /* Next chr not a digit? */ { big_end_string = src; big_errno = BIG_ARGERR; return big_errno; } while (*src == '0') /* Next chr a '0'? */ { src++; /* Read past all '0'es */ } chrptr = strchr(big_base_digits, *src); if ((chrptr == NULL) || (chrptr >= maxdigit)) /* Next char not a digit */ { big_end_string = src; a->sign = BIG_SIGN_0; /* It was just a plain 0 */ return big_errno; } dig_sum = 0; cnt = 0; while ((chrptr = strchr(big_base_digits, *src)), (chrptr != NULL) && (chrptr < maxdigit)) { dig_sum = dig_sum * base + ((chrptr - big_base_digits) >> 1); if (++cnt >= maxcnt) { umul_digit(a, dig_base); add_digit(a, dig_sum); dig_sum = 0; cnt = 0; } src++; } if (cnt > 0) { last_base = base; while (--cnt > 0) { last_base *= base; } umul_digit(a, last_base); add_digit(a, dig_sum); } big_end_string = src; return big_errno; } int big_ulong(a, n) bignum *a; ulong *n; { ulong old_n; DIGIT *dp; if (big_errno != BIG_OK) { return FALSE; } if (a->dgs_used > sizeof(ulong) / sizeof(DIGIT)) { return FALSE; } dp = a->dp + a->dgs_used - 1; old_n = *n = *dp--; while ((dp >= a->dp) && (old_n < *n)) { old_n = *n; *n = SHIFT_DIGIT_UP(*n) + *dp--; } if (old_n >= *n) { return FALSE; } return FALSE; } int big_long(a, n) bignum *a; long *n; { long old_n; DIGIT *dp; if (a->dgs_used > sizeof(ulong) / sizeof(DIGIT)) { return FALSE; } dp = a->dp + a->dgs_used - 1; old_n = *n = *dp--; while ((dp >= a->dp) && (old_n <= *n)) { old_n = *n; *n = SHIFT_DIGIT_UP(*n) + *dp--; } if (old_n >= *n) { return FALSE; } if (a->sign == BIG_SIGN_MINUS) { *n = -*n; } return FALSE; } char * big_string(a, base) bignum *a; int base; { ulong bit_length, str_length; char *dst; DIGIT dig_base, rem; int cnt, maxcnt; if (big_errno != BIG_OK) { return _(error_string); } big_set_big(a, &tmp_string); /* Need more room than minimum here... */ bit_length = tmp_string.dgs_used * DIGIT_BITS; /* bit_length = big_bitcount(&tmp_string); */ str_length = (ulong)(bit_length / log2tbl[base] + 4); if (str_length > length_last_string) { if (last_string != NULL) free(last_string); if ((last_string = calloc(1,str_length)) == NULL) { big_errno = BIG_MEMERR; return NULL; } length_last_string = str_length; } dst = last_string + length_last_string - 1; *dst-- = '\0'; maxcnt = big_block[base].digCnt; dig_base = big_block[base].dig_base; while (tmp_string.dgs_used > 1) { rem = udiv_digit(&tmp_string, dig_base); for (cnt = 0; cnt < maxcnt; cnt++) { *dst-- = big_base_digits[(rem % base) << 1]; rem /= base; } } rem = *tmp_string.dp; do { *dst-- = big_base_digits[(rem % base) << 1]; rem /= base; } while (rem != 0); if (a->sign == BIG_SIGN_MINUS) { *dst = '-'; } else { dst++; } return dst; } bigerr_t big_negate(a, b) bignum *a; bignum *b; { big_set_big(a, b); b->sign = NEGATE_SIGN(a->sign); return big_errno; } int big_sign(a) bignum *a; { return a->sign; } bigerr_t big_abs(a, b) bignum *a; bignum *b; { big_set_big(a, b); if (a->sign == BIG_SIGN_MINUS) { b->sign = NEGATE_SIGN(a->sign); } return big_errno; } int big_compare(a, b) bignum *a; bignum *b; { if (a->sign == b->sign) { if (a->sign == 0) { return 0; } return (a->sign == BIG_SIGN_MINUS) ? -ucompare_digits(a, b) : ucompare_digits(a, b); } return b->sign - a->sign; } int big_lessp(a, b) bignum *a; bignum *b; { return big_compare(a, b) < 0; } int big_leqp(a, b) bignum *a; bignum *b; { return !(big_compare(a, b) > 0); } int big_equalp(a, b) bignum *a; bignum *b; { return big_compare(a, b) == 0; } int big_geqp(a, b) bignum *a; bignum *b; { return !(big_compare(a, b) < 0); } int big_greaterp(a, b) bignum *a; bignum *b; { return big_compare(a, b) > 0; } int big_zerop(a) bignum *a; { return a->sign == BIG_SIGN_0; } int big_evenp(a) bignum *a; { return ((*a->dp & 0x01) == 0); } int big_oddp(a) bignum *a; { return ((*a->dp & 0x01) == 1); } bigerr_t big_add(a, b, c) bignum *a; bignum *b; bignum *c; { int cmp; if (big_errno != BIG_OK) { return big_errno; } print_digits("add:\ta = ", a); print_digits("\tb = ", b); if (a->sign == b->sign) { uadd_digits(a, b, c); c->sign = a->sign; } else { cmp = ucompare_digits(a, b); if (cmp < 0) { usub_digits(b, a, c); if (zerop(c)) { c->sign = BIG_SIGN_0; } else { c->sign = b->sign; } } else if (cmp > 0) { usub_digits(a, b, c); if (zerop(c)) { c->sign = BIG_SIGN_0; } else { c->sign = a->sign; } } else { c->dgs_used = 1; *c->dp = 0; c->sign = BIG_SIGN_0; } } print_digits("\tc = ", c); return big_errno; } bigerr_t big_sub(a, b, c) bignum *a; bignum *b; bignum *c; { int cmp; if (big_errno != BIG_OK) { return big_errno; } print_digits("sub:\ta = ", a); print_digits("\tb = ", b); if (a->sign == BIG_SIGN_0) { big_set_big(b, c); big_negate(c, c); print_digits("\tc = ", c); return big_errno; } if (b->sign == BIG_SIGN_0) { big_set_big(a, c); print_digits("\tc = ", c); return big_errno; } cmp = ucompare_digits(a, b); if (cmp <= 0) { if (a->sign != b->sign) { uadd_digits(a, b, c); c->sign = a->sign; } else { usub_digits(b, a, c); c->sign = (zerop(c) ? BIG_SIGN_0 : NEGATE_SIGN(a->sign)); } } else if (cmp > 0) { if (a->sign != b->sign) { uadd_digits(a, b, c); c->sign = a->sign; } else { usub_digits(a, b, c); c->sign = (zerop(c) ? BIG_SIGN_0 : b->sign); } } else { c->dgs_used = 1; *c->dp = 0; c->sign = BIG_SIGN_0; } print_digits("\tc = ", c); return big_errno; } bigerr_t big_mul(a, b, c) bignum *a; bignum *b; bignum *c; { DIGIT *dp_x, *dp_xstart, *dp_xend; DIGIT *dp_y, *dp_ystart, *dp_yend; DIGIT *dp_z, *dp_zstart, *dp_zend, *dp_zsumstart; ulong len_x, len_y /* , len_z */; DIGIT2 res; DIGIT tmp_res; /* Without use of this, gcc (v 1.39) generates */ /* erroneous code on a sun386 machine */ /* Should be removed with #ifdef's, when */ /* not on a sun386, but... */ if (big_errno != BIG_OK) { return big_errno; } print_digits("mul:\ta = ", a); print_digits("\tb = ", b); if (zerop(a) || zerop(b)) { c->sign = BIG_SIGN_0; c->dgs_used = 1; *c->dp = 0; print_digits("(a=0 || b=0)c = ", c); return big_errno; } if (uonep(a)) { big_set_big(b, c); c->sign = (a->sign == b->sign) ? BIG_SIGN_PLUS : BIG_SIGN_MINUS; print_digits("(abs(a)=1)c = ", c); return big_errno; } if (uonep(b)) { big_set_big(a, c); c->sign = (a->sign == b->sign) ? BIG_SIGN_PLUS : BIG_SIGN_MINUS; print_digits("(abs(b)=1)c = ", c); return big_errno; } if (a->dgs_used < b->dgs_used) { dp_xstart = a->dp; len_x = a->dgs_used; dp_ystart = b->dp; len_y = b->dgs_used; } else { dp_xstart = b->dp; len_x = b->dgs_used; dp_ystart = a->dp; len_y = a->dgs_used; } if ((c == a) || (c == b)) { if (newsize(&tmp_mul.dp, &tmp_mul.dgs_alloc, len_x + len_y, len_x + len_y + 2) != BIG_OK) { return big_errno; } dp_zsumstart = tmp_mul.dp; /* len_z = tmp_mul.dgs_alloc; */ } else { if (newsize(&c->dp, &c->dgs_alloc, len_x + len_y, len_x + len_y + 2) != BIG_OK) { return big_errno; } dp_zsumstart = c->dp; /* len_z = c->dgs_alloc; */ } dp_xend = dp_xstart + len_x; dp_yend = dp_ystart + len_y; dp_zend = dp_zsumstart + len_y; for (dp_z = dp_zsumstart; dp_z < dp_zend; dp_z++) { *dp_z = 0; /* Zero out rightmost digits */ } dp_zstart = dp_zsumstart; for (dp_x = dp_xstart; dp_x < dp_xend; dp_x++) { dp_z = dp_zstart; tmp_res = 0; for (dp_y = dp_ystart; dp_y < dp_yend; dp_y++) { res = (DIGIT2)(*dp_x) * (*dp_y) + (*dp_z) + tmp_res; *dp_z++ = DIGIT_PART(res); tmp_res = SHIFT_DIGIT_DOWN(res); } *dp_z = tmp_res; dp_zstart++; } if (dp_zsumstart != c->dp) { tmp_digits = c->dp; c->dp = tmp_mul.dp; tmp_mul.dp = tmp_digits; tmp_ulong = c->dgs_alloc; c->dgs_alloc = tmp_mul.dgs_alloc; tmp_mul.dgs_alloc = tmp_ulong; } if (*dp_z == 0) { dp_z--; } c->dgs_used = dp_z - dp_zsumstart + 1; c->sign = a->sign * b->sign; print_digits("\tc = ", c); return big_errno; } bigerr_t big_trunc(a, b, q, r) bignum *a; bignum *b; bignum *q; bignum *r; { DIGIT *v_end, *q_end, *r_end, *src, *dst; DIGIT norm, qhat, t1, t2, t3, v1, v2, u1, u2, u3; DIGIT2 temp, res, carry; long a_l, b_l, q_l, r_l; #if 0 ulong i; #endif ulong j, m, n; int cmp, q_eq_a, q_eq_b, r_eq_a, r_eq_b; if (big_errno != BIG_OK) { return big_errno; } print_digits("div:\ta = ", a); print_digits("\tb = ", b); if (zerop(b)) { big_errno = BIG_DIV_ZERO; return big_errno; } if (q == r) { big_errno = BIG_ARGERR; return big_errno; } if (b->dgs_used == 1) { big_set_big(a, q); q->sign = ((a->sign == b->sign) ? BIG_SIGN_PLUS : BIG_SIGN_MINUS); *r->dp = udiv_digit(q, *b->dp); r->dgs_used = 1; r->sign = (zerop(r) ? BIG_SIGN_0 : a->sign); print_digits("\t3:q = ", q); print_digits("\t r = ", r); return big_errno; } if (big_long(a, &a_l)) /* Pretend it is a signed value */ { big_long(b, &b_l); /* |a| < |b| so this will succeed */ q_l = a_l / b_l; /* Compute with unsigned operators */ r_l = a_l % b_l; big_set_long((long)q_l, q); big_set_long((long)r_l, r); print_digits("\t4:q = ", q); print_digits("\t r = ", r); return big_errno; } cmp = ucompare_digits(a, b); /* Unsigned compare, that is... */ if (cmp < 0) { big_set_big(a, r); /* r = a (don't care about big_errno here) */ q->sign = BIG_SIGN_0; /* q = 0 */ *q->dp = 0; q->dgs_used = 1; print_digits("\t1:q = ", q); print_digits("\t r = ", r); return big_errno; } else if (cmp == 0) { q->sign = ((a->sign == b->sign) ? BIG_SIGN_PLUS : BIG_SIGN_MINUS); *q->dp = 1; /* q = 1 */ q->dgs_used = 1; r->sign = BIG_SIGN_0; /* r = 0 */ *r->dp = 0; r->dgs_used = 1; print_digits("\t2:q = ", q); print_digits("\t r = ", r); return big_errno; } q_eq_a = (q == a); q_eq_b = (q == b); if (q_eq_a || q_eq_b) { q = &tmp_q; } r_eq_a = (r == a); r_eq_b = (r == b); if (r_eq_a || r_eq_b) { r = &tmp_r; } if (newsize(&r->dp, &r->dgs_alloc, /* At least one more dig in r */ a->dgs_used + 1, a->dgs_used + 2) != BIG_OK) { return big_errno; } big_set_big(a, r); r->dp[a->dgs_used] = 0; /* In case no overflow in mult. */ n = b->dgs_used; v_end = b->dp + n - 1; norm = DIGIT_PART((DIGIT_BASE / ((DIGIT2)*v_end + 1))); if (norm != 1) { umul_digit(r, norm); umul_digit(b, norm); print_digits("r = ", r); print_digits("b = ", b); } m = a->dgs_used + 1 - b->dgs_used; r_end = r->dp + a->dgs_used; if (newsize(&q->dp, &q->dgs_alloc, m, m + 2) != BIG_OK) { return big_errno; } q_end = q->dp + m - 1; v1 = *v_end; v2 = *(v_end - 1); for (j = 0; j < m; j++) /* m steps through division */ { u1 = *r_end; /* routine */ u2 = *(r_end - 1); u3 = *(r_end - 2); qhat = ((u1 == v1) ? MAX_DIGIT : DIGIT_PART(((DIGIT2)u1 * DIGIT_BASE + u2) / v1)); while (1) { t3 = DIGIT_PART(temp = (DIGIT2)qhat * v2); t2 = DIGIT_PART(temp = SHIFT_DIGIT_DOWN(temp) + v1 * (DIGIT2)qhat); t1 = DIGIT_PART(SHIFT_DIGIT_DOWN(temp)); #if 0 printf("t1 = %lu, ", (ulong)t1); printf("t2 = %lu, ", (ulong)t2); printf("t3 = %lu\n", (ulong)t3); #endif if (t1 < u1) break; if (t1 > u1) {--qhat; continue; } if (t2 < u2) break; if (t2 > u2) { --qhat; continue; } if (t3 <= u3) break; qhat--; } /* This is a tricky one - multiply and subtract simultaneously */ carry = 1; res = 0; src = b->dp; dst = r->dp + m - j - 1; while (src <= v_end) { res = (DIGIT2)qhat * *(src++) + SHIFT_DIGIT_DOWN(res); carry += (DIGIT2)(*dst) + MAX_DIGIT - DIGIT_PART(res); *(dst++) = DIGIT_PART(carry); carry = DIGIT_PART(SHIFT_DIGIT_DOWN(carry)); } carry += (DIGIT2)(*dst) + MAX_DIGIT - SHIFT_DIGIT_DOWN(res); *dst = DIGIT_PART(carry); carry = DIGIT_PART(SHIFT_DIGIT_DOWN(carry)); if (carry == 0) { qhat--; src = b->dp; dst = r->dp + m - j - 1; while (dst <= r_end) { carry = (DIGIT2)(*dst) + *src++ + SHIFT_DIGIT_DOWN(carry); *dst++ = DIGIT_PART(carry); } *dst = 0; } *(q_end - j) = DIGIT_PART(qhat); r_end--; } r->sign = a->sign; #if 0 i = r->dgs_used; #endif while ((*r_end == 0) && (r_end > r->dp)) { r_end--; } if (r_end == r->dp) { r->dgs_used = 1; r->sign = BIG_SIGN_0; } else { r->dgs_used = r_end - r->dp + 1; r->sign = a->sign; } if (norm != 1) { udiv_digit(b, norm); udiv_digit(r, norm); } while ((*q_end == 0) && (q_end > q->dp)) /* This is not needed!(?) */ { q_end--; } /* was ifdef'd out already in original */ #if 0 i = m - 1; while ((i > 0) && (q->dp[i--] == 0)) { /* Loop through all zeroes */ } #endif q->dgs_used = q_end - q->dp + 1; q->sign = ((a->sign == b->sign) ? BIG_SIGN_PLUS : BIG_SIGN_MINUS); if (q_eq_a) { big_set_big(q, a); } else if (q_eq_b) { big_set_big(q, b); } if (r_eq_b) { big_set_big(r, b); } else if (r_eq_a) { big_set_big(r, a); } print_digits("\t5:q = ", q); print_digits("\t r = ", r); return big_errno; } bigerr_t big_floor(a, b, q, r) bignum *a; bignum *b; bignum *q; bignum *r; { int b_eq_qr, sign_eq; if ( ( b_eq_qr = ((b == q) || (b == r)) ) ) { big_set_big(b, &tmp_mul); } sign_eq = a->sign == b->sign; big_trunc(a, b, q, r); if (sign_eq) { return big_errno; } if (!zerop(r)) { if (b_eq_qr) { big_add(r, &tmp_mul, r); } else { big_add(r, b, r); } print_digits("big_one = ", &big_one); big_sub(q, &big_one, q); } return big_errno; } bigerr_t big_ceil(a, b, q, r) bignum *a; bignum *b; bignum *q; bignum *r; { int b_eq_qr, sign_diff; if ( ( b_eq_qr = ((b == q) || (b == r)) ) ) { big_set_big(b, &tmp_mul); } sign_diff = a->sign != b->sign; big_trunc(a, b, q, r); if (sign_diff) { return big_errno; } if (!zerop(r)) { if (b_eq_qr) { big_sub(r, &tmp_mul, r); } else { big_sub(r, b, r); } big_add(q, &big_one, q); } return big_errno; } #ifdef UNUSED_CODE /* This one doesn't work to 100%. I was a little braindamaged when I wrote * this, but I'll eventually fix it. Zzzzzzz. */ bigerr_t big_round(a, b, q, r) bignum *a; bignum *b; bignum *q; bignum *r; { int b_eq_qr, b_neg_p, a_sgn_neq_b_sgn; if ( ( b_eq_qr = ((b == q) || (b == r)) ) ) { big_set_big(b, &tmp_round); } b_neg_p = b->sign == BIG_SIGN_MINUS; a_sgn_neq_b_sgn = a->sign != b->sign; big_trunc(a, b, q, r); big_set_big(r, &tmp_add); umul_digit(&tmp_add, 2); if (ucompare_digits(&tmp_add, b) > 0) /* |2 * r| > |b| */ { if (q->sign == BIG_SIGN_0) { if (a_sgn_neq_b_sgn) { big_sub(q, &big_one, q); } else { big_add(q, &big_one, q); } } else { if (q->sign == BIG_SIGN_MINUS) { big_sub(q, &big_one, q); } else { big_add(q, &big_one, q); } } if (b_eq_qr) { if (q->sign == BIG_SIGN_PLUS) { big_sub(r, &tmp_round, r); } else { big_add(r, &tmp_round, r); } } else { if (q->sign == BIG_SIGN_PLUS) { big_sub(r, b, r); } else { big_add(r, b, r); } } } return big_errno; } bigerr_t big_random(a, b) bignum *a; bignum *b; { unsigned long i; int a_sgn = a->sign; if (big_errno != BIG_OK) { return big_errno; } if (zerop(a)) /* a = 0 -> big_random => 0 (special case) */ { *b->dp = 0; b->dgs_used = 1; b->sign = a_sgn; return big_errno; } if (newsize(&tmp_rand.dp, &tmp_rand.dgs_alloc, a->dgs_used + 1, a->dgs_used + 1) != BIG_OK) { return big_errno; } for (i = 0; i <= a->dgs_used; i++) { tmp_rand.dp[i] = rand(); } while (tmp_rand.dp[a->dgs_used] == 0) /* Make sure high digit is non-0 */ { tmp_rand.dp[a->dgs_used] = rand(); } tmp_rand.dgs_used = a->dgs_used + 1; tmp_rand.sign = BIG_SIGN_PLUS; a->sign = BIG_SIGN_PLUS; big_trunc(&tmp_rand, a, &tmp_q, b); /* Dangerous to use tmp_q here... */ a->sign = a_sgn; b->sign = zerop(&tmp_rand) ? BIG_SIGN_0 : a_sgn; return big_errno; } /* UNUSED_CODE */ #endif /* ---------------------------------------------------------------------- * External functions that do not need to know anything about the internal * representation of a bignum. */ #ifdef UNUSED_CODE int big_expt(a, z, x) bignum *a; unsigned long z; bignum *x; { bignum b; big_create(&b); big_set_big(a, &b); big_set_long((long)1, x); while ((z != 0) && (big_errno == BIG_OK)) { while ((z & 0x01) == 0) { z >>= 1; big_mul(&b, &b, &b); } z -= 1; big_mul(x, &b, x); } big_destroy(&b); return big_errno; } /* UNUSED_CODE */ #endif int big_exptmod(a_in, z_in, n, x) bignum *a_in; bignum *z_in; bignum *n; bignum *x; { bignum a, z, b0, b1, b2, dummy; big_create(&a); big_create(&z); big_create(&b0); big_create(&b1); big_create(&b2); big_create(&dummy); big_set_big(a_in, &a); big_set_big(z_in, &z); big_set_long((long)1, x); big_set_long((long)0, &b0); big_set_long((long)1, &b1); big_set_long((long)2, &b2); /* No foolproof testing on big_errno - it really ought to be done */ while ((big_compare(&z, &b0) != 0) && (big_errno == BIG_OK)) { while (big_evenp(&z) && (big_errno == BIG_OK)) { big_trunc(&z, &b2, &z, &dummy); big_mul(&a, &a, &a); big_trunc(&a, n, &dummy, &a); } big_sub(&z, &b1, &z); big_mul(x, &a, x); big_trunc(x, n, &dummy, x); } big_destroy(&dummy); big_destroy(&b2); big_destroy(&b1); big_destroy(&b0); big_destroy(&z); big_destroy(&a); return big_errno; } #ifdef UNUSED_CODE int big_gcd(a, b, g) bignum *a; bignum *b; bignum *g; { bignum a1, b1, tmp; if (big_zerop(b)) { big_abs(a, g); return big_errno; } big_create(&a1); big_create(&b1); big_create(&tmp); big_abs(a, &a1); big_abs(b, &b1); while (!big_zerop(&b1) && (big_errno == BIG_OK)) { big_floor(&a1, &b1, &tmp, &a1); if (big_zerop(&a1)) { break; } big_floor(&b1, &a1, &tmp, &b1); } if (big_zerop(&a1)) { big_set_big(&b1, g); } else { big_set_big(&a1, g); } big_destroy(&tmp); big_destroy(&b1); big_destroy(&a1); return big_errno; } /* UNUSED_CODE */ #endif /* #if (defined (SH_WITH_CLIENT) || defined (SH_WITH_SERVER)) */ #endif