diff --git a/bn.c b/bn.c index 0debe07..33d027f 100644 --- a/bn.c +++ b/bn.c @@ -2888,11 +2888,6 @@ int mp_exptmod(mp_int *G, mp_int *X, mp_int *P, mp_int *Y) * * The M table contains powers of the input base, e.g. M[x] = G^x mod P * - * This table is not made in the straight forward manner of a for loop with only - * multiplications. Since squaring is faster than multiplication we use as many - * squarings as possible. As a result about half of the steps to make the M - * table are squarings. - * * The first half of the table is not computed though accept for M[0] and M[1] */ mp_set(&M[0], 1); @@ -2914,7 +2909,6 @@ int mp_exptmod(mp_int *G, mp_int *X, mp_int *P, mp_int *Y) } } - /* create upper table */ for (x = (1<<(winsize-1))+1; x < (1 << winsize); x++) { if ((err = mp_mul(&M[x-1], &M[1], &M[x])) != MP_OKAY) { @@ -3132,6 +3126,104 @@ __T1: mp_clear(&t1); return res; } +/* computes the jacobi c = (a | n) (or Legendre if b is prime) + * HAC pp. 73 Algorithm 2.149 + */ +int mp_jacobi(mp_int *a, mp_int *n, int *c) +{ + mp_int a1, n1, e; + int s, r, res; + mp_digit residue; + + /* step 1. if a == 0, return 0 */ + if (mp_iszero(a) == 1) { + *c = 0; + return MP_OKAY; + } + + /* step 2. if a == 1, return 1 */ + if (mp_cmp_d(a, 1) == MP_EQ) { + *c = 1; + return MP_OKAY; + } + + /* default */ + s = 0; + + /* step 3. write a = a1 * 2^e */ + if ((res = mp_init_copy(&a1, a)) != MP_OKAY) { + return res; + } + + if ((res = mp_init(&n1)) != MP_OKAY) { + goto __A1; + } + + if ((res = mp_init(&e)) != MP_OKAY) { + goto __N1; + } + + while (mp_iseven(&a1) == 1) { + if ((res = mp_add_d(&e, 1, &e)) != MP_OKAY) { + goto __E; + } + + if ((res = mp_div_2(&a1, &a1)) != MP_OKAY) { + goto __E; + } + } + + /* step 4. if e is even set s=1 */ + if (mp_iseven(&e) == 1) { + s = 1; + } else { + /* else set s=1 if n = 1/7 (mod 8) or s=-1 if n = 3/5 (mod 8) */ + if ((res = mp_mod_d(n, 8, &residue)) != MP_OKAY) { + goto __E; + } + + if (residue == 1 || residue == 7) { + s = 1; + } else if (residue == 3 || residue == 5) { + s = -1; + } + } + + /* step 5. if n == 3 (mod 4) *and* a1 == 3 (mod 4) then s = -s */ + if ((res = mp_mod_d(n, 4, &residue)) != MP_OKAY) { + goto __E; + } + if (residue == 3) { + if ((res = mp_mod_d(&a1, 4, &residue)) != MP_OKAY) { + goto __E; + } + if (residue == 3) { + s = -s; + } + } + + /* if a1 == 1 we're done */ + if (mp_cmp_d(&a1, 1) == MP_EQ) { + *c = s; + } else { + /* n1 = n mod a1 */ + if ((res = mp_mod(n, &a1, &n1)) != MP_OKAY) { + goto __E; + } + if ((res = mp_jacobi(&n1, &a1, &r)) != MP_OKAY) { + goto __E; + } + *c = s * r; + } + + /* done */ + res = MP_OKAY; +__E: mp_clear(&e); +__N1: mp_clear(&n1); +__A1: mp_clear(&a1); + return res; +} + /* --> radix conversion <-- */ /* reverse an array, used for radix code */ static void reverse(unsigned char *s, int len) diff --git a/bn.h b/bn.h index 4493c65..903e7d6 100644 --- a/bn.h +++ b/bn.h @@ -21,6 +21,11 @@ #include #include +#ifdef __cplusplus +extern "C" { +#endif + + /* some default configurations. * * A "mp_digit" must be able to hold DIGIT_BIT + 1 bits @@ -239,6 +244,9 @@ int mp_n_root(mp_int *a, mp_digit b, mp_int *c); /* shortcut for square root */ #define mp_sqrt(a, b) mp_n_root(a, 2, b) +/* computes the jacobi c = (a | n) (or Legendre if b is prime) */ +int mp_jacobi(mp_int *a, mp_int *n, int *c); + /* used to setup the Barrett reduction for a given modulus b */ int mp_reduce_setup(mp_int *a, mp_int *b); @@ -280,5 +288,10 @@ int mp_radix_size(mp_int *a, int radix); #define mp_todecimal(M, S) mp_toradix((M), (S), 10) #define mp_tohex(M, S) mp_toradix((M), (S), 16) +#ifdef __cplusplus + } +#endif + + #endif diff --git a/bn.pdf b/bn.pdf index 38011e2..3bea8f2 100644 Binary files a/bn.pdf and b/bn.pdf differ diff --git a/bn.tex b/bn.tex index cdf0213..ed2a46d 100644 --- a/bn.tex +++ b/bn.tex @@ -1,7 +1,7 @@ \documentclass{article} \begin{document} -\title{LibTomMath v0.08 \\ A Free Multiple Precision Integer Library} +\title{LibTomMath v0.09 \\ A Free Multiple Precision Integer Library} \author{Tom St Denis \\ tomstdenis@iahu.ca} \maketitle \newpage @@ -23,8 +23,8 @@ LibTomMath was designed with the following goals in mind: \item Be written entirely in portable C. \end{enumerate} -All three goals have been achieved. Particularly the speed increase goal. For example, a 512-bit modular exponentiation is -four times faster\footnote{On an Athlon XP with GCC 3.2} with LibTomMath compared to MPI. +All three goals have been achieved. Particularly the speed increase goal. For example, a 512-bit modular exponentiation +is eight times faster\footnote{On an Athlon XP with GCC 3.2} with LibTomMath compared to MPI. Being compatible with MPI means that applications that already use it can be ported fairly quickly. Currently there are a few differences but there are many similarities. In fact the average MPI based application can be ported in under 15 @@ -51,9 +51,7 @@ with #include "bn.h" \end{verbatim} -Remove ``mpi.c'' from your project and replace it with ``bn.c''. Note that currently MPI has a few more functions than -LibTomMath has (e.g. no square-root code and a few others). Those are planned for future releases. In the interim work -arounds can be sought. Note that LibTomMath doesn't lack any functions required to build a cryptosystem. +Remove ``mpi.c'' from your project and replace it with ``bn.c''. \section{Programming with LibTomMath} @@ -278,6 +276,9 @@ int mp_lcm(mp_int *a, mp_int *b, mp_int *c); /* find the b'th root of a */ int mp_n_root(mp_int *a, mp_digit b, mp_int *c); +/* computes the jacobi c = (a | n) (or Legendre if b is prime) */ +int mp_jacobi(mp_int *a, mp_int *n, int *c); + /* d = a^b (mod c) */ int mp_exptmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d); \end{verbatim} @@ -444,6 +445,14 @@ requires $b$ multiplications and one division for a total work of $O(6N^2 \cdot If the input $a$ is negative and $b$ is even the function returns an error. Otherwise the function will return a root that has a sign that agrees with the sign of $a$. +\subsubsection{mp\_jacobi(mp\_int *a, mp\_int *n, int *c)} +Computes $c = \left ( {a \over n} \right )$ or the Jacobi function of $(a, n)$ and stores the result in an integer addressed +by $c$. Since the result of the Jacobi function $\left ( {a \over n} \right ) \in \lbrace -1, 0, 1 \rbrace$ it seemed +natural to store the result in a simple C style \textbf{int}. If $n$ is prime then the Jacobi function produces +the same results as the Legendre function\footnote{Source: Handbook of Applied Cryptography, pp. 73}. This means if +$n$ is prime then $\left ( {a \over n} \right )$ is equal to $1$ if $a$ is a quadratic residue modulo $n$ or $-1$ if +it is not. + \subsubsection{mp\_exptmod(mp\_int *a, mp\_int *b, mp\_int *c, mp\_int *d)} Computes $d = a^b \mbox{ (mod }c\mbox{)}$ using a sliding window $k$-ary exponentiation algorithm. For an $\alpha$-bit exponent it performs $\alpha$ squarings and at most $\lfloor \alpha/k \rfloor + k$ multiplications. The value of $k$ is diff --git a/changes.txt b/changes.txt index 87773a2..2d84db9 100644 --- a/changes.txt +++ b/changes.txt @@ -1,3 +1,8 @@ +Jan 6th, 2003 +v0.09 -- Updated the manual to reflect recent changes. :-) + -- Added Jacobi function (mp_jacobi) to supplement the number theory side of the lib + -- Added a Mersenne prime finder demo in ./etc/mersenne.c + Jan 2nd, 2003 v0.08 -- Sped up the multipliers by moving the inner loop variables into a smaller scope -- Corrected a bunch of small "warnings" diff --git a/demo.c b/demo.c index 0bf5aac..f671758 100644 --- a/demo.c +++ b/demo.c @@ -94,7 +94,6 @@ int main(void) mp_init(&d); mp_init(&e); mp_init(&f); - mp_read_radix(&a, "V//////////////////////////////////////////////////////////////////////////////////////", 64); mp_reduce_setup(&b, &a); diff --git a/etc/makefile b/etc/makefile index f38ed47..ed8f915 100644 --- a/etc/makefile +++ b/etc/makefile @@ -1 +1 @@ -CFLAGS += -I../ -Wall -W -O3 -fomit-frame-pointer -funroll-loops ../bn.c \ No newline at end of file +CFLAGS += -I../ -Wall -W -Wshadow -ansi -O3 -fomit-frame-pointer -funroll-loops ../bn.c \ No newline at end of file diff --git a/etc/mersenne.c b/etc/mersenne.c new file mode 100644 index 0000000..b6bbd51 --- /dev/null +++ b/etc/mersenne.c @@ -0,0 +1,150 @@ +/* Finds Mersenne primes using the Lucas-Lehmer test + * + * Tom St Denis, tomstdenis@iahu.ca + */ +#include +#include + +int is_mersenne(long s, int *pp) +{ + mp_int n, u, mu; + int res, k; + long ss; + + *pp = 0; + + if ((res = mp_init(&n)) != MP_OKAY) { + return res; + } + + if ((res = mp_init(&u)) != MP_OKAY) { + goto __N; + } + + if ((res = mp_init(&mu)) != MP_OKAY) { + goto __U; + } + + /* n = 2^s - 1 */ + mp_set(&n, 1); + ss = s; + while (ss--) { + if ((res = mp_mul_2(&n, &n)) != MP_OKAY) { + goto __MU; + } + } + if ((res = mp_sub_d(&n, 1, &n)) != MP_OKAY) { + goto __MU; + } + + /* setup mu */ + if ((res = mp_reduce_setup(&mu, &n)) != MP_OKAY) { + goto __MU; + } + + /* set u=4 */ + mp_set(&u, 4); + + /* for k=1 to s-2 do */ + for (k = 1; k <= s - 2; k++) { + /* u = u^2 - 2 mod n */ + if ((res = mp_sqr(&u, &u)) != MP_OKAY) { + goto __MU; + } + if ((res = mp_sub_d(&u, 2, &u)) != MP_OKAY) { + goto __MU; + } + + /* make sure u is positive */ + if (u.sign == MP_NEG) { + if ((res = mp_add(&u, &n, &u)) != MP_OKAY) { + goto __MU; + } + } + + /* reduce */ + if ((res = mp_reduce(&u, &n, &mu)) != MP_OKAY) { + goto __MU; + } + } + + /* if u == 0 then its prime */ + if (mp_iszero(&u) == 1) { + *pp = 1; + } + + res = MP_OKAY; +__MU: mp_clear(&mu); +__U: mp_clear(&u); +__N: mp_clear(&n); + return res; +} + +/* square root of a long < 65536 */ +long i_sqrt(long x) +{ + long x1, x2; + + x2 = 16; + do { + x1 = x2; + x2 = x1 - ((x1 * x1) - x)/(2*x1); + } while (x1 != x2); + + if (x1*x1 > x) { + --x1; + } + + return x1; +} + +/* is the long prime by brute force */ +int isprime(long k) +{ + long y, z; + + y = i_sqrt(k); + for (z = 2; z <= y; z++) { + if ((k % z) == 0) return 0; + } + return 1; +} + + +int main(void) +{ + int pp; + long k; + clock_t tt; + + k = 3; + + for (;;) { + /* start time */ + tt = clock(); + + /* test if 2^k - 1 is prime */ + if (is_mersenne(k, &pp) != MP_OKAY) { + printf("Whoa error\n"); + return -1; + } + + if (pp == 1) { + /* count time */ + tt = clock() - tt; + + /* display if prime */ + printf("2^%-5ld - 1 is prime, test took %ld ticks\n", k, tt); + } + + /* goto next odd exponent */ + k += 2; + + /* but make sure its prime */ + while (isprime(k) == 0) { + k += 2; + } + } + return 0; +} + diff --git a/etc/pprime.c b/etc/pprime.c index 84cf79c..fb987e3 100644 --- a/etc/pprime.c +++ b/etc/pprime.c @@ -56,7 +56,7 @@ static mp_digit prime_digit() ++y; next = (y+1)*(y+1); } - + /* loop if divisible by 3,5,7,11,13,17,19,23,29 */ if ((r % 3) == 0) { x = 0; continue; } if ((r % 5) == 0) { x = 0; continue; } @@ -138,7 +138,7 @@ int pprime(int k, mp_int *p, mp_int *q) /* now loop making the single digit */ while (mp_count_bits(&a) < k) { - printf("prime is %4d bits left\r", k - mp_count_bits(&a)); fflush(stdout); + printf("prime has %4d bits left\r", k - mp_count_bits(&a)); fflush(stdout); top: mp_set(&b, prime_digit()); diff --git a/makefile b/makefile index cbb5ac7..d448933 100644 --- a/makefile +++ b/makefile @@ -1,13 +1,13 @@ CC = gcc -CFLAGS += -Wall -W -O3 -fomit-frame-pointer -funroll-loops +CFLAGS += -Wall -W -Wshadow -ansi -O3 -fomit-frame-pointer -funroll-loops -VERSION=0.08 +VERSION=0.09 default: test test: bn.o demo.o $(CC) bn.o demo.o -o demo - cd mtest ; gcc -O3 -fomit-frame-pointer -funroll-loops mtest.c -o mtest.exe -s + cd mtest ; gcc $(CFLAGS) mtest.c -o mtest.exe -s # builds the x86 demo test86: diff --git a/mtest/mtest.c b/mtest/mtest.c index de04e2b..df5422f 100644 --- a/mtest/mtest.c +++ b/mtest/mtest.c @@ -41,7 +41,7 @@ void rand_num(mp_int *a) unsigned char buf[512]; top: - size = 1 + ((fgetc(rng)*fgetc(rng)) % 32); + size = 1 + ((fgetc(rng)*fgetc(rng)) % 96); buf[0] = (fgetc(rng)&1)?1:0; fread(buf+1, 1, size, rng); for (n = 0; n < size; n++) { @@ -57,7 +57,7 @@ void rand_num2(mp_int *a) unsigned char buf[512]; top: - size = 1 + ((fgetc(rng)*fgetc(rng)) % 32); + size = 1 + ((fgetc(rng)*fgetc(rng)) % 96); buf[0] = (fgetc(rng)&1)?1:0; fread(buf+1, 1, size, rng); for (n = 0; n < size; n++) {