AuroraMiddleware/libtommath
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

added libtommath-0.28

Browse Source
This commit is contained in:
Tom St Denis 2003-12-24 18:59:22 +00:00 committed by Steffen Jaeckel
parent c343371bb2
commit 455bb4db20
83 changed files with 40276 additions and 1428 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
bn.ilg Normal file
View File

@ -0,0 +1,6 @@
This is makeindex, version 2.14 [02-Oct-2002] (kpathsea + Thai support).
Scanning input file bn.idx....done (53 entries accepted, 0 rejected).
Sorting entries....done (317 comparisons).
Generating output file bn.ind....done (56 lines written, 0 warnings).
Output written in bn.ind.
Transcript written in bn.ilg.

56
bn.ind Normal file
View File

@ -0,0 +1,56 @@
\begin{theindex}
\item mp\_add, \hyperpage{23}
\item mp\_and, \hyperpage{23}
\item mp\_clear, \hyperpage{7}
\item mp\_clear\_multi, \hyperpage{8}
\item mp\_cmp, \hyperpage{18}
\item mp\_cmp\_d, \hyperpage{20}
\item mp\_cmp\_mag, \hyperpage{17}
\item mp\_div, \hyperpage{29}
\item mp\_div\_2, \hyperpage{21}
\item mp\_div\_2d, \hyperpage{22}
\item MP\_EQ, \hyperpage{17}
\item mp\_error\_to\_string, \hyperpage{6}
\item mp\_expt\_d, \hyperpage{31}
\item mp\_exptmod, \hyperpage{31}
\item mp\_gcd, \hyperpage{39}
\item mp\_grow, \hyperpage{12}
\item MP\_GT, \hyperpage{17}
\item mp\_init, \hyperpage{7}
\item mp\_init\_copy, \hyperpage{9}
\item mp\_init\_multi, \hyperpage{8}
\item mp\_init\_size, \hyperpage{10}
\item mp\_int, \hyperpage{6}
\item mp\_invmod, \hyperpage{40}
\item mp\_jacobi, \hyperpage{39}
\item mp\_lcm, \hyperpage{39}
\item mp\_lshd, \hyperpage{23}
\item MP\_LT, \hyperpage{17}
\item MP\_MEM, \hyperpage{5}
\item mp\_mul, \hyperpage{25}
\item mp\_mul\_2, \hyperpage{21}
\item mp\_mul\_2d, \hyperpage{22}
\item mp\_n\_root, \hyperpage{31}
\item mp\_neg, \hyperpage{24}
\item MP\_NO, \hyperpage{5}
\item MP\_OKAY, \hyperpage{5}
\item mp\_or, \hyperpage{23}
\item mp\_prime\_fermat, \hyperpage{33}
\item mp\_prime\_is\_divisible, \hyperpage{33}
\item mp\_prime\_is\_prime, \hyperpage{34}
\item mp\_prime\_miller\_rabin, \hyperpage{33}
\item mp\_prime\_next\_prime, \hyperpage{34}
\item mp\_prime\_rabin\_miller\_trials, \hyperpage{34}
\item mp\_prime\_random, \hyperpage{35}
\item mp\_rshd, \hyperpage{23}
\item mp\_set, \hyperpage{15}
\item mp\_set\_int, \hyperpage{16}
\item mp\_shrink, \hyperpage{11}
\item mp\_sqr, \hyperpage{25}
\item mp\_sub, \hyperpage{23}
\item MP\_VAL, \hyperpage{5}
\item mp\_xor, \hyperpage{23}
\item MP\_YES, \hyperpage{5}
\end{theindex}

BIN
bn.pdf
View File

Binary file not shown.

2076
bn.tex
View File

File diff suppressed because it is too large Load Diff

7
bn_fast_mp_invmod.c
View File

