# moritz/libtommath forked from libtom/libtommath

1 parent 2cfbb89 commit 40c00add00c419fa7ee7525ad0aea3ab6935e81f Tom St Denis committed with sjaeckel Feb 28, 2003
Showing with 289 additions and 21 deletions.
1. +98 −6 bn.c
2. +13 −0 bn.h
3. bn.pdf
4. +15 −6 bn.tex
5. +5 −0 changes.txt
6. +0 −1 demo.c
7. +1 −1 etc/makefile
8. +150 −0 etc/mersenne.c
9. +2 −2 etc/pprime.c
10. +3 −3 makefile
11. +2 −2 mtest/mtest.c
104 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)
13 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
BIN bn.pdf
Binary file not shown.
21 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 @@ \section{Introduction} \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 @@ \section{Building Against LibTomMath} #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 @@ \subsection{Modular Arithmetic} /* 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 @@ \subsubsection{mp\_n\_root(mp\_int *a, mp\_digit b, mp\_int c)} 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
 @@ -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"
 @@ -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);
 @@ -1 +1 @@ -CFLAGS += -I../ -Wall -W -O3 -fomit-frame-pointer -funroll-loops ../bn.c +CFLAGS += -I../ -Wall -W -Wshadow -ansi -O3 -fomit-frame-pointer -funroll-loops ../bn.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; +} +