?? bigdigits.c
字號:
/* Create some temp storage for int values */
n = strlen(s);
if (0 == n) return 0;
newlen = uiceil(n * 0.5); /* log(16)/log(256)=0.5 */
newdigits = calloc(newlen, 1);
if (!newdigits) mpFail("mpConvFromHex: Not enough memory: " __FILE__);
/* Work through zero-terminated string */
for (i = 0; s[i]; i++)
{
t = s[i];
if ((t >= '0') && (t <= '9')) t = (t - '0');
else if ((t >= 'a') && (t <= 'f')) t = (t - 'a' + 10);
else if ((t >= 'A') && (t <= 'F')) t = (t - 'A' + 10);
else continue;
for (j = newlen; j > 0; j--)
{
t += (unsigned long)newdigits[j-1] << 4;
newdigits[j-1] = (unsigned char)(t & 0xFF);
t >>= 8;
}
}
/* Convert bytes to big digits */
n = mpConvFromOctets(a, ndigits, newdigits, newlen);
/* Clean up */
free(newdigits);
return n;
}
int mpModulo(DIGIT_T r[], const DIGIT_T u[], size_t udigits,
DIGIT_T v[], size_t vdigits)
{
/* Computes r = u mod v
where r, v are multiprecision integers of length vdigits
and u is a multiprecision integer of length udigits.
r may overlap v.
Note that r here is only vdigits long,
whereas in mpDivide it is udigits long.
Use remainder from mpDivide function.
Updated Version 2 to use malloc.
*/
DIGIT_T *qq, *rr;
size_t nn = max(udigits, vdigits);
qq = mpAlloc(udigits);
rr = mpAlloc(nn);
/* rr[nn] = u mod v */
mpDivide(qq, rr, u, udigits, v, vdigits);
/* Final r is only vdigits long */
mpSetEqual(r, rr, vdigits);
mpSetZero(rr, udigits);
mpSetZero(qq, udigits);
mpFree(&rr);
mpFree(&qq);
return 0;
}
int mpModMult(DIGIT_T a[], const DIGIT_T x[], const DIGIT_T y[],
const DIGIT_T m[], size_t ndigits)
{ /* Computes a = (x * y) mod m */
/* Updated Version 2 to use malloc. */
/* Double-length temp variable p */
DIGIT_T *p;
DIGIT_T *tm;
p = mpAlloc(ndigits * 2);
tm = mpAlloc(ndigits);
mpSetEqual(tm, m, ndigits);
/* Calc p[2n] = x * y */
mpMultiply(p, x, y, ndigits);
/* Then modulo */
mpModulo(a, p, ndigits * 2, tm, ndigits);
mpSetZero(p, ndigits * 2);
mpSetZero(tm, ndigits);
mpFree(&p);
mpFree(&tm);
return 0;
}
/* mpModExp */
/* Version 2: added new funcs with temps
to avoid having to alloc and free repeatedly
and added mpSquare function for slight speed improvement
*/
static int moduloTemp(DIGIT_T r[], const DIGIT_T u[], size_t udigits,
DIGIT_T v[], size_t vdigits, DIGIT_T tqq[], DIGIT_T trr[])
{
/* Calculates r = u mod v
where r, v are multiprecision integers of length vdigits
and u is a multiprecision integer of length udigits.
r may overlap v.
Same as mpModulo without allocs & free.
Requires temp mp's tqq and trr of length udigits each
*/
mpDivide(tqq, trr, u, udigits, v, vdigits);
/* Final r is only vdigits long */
mpSetEqual(r, trr, vdigits);
return 0;
}
static int modMultTemp(DIGIT_T a[], const DIGIT_T x[], const DIGIT_T y[],
DIGIT_T m[], size_t ndigits,
DIGIT_T temp[], DIGIT_T tqq[], DIGIT_T trr[])
{ /* Computes a = (x * y) mod m */
/* Requires 3 x temp mp's of length 2 * ndigits each */
/* Calc p[2n] = x * y */
mpMultiply(temp, x, y, ndigits);
/* Then modulo m */
moduloTemp(a, temp, ndigits * 2, m, ndigits, tqq, trr);
return 0;
}
static int modSquareTemp(DIGIT_T a[], const DIGIT_T x[],
DIGIT_T m[], size_t ndigits,
DIGIT_T temp[], DIGIT_T tqq[], DIGIT_T trr[])
{ /* Computes a = (x * x) mod m */
/* Requires 3 x temp mp's of length 2 * ndigits each */
/* Calc p[2n] = x^2 */
mpSquare(temp, x, ndigits);
/* Then modulo m */
moduloTemp(a, temp, ndigits * 2, m, ndigits, tqq, trr);
return 0;
}
int mpModExp(DIGIT_T yout[], const DIGIT_T x[],
const DIGIT_T e[], const DIGIT_T m[], size_t ndigits)
{ /* Computes y = x^e mod m */
/* Binary left-to-right method */
DIGIT_T mask;
size_t n, nn;
DIGIT_T *t1, *t2, *t3, *tm, *y;
if (ndigits == 0) return -1;
/* Create some temps */
nn = ndigits * 2;
t1 = mpAlloc(nn);
t2 = mpAlloc(nn);
t3 = mpAlloc(nn);
tm = mpAlloc(ndigits);
y = mpAlloc(ndigits);
mpSetEqual(tm, m, ndigits);
/* Find second-most significant bit in e */
n = mpSizeof(e, ndigits);
for (mask = HIBITMASK; mask > 0; mask >>= 1)
{
if (e[n-1] & mask)
break;
}
mpNEXTBITMASK(mask, n);
/* Set y = x */
mpSetEqual(y, x, ndigits);
/* For bit j = k-2 downto 0 */
while (n)
{
/* Square y = y * y mod n */
modSquareTemp(y, y, tm, ndigits, t1, t2, t3);
if (e[n-1] & mask)
{ /* if e(j) == 1 then multiply
y = y * x mod n */
modMultTemp(y, y, x, tm, ndigits, t1, t2, t3);
}
/* Move to next bit */
mpNEXTBITMASK(mask, n);
}
/* return y */
mpSetEqual(yout, y, ndigits);
mpSetZero(t1, nn);
mpSetZero(t2, nn);
mpSetZero(t3, nn);
mpSetZero(tm, ndigits);
mpSetZero(y, ndigits);
mpFree(&t1);
mpFree(&t2);
mpFree(&t3);
mpFree(&tm);
mpFree(&y);
return 0;
}
int mpModInv(DIGIT_T inv[], const DIGIT_T u[], const DIGIT_T v[], size_t ndigits)
{ /* Computes inv = u^(-1) mod v */
/* Ref: Knuth Algorithm X Vol 2 p 342
ignoring u2, v2, t2
and avoiding negative numbers.
Updated Version 2 to use malloc
and to return non-zero if inverse undefined.
*/
/* Hard-coded arrays removed in version 2:
//DIGIT_T u1[MAX_DIG_LEN], u3[MAX_DIG_LEN], v1[MAX_DIG_LEN], v3[MAX_DIG_LEN];
//DIGIT_T t1[MAX_DIG_LEN], t3[MAX_DIG_LEN], q[MAX_DIG_LEN];
//DIGIT_T w[2*MAX_DIG_LEN];
*/
DIGIT_T *u1, *u3, *v1, *v3, *t1, *t3, *q, *w;
int bIterations;
int result;
/* Allocate temp storage */
u1 = mpAlloc(ndigits);
u3 = mpAlloc(ndigits);
v1 = mpAlloc(ndigits);
v3 = mpAlloc(ndigits);
t1 = mpAlloc(ndigits);
t3 = mpAlloc(ndigits);
q = mpAlloc(ndigits);
w = mpAlloc(2 * ndigits);
/* Step X1. Initialise */
mpSetDigit(u1, 1, ndigits); /* u1 = 1 */
mpSetEqual(u3, u, ndigits); /* u3 = u */
mpSetZero(v1, ndigits); /* v1 = 0 */
mpSetEqual(v3, v, ndigits); /* v3 = v */
bIterations = 1; /* Remember odd/even iterations */
while (!mpIsZero(v3, ndigits)) /* Step X2. Loop while v3 != 0 */
{ /* Step X3. Divide and "Subtract" */
mpDivide(q, t3, u3, ndigits, v3, ndigits);
/* q = u3 / v3, t3 = u3 % v3 */
mpMultiply(w, q, v1, ndigits); /* w = q * v1 */
mpAdd(t1, u1, w, ndigits); /* t1 = u1 + w */
/* Swap u1 = v1; v1 = t1; u3 = v3; v3 = t3 */
mpSetEqual(u1, v1, ndigits);
mpSetEqual(v1, t1, ndigits);
mpSetEqual(u3, v3, ndigits);
mpSetEqual(v3, t3, ndigits);
bIterations = -bIterations;
}
if (bIterations < 0)
mpSubtract(inv, v, u1, ndigits); /* inv = v - u1 */
else
mpSetEqual(inv, u1, ndigits); /* inv = u1 */
/* Make sure u3 = gcd(u,v) == 1 */
if (mpShortCmp(u3, 1, ndigits) != 0)
{
result = 1;
mpSetZero(inv, ndigits);
}
else
result = 0;
/* Clear up */
mpSetZero(u1, ndigits);
mpSetZero(v1, ndigits);
mpSetZero(t1, ndigits);
mpSetZero(u3, ndigits);
mpSetZero(v3, ndigits);
mpSetZero(t3, ndigits);
mpSetZero(q, ndigits);
mpSetZero(w, 2*ndigits);
mpFree(&u1);
mpFree(&v1);
mpFree(&t1);
mpFree(&u3);
mpFree(&v3);
mpFree(&t3);
mpFree(&q);
mpFree(&w);
return result;
}
int mpGcd(DIGIT_T g[], const DIGIT_T x[], const DIGIT_T y[], size_t ndigits)
{
/* Computes g = gcd(x, y) */
/* Ref: Schneier */
/* NB This function requires temp storage */
/* Updated Version 2 to use malloc.
//DIGIT_T yy[MAX_DIG_LEN], xx[MAX_DIG_LEN];
*/
DIGIT_T *yy, *xx;
yy = mpAlloc(ndigits);
xx = mpAlloc(ndigits);
mpSetZero(yy, ndigits);
mpSetZero(xx, ndigits);
mpSetEqual(xx, x, ndigits);
mpSetEqual(yy, y, ndigits);
mpSetEqual(g, yy, ndigits); /* g = y */
while (!mpIsZero(xx, ndigits)) /* while (x > 0) */
{
mpSetEqual(g, xx, ndigits); /* g = x */
mpModulo(xx, yy, ndigits, xx, ndigits); /* x = y mod x */
mpSetEqual(yy, g, ndigits); /* y = g; */
}
mpSetZero(xx, ndigits);
mpSetZero(yy, ndigits);
mpFree(&xx);
mpFree(&yy);
return 0; /* gcd is in g */
}
int mpSqrt(DIGIT_T s[], const DIGIT_T x[], size_t ndigits)
/* Computes integer square root s = floor(sqrt(x)) */
/* [Added v2.1] Based on lsqrt() function */
{
DIGIT_T *v0, *q0, *x0, *x1, *t;
/* if (x <= 1) return x */
if (mpShortCmp(x, 1, ndigits) <= 0)
{
mpSetEqual(s, x, ndigits);
return 0;
}
/* Allocate temp storage */
v0 = mpAlloc(ndigits);
q0 = mpAlloc(ndigits);
x0 = mpAlloc(ndigits);
x1 = mpAlloc(ndigits);
t = mpAlloc(ndigits);
/* v0 = x - 1 */
mpShortSub(v0, x, 1, ndigits);
/* x0 = x/2 */
mpShortDiv(x0, x, 2, ndigits);
while (1)
{
/* q0 = v0/x0 */
mpDivide(q0, t, v0, ndigits, x0, ndigits);
/* x1 = (x0 + q0)/2 */
mpAdd(t, x0, q0, ndigits);
mpShortDiv(x1, t, 2, ndigits);
/* if (q0 >= x0) break */
if (mpCompare(q0, x0, ndigits) >= 0)
break;
/* x0 = x1 */
mpSetEqual(x0, x1, ndigits);
}
/* return x1 */
mpSetEqual(s, x1, ndigits);
mpFree(&v0);
mpFree(&q0);
mpFree(&x0);
mpFree(&x1);
mpFree(&t);
return 0;
}
/* mpIsPrime: Changes in Version 2:
Added mpAlloc for dynamic allocation
Increased no of small primes
Broke out mpRabinMiller() as a separate function
*/
static DIGIT_T SMALL_PRIMES[] = {
3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43,
47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101,
103, 107, 109, 113,
127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
467, 479, 487, 491, 499, 503, 509, 521, 523, 541,
547, 557, 563, 569, 571, 577, 587, 593, 599, 601,
607, 613, 617, 619, 631, 641, 643, 647, 653, 659,
661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
739, 743, 751, 757, 761, 769, 773, 787, 797, 809,
811, 821, 823, 827, 829, 839, 853, 857, 859, 863,
877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
947, 953, 967, 971, 977, 983, 991, 997,
};
#define N_SMALL_PRIMES (sizeof(SMALL_PRIMES)/sizeof(DIGIT_T))
int mpIsPrime(const DIGIT_T w[], size_t ndigits, size_t t)
{ /* Returns true if w > 2 is a probable prime */
/* Version 2: split out mpRabinMiller. */
size_t i;
/* Check the obvious */
if (mpISEVEN(w, ndigits))
return 0;
/* First check for small primes, unless we could be one ourself */
if (mpShortCmp(w, SMALL_PRIMES[N_SMALL_PRIMES-1], ndigits) > 0)
{
for (i = 0; i < N_SMALL_PRIMES; i++)
{
if (mpShortMod(w, SMALL_PRIMES[i], ndigits) == 0)
return 0; /* Failed */
}
}
else
{ /* w is a small number, so check directly */
for (i = 0; i < N_SMALL_PRIMES; i++)
{
if (mpShortCmp(w, SMALL_PRIMES[i], ndigits) == 0)
return 1; /* w is a small prime */
}
return 0; /* w is not a small prime */
}
return mpRabinMiller(w, ndigits, t);
}
/* Local, simple rng functions used in Rabin-Miller */
static void rand_seed(DIGIT_T seed);
static DIGIT_T rand_between(DIGIT_T lower, DIGIT_T upper);
int mpRabinMiller(const DIGIT_T w[], size_t ndigits, size_t t)
{
/* Returns true (1) if w is a probable prime using the
Rabin-Miller Probabilistic Primality Test.
Carries out t iterations specified by user.
Ref: FIPS-186-2 Appendix 2 Section 2.1.
Also Schneier 2nd ed p 260 & Knuth Vol 2, p 379
and ANSI 9.42-2003 Annex B.1.1.
DSS Standard and ANSI 9.42 recommend using t >= 50
for probability of error less than or equal to 2^-100.
Ferguson & Schneier recommend t = 64 for prob error < 2^-128
In practice, most random composites are caught in the first
round or two and so specifying a large t will only affect
the final check.
[v2.1] Updated range of bases from [2, N-1] to [2, N-2]
See ANSI 9.42-2003 Annex F.1.1 `Range of bases in Miller-Rabin test'
(NB this does not impact existing implementations because N-1
is unlikely to be chosen as a base).
*/
/* Temp big digits */
DIGIT_T *m, *a, *b, *z, *w1, *j;
DIGIT_T maxrand;
int failed, isprime;
size_t i;
/* Catch w <= 1 */
if (mpShortCmp(w, 1, ndigits) <= 0) return 0;
/* Allocate temp storage */
m = mpAlloc(ndigits);
a = mpAlloc(ndigits);
b = mpAlloc(ndigits);
z = mpAlloc(ndigits);
w1 = mpAlloc(ndigits);
j = mpAlloc(ndigits);
/* Seed the simple RNG for later on */
rand_seed((DIGIT_T)time(NULL));
/* Rabin-Miller from FIPS-186-2 Appendix 2.
Step 1. Set i = 1 [but do # tests requested by user].
Step 2. Find a and m where w = 1 + (2^a)m
m is odd and 2^a is largest power of 2 dividing w - 1
*/
mpShortSub(w1, w, 1, ndigits); /* Store w1 = w - 1 */
mpSetEqual(m, w1, ndigits); /* Set m = w - 1 */
/* for (a = 0; iseven(m); a++) */
for (mpSetZero(a, ndigits); mpISEVEN(m, ndigits);
mpShortAdd(a, a, 1, ndigits))
{ /* Divide by 2 until m is odd */
mpShiftRight(m, m, 1, ndigits);
}
/* assert((1 << a) * m + 1 == w); */
/* Catch a small w */
if (mpSizeof(w, ndigits) == 1)
maxrand = w[0] - 2; /* [v2.1] changed 1 to 2 */
else
maxrand = MAX_DIGIT;
isprime = 1;
for (i = 0; i < t; i++)
{
failed = 1; /* Assume fail unless passed in loop */
/* Step 3. Generate random integer b, 1 < b < w */
/* [v2.1] changed to 1 < b < w-1 (see ANSI X9.42-2003 Annex B.1.1) */
mpSetZero(b, ndigits);
do
{
b[0] = rand_between(2, maxrand);
} while (mpCompare(b, w, ndigits) >= 0);
/* assert(1 < b && b < w); */
/* Step 4. Set j = 0 and z = b^m mod w */
mpSetZero(j, ndigits);
mpModExp(z, b, m, w, ndigits);
do
{
/* Step 5. If j = 0 and z = 1, or if z = w - 1 */
/* i.e. if ((j == 0 && z == 1) || (z == w - 1)) */
if ((mpIsZero(j, ndigits)
&& mpShortCmp(z, 1, ndigits) == 0)
|| (mpCompare(z, w1, ndigits) == 0))
{ /* Passes on this loop - go to Step 9 */
failed = 0;
break;
}
/* Step 6. If j > 0 and z = 1 */
if (!mpIsZero(j, ndigits)
&& (mpShortCmp(z, 1, ndigits) == 0))
{ /* Fails - go to Step 8 */
failed = 1;
break;
}
/* Step 7. j = j + 1. If j < a set z = z^2 mod w */
mpShortAdd(j, j, 1, ndigits);
if (mpCompare(j, a, ndigits) < 0)
mpModMult(z, z, z, w, ndigits);
/* Loop: if j < a go to Step 5 */
} while (mpCompare(j, a, ndigits) < 0);
if (failed)
{ /* Step 8. Not a prime - stop */
isprime = 0;
break;
}
} /* Step 9. Go to Step 3 until i >= n */
/* Else, if i = n, w is probably prime => success */
/* Clean up */
mpSetZero(m, ndigits);
mpSetZero(a, ndigits);
mpSetZero(b, ndigits);
mpSetZero(z, ndigits);
mpSetZero(w1, ndigits);
mpSetZero(j, ndigits);
mpFree(&m);
mpFree(&a);
mpFree(&b);
mpFree(&z);
mpFree(&w1);
mpFree(&j);
return isprime;
}
/* Internal functions used for "simple" random numbers */
static void rand_seed(DIGIT_T seed)
{
srand(seed);
}
static DIGIT_T rand_between(DIGIT_T lower, DIGIT_T upper)
/* Returns a single pseudo-random digit between lower and upper.
Uses rand(). Assumes srand() already called. */
{
DIGIT_T d, range;
unsigned char *bp;
int i, nbits;
DIGIT_T mask;
if (upper <= lower) return lower;
range = upper - lower;
do
{
/* Generate a random DIGIT byte-by-byte using rand() */
bp = (unsigned char *)&d;
for (i = 0; i < sizeof(DIGIT_T); i++)
{
bp[i] = rand() & 0xFF;
}
/* Trim to next highest bit above required range */
mask = HIBITMASK;
for (nbits = BITS_PER_DIGIT; nbits > 0; nbits--, mask >>= 1)
{
if (range & mask)
break;
}
if (nbits < BITS_PER_DIGIT)
{
mask <<= 1;
mask--;
}
else
mask = MAX_DIGIT;
d &= mask;
} while (d > range);
return (lower + d);
}
DIGIT_T spSimpleRand(DIGIT_T lower, DIGIT_T upper)
{ /* Returns a pseudo-random digit.
Handles own seeding using time.
NOT for cryptographically-secure random numbers.
NOT thread-safe because of static variable.
Changed in Version 2 to use internal funcs.
*/
static unsigned seeded = 0;
if (!seeded)
{
/* seed with system time */
rand_seed((DIGIT_T)time(NULL));
seeded = 1;
}
return rand_between(lower, upper);
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -