Import imath 1.14

This commit is contained in:
Love Hornquist Astrand
2009-08-17 16:08:12 +02:00
parent 62433c844c
commit dd673af0b6
2 changed files with 223 additions and 224 deletions

View File

@@ -2,7 +2,7 @@
Name: imath.c
Purpose: Arbitrary precision integer arithmetic routines.
Author: M. J. Fromberger <http://spinning-yarns.org/michael/>
Info: $Id: imath.c 645 2008-08-03 04:00:30Z sting $
Info: $Id: imath.c 826 2009-02-11 16:21:04Z sting $
Copyright (C) 2002-2007 Michael J. Fromberger, All Rights Reserved.
@@ -40,7 +40,9 @@
#include <assert.h>
#if DEBUG
#define static
#define STATIC /* public */
#else
#define STATIC static
#endif
/* {{{ Constants */
@@ -58,8 +60,8 @@ const mp_result MP_MINERR = -6;
const mp_sign MP_NEG = 1; /* value is strictly negative */
const mp_sign MP_ZPOS = 0; /* value is non-negative */
static const char *s_unknown_err = "unknown result code";
static const char *s_error_msg[] = {
STATIC const char *s_unknown_err = "unknown result code";
STATIC const char *s_error_msg[] = {
"error code 0",
"boolean true",
"out of memory",
@@ -88,7 +90,7 @@ static const char *s_error_msg[] = {
The use of this table eliminates a dependency upon linkage against
the standard math libraries.
*/
static const double s_log2[] = {
STATIC const double s_log2[] = {
0.000000000, 0.000000000, 1.000000000, 0.630929754, /* 0 1 2 3 */
0.500000000, 0.430676558, 0.386852807, 0.356207187, /* 4 5 6 7 */
0.333333333, 0.315464877, 0.301029996, 0.289064826, /* 8 9 10 11 */
@@ -172,144 +174,144 @@ do{mp_size ua_=MP_USED(X),o_=ua_+ua_;ZERO(MP_DIGITS(Z),o_);\
#if IMATH_TEST
mp_size default_precision = MP_DEFAULT_PREC;
#else
static const mp_size default_precision = MP_DEFAULT_PREC;
STATIC const mp_size default_precision = MP_DEFAULT_PREC;
#endif
/* Minimum number of digits to invoke recursive multiply */
#if IMATH_TEST
mp_size multiply_threshold = MP_MULT_THRESH;
#else
static const mp_size multiply_threshold = MP_MULT_THRESH;
STATIC const mp_size multiply_threshold = MP_MULT_THRESH;
#endif
/* }}} */
/* Allocate a buffer of (at least) num digits, or return
NULL if that couldn't be done. */
static mp_digit *s_alloc(mp_size num);
STATIC mp_digit *s_alloc(mp_size num);
/* Release a buffer of digits allocated by s_alloc(). */
static void s_free(void *ptr);
STATIC void s_free(void *ptr);
/* Insure that z has at least min digits allocated, resizing if
necessary. Returns true if successful, false if out of memory. */
static int s_pad(mp_int z, mp_size min);
STATIC int s_pad(mp_int z, mp_size min);
/* Fill in a "fake" mp_int on the stack with a given value */
static void s_fake(mp_int z, mp_small value, mp_digit vbuf[]);
STATIC void s_fake(mp_int z, mp_small value, mp_digit vbuf[]);
/* Compare two runs of digits of given length, returns <0, 0, >0 */
static int s_cdig(mp_digit *da, mp_digit *db, mp_size len);
STATIC int s_cdig(mp_digit *da, mp_digit *db, mp_size len);
/* Pack the unsigned digits of v into array t */
static int s_vpack(mp_small v, mp_digit t[]);
STATIC int s_vpack(mp_small v, mp_digit t[]);
/* Compare magnitudes of a and b, returns <0, 0, >0 */
static int s_ucmp(mp_int a, mp_int b);
STATIC int s_ucmp(mp_int a, mp_int b);
/* Compare magnitudes of a and v, returns <0, 0, >0 */
static int s_vcmp(mp_int a, mp_small v);
STATIC int s_vcmp(mp_int a, mp_small v);
/* Unsigned magnitude addition; assumes dc is big enough.
Carry out is returned (no memory allocated). */
static mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc,
STATIC mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc,
mp_size size_a, mp_size size_b);
/* Unsigned magnitude subtraction. Assumes dc is big enough. */
static void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc,
STATIC void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc,
mp_size size_a, mp_size size_b);
/* Unsigned recursive multiplication. Assumes dc is big enough. */
static int s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc,
STATIC int s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc,
mp_size size_a, mp_size size_b);
/* Unsigned magnitude multiplication. Assumes dc is big enough. */
static void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc,
STATIC void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc,
mp_size size_a, mp_size size_b);
/* Unsigned recursive squaring. Assumes dc is big enough. */
static int s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a);
STATIC int s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a);
/* Unsigned magnitude squaring. Assumes dc is big enough. */
static void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a);
STATIC void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a);
/* Single digit addition. Assumes a is big enough. */
static void s_dadd(mp_int a, mp_digit b);
STATIC void s_dadd(mp_int a, mp_digit b);
/* Single digit multiplication. Assumes a is big enough. */
static void s_dmul(mp_int a, mp_digit b);
STATIC void s_dmul(mp_int a, mp_digit b);
/* Single digit multiplication on buffers; assumes dc is big enough. */
static void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc,
STATIC void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc,
mp_size size_a);
/* Single digit division. Replaces a with the quotient,
returns the remainder. */
static mp_digit s_ddiv(mp_int a, mp_digit b);
STATIC mp_digit s_ddiv(mp_int a, mp_digit b);
/* Quick division by a power of 2, replaces z (no allocation) */
static void s_qdiv(mp_int z, mp_size p2);
STATIC void s_qdiv(mp_int z, mp_size p2);
/* Quick remainder by a power of 2, replaces z (no allocation) */
static void s_qmod(mp_int z, mp_size p2);
STATIC void s_qmod(mp_int z, mp_size p2);
/* Quick multiplication by a power of 2, replaces z.
Allocates if necessary; returns false in case this fails. */
static int s_qmul(mp_int z, mp_size p2);
STATIC int s_qmul(mp_int z, mp_size p2);
/* Quick subtraction from a power of 2, replaces z.
Allocates if necessary; returns false in case this fails. */
static int s_qsub(mp_int z, mp_size p2);
STATIC int s_qsub(mp_int z, mp_size p2);
/* Return maximum k such that 2^k divides z. */
static int s_dp2k(mp_int z);
STATIC int s_dp2k(mp_int z);
/* Return k >= 0 such that z = 2^k, or -1 if there is no such k. */
static int s_isp2(mp_int z);
STATIC int s_isp2(mp_int z);
/* Set z to 2^k. May allocate; returns false in case this fails. */
static int s_2expt(mp_int z, mp_small k);
STATIC int s_2expt(mp_int z, mp_small k);
/* Normalize a and b for division, returns normalization constant */
static int s_norm(mp_int a, mp_int b);
STATIC int s_norm(mp_int a, mp_int b);
/* Compute constant mu for Barrett reduction, given modulus m, result
replaces z, m is untouched. */
static mp_result s_brmu(mp_int z, mp_int m);
STATIC mp_result s_brmu(mp_int z, mp_int m);
/* Reduce a modulo m, using Barrett's algorithm. */
static int s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2);
STATIC int s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2);
/* Modular exponentiation, using Barrett reduction */
static mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c);
STATIC mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c);
/* Unsigned magnitude division. Assumes |a| > |b|. Allocates
temporaries; overwrites a with quotient, b with remainder. */
static mp_result s_udiv(mp_int a, mp_int b);
STATIC mp_result s_udiv(mp_int a, mp_int b);
/* Compute the number of digits in radix r required to represent the
given value. Does not account for sign flags, terminators, etc. */
static int s_outlen(mp_int z, mp_size r);
STATIC int s_outlen(mp_int z, mp_size r);
/* Guess how many digits of precision will be needed to represent a
radix r value of the specified number of digits. Returns a value
guaranteed to be no smaller than the actual number required. */
static mp_size s_inlen(int len, mp_size r);
STATIC mp_size s_inlen(int len, mp_size r);
/* Convert a character to a digit value in radix r, or
-1 if out of range */
static int s_ch2val(char c, int r);
STATIC int s_ch2val(char c, int r);
/* Convert a digit value to a character */
static char s_val2ch(int v, int caps);
STATIC char s_val2ch(int v, int caps);
/* Take 2's complement of a buffer in place */
static void s_2comp(unsigned char *buf, int len);
STATIC void s_2comp(unsigned char *buf, int len);
/* Convert a value to binary, ignoring sign. On input, *limpos is the
bound on how many bytes should be written to buf; on output, *limpos
is set to the number of bytes actually written. */
static mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad);
STATIC mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad);
#if DEBUG
/* Dump a representation of the mp_int to standard output */
@@ -2115,7 +2117,7 @@ const char *mp_error_string(mp_result res)
/* {{{ s_alloc(num) */
static mp_digit *s_alloc(mp_size num)
STATIC mp_digit *s_alloc(mp_size num)
{
mp_digit *out = malloc(num * sizeof(mp_digit));
@@ -2137,7 +2139,7 @@ static mp_digit *s_alloc(mp_size num)
/* {{{ s_realloc(old, osize, nsize) */
static mp_digit *s_realloc(mp_digit *old, mp_size osize, mp_size nsize)
STATIC mp_digit *s_realloc(mp_digit *old, mp_size osize, mp_size nsize)
{
#if DEBUG > 1
mp_digit *new = s_alloc(nsize);
@@ -2159,7 +2161,7 @@ static mp_digit *s_realloc(mp_digit *old, mp_size osize, mp_size nsize)
/* {{{ s_free(ptr) */
static void s_free(void *ptr)
STATIC void s_free(void *ptr)
{
free(ptr);
}
@@ -2168,7 +2170,7 @@ static void s_free(void *ptr)
/* {{{ s_pad(z, min) */
static int s_pad(mp_int z, mp_size min)
STATIC int s_pad(mp_int z, mp_size min)
{
if(MP_ALLOC(z) < min) {
mp_size nsize = ROUND_PREC(min);
@@ -2194,7 +2196,7 @@ static int s_pad(mp_int z, mp_size min)
/* {{{ s_fake(z, value, vbuf) */
static void s_fake(mp_int z, mp_small value, mp_digit vbuf[])
STATIC void s_fake(mp_int z, mp_small value, mp_digit vbuf[])
{
mp_size uv = (mp_size) s_vpack(value, vbuf);
@@ -2208,7 +2210,7 @@ static void s_fake(mp_int z, mp_small value, mp_digit vbuf[])
/* {{{ s_cdig(da, db, len) */
static int s_cdig(mp_digit *da, mp_digit *db, mp_size len)
STATIC int s_cdig(mp_digit *da, mp_digit *db, mp_size len)
{
mp_digit *dat = da + len - 1, *dbt = db + len - 1;
@@ -2226,7 +2228,7 @@ static int s_cdig(mp_digit *da, mp_digit *db, mp_size len)
/* {{{ s_vpack(v, t[]) */
static int s_vpack(mp_small v, mp_digit t[])
STATIC int s_vpack(mp_small v, mp_digit t[])
{
mp_usmall uv = (mp_usmall) ((v < 0) ? -v : v);
int ndig = 0;
@@ -2248,7 +2250,7 @@ static int s_vpack(mp_small v, mp_digit t[])
/* {{{ s_ucmp(a, b) */
static int s_ucmp(mp_int a, mp_int b)
STATIC int s_ucmp(mp_int a, mp_int b)
{
mp_size ua = MP_USED(a), ub = MP_USED(b);
@@ -2264,7 +2266,7 @@ static int s_ucmp(mp_int a, mp_int b)
/* {{{ s_vcmp(a, v) */
static int s_vcmp(mp_int a, mp_small v)
STATIC int s_vcmp(mp_int a, mp_small v)
{
mp_digit vdig[MP_VALUE_DIGITS(v)];
int ndig = 0;
@@ -2284,7 +2286,7 @@ static int s_vcmp(mp_int a, mp_small v)
/* {{{ s_uadd(da, db, dc, size_a, size_b) */
static mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc,
STATIC mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc,
mp_size size_a, mp_size size_b)
{
mp_size pos;
@@ -2319,7 +2321,7 @@ static mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc,
/* {{{ s_usub(da, db, dc, size_a, size_b) */
static void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc,
STATIC void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc,
mp_size size_a, mp_size size_b)
{
mp_size pos;
@@ -2354,7 +2356,7 @@ static void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc,
/* {{{ s_kmul(da, db, dc, size_a, size_b) */
static int s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc,
STATIC int s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc,
mp_size size_a, mp_size size_b)
{
mp_size bot_size;
@@ -2443,7 +2445,7 @@ static int s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc,
/* {{{ s_umul(da, db, dc, size_a, size_b) */
static void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc,
STATIC void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc,
mp_size size_a, mp_size size_b)
{
mp_size a, b;
@@ -2472,7 +2474,7 @@ static void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc,
/* {{{ s_ksqr(da, dc, size_a) */
static int s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a)
STATIC int s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a)
{
if(multiply_threshold && size_a > multiply_threshold) {
mp_size bot_size = (size_a + 1) / 2;
@@ -2529,7 +2531,7 @@ static int s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a)
/* {{{ s_usqr(da, dc, size_a) */
static void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a)
STATIC void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a)
{
mp_size i, j;
mp_word w;
@@ -2585,7 +2587,7 @@ static void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a)
/* {{{ s_dadd(a, b) */
static void s_dadd(mp_int a, mp_digit b)
STATIC void s_dadd(mp_int a, mp_digit b)
{
mp_word w = 0;
mp_digit *da = MP_DIGITS(a);
@@ -2612,7 +2614,7 @@ static void s_dadd(mp_int a, mp_digit b)
/* {{{ s_dmul(a, b) */
static void s_dmul(mp_int a, mp_digit b)
STATIC void s_dmul(mp_int a, mp_digit b)
{
mp_word w = 0;
mp_digit *da = MP_DIGITS(a);
@@ -2635,7 +2637,7 @@ static void s_dmul(mp_int a, mp_digit b)
/* {{{ s_dbmul(da, b, dc, size_a) */
static void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, mp_size size_a)
STATIC void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, mp_size size_a)
{
mp_word w = 0;
@@ -2655,7 +2657,7 @@ static void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, mp_size size_a)
/* {{{ s_ddiv(da, d, dc, size_a) */
static mp_digit s_ddiv(mp_int a, mp_digit b)
STATIC mp_digit s_ddiv(mp_int a, mp_digit b)
{
mp_word w = 0, qdigit;
mp_size ua = MP_USED(a);
@@ -2683,7 +2685,7 @@ static mp_digit s_ddiv(mp_int a, mp_digit b)
/* {{{ s_qdiv(z, p2) */
static void s_qdiv(mp_int z, mp_size p2)
STATIC void s_qdiv(mp_int z, mp_size p2)
{
mp_size ndig = p2 / MP_DIGIT_BIT, nbits = p2 % MP_DIGIT_BIT;
mp_size uz = MP_USED(z);
@@ -2730,7 +2732,7 @@ static void s_qdiv(mp_int z, mp_size p2)
/* {{{ s_qmod(z, p2) */
static void s_qmod(mp_int z, mp_size p2)
STATIC void s_qmod(mp_int z, mp_size p2)
{
mp_size start = p2 / MP_DIGIT_BIT + 1, rest = p2 % MP_DIGIT_BIT;
mp_size uz = MP_USED(z);
@@ -2747,7 +2749,7 @@ static void s_qmod(mp_int z, mp_size p2)
/* {{{ s_qmul(z, p2) */
static int s_qmul(mp_int z, mp_size p2)
STATIC int s_qmul(mp_int z, mp_size p2)
{
mp_size uz, need, rest, extra, i;
mp_digit *from, *to, d;
@@ -2815,7 +2817,7 @@ static int s_qmul(mp_int z, mp_size p2)
/* Compute z = 2^p2 - |z|; requires that 2^p2 >= |z|
The sign of the result is always zero/positive.
*/
static int s_qsub(mp_int z, mp_size p2)
STATIC int s_qsub(mp_int z, mp_size p2)
{
mp_digit hi = (1 << (p2 % MP_DIGIT_BIT)), *zp;
mp_size tdig = (p2 / MP_DIGIT_BIT), pos;
@@ -2846,7 +2848,7 @@ static int s_qsub(mp_int z, mp_size p2)
/* {{{ s_dp2k(z) */
static int s_dp2k(mp_int z)
STATIC int s_dp2k(mp_int z)
{
int k = 0;
mp_digit *dp = MP_DIGITS(z), d;
@@ -2872,7 +2874,7 @@ static int s_dp2k(mp_int z)
/* {{{ s_isp2(z) */
static int s_isp2(mp_int z)
STATIC int s_isp2(mp_int z)
{
mp_size uz = MP_USED(z), k = 0;
mp_digit *dz = MP_DIGITS(z), d;
@@ -2898,7 +2900,7 @@ static int s_isp2(mp_int z)
/* {{{ s_2expt(z, k) */
static int s_2expt(mp_int z, mp_small k)
STATIC int s_2expt(mp_int z, mp_small k)
{
mp_size ndig, rest;
mp_digit *dz;
@@ -2921,7 +2923,7 @@ static int s_2expt(mp_int z, mp_small k)
/* {{{ s_norm(a, b) */
static int s_norm(mp_int a, mp_int b)
STATIC int s_norm(mp_int a, mp_int b)
{
mp_digit d = b->digits[MP_USED(b) - 1];
int k = 0;
@@ -2944,7 +2946,7 @@ static int s_norm(mp_int a, mp_int b)
/* {{{ s_brmu(z, m) */
static mp_result s_brmu(mp_int z, mp_int m)
STATIC mp_result s_brmu(mp_int z, mp_int m)
{
mp_size um = MP_USED(m) * 2;
@@ -2959,7 +2961,7 @@ static mp_result s_brmu(mp_int z, mp_int m)
/* {{{ s_reduce(x, m, mu, q1, q2) */
static int s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2)
STATIC int s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2)
{
mp_size um = MP_USED(m), umb_p1, umb_m1;
@@ -3008,7 +3010,7 @@ static int s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2)
/* Perform modular exponentiation using Barrett's method, where mu is
the reduction constant for m. Assumes a < m, b > 0. */
static mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c)
STATIC mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c)
{
mp_digit *db, *dbt, umu, d;
mpz_t temp[3];
@@ -3088,7 +3090,7 @@ static mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c)
/* Precondition: a >= b and b > 0
Postcondition: a' = a / b, b' = a % b
*/
static mp_result s_udiv(mp_int a, mp_int b)
STATIC mp_result s_udiv(mp_int a, mp_int b)
{
mpz_t q, r, t;
mp_size ua, ub, qpos = 0;
@@ -3139,10 +3141,7 @@ static mp_result s_udiv(mp_int a, mp_int b)
qdigit = pfx / btop;
if(qdigit > MP_DIGIT_MAX) {
if(qdigit & MP_DIGIT_MAX)
qdigit = MP_DIGIT_MAX;
else
qdigit = 1;
}
s_dbmul(MP_DIGITS(b), (mp_digit) qdigit, t.digits, ub);
@@ -3184,7 +3183,7 @@ static mp_result s_udiv(mp_int a, mp_int b)
/* {{{ s_outlen(z, r) */
static int s_outlen(mp_int z, mp_size r)
STATIC int s_outlen(mp_int z, mp_size r)
{
mp_result bits;
double raw;
@@ -3201,7 +3200,7 @@ static int s_outlen(mp_int z, mp_size r)
/* {{{ s_inlen(len, r) */
static mp_size s_inlen(int len, mp_size r)
STATIC mp_size s_inlen(int len, mp_size r)
{
double raw = (double)len / s_log2[r];
mp_size bits = (mp_size)(raw + 0.5);
@@ -3213,7 +3212,7 @@ static mp_size s_inlen(int len, mp_size r)
/* {{{ s_ch2val(c, r) */
static int s_ch2val(char c, int r)
STATIC int s_ch2val(char c, int r)
{
int out;
@@ -3231,7 +3230,7 @@ static int s_ch2val(char c, int r)
/* {{{ s_val2ch(v, caps) */
static char s_val2ch(int v, int caps)
STATIC char s_val2ch(int v, int caps)
{
assert(v >= 0);
@@ -3251,7 +3250,7 @@ static char s_val2ch(int v, int caps)
/* {{{ s_2comp(buf, len) */
static void s_2comp(unsigned char *buf, int len)
STATIC void s_2comp(unsigned char *buf, int len)
{
int i;
unsigned short s = 1;
@@ -3273,7 +3272,7 @@ static void s_2comp(unsigned char *buf, int len)
/* {{{ s_tobin(z, buf, *limpos) */
static mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad)
STATIC mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad)
{
mp_size uz;
mp_digit *dz;