@ -26,11 +26,8 @@ fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
mp_int x, y, u, v, B, D;
int res, neg;
/* 2. [modified] if a,b are both even then return an error!
*
* That is if gcd(a,b) = 2**k * q then obviously there is no inverse.
*/
if (mp_iseven (a) == 1 && mp_iseven (b) == 1) {
/* 2. [modified] b must be odd */
if (mp_iseven (b) == 1) {
return MP_VAL;
}

2
bn_fast_mp_montgomery_reduce.c
View File

@ -73,7 +73,7 @@ fast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
* that W[ix-1] have the carry cleared (see after the inner loop)
*/
register mp_digit mu;
mu = ((W[ix] & MP_MASK) * rho) & MP_MASK;
mu = (mp_digit) (((W[ix] & MP_MASK) * rho) & MP_MASK);
/* a = a + mu * m * b**i
*

3
bn_mp_add.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* high level addition (handles signs) */
int
mp_add (mp_int * a, mp_int * b, mp_int * c)
int mp_add (mp_int * a, mp_int * b, mp_int * c)
{
int sa, sb, res;

2
bn_mp_clear.c
View File

@ -24,7 +24,7 @@ mp_clear (mp_int * a)
memset (a->dp, 0, sizeof (mp_digit) * a->used);
/* free ram */
free (a->dp);
XFREE(a->dp);
/* reset members to make debugging easier */
a->dp = NULL;

3
bn_mp_cmp_d.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* compare a digit */
int
mp_cmp_d (mp_int * a, mp_digit b)
int mp_cmp_d(mp_int * a, mp_digit b)
{
/* compare based on sign */
if (a->sign == MP_NEG) {

3
bn_mp_cmp_mag.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* compare maginitude of two ints (unsigned) */
int
mp_cmp_mag (mp_int * a, mp_int * b)
int mp_cmp_mag (mp_int * a, mp_int * b)
{
int n;
mp_digit *tmpa, *tmpb;

3
bn_mp_div.c
View File

@ -27,8 +27,7 @@
* The overall algorithm is as described as
* 14.20 from HAC but fixed to treat these cases.
*/
int
mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
int mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
{
mp_int q, x, y, t1, t2;
int res, n, t, i, norm, neg;

3
bn_mp_div_2.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* b = a/2 */
int
mp_div_2 (mp_int * a, mp_int * b)
int mp_div_2(mp_int * a, mp_int * b)
{
int x, res, oldused;

3
bn_mp_div_2d.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* shift right by a certain bit count (store quotient in c, optional remainder in d) */
int
mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d)
int mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d)
{
mp_digit D, r, rr;
int x, res;

2
bn_mp_div_d.c
View File

@ -81,7 +81,7 @@ mp_div_d (mp_int * a, mp_digit b, mp_int * c, mp_digit * d)
if (w >= b) {
t = (mp_digit)(w / b);
w = w % b;
w -= ((mp_word)t) * ((mp_word)b);
} else {
t = 0;
}

3
bn_mp_expt_d.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* calculate c = a**b using a square-multiply algorithm */
int
mp_expt_d (mp_int * a, mp_digit b, mp_int * c)
int mp_expt_d (mp_int * a, mp_digit b, mp_int * c)
{
int res, x;
mp_int g;

3
bn_mp_exptmod.c
View File

@ -20,8 +20,7 @@
* embedded in the normal function but that wasted alot of stack space
* for nothing (since 99% of the time the Montgomery code would be called)
*/
int
mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
{
int dr;

8
bn_mp_fwrite.c
View File

@ -24,24 +24,24 @@ int mp_fwrite(mp_int *a, int radix, FILE *stream)
return MP_VAL;
}
buf = malloc(len);
buf = XMALLOC (len);
if (buf == NULL) {
return MP_MEM;
}
if ((err = mp_toradix(a, buf, radix)) != MP_OKAY) {
free(buf);
XFREE (buf);
return err;
}
for (x = 0; x < len; x++) {
if (fputc(buf[x], stream) == EOF) {
free(buf);
XFREE (buf);
return MP_VAL;
}
}
free(buf);
XFREE (buf);
return MP_OKAY;
}

3
bn_mp_gcd.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* Greatest Common Divisor using the binary method */
int
mp_gcd (mp_int * a, mp_int * b, mp_int * c)
int mp_gcd (mp_int * a, mp_int * b, mp_int * c)
{
mp_int u, v;
int k, u_lsb, v_lsb, res;

6
bn_mp_grow.c
View File

@ -15,13 +15,11 @@
#include <tommath.h>
/* grow as required */
int
mp_grow (mp_int * a, int size)
int mp_grow (mp_int * a, int size)
{
int i;
mp_digit *tmp;
/* if the alloc size is smaller alloc more ram */
if (a->alloc < size) {
/* ensure there are always at least MP_PREC digits extra on top */
@ -33,7 +31,7 @@ mp_grow (mp_int * a, int size)
* in case the operation failed we don't want
* to overwrite the dp member of a.
*/
tmp = OPT_CAST realloc (a->dp, sizeof (mp_digit) * size);
tmp = OPT_CAST XREALLOC (a->dp, sizeof (mp_digit) * size);
if (tmp == NULL) {
/* reallocation failed but "a" is still valid [can be freed] */
return MP_MEM;

5
bn_mp_init.c
View File

@ -15,11 +15,10 @@
#include <tommath.h>
/* init a new bigint */
int
mp_init (mp_int * a)
int mp_init (mp_int * a)
{
/* allocate memory required and clear it */
a->dp = OPT_CAST calloc (sizeof (mp_digit), MP_PREC);
a->dp = OPT_CAST XCALLOC (sizeof (mp_digit), MP_PREC);
if (a->dp == NULL) {
return MP_MEM;
}

3
bn_mp_init_copy.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* creates "a" then copies b into it */
int
mp_init_copy (mp_int * a, mp_int * b)
int mp_init_copy (mp_int * a, mp_int * b)
{
int res;

5
bn_mp_init_size.c
View File

@ -15,14 +15,13 @@
#include <tommath.h>
/* init an mp_init for a given size */
int
mp_init_size (mp_int * a, int size)
int mp_init_size (mp_int * a, int size)
{
/* pad size so there are always extra digits */
size += (MP_PREC * 2) - (size % MP_PREC);
/* alloc mem */
a->dp = OPT_CAST calloc (sizeof (mp_digit), size);
a->dp = OPT_CAST XCALLOC (sizeof (mp_digit), size);
if (a->dp == NULL) {
return MP_MEM;
}

3
bn_mp_invmod.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* hac 14.61, pp608 */
int
mp_invmod (mp_int * a, mp_int * b, mp_int * c)
int mp_invmod (mp_int * a, mp_int * b, mp_int * c)
{
mp_int x, y, u, v, A, B, C, D;
int res;

3
bn_mp_jacobi.c
View File

@ -17,8 +17,7 @@
/* computes the jacobi c = (a | n) (or Legendre if n is prime)
* HAC pp. 73 Algorithm 2.149
*/
int
mp_jacobi (mp_int * a, mp_int * p, int *c)
int mp_jacobi (mp_int * a, mp_int * p, int *c)
{
mp_int a1, p1;
int k, s, r, res;

3
bn_mp_lcm.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* computes least common multiple as |a*b|/(a, b) */
int
mp_lcm (mp_int * a, mp_int * b, mp_int * c)
int mp_lcm (mp_int * a, mp_int * b, mp_int * c)
{
int res;
mp_int t1, t2;

3
bn_mp_lshd.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* shift left a certain amount of digits */
int
mp_lshd (mp_int * a, int b)
int mp_lshd (mp_int * a, int b)
{
int x, res;

2
bn_mp_montgomery_reduce.c
View File

@ -51,7 +51,7 @@ mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
* following inner loop to reduce the
* input one digit at a time
*/
mu = ((mp_word)x->dp[ix]) * ((mp_word)rho) & MP_MASK;
mu = (mp_digit) (((mp_word)x->dp[ix]) * ((mp_word)rho) & MP_MASK);
/* a = a + mu * m * b**i */
{

3
bn_mp_mul.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* high level multiplication (handles sign) */
int
mp_mul (mp_int * a, mp_int * b, mp_int * c)
int mp_mul (mp_int * a, mp_int * b, mp_int * c)
{
int res, neg;
neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;

3
bn_mp_mul_2.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* b = a*2 */
int
mp_mul_2 (mp_int * a, mp_int * b)
int mp_mul_2(mp_int * a, mp_int * b)
{
int x, res, oldused;

3
bn_mp_mul_2d.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* shift left by a certain bit count */
int
mp_mul_2d (mp_int * a, int b, mp_int * c)
int mp_mul_2d (mp_int * a, int b, mp_int * c)
{
mp_digit d;
int res;

1
bn_mp_mulmod.c
View File

@ -21,7 +21,6 @@ mp_mulmod (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
int res;
mp_int t;
if ((res = mp_init (&t)) != MP_OKAY) {
return res;
}

3
bn_mp_n_root.c
View File

@ -24,8 +24,7 @@
* each step involves a fair bit. This is not meant to
* find huge roots [square and cube, etc].
*/
int
mp_n_root (mp_int * a, mp_digit b, mp_int * c)
int mp_n_root (mp_int * a, mp_digit b, mp_int * c)
{
mp_int t1, t2, t3;
int res, neg;

3
bn_mp_neg.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* b = -a */
int
mp_neg (mp_int * a, mp_int * b)
int mp_neg (mp_int * a, mp_int * b)
{
int res;
if ((res = mp_copy (a, b)) != MP_OKAY) {

3
bn_mp_or.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* OR two ints together */
int
mp_or (mp_int * a, mp_int * b, mp_int * c)
int mp_or (mp_int * a, mp_int * b, mp_int * c)
{
int res, ix, px;
mp_int t, *x;

7
bn_mp_prime_fermat.c
View File

@ -22,14 +22,13 @@
*
* Sets result to 1 if the congruence holds, or zero otherwise.
*/
int
mp_prime_fermat (mp_int * a, mp_int * b, int *result)
int mp_prime_fermat (mp_int * a, mp_int * b, int *result)
{
mp_int t;
int err;
/* default to composite */
*result = 0;
*result = MP_NO;
/* ensure b > 1 */
if (mp_cmp_d(b, 1) != MP_GT) {
@ -48,7 +47,7 @@ mp_prime_fermat (mp_int * a, mp_int * b, int *result)
/* is it equal to b? */
if (mp_cmp (&t, b) == MP_EQ) {
*result = 1;
*result = MP_YES;
}
err = MP_OKAY;

7
bn_mp_prime_is_divisible.c
View File

@ -19,14 +19,13 @@
*
* sets result to 0 if not, 1 if yes
*/
int
mp_prime_is_divisible (mp_int * a, int *result)
int mp_prime_is_divisible (mp_int * a, int *result)
{
int err, ix;
mp_digit res;
/* default to not */
*result = 0;
*result = MP_NO;
for (ix = 0; ix < PRIME_SIZE; ix++) {
/* what is a mod __prime_tab[ix] */
@ -36,7 +35,7 @@ mp_prime_is_divisible (mp_int * a, int *result)
/* is the residue zero? */
if (res == 0) {
*result = 1;
*result = MP_YES;
return MP_OKAY;
}
}

11
bn_mp_prime_is_prime.c
View File

@ -21,14 +21,13 @@
*
* Sets result to 1 if probably prime, 0 otherwise
*/
int
mp_prime_is_prime (mp_int * a, int t, int *result)
int mp_prime_is_prime (mp_int * a, int t, int *result)
{
mp_int b;
int ix, err, res;
/* default to no */
*result = 0;
*result = MP_NO;
/* valid value of t? */
if (t <= 0 || t > PRIME_SIZE) {
@ -49,7 +48,7 @@ mp_prime_is_prime (mp_int * a, int t, int *result)
}
/* return if it was trivially divisible */
if (res == 1) {
if (res == MP_YES) {
return MP_OKAY;
}
@ -66,13 +65,13 @@ mp_prime_is_prime (mp_int * a, int t, int *result)
goto __B;
}
if (res == 0) {
if (res == MP_NO) {
goto __B;
}
}
/* passed the test */
*result = 1;
*result = MP_YES;
__B:mp_clear (&b);
return err;
}

7
bn_mp_prime_miller_rabin.c
View File

@ -21,14 +21,13 @@
* Randomly the chance of error is no more than 1/4 and often
* very much lower.
*/
int
mp_prime_miller_rabin (mp_int * a, mp_int * b, int *result)
int mp_prime_miller_rabin (mp_int * a, mp_int * b, int *result)
{
mp_int n1, y, r;
int s, j, err;
/* default */
*result = 0;
*result = MP_NO;
/* ensure b > 1 */
if (mp_cmp_d(b, 1) != MP_GT) {
@ -90,7 +89,7 @@ mp_prime_miller_rabin (mp_int * a, mp_int * b, int *result)
}
/* probably prime now */
*result = 1;
*result = MP_YES;
__Y:mp_clear (&y);
__R:mp_clear (&r);
__N1:mp_clear (&n1);

4
bn_mp_prime_next_prime.c
View File

@ -146,12 +146,12 @@ int mp_prime_next_prime(mp_int *a, int t, int bbs_style)
if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
goto __ERR;
}
if (res == 0) {
if (res == MP_NO) {
break;
}
}
if (res == 1) {
if (res == MP_YES) {
break;
}
}

74
bn_mp_prime_random.c Normal file
View File

@ -0,0 +1,74 @@
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
* LibTomMath is a library that provides multiple-precision
* integer arithmetic as well as number theoretic functionality.
*
* The library was designed directly after the MPI library by
* Michael Fromberger but has been written from scratch with
* additional optimizations in place.
*
* The library is free for all purposes without any express
* guarantee it works.
*
* Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
*/
#include <tommath.h>
/* makes a truly random prime of a given size (bytes),
* call with bbs = 1 if you want it to be congruent to 3 mod 4
*
* You have to supply a callback which fills in a buffer with random bytes. "dat" is a parameter you can
* have passed to the callback (e.g. a state or something). This function doesn't use "dat" itself
* so it can be NULL
*
* The prime generated will be larger than 2^(8*size).
*/
/* this sole function may hold the key to enslaving all mankind! */
int mp_prime_random(mp_int *a, int t, int size, int bbs, ltm_prime_callback cb, void *dat)
{
unsigned char *tmp;
int res, err;
/* sanity check the input */
if (size <= 0) {
return MP_VAL;
}
/* we need a buffer of size+1 bytes */
tmp = XMALLOC(size+1);
if (tmp == NULL) {
return MP_MEM;
}
/* fix MSB */
tmp[0] = 1;
do {
/* read the bytes */
if (cb(tmp+1, size, dat) != size) {
err = MP_VAL;
goto error;
}
/* fix the LSB */
tmp[size] |= (bbs ? 3 : 1);
/* read it in */
if ((err = mp_read_unsigned_bin(a, tmp, size+1)) != MP_OKAY) {
goto error;
}
/* is it prime? */
if ((err = mp_prime_is_prime(a, t, &res)) != MP_OKAY) {
goto error;
}
} while (res == MP_NO);
err = MP_OKAY;
error:
XFREE(tmp);
return err;
}

10
bn_mp_radix_size.c
View File

@ -27,28 +27,36 @@ mp_radix_size (mp_int * a, int radix)
return mp_count_bits (a) + (a->sign == MP_NEG ? 1 : 0) + 1;
}
/* make sure the radix is in range */
if (radix < 2 || radix > 64) {
return 0;
}
/* init a copy of the input */
if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
return 0;
}
/* digs is the digit count */
digs = 0;
/* if it's negative add one for the sign */
if (t.sign == MP_NEG) {
++digs;
t.sign = MP_ZPOS;
}
/* fetch out all of the digits */
while (mp_iszero (&t) == 0) {
if ((res = mp_div_d (&t, (mp_digit) radix, &t, &d)) != MP_OKAY) {
mp_clear (&t);
return 0;
return res;
}
++digs;
}
mp_clear (&t);
/* return digs + 1, the 1 is for the NULL byte that would be required. */
return digs + 1;
}

17
bn_mp_reduce_is_2k.c
View File

@ -15,10 +15,9 @@
#include <tommath.h>
/* determines if mp_reduce_2k can be used */
int
mp_reduce_is_2k(mp_int *a)
int mp_reduce_is_2k(mp_int *a)
{
int ix, iy;
int ix, iy, iz, iw;
if (a->used == 0) {
return 0;
@ -26,11 +25,19 @@ mp_reduce_is_2k(mp_int *a)
return 1;
} else if (a->used > 1) {
iy = mp_count_bits(a);
iz = 1;
iw = 1;
/* Test every bit from the second digit up, must be 1 */
for (ix = DIGIT_BIT; ix < iy; ix++) {
if ((a->dp[ix/DIGIT_BIT] &
((mp_digit)1 << (mp_digit)(ix % DIGIT_BIT))) == 0) {
if ((a->dp[iw] & iz) == 0) {
return 0;
}
iz <<= 1;
if (iz > (int)MP_MASK) {
++iw;
iz = 1;
}
}
}
return 1;

3
bn_mp_rshd.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* shift right a certain amount of digits */
void
mp_rshd (mp_int * a, int b)
void mp_rshd (mp_int * a, int b)
{
int x;

3
bn_mp_set.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* set to a digit */
void
mp_set (mp_int * a, mp_digit b)
void mp_set (mp_int * a, mp_digit b)
{
mp_zero (a);
a->dp[0] = b & MP_MASK;

3
bn_mp_set_int.c
View File

@ -15,8 +15,7 @@
#include <tommath.h>
/* set a 32-bit const */
int
mp_set_int (mp_int * a, unsigned int b)
int mp_set_int (mp_int * a, unsigned long b)
{
int x, res;

7
bn_mp_shrink.c
View File

@ -15,13 +15,14 @@
#include <tommath.h>
/* shrink a bignum */
int
mp_shrink (mp_int * a)
int mp_shrink (mp_int * a)
{
mp_digit *tmp;
if (a->alloc != a->used) {
if ((a->dp = OPT_CAST realloc (a->dp, sizeof (mp_digit) * a->used)) == NULL) {
if ((tmp = OPT_CAST XREALLOC (a->dp, sizeof (mp_digit) * a->used)) == NULL) {
return MP_MEM;
}
a->dp = tmp;
a->alloc = a->used;
}
return MP_OKAY;

4
bn_mp_toradix.c
View File

@ -23,6 +23,7 @@ mp_toradix (mp_int * a, char *str, int radix)
mp_digit d;
char *_s = str;
/* check range of the radix */
if (radix < 2 || radix > 64) {
return MP_VAL;
}
@ -61,8 +62,7 @@ mp_toradix (mp_int * a, char *str, int radix)
bn_reverse ((unsigned char *)_s, digs);
/* append a NULL so the string is properly terminated */
*str++ = '\0';
*str = '\0';
mp_clear (&t);
return MP_OKAY;

69
bn_prime_sizes_tab.c Normal file
View File

@ -0,0 +1,69 @@
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
* LibTomMath is a library that provides multiple-precision
* integer arithmetic as well as number theoretic functionality.
*
* The library was designed directly after the MPI library by
* Michael Fromberger but has been written from scratch with
* additional optimizations in place.
*
* The library is free for all purposes without any express
* guarantee it works.
*
* Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
*/
#include <tommath.h>
/* this table gives the # of rabin miller trials for a prob of failure lower than 2^-96 */
static const struct {
int k, t;
} sizes[] = {
{ 128, 28 },
{ 256, 16 },
{ 384, 10 },
{ 512, 7 },
{ 640, 6 },
{ 768, 5 },
{ 896, 4 },
{ 1024, 4 },
{ 1152, 3 },
{ 1280, 3 },
{ 1408, 3 },
{ 1536, 3 },
{ 1664, 3 },
{ 1792, 2 },
{ 1920, 2 },
{ 2048, 2 },
{ 2176, 2 },
{ 2304, 2 },
{ 2432, 2 },
{ 2560, 2 },
{ 2688, 2 },
{ 2816, 2 },
{ 2944, 2 },
{ 3072, 2 },
{ 3200, 2 },
{ 3328, 2 },
{ 3456, 2 },
{ 3584, 2 },
{ 3712, 2 },
{ 3840, 1 },
{ 3968, 1 },
{ 4096, 1 } };
/* returns # of RM trials required for a given bit size */
int mp_prime_rabin_miller_trials(int size)
{
int x;
for (x = 0; x < (int)(sizeof(sizes)/(sizeof(sizes[0]))); x++) {
if (sizes[x].k == size) {
return sizes[x].t;
} else if (sizes[x].k > size) {
return (x == 0) ? sizes[0].t : sizes[x - 1].t;
}
}
return 1;
}

12
bn_s_mp_sqr.c
View File

@ -27,6 +27,8 @@ s_mp_sqr (mp_int * a, mp_int * b)
if ((res = mp_init_size (&t, 2*pa + 1)) != MP_OKAY) {
return res;
}
/* default used is maximum possible size */
t.used = 2*pa + 1;
for (ix = 0; ix < pa; ix++) {
@ -36,20 +38,20 @@ s_mp_sqr (mp_int * a, mp_int * b)
((mp_word)a->dp[ix])*((mp_word)a->dp[ix]);
/* store lower part in result */
t.dp[2*ix] = (mp_digit) (r & ((mp_word) MP_MASK));
t.dp[ix+ix] = (mp_digit) (r & ((mp_word) MP_MASK));
/* get the carry */
u = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
u = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
/* left hand side of A[ix] * A[iy] */
tmpx = a->dp[ix];
tmpx = a->dp[ix];
/* alias for where to store the results */
tmpt = t.dp + (2*ix + 1);
tmpt = t.dp + (2*ix + 1);
for (iy = ix + 1; iy < pa; iy++) {
/* first calculate the product */
r = ((mp_word)tmpx) * ((mp_word)a->dp[iy]);
r = ((mp_word)tmpx) * ((mp_word)a->dp[iy]);
/* now calculate the double precision result, note we use
* addition instead of *2 since it's easier to optimize

21
changes.txt
View File

@ -1,3 +1,24 @@
Dec 24th, 2003
v0.28 -- Henrik suggested I add casts to the montomgery code [stores into mu...] so compilers wouldn't
spew [erroneous] diagnostics... fixed.
-- Henrik also spotted two typos. One in mp_radix_size() and another in mp_toradix().
-- Added fix to mp_shrink() to avoid a memory leak.
-- Added mp_prime_random() which requires a callback to make truly random primes of a given nature
(idea from chat with Niels Ferguson at Crypto'03)
-- Picked up a second wind. I'm filled with Gooo. Mission Gooo!
-- Removed divisions from mp_reduce_is_2k()
-- Sped up mp_div_d() [general case] to use only one division per digit instead of two.
-- Added the heap macros from LTC to LTM. Now you can easily [by editing four lines of tommath.h]
change the name of the heap functions used in LTM [also compatible with LTC via MPI mode]
-- Added bn_prime_rabin_miller_trials() which gives the number of Rabin-Miller trials to achieve
a failure rate of less than 2^-96
-- fixed bug in fast_mp_invmod(). The initial testing logic was wrong. An invalid input is not when
"a" and "b" are even it's when "b" is even [the algo is for odd moduli only].
-- Started a new manual [finally]. It is incomplete and will be finished as time goes on. I had to stop
adding full demos around half way in chapter three so I could at least get a good portion of the
manual done. If you really need help using the library you can always email me!
-- My Textbook is now included as part of the package [all Public Domain]
Sept 19th, 2003
v0.27 -- Removed changes.txt~ which was made by accident since "kate" decided it was
a good time to re-enable backups... [kde is fun!]

39
demo/demo.c
View File

@ -111,10 +111,10 @@ int main(void)
/* test mp_reduce_2k */
#if 0
for (cnt = 3; cnt <= 256; ++cnt) {
for (cnt = 3; cnt <= 384; ++cnt) {
mp_digit tmp;
mp_2expt(&a, cnt);
mp_sub_d(&a, 1, &a); /* a = 2**cnt - 1 */
mp_sub_d(&a, 2, &a); /* a = 2**cnt - 2 */
printf("\nTesting %4d bits", cnt);
@ -138,11 +138,11 @@ int main(void)
/* test mp_div_3 */
#if 0
for (cnt = 0; cnt < 10000; ) {
for (cnt = 0; cnt < 1000000; ) {
mp_digit r1, r2;
if (!(++cnt & 127)) printf("%9d\r", cnt);
mp_rand(&a, abs(rand()) % 32 + 1);
mp_rand(&a, abs(rand()) % 128 + 1);
mp_div_d(&a, 3, &b, &r1);
mp_div_3(&a, &c, &r2);
@ -155,7 +155,7 @@ int main(void)
/* test the DR reduction */
#if 0
for (cnt = 2; cnt < 32; cnt++) {
for (cnt = 2; cnt < 128; cnt++) {
printf("%d digit modulus\n", cnt);
mp_grow(&a, cnt);
mp_zero(&a);
@ -181,7 +181,7 @@ int main(void)
printf("Failed on trial %lu\n", rr); exit(-1);
}
} while (++rr < 10000);
} while (++rr < 100000);
printf("Passed DR test for %d digits\n", cnt);
}
#endif
@ -203,8 +203,8 @@ int main(void)
rr += 16;
} while (rdtsc() < (CLOCKS_PER_SEC * 2));
tt = rdtsc();
printf("Adding\t\t%4d-bit => %9llu/sec, %9llu ticks\n", mp_count_bits(&a), (((unsigned long long)rr)*CLOCKS_PER_SEC)/tt, tt);
fprintf(log, "%d %9llu\n", cnt*DIGIT_BIT, (((unsigned long long)rr)*CLOCKS_PER_SEC)/tt);
printf("Adding\t\t%4d-bit => %9llu/sec, %9llu ticks\n", mp_count_bits(&a), (((ulong64)rr)*CLOCKS_PER_SEC)/tt, tt);
fprintf(log, "%d %9llu\n", cnt*DIGIT_BIT, (((ulong64)rr)*CLOCKS_PER_SEC)/tt);
}
fclose(log);
@ -220,12 +220,13 @@ int main(void)
rr += 16;
} while (rdtsc() < (CLOCKS_PER_SEC * 2));
tt = rdtsc();
printf("Subtracting\t\t%4d-bit => %9llu/sec, %9llu ticks\n", mp_count_bits(&a), (((unsigned long long)rr)*CLOCKS_PER_SEC)/tt, tt);
fprintf(log, "%d %9llu\n", cnt*DIGIT_BIT, (((unsigned long long)rr)*CLOCKS_PER_SEC)/tt);
printf("Subtracting\t\t%4d-bit => %9llu/sec, %9llu ticks\n", mp_count_bits(&a), (((ulong64)rr)*CLOCKS_PER_SEC)/tt, tt);
fprintf(log, "%d %9llu\n", cnt*DIGIT_BIT, (((ulong64)rr)*CLOCKS_PER_SEC)/tt);
}
fclose(log);
/* do mult/square twice, first without karatsuba and second with */
mult_test:
old_kara_m = KARATSUBA_MUL_CUTOFF;
old_kara_s = KARATSUBA_SQR_CUTOFF;
for (ix = 0; ix < 2; ix++) {
@ -246,8 +247,8 @@ int main(void)
rr += 16;
} while (rdtsc() < (CLOCKS_PER_SEC * 2));
tt = rdtsc();
printf("Multiplying\t%4d-bit => %9llu/sec, %9llu ticks\n", mp_count_bits(&a), (((unsigned long long)rr)*CLOCKS_PER_SEC)/tt, tt);
fprintf(log, "%d %9llu\n", cnt*DIGIT_BIT, (((unsigned long long)rr)*CLOCKS_PER_SEC)/tt);
printf("Multiplying\t%4d-bit => %9llu/sec, %9llu ticks\n", mp_count_bits(&a), (((ulong64)rr)*CLOCKS_PER_SEC)/tt, tt);
fprintf(log, "%d %9llu\n", mp_count_bits(&a), (((ulong64)rr)*CLOCKS_PER_SEC)/tt);
}
fclose(log);
@ -262,8 +263,8 @@ int main(void)
rr += 16;
} while (rdtsc() < (CLOCKS_PER_SEC * 2));
tt = rdtsc();
printf("Squaring\t%4d-bit => %9llu/sec, %9llu ticks\n", mp_count_bits(&a), (((unsigned long long)rr)*CLOCKS_PER_SEC)/tt, tt);
fprintf(log, "%d %9llu\n", cnt*DIGIT_BIT, (((unsigned long long)rr)*CLOCKS_PER_SEC)/tt);
printf("Squaring\t%4d-bit => %9llu/sec, %9llu ticks\n", mp_count_bits(&a), (((ulong64)rr)*CLOCKS_PER_SEC)/tt, tt);
fprintf(log, "%d %9llu\n", mp_count_bits(&a), (((ulong64)rr)*CLOCKS_PER_SEC)/tt);
}
fclose(log);
@ -297,12 +298,14 @@ int main(void)
"1214855636816562637502584060163403830270705000634713483015101384881871978446801224798536155406895823305035467591632531067547890948695117172076954220727075688048751022421198712032848890056357845974246560748347918630050853933697792254955890439720297560693579400297062396904306270145886830719309296352765295712183040773146419022875165382778007040109957609739589875590885701126197906063620133954893216612678838507540777138437797705602453719559017633986486649523611975865005712371194067612263330335590526176087004421363598470302731349138773205901447704682181517904064735636518462452242791676541725292378925568296858010151852326316777511935037531017413910506921922450666933202278489024521263798482237150056835746454842662048692127173834433089016107854491097456725016327709663199738238442164843147132789153725513257167915555162094970853584447993125488607696008169807374736711297007473812256272245489405898470297178738029484459690836250560495461579533254473316340608217876781986188705928270735695752830825527963838355419762516246028680280988020401914551825487349990306976304093109384451438813251211051597392127491464898797406789175453067960072008590614886532333015881171367104445044718144312416815712216611576221546455968770801413440778423979",
NULL
};
expt_test:
log = fopen("logs/expt.log", "w");
logb = fopen("logs/expt_dr.log", "w");
logc = fopen("logs/expt_2k.log", "w");
for (n = 0; primes[n]; n++) {
SLEEP;
mp_read_radix(&a, primes[n], 10);
printf("Different (%d)!!!\n", mp_count_bits(&a));
mp_zero(&b);
for (rr = 0; rr < mp_count_bits(&a); rr++) {
mp_mul_2(&b, &b);
@ -328,8 +331,8 @@ int main(void)
draw(&d);
exit(0);
}
printf("Exponentiating\t%4d-bit => %9llu/sec, %9llu ticks\n", mp_count_bits(&a), (((unsigned long long)rr)*CLOCKS_PER_SEC)/tt, tt);
fprintf((n < 6) ? logc : (n < 13) ? logb : log, "%d %9llu\n", mp_count_bits(&a), (((unsigned long long)rr)*CLOCKS_PER_SEC)/tt);
printf("Exponentiating\t%4d-bit => %9llu/sec, %9llu ticks\n", mp_count_bits(&a), (((ulong64)rr)*CLOCKS_PER_SEC)/tt, tt);
fprintf((n < 6) ? logc : (n < 13) ? logb : log, "%d %9llu\n", mp_count_bits(&a), (((ulong64)rr)*CLOCKS_PER_SEC)/tt);
}
}
fclose(log);
@ -359,8 +362,8 @@ int main(void)
printf("Failed to invert\n");
return 0;
}
printf("Inverting mod\t%4d-bit => %9llu/sec, %9llu ticks\n", mp_count_bits(&a), (((unsigned long long)rr)*CLOCKS_PER_SEC)/tt, tt);
fprintf(log, "%d %9llu\n", cnt*DIGIT_BIT, (((unsigned long long)rr)*CLOCKS_PER_SEC)/tt);
printf("Inverting mod\t%4d-bit => %9llu/sec, %9llu ticks\n", mp_count_bits(&a), (((ulong64)rr)*CLOCKS_PER_SEC)/tt, tt);
fprintf(log, "%d %9llu\n", cnt*DIGIT_BIT, (((ulong64)rr)*CLOCKS_PER_SEC)/tt);
}
fclose(log);

2
etc/2kprime.1
View File

@ -0,0 +1,2 @@
256-bits (k = 36113) = 115792089237316195423570985008687907853269984665640564039457584007913129603823
512-bits (k = 38117) = 13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006045979

3
etc/drprimes.txt
View File

@ -1,3 +1,6 @@
280-bit prime:
p == 1942668892225729070919461906823518906642406839052139521251812409738904285204940164839
532-bit prime:
p == 14059105607947488696282932836518693308967803494693489478439861164411992439598399594747002144074658928593502845729752797260025831423419686528151609940203368691747

5
etc/makefile
View File

@ -19,6 +19,11 @@ tune86: tune.c
nasm -f coff timer.asm
$(CC) -DX86_TIMER $(CFLAGS) tune.c timer.o $(LIBNAME) -o tune86
# for cygwin
tune86c: tune.c
nasm -f gnuwin32 timer.asm
$(CC) -DX86_TIMER $(CFLAGS) tune.c timer.o $(LIBNAME) -o tune86
#make tune86 for linux or any ELF format
tune86l: tune.c
nasm -f elf -DUSE_ELF timer.asm

32
logs/add.log
View File

@ -1,16 +1,16 @@
224 12849536
448 9956080
672 8372000
896 7065464
1120 5824864
1344 5141728
1568 4511808
1792 4170480
2016 3802536
2240 3500936
2464 3193352
2688 2991976
2912 2818672
3136 2661448
3360 2506560
3584 2343304
224 8622881
448 7241417
672 5844712
896 4938016
1120 4256543
1344 3728000
1568 3328263
1792 3012161
2016 2743688
2240 2512095
2464 2234464
2688 1960139
2912 2013395
3136 1879636
3360 1756301
3584 1680982

7
logs/expt.log
View File

@ -1,7 +0,0 @@
513 600
769 221
1025 103
2049 15
2561 8
3073 4
4097 2

6
logs/expt_2k.log
View File

@ -1,6 +0,0 @@
521 728
607 549
1279 100
2203 29
3217 11
4253 5

7
logs/expt_dr.log
View File

@ -1,7 +0,0 @@
532 1032
784 424
1036 214
1540 81
2072 38
3080 13
4116 5

34
logs/mult.log
View File

@ -1,17 +1,17 @@
896 301008
1344 141872
1792 84424
2240 55864
2688 39784
3136 29624
3584 22952
4032 18304
4480 14944
4928 12432
5376 10496
5824 8976
6272 7776
6720 6792
7168 1656
7616 1472
8064 1312
896 348504
1344 165040
1792 98696
2240 65400
2688 46672
3136 34968
3584 27144
4032 21648
4480 17672
4928 14768
5376 12416
5824 10696
6272 9184
6720 8064
7168 1896
7616 1680
8064 1504

34
logs/sqr.log
View File

@ -1,17 +1,17 @@
896 371856
1344 196352
1792 122312
2240 83144
2688 60304
3136 45832
3584 12760
4032 10160
4480 8352
4928 6944
5376 5824
5824 5008
6272 4336
6720 3768
7168 3280
7616 2952
8064 2640
911 167013
1359 83796
1807 50308
2254 33637
2703 24067
3151 17997
3599 5751
4047 4561
4490 3714
4943 3067
5391 2597
5839 2204
6286 1909
6735 1637
7183 1461
7631 1302
8078 1158

17
logs/sqr.old Normal file
View File

@ -0,0 +1,17 @@
896 382617
1344 207161
1792 131522
2240 90775
2688 66652
3136 50955
3584 11678
4032 9342
4480 7684
4928 6382
5376 5399
5824 4545
6272 3994
6720 3490
7168 3075
7616 2733
8064 2428

34
logs/sqr_kara.log
View File

@ -1,17 +1,17 @@
896 372256
1344 196368
1792 122272
2240 82976
2688 60480
3136 45808
3584 33296
4032 27888
4480 23608
4928 20296
5376 17576
5824 15416
6272 13600
6720 12104
7168 10080
7616 9232
8064 8008
910 165312
1358 84355
1806 50316
2255 33661
2702 24027
3151 18068
3599 14721
4046 12101
4493 10112
4942 8591
5390 7364
5839 6398
6285 5607
6735 4952
7182 4625
7631 4193
8079 3810

32
logs/sub.log
View File

@ -1,16 +1,16 @@
224 9325944
448 8075808
672 7054912
896 5757992
1120 5081768
1344 4669384
1568 4422384
1792 3900416
2016 3548872
2240 3428912
2464 3216968
2688 2905280
2912 2782664
3136 2591440
3360 2475728
3584 2282216
224 10295756
448 7577910
672 6279588
896 5345182
1120 4646989
1344 4101759
1568 3685447
1792 3337497
2016 3051095
2240 2811900
2464 2605371
2688 2420561
2912 2273174
3136 2134662
3360 2014354
3584 1901723

34
makefile
View File

@ -6,7 +6,7 @@ CFLAGS += -I./ -Wall -W -Wshadow -O3 -funroll-loops
#x86 optimizations [should be valid for any GCC install though]
CFLAGS += -fomit-frame-pointer
VERSION=0.27
VERSION=0.28
default: libtommath.a
@ -44,7 +44,7 @@ bn_mp_toom_mul.o bn_mp_toom_sqr.o bn_mp_div_3.o bn_s_mp_exptmod.o \
bn_mp_reduce_2k.o bn_mp_reduce_is_2k.o bn_mp_reduce_2k_setup.o \
bn_mp_radix_smap.o bn_mp_read_radix.o bn_mp_toradix.o bn_mp_radix_size.o \
bn_mp_fread.o bn_mp_fwrite.o bn_mp_cnt_lsb.o bn_error.o \
bn_mp_init_multi.o bn_mp_clear_multi.o
bn_mp_init_multi.o bn_mp_clear_multi.o bn_mp_prime_random.o bn_prime_sizes_tab.o
libtommath.a: $(OBJECTS)
$(AR) $(ARFLAGS) libtommath.a $(OBJECTS)
@ -87,23 +87,29 @@ docs:
latex tommath > /dev/null
pdflatex tommath
rm -f tommath.log tommath.aux tommath.dvi tommath.idx tommath.toc tommath.lof tommath.ind tommath.ilg
cd pics ; make clean
#the old manual being phased out
manual:
latex bn
pdflatex bn
rm -f bn.aux bn.dvi bn.log
#LTM user manual
mandvi: bn.tex
echo "hello" > bn.ind
latex bn > /dev/null
makeindex bn
latex bn > /dev/null
#LTM user manual [pdf]
manual: mandvi
pdflatex bn >/dev/null
rm -f bn.aux bn.dvi bn.log bn.idx bn.lof bn.out bn.toc
clean:
rm -f *.bat *.pdf *.o *.a *.obj *.lib *.exe etclib/*.o demo/demo.o test ltmtest mpitest mtest/mtest mtest/mtest.exe \
tommath.idx tommath.toc tommath.log tommath.aux tommath.dvi tommath.lof tommath.ind tommath.ilg *.ps *.pdf *.log *.s mpi.c \
poster.aux poster.dvi poster.log
rm -f *.bat *.pdf *.o *.a *.obj *.lib *.exe *.dll etclib/*.o demo/demo.o test ltmtest mpitest mtest/mtest mtest/mtest.exe \
*.idx *.toc *.log *.aux *.dvi *.lof *.ind *.ilg *.ps *.log *.s mpi.c
cd etc ; make clean
cd pics ; make clean
zipup: clean manual poster
zipup: clean manual poster docs
perl gen.pl ; mv mpi.c pre_gen/ ; \
cd .. ; rm -rf ltm* libtommath-$(VERSION) ; mkdir libtommath-$(VERSION) ; \
cp -R ./libtommath/* ./libtommath-$(VERSION)/ ; cd ./libtommath-$(VERSION) ; rm -f tommath.src tommath.tex tommath.out ; cd pics ; rm -f *.tif *.ps *.pdf ; cd .. ; cd .. ; ls ; \
tar -c libtommath-$(VERSION)/* > ltm-$(VERSION).tar ; \
bzip2 -9vv ltm-$(VERSION).tar ; zip -9 -r ltm-$(VERSION).zip libtommath-$(VERSION)/*
cp -R ./libtommath/* ./libtommath-$(VERSION)/ ; \
tar -c libtommath-$(VERSION)/* | bzip2 -9vvc > ltm-$(VERSION).tar.bz2 ; \
zip -9 -r ltm-$(VERSION).zip libtommath-$(VERSION)/*

2
makefile.bcc
View File

@ -29,7 +29,7 @@ bn_mp_toom_mul.obj bn_mp_toom_sqr.obj bn_mp_div_3.obj bn_s_mp_exptmod.obj \
bn_mp_reduce_2k.obj bn_mp_reduce_is_2k.obj bn_mp_reduce_2k_setup.obj \
bn_mp_radix_smap.obj bn_mp_read_radix.obj bn_mp_toradix.obj bn_mp_radix_size.obj \
bn_mp_fread.obj bn_mp_fwrite.obj bn_mp_cnt_lsb.obj bn_error.obj \
bn_mp_init_multi.obj bn_mp_clear_multi.obj
bn_mp_init_multi.obj bn_mp_clear_multi.obj bn_mp_prime_random.obj bn_prime_sizes_tab.obj
TARGET = libtommath.lib

47
makefile.cygwin_dll Normal file
View File

@ -0,0 +1,47 @@
#Makefile for Cygwin-GCC
#
#This makefile will build a Windows DLL [doesn't require cygwin to run] in the file
#libtommath.dll. The import library is in libtommath.dll.a. Remember to add
#"-Wl,--enable-auto-import" to your client build to avoid the auto-import warnings
#
#Tom St Denis
CFLAGS += -I./ -Wall -W -Wshadow -O3 -funroll-loops -mno-cygwin
#x86 optimizations [should be valid for any GCC install though]
CFLAGS += -fomit-frame-pointer
default: windll
OBJECTS=bncore.o bn_mp_init.o bn_mp_clear.o bn_mp_exch.o bn_mp_grow.o bn_mp_shrink.o \
bn_mp_clamp.o bn_mp_zero.o bn_mp_set.o bn_mp_set_int.o bn_mp_init_size.o bn_mp_copy.o \
bn_mp_init_copy.o bn_mp_abs.o bn_mp_neg.o bn_mp_cmp_mag.o bn_mp_cmp.o bn_mp_cmp_d.o \
bn_mp_rshd.o bn_mp_lshd.o bn_mp_mod_2d.o bn_mp_div_2d.o bn_mp_mul_2d.o bn_mp_div_2.o \
bn_mp_mul_2.o bn_s_mp_add.o bn_s_mp_sub.o bn_fast_s_mp_mul_digs.o bn_s_mp_mul_digs.o \
bn_fast_s_mp_mul_high_digs.o bn_s_mp_mul_high_digs.o bn_fast_s_mp_sqr.o bn_s_mp_sqr.o \
bn_mp_add.o bn_mp_sub.o bn_mp_karatsuba_mul.o bn_mp_mul.o bn_mp_karatsuba_sqr.o \
bn_mp_sqr.o bn_mp_div.o bn_mp_mod.o bn_mp_add_d.o bn_mp_sub_d.o bn_mp_mul_d.o \
bn_mp_div_d.o bn_mp_mod_d.o bn_mp_expt_d.o bn_mp_addmod.o bn_mp_submod.o \
bn_mp_mulmod.o bn_mp_sqrmod.o bn_mp_gcd.o bn_mp_lcm.o bn_fast_mp_invmod.o bn_mp_invmod.o \
bn_mp_reduce.o bn_mp_montgomery_setup.o bn_fast_mp_montgomery_reduce.o bn_mp_montgomery_reduce.o \
bn_mp_exptmod_fast.o bn_mp_exptmod.o bn_mp_2expt.o bn_mp_n_root.o bn_mp_jacobi.o bn_reverse.o \
bn_mp_count_bits.o bn_mp_read_unsigned_bin.o bn_mp_read_signed_bin.o bn_mp_to_unsigned_bin.o \
bn_mp_to_signed_bin.o bn_mp_unsigned_bin_size.o bn_mp_signed_bin_size.o \
bn_mp_xor.o bn_mp_and.o bn_mp_or.o bn_mp_rand.o bn_mp_montgomery_calc_normalization.o \
bn_mp_prime_is_divisible.o bn_prime_tab.o bn_mp_prime_fermat.o bn_mp_prime_miller_rabin.o \
bn_mp_prime_is_prime.o bn_mp_prime_next_prime.o bn_mp_dr_reduce.o \
bn_mp_dr_is_modulus.o bn_mp_dr_setup.o bn_mp_reduce_setup.o \
bn_mp_toom_mul.o bn_mp_toom_sqr.o bn_mp_div_3.o bn_s_mp_exptmod.o \
bn_mp_reduce_2k.o bn_mp_reduce_is_2k.o bn_mp_reduce_2k_setup.o \
bn_mp_radix_smap.o bn_mp_read_radix.o bn_mp_toradix.o bn_mp_radix_size.o \
bn_mp_fread.o bn_mp_fwrite.o bn_mp_cnt_lsb.o bn_error.o \
bn_mp_init_multi.o bn_mp_clear_multi.o bn_mp_prime_random.o bn_prime_sizes_tab.o
# make a Windows DLL via Cygwin
windll: $(OBJECTS)
gcc -mno-cygwin -mdll -o libtommath.dll -Wl,--out-implib=libtommath.dll.a -Wl,--export-all-symbols *.o
ranlib libtommath.dll.a
# build the test program using the windows DLL
test: $(OBJECTS) windll
gcc $(CFLAGS) demo/demo.c libtommath.dll.a -Wl,--enable-auto-import -o test -s
cd mtest ; $(CC) -O3 -fomit-frame-pointer -funroll-loops mtest.c -o mtest -s

4
makefile.msvc
View File

@ -2,7 +2,7 @@
#
#Tom St Denis
CFLAGS = /I. /Ox /DWIN32 /W3
CFLAGS = /I. /Ox /DWIN32 /W4
default: library
@ -28,7 +28,7 @@ bn_mp_toom_mul.obj bn_mp_toom_sqr.obj bn_mp_div_3.obj bn_s_mp_exptmod.obj \
bn_mp_reduce_2k.obj bn_mp_reduce_is_2k.obj bn_mp_reduce_2k_setup.obj \
bn_mp_radix_smap.obj bn_mp_read_radix.obj bn_mp_toradix.obj bn_mp_radix_size.obj \
bn_mp_fread.obj bn_mp_fwrite.obj bn_mp_cnt_lsb.obj bn_error.obj \
bn_mp_init_multi.obj bn_mp_clear_multi.obj
bn_mp_init_multi.obj bn_mp_clear_multi.obj bn_mp_prime_random.obj bn_prime_sizes_tab.obj
library: $(OBJECTS)
lib /out:tommath.lib $(OBJECTS)

23
mtest/test.c~
View File

@ -1,23 +0,0 @@
#include <stdio.h>
#include "mpi.c"
int main(void)
{
mp_int a, b;
int ix;
char buf[1024];
mp_init(&a);
mp_init(&b);
mp_set(&a, 0x1B);
mp_neg(&a, &a);
ix = 0;
mp_add_d(&a, ix, &b);
mp_toradix(&b, buf, 64);
printf("b == %s\n", buf);
return 0;
}

BIN
pics/design_process.tif Normal file
View File

Binary file not shown.

BIN
pics/expt_state.tif Normal file
View File

Binary file not shown.

BIN
pics/primality.tif Normal file
View File

Binary file not shown.

BIN
pics/radix.sxd Normal file
View File

Binary file not shown.

BIN
pics/sliding_window.tif Normal file
View File

Binary file not shown.

BIN
poster.pdf
View File

Binary file not shown.

1
poster.tex
View File

@ -2,7 +2,6 @@
\usepackage{amsmath, amssymb}
\usepackage{hyperref}
\begin{document}
\hspace*{-3in}
\begin{tabular}{llllll}
$c = a + b$ & {\tt mp\_add(\&a, \&b, \&c)} & $b = 2a$ & {\tt mp\_mul\_2(\&a, \&b)} & \\

346
pre_gen/mpi.c
View File

@ -72,11 +72,8 @@ fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
mp_int x, y, u, v, B, D;
int res, neg;
/* 2. [modified] if a,b are both even then return an error!
*
* That is if gcd(a,b) = 2**k * q then obviously there is no inverse.
*/
if (mp_iseven (a) == 1 && mp_iseven (b) == 1) {
/* 2. [modified] b must be odd */
if (mp_iseven (b) == 1) {
return MP_VAL;
}
@ -269,7 +266,7 @@ fast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
* that W[ix-1] have the carry cleared (see after the inner loop)
*/
register mp_digit mu;
mu = ((W[ix] & MP_MASK) * rho) & MP_MASK;
mu = (mp_digit) (((W[ix] & MP_MASK) * rho) & MP_MASK);
/* a = a + mu * m * b**i
*
@ -849,8 +846,7 @@ mp_abs (mp_int * a, mp_int * b)
#include <tommath.h>
/* high level addition (handles signs) */
int
mp_add (mp_int * a, mp_int * b, mp_int * c)
int mp_add (mp_int * a, mp_int * b, mp_int * c)
{
int sa, sb, res;
@ -1153,7 +1149,7 @@ mp_clear (mp_int * a)
memset (a->dp, 0, sizeof (mp_digit) * a->used);
/* free ram */
free (a->dp);
XFREE(a->dp);
/* reset members to make debugging easier */
a->dp = NULL;
@ -1255,8 +1251,7 @@ mp_cmp (mp_int * a, mp_int * b)
#include <tommath.h>
/* compare a digit */
int
mp_cmp_d (mp_int * a, mp_digit b)
int mp_cmp_d(mp_int * a, mp_digit b)
{
/* compare based on sign */
if (a->sign == MP_NEG) {
@ -1298,8 +1293,7 @@ mp_cmp_d (mp_int * a, mp_digit b)
#include <tommath.h>
/* compare maginitude of two ints (unsigned) */
int
mp_cmp_mag (mp_int * a, mp_int * b)
int mp_cmp_mag (mp_int * a, mp_int * b)
{
int n;
mp_digit *tmpa, *tmpb;
@ -1518,8 +1512,7 @@ mp_count_bits (mp_int * a)
* The overall algorithm is as described as
* 14.20 from HAC but fixed to treat these cases.
*/
int
mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
int mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
{
mp_int q, x, y, t1, t2;
int res, n, t, i, norm, neg;
@ -1722,8 +1715,7 @@ __Q:mp_clear (&q);
#include <tommath.h>
/* b = a/2 */
int
mp_div_2 (mp_int * a, mp_int * b)
int mp_div_2(mp_int * a, mp_int * b)
{
int x, res, oldused;
@ -1789,8 +1781,7 @@ mp_div_2 (mp_int * a, mp_int * b)
#include <tommath.h>
/* shift right by a certain bit count (store quotient in c, optional remainder in d) */
int
mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d)
int mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d)
{
mp_digit D, r, rr;
int x, res;
@ -2028,7 +2019,7 @@ mp_div_d (mp_int * a, mp_digit b, mp_int * c, mp_digit * d)
if (w >= b) {
t = (mp_digit)(w / b);
w = w % b;
w -= ((mp_word)t) * ((mp_word)b);
} else {
t = 0;
}
@ -2264,8 +2255,7 @@ mp_exch (mp_int * a, mp_int * b)
#include <tommath.h>
/* calculate c = a**b using a square-multiply algorithm */
int
mp_expt_d (mp_int * a, mp_digit b, mp_int * c)
int mp_expt_d (mp_int * a, mp_digit b, mp_int * c)
{
int res, x;
mp_int g;
@ -2325,8 +2315,7 @@ mp_expt_d (mp_int * a, mp_digit b, mp_int * c)
* embedded in the normal function but that wasted alot of stack space
* for nothing (since 99% of the time the Montgomery code would be called)
*/
int
mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
{
int dr;
@ -2768,24 +2757,24 @@ int mp_fwrite(mp_int *a, int radix, FILE *stream)
return MP_VAL;
}
buf = malloc(len);
buf = XMALLOC (len);
if (buf == NULL) {
return MP_MEM;
}
if ((err = mp_toradix(a, buf, radix)) != MP_OKAY) {
free(buf);
XFREE (buf);
return err;
}
for (x = 0; x < len; x++) {
if (fputc(buf[x], stream) == EOF) {
free(buf);
XFREE (buf);
return MP_VAL;
}
}
free(buf);
XFREE (buf);
return MP_OKAY;
}
@ -2810,8 +2799,7 @@ int mp_fwrite(mp_int *a, int radix, FILE *stream)
#include <tommath.h>
/* Greatest Common Divisor using the binary method */
int
mp_gcd (mp_int * a, mp_int * b, mp_int * c)
int mp_gcd (mp_int * a, mp_int * b, mp_int * c)
{
mp_int u, v;
int k, u_lsb, v_lsb, res;
@ -2922,13 +2910,11 @@ __U:mp_clear (&v);
#include <tommath.h>
/* grow as required */
int
mp_grow (mp_int * a, int size)
int mp_grow (mp_int * a, int size)
{
int i;
mp_digit *tmp;
/* if the alloc size is smaller alloc more ram */
if (a->alloc < size) {
/* ensure there are always at least MP_PREC digits extra on top */
@ -2940,7 +2926,7 @@ mp_grow (mp_int * a, int size)
* in case the operation failed we don't want
* to overwrite the dp member of a.
*/
tmp = OPT_CAST realloc (a->dp, sizeof (mp_digit) * size);
tmp = OPT_CAST XREALLOC (a->dp, sizeof (mp_digit) * size);
if (tmp == NULL) {
/* reallocation failed but "a" is still valid [can be freed] */
return MP_MEM;
@ -2979,11 +2965,10 @@ mp_grow (mp_int * a, int size)
#include <tommath.h>
/* init a new bigint */
int
mp_init (mp_int * a)
int mp_init (mp_int * a)
{
/* allocate memory required and clear it */
a->dp = OPT_CAST calloc (sizeof (mp_digit), MP_PREC);
a->dp = OPT_CAST XCALLOC (sizeof (mp_digit), MP_PREC);
if (a->dp == NULL) {
return MP_MEM;
}
@ -3017,8 +3002,7 @@ mp_init (mp_int * a)
#include <tommath.h>
/* creates "a" then copies b into it */
int
mp_init_copy (mp_int * a, mp_int * b)
int mp_init_copy (mp_int * a, mp_int * b)
{
int res;
@ -3105,14 +3089,13 @@ int mp_init_multi(mp_int *mp, ...)
#include <tommath.h>
/* init an mp_init for a given size */
int
mp_init_size (mp_int * a, int size)
int mp_init_size (mp_int * a, int size)
{
/* pad size so there are always extra digits */
size += (MP_PREC * 2) - (size % MP_PREC);
/* alloc mem */
a->dp = OPT_CAST calloc (sizeof (mp_digit), size);
a->dp = OPT_CAST XCALLOC (sizeof (mp_digit), size);
if (a->dp == NULL) {
return MP_MEM;
}
@ -3143,8 +3126,7 @@ mp_init_size (mp_int * a, int size)
#include <tommath.h>
/* hac 14.61, pp608 */
int
mp_invmod (mp_int * a, mp_int * b, mp_int * c)
int mp_invmod (mp_int * a, mp_int * b, mp_int * c)
{
mp_int x, y, u, v, A, B, C, D;
int res;
@ -3324,8 +3306,7 @@ __ERR:mp_clear_multi (&x, &y, &u, &v, &A, &B, &C, &D, NULL);
/* computes the jacobi c = (a | n) (or Legendre if n is prime)
* HAC pp. 73 Algorithm 2.149
*/
int
mp_jacobi (mp_int * a, mp_int * p, int *c)
int mp_jacobi (mp_int * a, mp_int * p, int *c)
{
mp_int a1, p1;
int k, s, r, res;
@ -3715,8 +3696,7 @@ ERR:
#include <tommath.h>
/* computes least common multiple as |a*b|/(a, b) */
int
mp_lcm (mp_int * a, mp_int * b, mp_int * c)
int mp_lcm (mp_int * a, mp_int * b, mp_int * c)
{
int res;
mp_int t1, t2;
@ -3774,8 +3754,7 @@ __T:
#include <tommath.h>
/* shift left a certain amount of digits */
int
mp_lshd (mp_int * a, int b)
int mp_lshd (mp_int * a, int b)
{
int x, res;
@ -4058,7 +4037,7 @@ mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
* following inner loop to reduce the
* input one digit at a time
*/
mu = ((mp_word)x->dp[ix]) * ((mp_word)rho) & MP_MASK;
mu = (mp_digit) (((mp_word)x->dp[ix]) * ((mp_word)rho) & MP_MASK);
/* a = a + mu * m * b**i */
{
@ -4195,8 +4174,7 @@ mp_montgomery_setup (mp_int * n, mp_digit * rho)
#include <tommath.h>
/* high level multiplication (handles sign) */
int
mp_mul (mp_int * a, mp_int * b, mp_int * c)
int mp_mul (mp_int * a, mp_int * b, mp_int * c)
{
int res, neg;
neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
@ -4249,8 +4227,7 @@ mp_mul (mp_int * a, mp_int * b, mp_int * c)
#include <tommath.h>
/* b = a*2 */
int
mp_mul_2 (mp_int * a, mp_int * b)
int mp_mul_2(mp_int * a, mp_int * b)
{
int x, res, oldused;
@ -4330,8 +4307,7 @@ mp_mul_2 (mp_int * a, mp_int * b)
#include <tommath.h>
/* shift left by a certain bit count */
int
mp_mul_2d (mp_int * a, int b, mp_int * c)
int mp_mul_2d (mp_int * a, int b, mp_int * c)
{
mp_digit d;
int res;
@ -4496,7 +4472,6 @@ mp_mulmod (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
int res;
mp_int t;
if ((res = mp_init (&t)) != MP_OKAY) {
return res;
}
@ -4539,8 +4514,7 @@ mp_mulmod (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
* each step involves a fair bit. This is not meant to
* find huge roots [square and cube, etc].
*/
int
mp_n_root (mp_int * a, mp_digit b, mp_int * c)
int mp_n_root (mp_int * a, mp_digit b, mp_int * c)
{
mp_int t1, t2, t3;
int res, neg;
@ -4661,8 +4635,7 @@ __T1:mp_clear (&t1);
#include <tommath.h>
/* b = -a */
int
mp_neg (mp_int * a, mp_int * b)
int mp_neg (mp_int * a, mp_int * b)
{
int res;
if ((res = mp_copy (a, b)) != MP_OKAY) {
@ -4694,8 +4667,7 @@ mp_neg (mp_int * a, mp_int * b)
#include <tommath.h>
/* OR two ints together */
int
mp_or (mp_int * a, mp_int * b, mp_int * c)
int mp_or (mp_int * a, mp_int * b, mp_int * c)
{
int res, ix, px;
mp_int t, *x;
@ -4750,14 +4722,13 @@ mp_or (mp_int * a, mp_int * b, mp_int * c)
*
* Sets result to 1 if the congruence holds, or zero otherwise.
*/
int
mp_prime_fermat (mp_int * a, mp_int * b, int *result)
int mp_prime_fermat (mp_int * a, mp_int * b, int *result)
{
mp_int t;
int err;
/* default to composite */
*result = 0;
*result = MP_NO;
/* ensure b > 1 */
if (mp_cmp_d(b, 1) != MP_GT) {
@ -4776,7 +4747,7 @@ mp_prime_fermat (mp_int * a, mp_int * b, int *result)
/* is it equal to b? */
if (mp_cmp (&t, b) == MP_EQ) {
*result = 1;
*result = MP_YES;
}
err = MP_OKAY;
@ -4808,14 +4779,13 @@ __T:mp_clear (&t);
*
* sets result to 0 if not, 1 if yes
*/
int
mp_prime_is_divisible (mp_int * a, int *result)
int mp_prime_is_divisible (mp_int * a, int *result)
{
int err, ix;
mp_digit res;
/* default to not */
*result = 0;
*result = MP_NO;
for (ix = 0; ix < PRIME_SIZE; ix++) {
/* what is a mod __prime_tab[ix] */
@ -4825,7 +4795,7 @@ mp_prime_is_divisible (mp_int * a, int *result)
/* is the residue zero? */
if (res == 0) {
*result = 1;
*result = MP_YES;
return MP_OKAY;
}
}
@ -4859,14 +4829,13 @@ mp_prime_is_divisible (mp_int * a, int *result)
*
* Sets result to 1 if probably prime, 0 otherwise
*/
int
mp_prime_is_prime (mp_int * a, int t, int *result)
int mp_prime_is_prime (mp_int * a, int t, int *result)
{
mp_int b;
int ix, err, res;
/* default to no */
*result = 0;
*result = MP_NO;
/* valid value of t? */
if (t <= 0 || t > PRIME_SIZE) {
@ -4887,7 +4856,7 @@ mp_prime_is_prime (mp_int * a, int t, int *result)
}
/* return if it was trivially divisible */
if (res == 1) {
if (res == MP_YES) {
return MP_OKAY;
}
@ -4904,13 +4873,13 @@ mp_prime_is_prime (mp_int * a, int t, int *result)
goto __B;
}
if (res == 0) {
if (res == MP_NO) {
goto __B;
}
}
/* passed the test */
*result = 1;
*result = MP_YES;
__B:mp_clear (&b);
return err;
}
@ -4941,14 +4910,13 @@ __B:mp_clear (&b);
* Randomly the chance of error is no more than 1/4 and often
* very much lower.
*/
int
mp_prime_miller_rabin (mp_int * a, mp_int * b, int *result)
int mp_prime_miller_rabin (mp_int * a, mp_int * b, int *result)
{
mp_int n1, y, r;
int s, j, err;
/* default */
*result = 0;
*result = MP_NO;
/* ensure b > 1 */
if (mp_cmp_d(b, 1) != MP_GT) {
@ -5010,7 +4978,7 @@ mp_prime_miller_rabin (mp_int * a, mp_int * b, int *result)
}
/* probably prime now */
*result = 1;
*result = MP_YES;
__Y:mp_clear (&y);
__R:mp_clear (&r);
__N1:mp_clear (&n1);
@ -5168,12 +5136,12 @@ int mp_prime_next_prime(mp_int *a, int t, int bbs_style)
if ((err = mp_prime_miller_rabin(a, &b, &res)) != MP_OKAY) {
goto __ERR;
}
if (res == 0) {
if (res == MP_NO) {
break;
}
}
if (res == 1) {
if (res == MP_YES) {
break;
}
}
@ -5187,6 +5155,84 @@ __ERR:
/* End: bn_mp_prime_next_prime.c */
/* Start: bn_mp_prime_random.c */
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
* LibTomMath is a library that provides multiple-precision
* integer arithmetic as well as number theoretic functionality.
*
* The library was designed directly after the MPI library by
* Michael Fromberger but has been written from scratch with
* additional optimizations in place.
*
* The library is free for all purposes without any express
* guarantee it works.
*
* Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
*/
#include <tommath.h>
/* makes a truly random prime of a given size (bytes),
* call with bbs = 1 if you want it to be congruent to 3 mod 4
*
* You have to supply a callback which fills in a buffer with random bytes. "dat" is a parameter you can
* have passed to the callback (e.g. a state or something). This function doesn't use "dat" itself
* so it can be NULL
*
* The prime generated will be larger than 2^(8*size).
*/
/* this sole function may hold the key to enslaving all mankind! */
int mp_prime_random(mp_int *a, int t, int size, int bbs, ltm_prime_callback cb, void *dat)
{
unsigned char *tmp;
int res, err;
/* sanity check the input */
if (size <= 0) {
return MP_VAL;
}
/* we need a buffer of size+1 bytes */
tmp = XMALLOC(size+1);
if (tmp == NULL) {
return MP_MEM;
}
/* fix MSB */
tmp[0] = 1;
do {
/* read the bytes */
if (cb(tmp+1, size, dat) != size) {
err = MP_VAL;
goto error;
}
/* fix the LSB */
tmp[size] |= (bbs ? 3 : 1);
/* read it in */
if ((err = mp_read_unsigned_bin(a, tmp, size+1)) != MP_OKAY) {
goto error;
}
/* is it prime? */
if ((err = mp_prime_is_prime(a, t, &res)) != MP_OKAY) {
goto error;
}
} while (res == MP_NO);
err = MP_OKAY;
error:
XFREE(tmp);
return err;
}
/* End: bn_mp_prime_random.c */
/* Start: bn_mp_radix_size.c */
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
@ -5217,28 +5263,36 @@ mp_radix_size (mp_int * a, int radix)
return mp_count_bits (a) + (a->sign == MP_NEG ? 1 : 0) + 1;
}
/* make sure the radix is in range */
if (radix < 2 || radix > 64) {
return 0;
}
/* init a copy of the input */
if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
return 0;
}
/* digs is the digit count */
digs = 0;
/* if it's negative add one for the sign */
if (t.sign == MP_NEG) {
++digs;
t.sign = MP_ZPOS;
}
/* fetch out all of the digits */
while (mp_iszero (&t) == 0) {
if ((res = mp_div_d (&t, (mp_digit) radix, &t, &d)) != MP_OKAY) {
mp_clear (&t);
return 0;
return res;
}
++digs;
}
mp_clear (&t);
/* return digs + 1, the 1 is for the NULL byte that would be required. */
return digs + 1;
}
@ -5707,10 +5761,9 @@ mp_reduce_2k_setup(mp_int *a, mp_digit *d)
#include <tommath.h>
/* determines if mp_reduce_2k can be used */
int
mp_reduce_is_2k(mp_int *a)
int mp_reduce_is_2k(mp_int *a)
{
int ix, iy;
int ix, iy, iz, iw;
if (a->used == 0) {
return 0;
@ -5718,11 +5771,19 @@ mp_reduce_is_2k(mp_int *a)
return 1;
} else if (a->used > 1) {
iy = mp_count_bits(a);
iz = 1;
iw = 1;
/* Test every bit from the second digit up, must be 1 */
for (ix = DIGIT_BIT; ix < iy; ix++) {
if ((a->dp[ix/DIGIT_BIT] &
((mp_digit)1 << (mp_digit)(ix % DIGIT_BIT))) == 0) {
if ((a->dp[iw] & iz) == 0) {
return 0;
}
iz <<= 1;
if (iz > (int)MP_MASK) {
++iw;
iz = 1;
}
}
}
return 1;
@ -5782,8 +5843,7 @@ mp_reduce_setup (mp_int * a, mp_int * b)
#include <tommath.h>
/* shift right a certain amount of digits */
void
mp_rshd (mp_int * a, int b)
void mp_rshd (mp_int * a, int b)
{
int x;
@ -5853,8 +5913,7 @@ mp_rshd (mp_int * a, int b)
#include <tommath.h>
/* set to a digit */
void
mp_set (mp_int * a, mp_digit b)
void mp_set (mp_int * a, mp_digit b)
{
mp_zero (a);
a->dp[0] = b & MP_MASK;
@ -5881,8 +5940,7 @@ mp_set (mp_int * a, mp_digit b)
#include <tommath.h>
/* set a 32-bit const */
int
mp_set_int (mp_int * a, unsigned int b)
int mp_set_int (mp_int * a, unsigned long b)
{
int x, res;
@ -5928,13 +5986,14 @@ mp_set_int (mp_int * a, unsigned int b)
#include <tommath.h>
/* shrink a bignum */
int
mp_shrink (mp_int * a)
int mp_shrink (mp_int * a)
{
mp_digit *tmp;
if (a->alloc != a->used) {
if ((a->dp = OPT_CAST realloc (a->dp, sizeof (mp_digit) * a->used)) == NULL) {
if ((tmp = OPT_CAST XREALLOC (a->dp, sizeof (mp_digit) * a->used)) == NULL) {
return MP_MEM;
}
a->dp = tmp;
a->alloc = a->used;
}
return MP_OKAY;
@ -6841,6 +6900,7 @@ mp_toradix (mp_int * a, char *str, int radix)
mp_digit d;
char *_s = str;
/* check range of the radix */
if (radix < 2 || radix > 64) {
return MP_VAL;
}
@ -6879,8 +6939,7 @@ mp_toradix (mp_int * a, char *str, int radix)
bn_reverse ((unsigned char *)_s, digs);
/* append a NULL so the string is properly terminated */
*str++ = '\0';
*str = '\0';
mp_clear (&t);
return MP_OKAY;
@ -6993,6 +7052,79 @@ mp_zero (mp_int * a)
/* End: bn_mp_zero.c */
/* Start: bn_prime_sizes_tab.c */
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
* LibTomMath is a library that provides multiple-precision
* integer arithmetic as well as number theoretic functionality.
*
* The library was designed directly after the MPI library by
* Michael Fromberger but has been written from scratch with
* additional optimizations in place.
*
* The library is free for all purposes without any express
* guarantee it works.
*
* Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
*/
#include <tommath.h>
/* this table gives the # of rabin miller trials for a prob of failure lower than 2^-96 */
static const struct {
int k, t;
} sizes[] = {
{ 128, 28 },
{ 256, 16 },
{ 384, 10 },
{ 512, 7 },
{ 640, 6 },
{ 768, 5 },
{ 896, 4 },
{ 1024, 4 },
{ 1152, 3 },
{ 1280, 3 },
{ 1408, 3 },
{ 1536, 3 },
{ 1664, 3 },
{ 1792, 2 },
{ 1920, 2 },
{ 2048, 2 },
{ 2176, 2 },
{ 2304, 2 },
{ 2432, 2 },
{ 2560, 2 },
{ 2688, 2 },
{ 2816, 2 },
{ 2944, 2 },
{ 3072, 2 },
{ 3200, 2 },
{ 3328, 2 },
{ 3456, 2 },
{ 3584, 2 },
{ 3712, 2 },
{ 3840, 1 },
{ 3968, 1 },
{ 4096, 1 } };
/* returns # of RM trials required for a given bit size */
int mp_prime_rabin_miller_trials(int size)
{
int x;
for (x = 0; x < (int)(sizeof(sizes)/(sizeof(sizes[0]))); x++) {
if (sizes[x].k == size) {
return sizes[x].t;
} else if (sizes[x].k > size) {
return (x == 0) ? sizes[0].t : sizes[x - 1].t;
}
}
return 1;
}
/* End: bn_prime_sizes_tab.c */
/* Start: bn_prime_tab.c */
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
@ -7631,6 +7763,8 @@ s_mp_sqr (mp_int * a, mp_int * b)
if ((res = mp_init_size (&t, 2*pa + 1)) != MP_OKAY) {
return res;
}
/* default used is maximum possible size */
t.used = 2*pa + 1;
for (ix = 0; ix < pa; ix++) {
@ -7640,20 +7774,20 @@ s_mp_sqr (mp_int * a, mp_int * b)
((mp_word)a->dp[ix])*((mp_word)a->dp[ix]);
/* store lower part in result */
t.dp[2*ix] = (mp_digit) (r & ((mp_word) MP_MASK));
t.dp[ix+ix] = (mp_digit) (r & ((mp_word) MP_MASK));
/* get the carry */
u = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
u = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
/* left hand side of A[ix] * A[iy] */
tmpx = a->dp[ix];
tmpx = a->dp[ix];
/* alias for where to store the results */
tmpt = t.dp + (2*ix + 1);
tmpt = t.dp + (2*ix + 1);
for (iy = ix + 1; iy < pa; iy++) {
/* first calculate the product */
r = ((mp_word)tmpx) * ((mp_word)a->dp[iy]);
r = ((mp_word)tmpx) * ((mp_word)a->dp[iy]);
/* now calculate the double precision result, note we use
* addition instead of *2 since it's easier to optimize

53
tommath.h
View File

@ -91,6 +91,24 @@ extern "C" {
#endif
#endif
/* define heap macros */
#ifndef CRYPT
/* default to libc stuff */
#ifndef XMALLOC
#define XMALLOC malloc
#define XFREE free
#define XREALLOC realloc
#define XCALLOC calloc
#endif
/* prototypes for our heap functions */
extern void *XMALLOC(size_t n);
extern void *REALLOC(void *p, size_t n);
extern void *XCALLOC(size_t n, size_t s);
extern void XFREE(void *p);
#endif
/* otherwise the bits per digit is calculated automatically from the size of a mp_digit */
#ifndef DIGIT_BIT
#define DIGIT_BIT ((int)((CHAR_BIT * sizeof(mp_digit) - 1))) /* bits per digit */
@ -113,6 +131,9 @@ extern "C" {
#define MP_VAL -3 /* invalid input */
#define MP_RANGE MP_VAL
#define MP_YES 1 /* yes response */
#define MP_NO 0 /* no response */
typedef int mp_err;
/* you'll have to tune these... */
@ -130,11 +151,16 @@ extern int KARATSUBA_MUL_CUTOFF,
/* size of comba arrays, should be at least 2 * 2**(BITS_PER_WORD - BITS_PER_DIGIT*2) */
#define MP_WARRAY (1 << (sizeof(mp_word) * CHAR_BIT - 2 * DIGIT_BIT + 1))
/* the infamous mp_int structure */
typedef struct {
int used, alloc, sign;
mp_digit *dp;
} mp_int;
/* callback for mp_prime_random, should fill dst with random bytes and return how many read [upto len] */
typedef int ltm_prime_callback(unsigned char *dst, int len, void *dat);
#define USED(m) ((m)->used)
#define DIGIT(m,k) ((m)->dp[(k)])
#define SIGN(m) ((m)->sign)
@ -168,9 +194,9 @@ int mp_grow(mp_int *a, int size);
int mp_init_size(mp_int *a, int size);
/* ---> Basic Manipulations <--- */
#define mp_iszero(a) (((a)->used == 0) ? 1 : 0)
#define mp_iseven(a) (((a)->used > 0 && (((a)->dp[0] & 1) == 0)) ? 1 : 0)
#define mp_isodd(a) (((a)->used > 0 && (((a)->dp[0] & 1) == 1)) ? 1 : 0)
#define mp_iszero(a) (((a)->used == 0) ? MP_YES : MP_NO)
#define mp_iseven(a) (((a)->used > 0 && (((a)->dp[0] & 1) == 0)) ? MP_YES : MP_NO)
#define mp_isodd(a) (((a)->used > 0 && (((a)->dp[0] & 1) == 1)) ? MP_YES : MP_NO)
/* set to zero */
void mp_zero(mp_int *a);
@ -179,7 +205,7 @@ void mp_zero(mp_int *a);
void mp_set(mp_int *a, mp_digit b);
/* set a 32-bit const */
int mp_set_int(mp_int *a, unsigned int b);
int mp_set_int(mp_int *a, unsigned long b);
/* copy, b = a */
int mp_copy(mp_int *a, mp_int *b);
@ -219,6 +245,8 @@ int mp_2expt(mp_int *a, int b);
/* Counts the number of lsbs which are zero before the first zero bit */
int mp_cnt_lsb(mp_int *a);
/* I Love Earth! */
/* makes a pseudo-random int of a given size */
int mp_rand(mp_int *a, int digits);
@ -392,6 +420,11 @@ int mp_prime_fermat(mp_int *a, mp_int *b, int *result);
*/
int mp_prime_miller_rabin(mp_int *a, mp_int *b, int *result);
/* This gives [for a given bit size] the number of trials required
* such that Miller-Rabin gives a prob of failure lower than 2^-96
*/
int mp_prime_rabin_miller_trials(int size);
/* performs t rounds of Miller-Rabin on "a" using the first
* t prime bases. Also performs an initial sieve of trial
* division. Determines if "a" is prime with probability
@ -408,6 +441,18 @@ int mp_prime_is_prime(mp_int *a, int t, int *result);
*/
int mp_prime_next_prime(mp_int *a, int t, int bbs_style);
/* makes a truly random prime of a given size (bytes),
* call with bbs = 1 if you want it to be congruent to 3 mod 4
*
* You have to supply a callback which fills in a buffer with random bytes. "dat" is a parameter you can
* have passed to the callback (e.g. a state or something). This function doesn't use "dat" itself
* so it can be NULL
*
* The prime generated will be larger than 2^(8*size).
*/
int mp_prime_random(mp_int *a, int t, int size, int bbs, ltm_prime_callback cb, void *dat);
/* ---> radix conversion <--- */
int mp_count_bits(mp_int *a);

139
tommath.out Normal file
View File

@ -0,0 +1,139 @@
\BOOKMARK [0][-]{chapter.1}{Introduction}{}
\BOOKMARK [1][-]{section.1.1}{Multiple Precision Arithmetic}{chapter.1}
\BOOKMARK [2][-]{subsection.1.1.1}{The Need for Multiple Precision Arithmetic}{section.1.1}
\BOOKMARK [2][-]{subsection.1.1.2}{What is Multiple Precision Arithmetic?}{section.1.1}
\BOOKMARK [2][-]{subsection.1.1.3}{Benefits of Multiple Precision Arithmetic}{section.1.1}
\BOOKMARK [1][-]{section.1.2}{Purpose of This Text}{chapter.1}
\BOOKMARK [1][-]{section.1.3}{Discussion and Notation}{chapter.1}
\BOOKMARK [2][-]{subsection.1.3.1}{Notation}{section.1.3}
\BOOKMARK [2][-]{subsection.1.3.2}{Precision Notation}{section.1.3}
\BOOKMARK [2][-]{subsection.1.3.3}{Algorithm Inputs and Outputs}{section.1.3}
\BOOKMARK [2][-]{subsection.1.3.4}{Mathematical Expressions}{section.1.3}
\BOOKMARK [2][-]{subsection.1.3.5}{Work Effort}{section.1.3}
\BOOKMARK [1][-]{section.1.4}{Exercises}{chapter.1}
\BOOKMARK [0][-]{chapter.2}{Introduction to LibTomMath}{}
\BOOKMARK [1][-]{section.2.1}{What is LibTomMath?}{chapter.2}
\BOOKMARK [1][-]{section.2.2}{Goals of LibTomMath}{chapter.2}
\BOOKMARK [1][-]{section.2.3}{Choice of LibTomMath}{chapter.2}
\BOOKMARK [2][-]{subsection.2.3.1}{Code Base}{section.2.3}
\BOOKMARK [2][-]{subsection.2.3.2}{API Simplicity}{section.2.3}
\BOOKMARK [2][-]{subsection.2.3.3}{Optimizations}{section.2.3}
\BOOKMARK [2][-]{subsection.2.3.4}{Portability and Stability}{section.2.3}
\BOOKMARK [2][-]{subsection.2.3.5}{Choice}{section.2.3}
\BOOKMARK [0][-]{chapter.3}{Getting Started}{}
\BOOKMARK [1][-]{section.3.1}{Library Basics}{chapter.3}
\BOOKMARK [1][-]{section.3.2}{What is a Multiple Precision Integer?}{chapter.3}
\BOOKMARK [2][-]{subsection.3.2.1}{The mp\137int Structure}{section.3.2}
\BOOKMARK [1][-]{section.3.3}{Argument Passing}{chapter.3}
\BOOKMARK [1][-]{section.3.4}{Return Values}{chapter.3}
\BOOKMARK [1][-]{section.3.5}{Initialization and Clearing}{chapter.3}
\BOOKMARK [2][-]{subsection.3.5.1}{Initializing an mp\137int}{section.3.5}
\BOOKMARK [2][-]{subsection.3.5.2}{Clearing an mp\137int}{section.3.5}
\BOOKMARK [1][-]{section.3.6}{Maintenance Algorithms}{chapter.3}
\BOOKMARK [2][-]{subsection.3.6.1}{Augmenting an mp\137int's Precision}{section.3.6}
\BOOKMARK [2][-]{subsection.3.6.2}{Initializing Variable Precision mp\137ints}{section.3.6}
\BOOKMARK [2][-]{subsection.3.6.3}{Multiple Integer Initializations and Clearings}{section.3.6}
\BOOKMARK [2][-]{subsection.3.6.4}{Clamping Excess Digits}{section.3.6}
\BOOKMARK [0][-]{chapter.4}{Basic Operations}{}
\BOOKMARK [1][-]{section.4.1}{Introduction}{chapter.4}
\BOOKMARK [1][-]{section.4.2}{Assigning Values to mp\137int Structures}{chapter.4}
\BOOKMARK [2][-]{subsection.4.2.1}{Copying an mp\137int}{section.4.2}
\BOOKMARK [2][-]{subsection.4.2.2}{Creating a Clone}{section.4.2}
\BOOKMARK [1][-]{section.4.3}{Zeroing an Integer}{chapter.4}
\BOOKMARK [1][-]{section.4.4}{Sign Manipulation}{chapter.4}
\BOOKMARK [2][-]{subsection.4.4.1}{Absolute Value}{section.4.4}
\BOOKMARK [2][-]{subsection.4.4.2}{Integer Negation}{section.4.4}
\BOOKMARK [1][-]{section.4.5}{Small Constants}{chapter.4}
\BOOKMARK [2][-]{subsection.4.5.1}{Setting Small Constants}{section.4.5}
\BOOKMARK [2][-]{subsection.4.5.2}{Setting Large Constants}{section.4.5}
\BOOKMARK [1][-]{section.4.6}{Comparisons}{chapter.4}
\BOOKMARK [2][-]{subsection.4.6.1}{Unsigned Comparisions}{section.4.6}
\BOOKMARK [2][-]{subsection.4.6.2}{Signed Comparisons}{section.4.6}
\BOOKMARK [0][-]{chapter.5}{Basic Arithmetic}{}
\BOOKMARK [1][-]{section.5.1}{Introduction}{chapter.5}
\BOOKMARK [1][-]{section.5.2}{Addition and Subtraction}{chapter.5}
\BOOKMARK [2][-]{subsection.5.2.1}{Low Level Addition}{section.5.2}
\BOOKMARK [2][-]{subsection.5.2.2}{Low Level Subtraction}{section.5.2}
\BOOKMARK [2][-]{subsection.5.2.3}{High Level Addition}{section.5.2}
\BOOKMARK [2][-]{subsection.5.2.4}{High Level Subtraction}{section.5.2}
\BOOKMARK [1][-]{section.5.3}{Bit and Digit Shifting}{chapter.5}
\BOOKMARK [2][-]{subsection.5.3.1}{Multiplication by Two}{section.5.3}
\BOOKMARK [2][-]{subsection.5.3.2}{Division by Two}{section.5.3}
\BOOKMARK [1][-]{section.5.4}{Polynomial Basis Operations}{chapter.5}
\BOOKMARK [2][-]{subsection.5.4.1}{Multiplication by x}{section.5.4}
\BOOKMARK [2][-]{subsection.5.4.2}{Division by x}{section.5.4}
\BOOKMARK [1][-]{section.5.5}{Powers of Two}{chapter.5}
\BOOKMARK [2][-]{subsection.5.5.1}{Multiplication by Power of Two}{section.5.5}
\BOOKMARK [2][-]{subsection.5.5.2}{Division by Power of Two}{section.5.5}
\BOOKMARK [2][-]{subsection.5.5.3}{Remainder of Division by Power of Two}{section.5.5}
\BOOKMARK [0][-]{chapter.6}{Multiplication and Squaring}{}
\BOOKMARK [1][-]{section.6.1}{The Multipliers}{chapter.6}
\BOOKMARK [1][-]{section.6.2}{Multiplication}{chapter.6}
\BOOKMARK [2][-]{subsection.6.2.1}{The Baseline Multiplication}{section.6.2}
\BOOKMARK [2][-]{subsection.6.2.2}{Faster Multiplication by the ``Comba'' Method}{section.6.2}
\BOOKMARK [2][-]{subsection.6.2.3}{Polynomial Basis Multiplication}{section.6.2}
\BOOKMARK [2][-]{subsection.6.2.4}{Karatsuba Multiplication}{section.6.2}
\BOOKMARK [2][-]{subsection.6.2.5}{Toom-Cook 3-Way Multiplication}{section.6.2}
\BOOKMARK [2][-]{subsection.6.2.6}{Signed Multiplication}{section.6.2}
\BOOKMARK [1][-]{section.6.3}{Squaring}{chapter.6}
\BOOKMARK [2][-]{subsection.6.3.1}{The Baseline Squaring Algorithm}{section.6.3}
\BOOKMARK [2][-]{subsection.6.3.2}{Faster Squaring by the ``Comba'' Method}{section.6.3}
\BOOKMARK [2][-]{subsection.6.3.3}{Polynomial Basis Squaring}{section.6.3}
\BOOKMARK [2][-]{subsection.6.3.4}{Karatsuba Squaring}{section.6.3}
\BOOKMARK [2][-]{subsection.6.3.5}{Toom-Cook Squaring}{section.6.3}
\BOOKMARK [2][-]{subsection.6.3.6}{High Level Squaring}{section.6.3}
\BOOKMARK [0][-]{chapter.7}{Modular Reduction}{}
\BOOKMARK [1][-]{section.7.1}{Basics of Modular Reduction}{chapter.7}
\BOOKMARK [1][-]{section.7.2}{The Barrett Reduction}{chapter.7}
\BOOKMARK [2][-]{subsection.7.2.1}{Fixed Point Arithmetic}{section.7.2}
\BOOKMARK [2][-]{subsection.7.2.2}{Choosing a Radix Point}{section.7.2}
\BOOKMARK [2][-]{subsection.7.2.3}{Trimming the Quotient}{section.7.2}
\BOOKMARK [2][-]{subsection.7.2.4}{Trimming the Residue}{section.7.2}
\BOOKMARK [2][-]{subsection.7.2.5}{The Barrett Algorithm}{section.7.2}
\BOOKMARK [2][-]{subsection.7.2.6}{The Barrett Setup Algorithm}{section.7.2}
\BOOKMARK [1][-]{section.7.3}{The Montgomery Reduction}{chapter.7}
\BOOKMARK [2][-]{subsection.7.3.1}{Digit Based Montgomery Reduction}{section.7.3}
\BOOKMARK [2][-]{subsection.7.3.2}{Baseline Montgomery Reduction}{section.7.3}
\BOOKMARK [2][-]{subsection.7.3.3}{Faster ``Comba'' Montgomery Reduction}{section.7.3}
\BOOKMARK [2][-]{subsection.7.3.4}{Montgomery Setup}{section.7.3}
\BOOKMARK [1][-]{section.7.4}{The Diminished Radix Algorithm}{chapter.7}
\BOOKMARK [2][-]{subsection.7.4.1}{Choice of Moduli}{section.7.4}
\BOOKMARK [2][-]{subsection.7.4.2}{Choice of k}{section.7.4}
\BOOKMARK [2][-]{subsection.7.4.3}{Restricted Diminished Radix Reduction}{section.7.4}
\BOOKMARK [2][-]{subsection.7.4.4}{Unrestricted Diminished Radix Reduction}{section.7.4}
\BOOKMARK [1][-]{section.7.5}{Algorithm Comparison}{chapter.7}
\BOOKMARK [0][-]{chapter.8}{Exponentiation}{}
\BOOKMARK [1][-]{section.8.1}{Exponentiation Basics}{chapter.8}
\BOOKMARK [2][-]{subsection.8.1.1}{Single Digit Exponentiation}{section.8.1}
\BOOKMARK [1][-]{section.8.2}{k-ary Exponentiation}{chapter.8}
\BOOKMARK [2][-]{subsection.8.2.1}{Optimal Values of k}{section.8.2}
\BOOKMARK [2][-]{subsection.8.2.2}{Sliding-Window Exponentiation}{section.8.2}
\BOOKMARK [1][-]{section.8.3}{Modular Exponentiation}{chapter.8}
\BOOKMARK [2][-]{subsection.8.3.1}{Barrett Modular Exponentiation}{section.8.3}
\BOOKMARK [1][-]{section.8.4}{Quick Power of Two}{chapter.8}
\BOOKMARK [0][-]{chapter.9}{Higher Level Algorithms}{}
\BOOKMARK [1][-]{section.9.1}{Integer Division with Remainder}{chapter.9}
\BOOKMARK [2][-]{subsection.9.1.1}{Quotient Estimation}{section.9.1}
\BOOKMARK [2][-]{subsection.9.1.2}{Normalized Integers}{section.9.1}
\BOOKMARK [2][-]{subsection.9.1.3}{Radix- Division with Remainder}{section.9.1}
\BOOKMARK [1][-]{section.9.2}{Single Digit Helpers}{chapter.9}
\BOOKMARK [2][-]{subsection.9.2.1}{Single Digit Addition and Subtraction}{section.9.2}
\BOOKMARK [2][-]{subsection.9.2.2}{Single Digit Multiplication}{section.9.2}
\BOOKMARK [2][-]{subsection.9.2.3}{Single Digit Division}{section.9.2}
\BOOKMARK [2][-]{subsection.9.2.4}{Single Digit Root Extraction}{section.9.2}
\BOOKMARK [1][-]{section.9.3}{Random Number Generation}{chapter.9}
\BOOKMARK [1][-]{section.9.4}{Formatted Representations}{chapter.9}
\BOOKMARK [2][-]{subsection.9.4.1}{Reading Radix-n Input}{section.9.4}
\BOOKMARK [2][-]{subsection.9.4.2}{Generating Radix-n Output}{section.9.4}
\BOOKMARK [0][-]{chapter.10}{Number Theoretic Algorithms}{}
\BOOKMARK [1][-]{section.10.1}{Greatest Common Divisor}{chapter.10}
\BOOKMARK [2][-]{subsection.10.1.1}{Complete Greatest Common Divisor}{section.10.1}
\BOOKMARK [1][-]{section.10.2}{Least Common Multiple}{chapter.10}
\BOOKMARK [1][-]{section.10.3}{Jacobi Symbol Computation}{chapter.10}
\BOOKMARK [2][-]{subsection.10.3.1}{Jacobi Symbol}{section.10.3}
\BOOKMARK [1][-]{section.10.4}{Modular Inverse}{chapter.10}
\BOOKMARK [2][-]{subsection.10.4.1}{General Case}{section.10.4}
\BOOKMARK [1][-]{section.10.5}{Primality Tests}{chapter.10}
\BOOKMARK [2][-]{subsection.10.5.1}{Trial Division}{section.10.5}
\BOOKMARK [2][-]{subsection.10.5.2}{The Fermat Test}{section.10.5}
\BOOKMARK [2][-]{subsection.10.5.3}{The Miller-Rabin Test}{section.10.5}

21175
tommath.pdf Normal file
View File

File diff suppressed because one or more lines are too long

6268
tommath.src Normal file
View File

File diff suppressed because it is too large Load Diff

10683
tommath.tex Normal file
View File

File diff suppressed because it is too large Load Diff