Browse files

added libtommath-0.02

  • Loading branch information...
1 parent 390fa39 commit 8c97c9e1814ad3b7c5e99a7c951521d6ea079834 Tom St Denis committed with sjaeckel Feb 28, 2003
Showing with 4,689 additions and 47 deletions.
  1. +2 −2 b.bat
  2. +71 −14 bn.c
  3. +17 −0 bn.h
  4. BIN bn.pdf
  5. +34 −24 bn.tex
  6. +10 −5 changes.txt
  7. +6 −1 demo.c
  8. +1 −1 makefile
  9. +20 −0 mtest/logtab.h
  10. +86 −0 mtest/mpi-config.h
  11. +16 −0 mtest/mpi-types.h
  12. +3,981 −0 mtest/mpi.c
  13. +227 −0 mtest/mpi.h
  14. +218 −0 mtest/mtest.c
View
4 b.bat
@@ -1,3 +1,3 @@
nasm -f coff timer.asm
-gcc -Wall -W -O3 -fomit-frame-pointer -funroll-loops -DTIMER demo.c bn.c timer.o -o demo
-gcc -I./mtest/ -DU_MPI -Wall -W -O3 -fomit-frame-pointer -funroll-loops -DTIMER demo.c mtest/mpi.c timer.o -o mpidemo
+gcc -Wall -W -O3 -fomit-frame-pointer -funroll-loops -DTIMER_X86 demo.c bn.c timer.o -o demo
+gcc -I./mtest/ -DU_MPI -Wall -W -O3 -fomit-frame-pointer -funroll-loops -DTIMER_X86 demo.c mtest/mpi.c timer.o -o mpidemo
View
85 bn.c
@@ -168,8 +168,31 @@ int mp_init_copy(mp_int *a, mp_int *b)
return mp_copy(b, a);
}
+/* b = |a| */
+int mp_abs(mp_int *a, mp_int *b)
+{
+ int res;
+ if ((res = mp_copy(a, b)) != MP_OKAY) {
+ return res;
+ }
+ b->sign = MP_ZPOS;
+ return MP_OKAY;
+}
+
+/* b = -a */
+int mp_neg(mp_int *a, mp_int *b)
+{
+ int res;
+ if ((res = mp_copy(a, b)) != MP_OKAY) {
+ return res;
+ }
+ b->sign = (a->sign == MP_ZPOS) ? MP_NEG : MP_ZPOS;
+ return MP_OKAY;
+}
+
+
/* compare maginitude of two ints (unsigned) */
-static int s_mp_cmp(mp_int *a, mp_int *b)
+int mp_cmp_mag(mp_int *a, mp_int *b)
{
int n;
@@ -200,7 +223,7 @@ int mp_cmp(mp_int *a, mp_int *b)
} else if (a->sign == MP_ZPOS && b->sign == MP_NEG) {
return MP_GT;
}
- return s_mp_cmp(a, b);
+ return mp_cmp_mag(a, b);
}
/* compare a digit */
@@ -500,7 +523,7 @@ static int fast_s_mp_mul_digs(mp_int *a, mp_int *b, mp_int *c, int digs)
{
mp_int t;
int res, pa, pb, ix, iy;
- mp_word W[512];
+ mp_word W[512], *_W;
mp_digit tmpx, *tmpt, *tmpy;
// printf("\nHOLA\n");
@@ -519,8 +542,9 @@ static int fast_s_mp_mul_digs(mp_int *a, mp_int *b, mp_int *c, int digs)
tmpx = a->dp[ix];
tmpt = &(t.dp[ix]);
tmpy = b->dp;
+ _W = &(W[ix]);
for (iy = 0; iy < pb; iy++) {
- W[iy+ix] += ((mp_word)tmpx) * ((mp_word)*tmpy++);
+ *_W++ += ((mp_word)tmpx) * ((mp_word)*tmpy++);
}
}
@@ -590,7 +614,7 @@ static int fast_s_mp_mul_high_digs(mp_int *a, mp_int *b, mp_int *c, int digs)
{
mp_int t;
int res, pa, pb, ix, iy;
- mp_word W[512];
+ mp_word W[512], *_W;
mp_digit tmpx, *tmpt, *tmpy;
if ((res = mp_init_size(&t, a->used + b->used + 1)) != MP_OKAY) {
@@ -605,8 +629,9 @@ static int fast_s_mp_mul_high_digs(mp_int *a, mp_int *b, mp_int *c, int digs)
tmpx = a->dp[ix];
tmpt = &(t.dp[digs]);
tmpy = b->dp + (digs - ix);
+ _W = &(W[digs]);
for (iy = digs - ix; iy < pb; iy++) {
- W[ix+iy] += ((mp_word)tmpx) * ((mp_word)*tmpy++);
+ *_W++ += ((mp_word)tmpx) * ((mp_word)*tmpy++);
}
}
@@ -678,7 +703,7 @@ static int fast_s_mp_sqr(mp_int *a, mp_int *b)
{
mp_int t;
int res, ix, iy, pa;
- mp_word r, W[512];
+ mp_word W[512], *_W;
mp_digit tmpx, *tmpy;
pa = a->used;
@@ -694,9 +719,9 @@ static int fast_s_mp_sqr(mp_int *a, mp_int *b)
W[ix+ix] += ((mp_word)a->dp[ix]) * ((mp_word)a->dp[ix]);
tmpx = a->dp[ix];
tmpy = &(a->dp[ix+1]);
+ _W = &(W[ix+ix+1]);
for (iy = ix + 1; iy < pa; iy++) {
- r = ((mp_word)tmpx) * ((mp_word)*tmpy++);
- W[ix+iy] += r + r;
+ *_W++ += ((mp_word)tmpx) * ((mp_word)*tmpy++) << ((mp_word)1);
}
}
for (ix = 1; ix < (pa+pa+1); ix++) {
@@ -781,7 +806,7 @@ int mp_add(mp_int *a, mp_int *b, mp_int *c)
c->sign = MP_ZPOS;
} else if (sa == MP_ZPOS && sb == MP_NEG) {
/* a + -b == a - b, but if b>a then we do it as -(b-a) */
- if (s_mp_cmp(a, b) == MP_LT) {
+ if (mp_cmp_mag(a, b) == MP_LT) {
res = s_mp_sub(b, a, c);
c->sign = MP_NEG;
} else {
@@ -790,7 +815,7 @@ int mp_add(mp_int *a, mp_int *b, mp_int *c)
}
} else if (sa == MP_NEG && sb == MP_ZPOS) {
/* -a + b == b - a, but if a>b then we do it as -(a-b) */
- if (s_mp_cmp(a, b) == MP_GT) {
+ if (mp_cmp_mag(a, b) == MP_GT) {
res = s_mp_sub(a, b, c);
c->sign = MP_NEG;
} else {
@@ -816,7 +841,7 @@ int mp_sub(mp_int *a, mp_int *b, mp_int *c)
/* handle four cases */
if (sa == MP_ZPOS && sb == MP_ZPOS) {
/* both positive, a - b, but if b>a then we do -(b - a) */
- if (s_mp_cmp(a, b) == MP_LT) {
+ if (mp_cmp_mag(a, b) == MP_LT) {
/* b>a */
res = s_mp_sub(b, a, c);
c->sign = MP_NEG;
@@ -834,7 +859,7 @@ int mp_sub(mp_int *a, mp_int *b, mp_int *c)
c->sign = MP_NEG;
} else {
/* -a - -b == b - a, but if a>b == -(a - b) */
- if (s_mp_cmp(a, b) == MP_GT) {
+ if (mp_cmp_mag(a, b) == MP_GT) {
res = s_mp_sub(a, b, c);
c->sign = MP_NEG;
} else {
@@ -1023,7 +1048,7 @@ int mp_div(mp_int *a, mp_int *b, mp_int *c, mp_int *d)
int res, n, t, i, norm, neg;
/* if a < b then q=0, r = a */
- if (s_mp_cmp(a, b) == MP_LT) {
+ if (mp_cmp_mag(a, b) == MP_LT) {
if (d != NULL) {
res = mp_copy(a, d);
d->sign = a->sign;
@@ -1981,3 +2006,35 @@ int mp_toradix(mp_int *a, unsigned char *str, int radix)
return MP_OKAY;
}
+/* returns size of ASCII reprensentation */
+int mp_radix_size(mp_int *a, int radix)
+{
+ int res, digs;
+ mp_int t;
+ mp_digit d;
+
+ digs = 0;
+
+ if (radix < 2 || radix > 64) {
+ return 0;
+ }
+
+ if ((res = mp_init_copy(&t, a)) != MP_OKAY) {
+ return 0;
+ }
+
+ if (t.sign == MP_NEG) {
+ ++digs;
+ t.sign = MP_ZPOS;
+ }
+
+ while (mp_iszero(&t) == 0) {
+ if ((res = mp_div_d(&t, (mp_digit)radix, &t, &d)) != MP_OKAY) {
+ return 0;
+ }
+ ++digs;
+ }
+ mp_clear(&t);
+ return digs + 1;
+}
+
View
17 bn.h
@@ -128,9 +128,18 @@ int mp_mod_2d(mp_int *a, int b, mp_int *c);
/* ---> Basic arithmetic <--- */
+/* b = -a */
+int mp_neg(mp_int *a, mp_int *b);
+
+/* b = |a| */
+int mp_abs(mp_int *a, mp_int *b);
+
/* compare a to b */
int mp_cmp(mp_int *a, mp_int *b);
+/* compare |a| to |b| */
+int mp_cmp_mag(mp_int *a, mp_int *b);
+
/* c = a + b */
int mp_add(mp_int *a, mp_int *b, mp_int *c);
@@ -215,6 +224,14 @@ int mp_to_signed_bin(mp_int *a, unsigned char *b);
int mp_read_radix(mp_int *a, unsigned char *str, int radix);
int mp_toradix(mp_int *a, unsigned char *str, int radix);
+int mp_radix_size(mp_int *a, int radix);
+
+#define mp_read_raw(mp, str, len) mp_read_signed_bin((mp), (str), (len))
+#define mp_raw_size(mp) mp_signed_bin_size(mp)
+#define mp_toraw(mp, str) mp_to_signed_bin((mp), (str))
+#define mp_read_mag(mp, str, len) mp_read_unsigned_bin((mp), (str), (len))
+#define mp_mag_size(mp) mp_unsigned_bin_size(mp)
+#define mp_tomag(mp, str) mp_to_unsigned_bin((mp), (str))
#define mp_tobinary(M, S) mp_toradix((M), (S), 2)
#define mp_tooctal(M, S) mp_toradix((M), (S), 8)
View
BIN bn.pdf
Binary file not shown.
View
58 bn.tex
@@ -1,7 +1,7 @@
\documentclass{article}
\begin{document}
-\title{LibTomMath v0.01 \\ A Free Multiple Precision Integer Library}
+\title{LibTomMath v0.02 \\ A Free Multiple Precision Integer Library}
\author{Tom St Denis \\ tomstdenis@iahu.ca}
\maketitle
\newpage
@@ -192,9 +192,18 @@ \subsection{Basic Arithmetic}
Next are the class of functions which provide basic arithmetic.
\begin{verbatim}
+/* b = -a */
+int mp_neg(mp_int *a, mp_int *b);
+
+/* b = |a| */
+int mp_abs(mp_int *a, mp_int *b);
+
/* compare a to b */
int mp_cmp(mp_int *a, mp_int *b);
+/* compare |a| to |b| */
+int mp_cmp_mag(mp_int *a, mp_int *b);
+
/* c = a + b */
int mp_add(mp_int *a, mp_int *b, mp_int *c);
@@ -298,6 +307,7 @@ \subsection{Radix Conversions}
int mp_read_radix(mp_int *a, unsigned char *str, int radix);
int mp_toradix(mp_int *a, unsigned char *str, int radix);
+int mp_radix_size(mp_int *a, int radix);
\end{verbatim}
The integers are stored in big endian format as most libraries (and MPI) expect. The \textbf{mp\_read\_radix} and
@@ -318,27 +328,27 @@ \subsection{Observed Timings}
\begin{tabular}{c|c|c|c}
\hline \textbf{Operation} & \textbf{Size (bits)} & \textbf{Time with MPI (cycles)} & \textbf{Time with LibTomMath (cycles)} \\
\hline
-Multiply & 128 & 1,394 & 915 \\
-Multiply & 256 & 2,559 & 1,893 \\
-Multiply & 512 & 7,919 & 3,770 \\
-Multiply & 1024 & 28,460 & 9,970 \\
-Multiply & 2048 & 109,637 & 32,264 \\
-Multiply & 4096 & 467,226 & 129,645 \\
+Multiply & 128 & 1,394 & 893 \\
+Multiply & 256 & 2,559 & 1,744 \\
+Multiply & 512 & 7,919 & 4,484 \\
+Multiply & 1024 & 28,460 & 9,326, \\
+Multiply & 2048 & 109,637 & 30,140 \\
+Multiply & 4096 & 467,226 & 122,290 \\
\hline
-Square & 128 & 1,288 & 1,147 \\
-Square & 256 & 1,705 & 2,129 \\
-Square & 512 & 5,365 & 3,755 \\
-Square & 1024 & 18,836 & 9,267 \\
-Square & 2048 & 72,334 & 28,387 \\
-Square & 4096 & 306,252 & 112,391 \\
+Square & 128 & 1,288 & 1,172 \\
+Square & 256 & 1,705 & 2,162 \\
+Square & 512 & 5,365 & 3,723 \\
+Square & 1024 & 18,836 & 9,063 \\
+Square & 2048 & 72,334 & 27,489 \\
+Square & 4096 & 306,252 & 110,372 \\
\hline
-Exptmod & 512 & 30,497,732 & 7,222,872 \\
-Exptmod & 768 & 98,943,020 & 16,474,567 \\
-Exptmod & 1024 & 221,123,749 & 30,070,883 \\
-Exptmod & 2048 & 1,694,796,907 & 154,697,320 \\
-Exptmod & 2560 & 3,262,360,107 & 318,998,183 \\
-Exptmod & 3072 & 5,647,243,373 & 494,313,122 \\
-Exptmod & 4096 & 13,345,194,048 & 1,036,254,558
+Exptmod & 512 & 30,497,732 & 6,898,504 \\
+Exptmod & 768 & 98,943,020 & 15,510,779 \\
+Exptmod & 1024 & 221,123,749 & 27,962,904 \\
+Exptmod & 2048 & 1,694,796,907 & 146,631,975 \\
+Exptmod & 2560 & 3,262,360,107 & 305,530,060 \\
+Exptmod & 3072 & 5,647,243,373 & 472,572,762 \\
+Exptmod & 4096 & 13,345,194,048 & 984,415,240
\end{tabular}
\end{center}
@@ -367,7 +377,7 @@ \subsection{Multiplication Algorithms}
will benefit from the algorithm.
MPI only implements the slower baseline multiplier where carries are dealt with in the inner loop. As a result even at
-smaller numbers (below the Karatsuba cutoff) the LibTomCrypt multipliers are faster.
+smaller numbers (below the Karatsuba cutoff) the LibTomMath multipliers are faster.
\subsection{Squaring Algorithms}
@@ -393,16 +403,16 @@ \subsection{Exponentiation Algorithms}
of $g(3)$, four squarings and and a multiplication by $g(1)$. In total there are 8 squarings and 3 multiplications.
MPI uses a binary square-multiply method. For the same exponent $e$ it would have had 8 squarings and 5 multiplications.
-There is a precomputation phase for the method LibTomCrypt uses but it generally cuts down considerably on the number
+There is a precomputation phase for the method LibTomMath uses but it generally cuts down considerably on the number
of multiplications. Consider a 512-bit exponent. The worst case for the LibTomMath method results in 512 squarings and
124 multiplications. The MPI method would have 512 squarings and 512 multiplications. Randomly every $2k$ bits another
multiplication is saved via the sliding-window technique on top of the savings the $k$-ary method provides.
Both LibTomMath and MPI use Barrett reduction instead of division to reduce the numbers modulo the modulus given.
However, LibTomMath can take advantage of the fact that the multiplications required within the Barrett reduction
-do not to give full precision. As a result the reduction step is much faster and just as accurate. The LibTomMath code
+do not have to give full precision. As a result the reduction step is much faster and just as accurate. The LibTomMath code
will automatically determine at run-time (e.g. when its called) whether the faster multiplier can be used. The
-faster multipliers have also been optimized into the two variants (baseline and fast baseline).
+faster multipliers have also been optimized into the two variants (baseline and comba baseline).
As a result of all these changes exponentiation in LibTomMath is much faster than compared to MPI.
View
15 changes.txt
@@ -1,6 +1,11 @@
+Dec 26th, 2002
+v0.02 -- Fixed a few "slips" in the manual. This is "LibTomMath" afterall :-)
+ -- Added mp_cmp_mag, mp_neg, mp_abs and mp_radix_size that were missing.
+ -- Sped up the fast [comba] multipliers more [yahoo!]
+
Dec 25th,2002
-v0.01 -- Initial release. Gimme a break.
- -- Todo list,
- add details to manual [e.g. algorithms]
- more comments in code
- example programs
+v0.01 -- Initial release. Gimme a break.
+ -- Todo list,
+ add details to manual [e.g. algorithms]
+ more comments in code
+ example programs
View
7 demo.c
@@ -11,9 +11,14 @@
#include "bn.h"
#endif
-#ifdef TIMER
+#ifdef TIMER_X86
+#define TIMER
extern unsigned long long rdtsc(void);
extern void reset(void);
+#else
+unsigned long long _tt;
+void reset(void) { _tt = clock(); }
+unsigned long long rdtsc(void) { return clock() - _tt; }
#endif
static void draw(mp_int *a)
View
2 makefile
@@ -1,7 +1,7 @@
CC = gcc
CFLAGS += -Wall -W -O3 -funroll-loops
-VERSION=0.01
+VERSION=0.02
default: test
View
20 mtest/logtab.h
@@ -0,0 +1,20 @@
+const float s_logv_2[] = {
+ 0.000000000, 0.000000000, 1.000000000, 0.630929754, /* 0 1 2 3 */
+ 0.500000000, 0.430676558, 0.386852807, 0.356207187, /* 4 5 6 7 */
+ 0.333333333, 0.315464877, 0.301029996, 0.289064826, /* 8 9 10 11 */
+ 0.278942946, 0.270238154, 0.262649535, 0.255958025, /* 12 13 14 15 */
+ 0.250000000, 0.244650542, 0.239812467, 0.235408913, /* 16 17 18 19 */
+ 0.231378213, 0.227670249, 0.224243824, 0.221064729, /* 20 21 22 23 */
+ 0.218104292, 0.215338279, 0.212746054, 0.210309918, /* 24 25 26 27 */
+ 0.208014598, 0.205846832, 0.203795047, 0.201849087, /* 28 29 30 31 */
+ 0.200000000, 0.198239863, 0.196561632, 0.194959022, /* 32 33 34 35 */
+ 0.193426404, 0.191958720, 0.190551412, 0.189200360, /* 36 37 38 39 */
+ 0.187901825, 0.186652411, 0.185449023, 0.184288833, /* 40 41 42 43 */
+ 0.183169251, 0.182087900, 0.181042597, 0.180031327, /* 44 45 46 47 */
+ 0.179052232, 0.178103594, 0.177183820, 0.176291434, /* 48 49 50 51 */
+ 0.175425064, 0.174583430, 0.173765343, 0.172969690, /* 52 53 54 55 */
+ 0.172195434, 0.171441601, 0.170707280, 0.169991616, /* 56 57 58 59 */
+ 0.169293808, 0.168613099, 0.167948779, 0.167300179, /* 60 61 62 63 */
+ 0.166666667
+};
+
View
86 mtest/mpi-config.h
@@ -0,0 +1,86 @@
+/* Default configuration for MPI library */
+/* $ID$ */
+
+#ifndef MPI_CONFIG_H_
+#define MPI_CONFIG_H_
+
+/*
+ For boolean options,
+ 0 = no
+ 1 = yes
+
+ Other options are documented individually.
+
+ */
+
+#ifndef MP_IOFUNC
+#define MP_IOFUNC 0 /* include mp_print() ? */
+#endif
+
+#ifndef MP_MODARITH
+#define MP_MODARITH 1 /* include modular arithmetic ? */
+#endif
+
+#ifndef MP_NUMTH
+#define MP_NUMTH 1 /* include number theoretic functions? */
+#endif
+
+#ifndef MP_LOGTAB
+#define MP_LOGTAB 1 /* use table of logs instead of log()? */
+#endif
+
+#ifndef MP_MEMSET
+#define MP_MEMSET 1 /* use memset() to zero buffers? */
+#endif
+
+#ifndef MP_MEMCPY
+#define MP_MEMCPY 1 /* use memcpy() to copy buffers? */
+#endif
+
+#ifndef MP_CRYPTO
+#define MP_CRYPTO 1 /* erase memory on free? */
+#endif
+
+#ifndef MP_ARGCHK
+/*
+ 0 = no parameter checks
+ 1 = runtime checks, continue execution and return an error to caller
+ 2 = assertions; dump core on parameter errors
+ */
+#define MP_ARGCHK 2 /* how to check input arguments */
+#endif
+
+#ifndef MP_DEBUG
+#define MP_DEBUG 0 /* print diagnostic output? */
+#endif
+
+#ifndef MP_DEFPREC
+#define MP_DEFPREC 64 /* default precision, in digits */
+#endif
+
+#ifndef MP_MACRO
+#define MP_MACRO 1 /* use macros for frequent calls? */
+#endif
+
+#ifndef MP_SQUARE
+#define MP_SQUARE 1 /* use separate squaring code? */
+#endif
+
+#ifndef MP_PTAB_SIZE
+/*
+ When building mpprime.c, we build in a table of small prime
+ values to use for primality testing. The more you include,
+ the more space they take up. See primes.c for the possible
+ values (currently 16, 32, 64, 128, 256, and 6542)
+ */
+#define MP_PTAB_SIZE 128 /* how many built-in primes? */
+#endif
+
+#ifndef MP_COMPAT_MACROS
+#define MP_COMPAT_MACROS 1 /* define compatibility macros? */
+#endif
+
+#endif /* ifndef MPI_CONFIG_H_ */
+
+
+/* crc==3287762869, version==2, Sat Feb 02 06:43:53 2002 */
View
16 mtest/mpi-types.h
@@ -0,0 +1,16 @@
+/* Type definitions generated by 'types.pl' */
+typedef char mp_sign;
+typedef unsigned short mp_digit; /* 2 byte type */
+typedef unsigned int mp_word; /* 4 byte type */
+typedef unsigned int mp_size;
+typedef int mp_err;
+
+#define MP_DIGIT_BIT (CHAR_BIT*sizeof(mp_digit))
+#define MP_DIGIT_MAX USHRT_MAX
+#define MP_WORD_BIT (CHAR_BIT*sizeof(mp_word))
+#define MP_WORD_MAX UINT_MAX
+
+#define MP_DIGIT_SIZE 2
+#define DIGIT_FMT "%04X"
+#define RADIX (MP_DIGIT_MAX+1)
+
View
3,981 mtest/mpi.c
3,981 additions, 0 deletions not shown because the diff is too large. Please use a local Git client to view these changes.
View
227 mtest/mpi.h
@@ -0,0 +1,227 @@
+/*
+ mpi.h
+
+ by Michael J. Fromberger <sting@linguist.dartmouth.edu>
+ Copyright (C) 1998 Michael J. Fromberger, All Rights Reserved
+
+ Arbitrary precision integer arithmetic library
+
+ $ID$
+ */
+
+#ifndef _H_MPI_
+#define _H_MPI_
+
+#include "mpi-config.h"
+
+#define MP_LT -1
+#define MP_EQ 0
+#define MP_GT 1
+
+#if MP_DEBUG
+#undef MP_IOFUNC
+#define MP_IOFUNC 1
+#endif
+
+#if MP_IOFUNC
+#include <stdio.h>
+#include <ctype.h>
+#endif
+
+#include <limits.h>
+
+#define MP_NEG 1
+#define MP_ZPOS 0
+
+/* Included for compatibility... */
+#define NEG MP_NEG
+#define ZPOS MP_ZPOS
+
+#define MP_OKAY 0 /* no error, all is well */
+#define MP_YES 0 /* yes (boolean result) */
+#define MP_NO -1 /* no (boolean result) */
+#define MP_MEM -2 /* out of memory */
+#define MP_RANGE -3 /* argument out of range */
+#define MP_BADARG -4 /* invalid parameter */
+#define MP_UNDEF -5 /* answer is undefined */
+#define MP_LAST_CODE MP_UNDEF
+
+#include "mpi-types.h"
+
+/* Included for compatibility... */
+#define DIGIT_BIT MP_DIGIT_BIT
+#define DIGIT_MAX MP_DIGIT_MAX
+
+/* Macros for accessing the mp_int internals */
+#define SIGN(MP) ((MP)->sign)
+#define USED(MP) ((MP)->used)
+#define ALLOC(MP) ((MP)->alloc)
+#define DIGITS(MP) ((MP)->dp)
+#define DIGIT(MP,N) (MP)->dp[(N)]
+
+#if MP_ARGCHK == 1
+#define ARGCHK(X,Y) {if(!(X)){return (Y);}}
+#elif MP_ARGCHK == 2
+#include <assert.h>
+#define ARGCHK(X,Y) assert(X)
+#else
+#define ARGCHK(X,Y) /* */
+#endif
+
+/* This defines the maximum I/O base (minimum is 2) */
+#define MAX_RADIX 64
+
+typedef struct {
+ mp_sign sign; /* sign of this quantity */
+ mp_size alloc; /* how many digits allocated */
+ mp_size used; /* how many digits used */
+ mp_digit *dp; /* the digits themselves */
+} mp_int;
+
+/*------------------------------------------------------------------------*/
+/* Default precision */
+
+unsigned int mp_get_prec(void);
+void mp_set_prec(unsigned int prec);
+
+/*------------------------------------------------------------------------*/
+/* Memory management */
+
+mp_err mp_init(mp_int *mp);
+mp_err mp_init_array(mp_int mp[], int count);
+mp_err mp_init_size(mp_int *mp, mp_size prec);
+mp_err mp_init_copy(mp_int *mp, mp_int *from);
+mp_err mp_copy(mp_int *from, mp_int *to);
+void mp_exch(mp_int *mp1, mp_int *mp2);
+void mp_clear(mp_int *mp);
+void mp_clear_array(mp_int mp[], int count);
+void mp_zero(mp_int *mp);
+void mp_set(mp_int *mp, mp_digit d);
+mp_err mp_set_int(mp_int *mp, long z);
+mp_err mp_shrink(mp_int *a);
+
+
+/*------------------------------------------------------------------------*/
+/* Single digit arithmetic */
+
+mp_err mp_add_d(mp_int *a, mp_digit d, mp_int *b);
+mp_err mp_sub_d(mp_int *a, mp_digit d, mp_int *b);
+mp_err mp_mul_d(mp_int *a, mp_digit d, mp_int *b);
+mp_err mp_mul_2(mp_int *a, mp_int *c);
+mp_err mp_div_d(mp_int *a, mp_digit d, mp_int *q, mp_digit *r);
+mp_err mp_div_2(mp_int *a, mp_int *c);
+mp_err mp_expt_d(mp_int *a, mp_digit d, mp_int *c);
+
+/*------------------------------------------------------------------------*/
+/* Sign manipulations */
+
+mp_err mp_abs(mp_int *a, mp_int *b);
+mp_err mp_neg(mp_int *a, mp_int *b);
+
+/*------------------------------------------------------------------------*/
+/* Full arithmetic */
+
+mp_err mp_add(mp_int *a, mp_int *b, mp_int *c);
+mp_err mp_sub(mp_int *a, mp_int *b, mp_int *c);
+mp_err mp_mul(mp_int *a, mp_int *b, mp_int *c);
+mp_err mp_mul_2d(mp_int *a, mp_digit d, mp_int *c);
+#if MP_SQUARE
+mp_err mp_sqr(mp_int *a, mp_int *b);
+#else
+#define mp_sqr(a, b) mp_mul(a, a, b)
+#endif
+mp_err mp_div(mp_int *a, mp_int *b, mp_int *q, mp_int *r);
+mp_err mp_div_2d(mp_int *a, mp_digit d, mp_int *q, mp_int *r);
+mp_err mp_expt(mp_int *a, mp_int *b, mp_int *c);
+mp_err mp_2expt(mp_int *a, mp_digit k);
+mp_err mp_sqrt(mp_int *a, mp_int *b);
+
+/*------------------------------------------------------------------------*/
+/* Modular arithmetic */
+
+#if MP_MODARITH
+mp_err mp_mod(mp_int *a, mp_int *m, mp_int *c);
+mp_err mp_mod_d(mp_int *a, mp_digit d, mp_digit *c);
+mp_err mp_addmod(mp_int *a, mp_int *b, mp_int *m, mp_int *c);
+mp_err mp_submod(mp_int *a, mp_int *b, mp_int *m, mp_int *c);
+mp_err mp_mulmod(mp_int *a, mp_int *b, mp_int *m, mp_int *c);
+#if MP_SQUARE
+mp_err mp_sqrmod(mp_int *a, mp_int *m, mp_int *c);
+#else
+#define mp_sqrmod(a, m, c) mp_mulmod(a, a, m, c)
+#endif
+mp_err mp_exptmod(mp_int *a, mp_int *b, mp_int *m, mp_int *c);
+mp_err mp_exptmod_d(mp_int *a, mp_digit d, mp_int *m, mp_int *c);
+#endif /* MP_MODARITH */
+
+/*------------------------------------------------------------------------*/
+/* Comparisons */
+
+int mp_cmp_z(mp_int *a);
+int mp_cmp_d(mp_int *a, mp_digit d);
+int mp_cmp(mp_int *a, mp_int *b);
+int mp_cmp_mag(mp_int *a, mp_int *b);
+int mp_cmp_int(mp_int *a, long z);
+int mp_isodd(mp_int *a);
+int mp_iseven(mp_int *a);
+
+/*------------------------------------------------------------------------*/
+/* Number theoretic */
+
+#if MP_NUMTH
+mp_err mp_gcd(mp_int *a, mp_int *b, mp_int *c);
+mp_err mp_lcm(mp_int *a, mp_int *b, mp_int *c);
+mp_err mp_xgcd(mp_int *a, mp_int *b, mp_int *g, mp_int *x, mp_int *y);
+mp_err mp_invmod(mp_int *a, mp_int *m, mp_int *c);
+#endif /* end MP_NUMTH */
+
+/*------------------------------------------------------------------------*/
+/* Input and output */
+
+#if MP_IOFUNC
+void mp_print(mp_int *mp, FILE *ofp);
+#endif /* end MP_IOFUNC */
+
+/*------------------------------------------------------------------------*/
+/* Base conversion */
+
+#define BITS 1
+#define BYTES CHAR_BIT
+
+mp_err mp_read_signed_bin(mp_int *mp, unsigned char *str, int len);
+int mp_signed_bin_size(mp_int *mp);
+mp_err mp_to_signed_bin(mp_int *mp, unsigned char *str);
+
+mp_err mp_read_unsigned_bin(mp_int *mp, unsigned char *str, int len);
+int mp_unsigned_bin_size(mp_int *mp);
+mp_err mp_to_unsigned_bin(mp_int *mp, unsigned char *str);
+
+int mp_count_bits(mp_int *mp);
+
+#if MP_COMPAT_MACROS
+#define mp_read_raw(mp, str, len) mp_read_signed_bin((mp), (str), (len))
+#define mp_raw_size(mp) mp_signed_bin_size(mp)
+#define mp_toraw(mp, str) mp_to_signed_bin((mp), (str))
+#define mp_read_mag(mp, str, len) mp_read_unsigned_bin((mp), (str), (len))
+#define mp_mag_size(mp) mp_unsigned_bin_size(mp)
+#define mp_tomag(mp, str) mp_to_unsigned_bin((mp), (str))
+#endif
+
+mp_err mp_read_radix(mp_int *mp, unsigned char *str, int radix);
+int mp_radix_size(mp_int *mp, int radix);
+int mp_value_radix_size(int num, int qty, int radix);
+mp_err mp_toradix(mp_int *mp, unsigned char *str, int radix);
+
+int mp_char2value(char ch, int r);
+
+#define mp_tobinary(M, S) mp_toradix((M), (S), 2)
+#define mp_tooctal(M, S) mp_toradix((M), (S), 8)
+#define mp_todecimal(M, S) mp_toradix((M), (S), 10)
+#define mp_tohex(M, S) mp_toradix((M), (S), 16)
+
+/*------------------------------------------------------------------------*/
+/* Error strings */
+
+const char *mp_strerror(mp_err ec);
+
+#endif /* end _H_MPI_ */
View
218 mtest/mtest.c
@@ -0,0 +1,218 @@
+/* makes a bignum test harness with NUM tests per operation
+ *
+ * the output is made in the following format [one parameter per line]
+
+operation
+operand1
+operand2
+[... operandN]
+result1
+result2
+[... resultN]
+
+So for example "a * b mod n" would be
+
+mulmod
+a
+b
+n
+a*b mod n
+
+e.g. if a=3, b=4 n=11 then
+
+mulmod
+3
+4
+11
+1
+
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include "mpi.c"
+
+FILE *rng;
+
+void rand_num(mp_int *a)
+{
+ int n, size;
+ unsigned char buf[512];
+
+top:
+ size = 1 + (fgetc(rng) % 96);
+ buf[0] = (fgetc(rng)&1)?1:0;
+ fread(buf+1, 1, size, rng);
+ for (n = 0; n < size; n++) {
+ if (buf[n+1]) break;
+ }
+ if (n == size) goto top;
+ mp_read_raw(a, buf, 1+size);
+}
+
+void rand_num2(mp_int *a)
+{
+ int n, size;
+ unsigned char buf[512];
+
+top:
+ size = 1 + (fgetc(rng) % 128);
+ buf[0] = (fgetc(rng)&1)?1:0;
+ fread(buf+1, 1, size, rng);
+ for (n = 0; n < size; n++) {
+ if (buf[n+1]) break;
+ }
+ if (n == size) goto top;
+ mp_read_raw(a, buf, 1+size);
+}
+
+int main(void)
+{
+ int n;
+ mp_int a, b, c, d, e;
+ char buf[4096];
+
+ mp_init(&a);
+ mp_init(&b);
+ mp_init(&c);
+ mp_init(&d);
+ mp_init(&e);
+
+ rng = fopen("/dev/urandom", "rb");
+
+ for (;;) {
+ n = fgetc(rng) % 10;
+
+ if (n == 0) {
+ /* add tests */
+ rand_num(&a);
+ rand_num(&b);
+ mp_add(&a, &b, &c);
+ printf("add\n");
+ mp_todecimal(&a, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&b, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&c, buf);
+ printf("%s\n", buf);
+ } else if (n == 1) {
+ /* sub tests */
+ rand_num(&a);
+ rand_num(&b);
+ mp_sub(&a, &b, &c);
+ printf("sub\n");
+ mp_todecimal(&a, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&b, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&c, buf);
+ printf("%s\n", buf);
+ } else if (n == 2) {
+ /* mul tests */
+ rand_num(&a);
+ rand_num(&b);
+ mp_mul(&a, &b, &c);
+ printf("mul\n");
+ mp_todecimal(&a, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&b, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&c, buf);
+ printf("%s\n", buf);
+ } else if (n == 3) {
+ /* div tests */
+ rand_num(&a);
+ rand_num(&b);
+ mp_div(&a, &b, &c, &d);
+ printf("div\n");
+ mp_todecimal(&a, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&b, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&c, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&d, buf);
+ printf("%s\n", buf);
+ } else if (n == 4) {
+ /* sqr tests */
+ rand_num(&a);
+ mp_sqr(&a, &b);
+ printf("sqr\n");
+ mp_todecimal(&a, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&b, buf);
+ printf("%s\n", buf);
+ } else if (n == 5) {
+ /* mul_2d test */
+ rand_num(&a);
+ mp_copy(&a, &b);
+ n = fgetc(rng) & 63;
+ mp_mul_2d(&b, n, &b);
+ mp_todecimal(&a, buf);
+ printf("mul2d\n");
+ printf("%s\n", buf);
+ printf("%d\n", n);
+ mp_todecimal(&b, buf);
+ printf("%s\n", buf);
+ } else if (n == 6) {
+ /* div_2d test */
+ rand_num(&a);
+ mp_copy(&a, &b);
+ n = fgetc(rng) & 63;
+ mp_div_2d(&b, n, &b, NULL);
+ mp_todecimal(&a, buf);
+ printf("div2d\n");
+ printf("%s\n", buf);
+ printf("%d\n", n);
+ mp_todecimal(&b, buf);
+ printf("%s\n", buf);
+ } else if (n == 7) {
+ /* gcd test */
+ rand_num(&a);
+ rand_num(&b);
+ a.sign = MP_ZPOS;
+ b.sign = MP_ZPOS;
+ mp_gcd(&a, &b, &c);
+ printf("gcd\n");
+ mp_todecimal(&a, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&b, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&c, buf);
+ printf("%s\n", buf);
+ } else if (n == 8) {
+ /* lcm test */
+ rand_num(&a);
+ rand_num(&b);
+ a.sign = MP_ZPOS;
+ b.sign = MP_ZPOS;
+ mp_lcm(&a, &b, &c);
+ printf("lcm\n");
+ mp_todecimal(&a, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&b, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&c, buf);
+ printf("%s\n", buf);
+ } else if (n == 9) {
+ /* lcm test */
+ rand_num2(&a);
+ rand_num2(&b);
+ rand_num2(&c);
+ a.sign = b.sign = c.sign = 0;
+ mp_exptmod(&a, &b, &c, &d);
+ printf("expt\n");
+ mp_todecimal(&a, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&b, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&c, buf);
+ printf("%s\n", buf);
+ mp_todecimal(&d, buf);
+ printf("%s\n", buf);
+ }
+ }
+ fclose(rng);
+ return 0;
+}

0 comments on commit 8c97c9e

Please sign in to comment.