TeaSpeak/tommath
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

added libtommath-0.15

Browse Source
This commit is contained in:
Tom St Denis 2003-03-22 15:10:20 +00:00 committed by Steffen Jaeckel
parent 82f4858291
commit b1756f2f98
37 changed files with 7154 additions and 90 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
bn.pdf
View File

Binary file not shown.

372
bn.tex
View File

@ -1,7 +1,7 @@
\documentclass{article} \documentclass{article}
\begin{document} \begin{document}
\title{LibTomMath v0.14 \\ A Free Multiple Precision Integer Library \\ http://math.libtomcrypt.org } \title{LibTomMath v0.15 \\ A Free Multiple Precision Integer Library \\ http://math.libtomcrypt.org }
\author{Tom St Denis \\ tomstdenis@iahu.ca} \author{Tom St Denis \\ tomstdenis@iahu.ca}
\maketitle \maketitle
\newpage \newpage
@ -100,6 +100,22 @@ in the order $x, y, z$. For example:
mp_div_2(&x, &y); /* y = x / 2 */ mp_div_2(&x, &y); /* y = x / 2 */
\end{verbatim} \end{verbatim}
\subsection{Various Optimizations}
Various routines come in several ``flavours'' which are optimized for particular cases of inputs. For instance
the multiplicative inverse function ``mp\_invmod()'' has a routine for odd and even moduli. Similarly the
``mp\_exptmod()'' function has several variants depending on the modulus as well. Several lower level
functions such as multiplication, squaring and reductions come in ``comba'' and ``baseline'' variants.
The design of LibTomMath is such that the end user does not have to concern themselves too much with these
details. This is why the functions provided will determine \textit{automatically} when an appropriate
optimal function can be used. For example, when you call ``mp\_mul()'' the routines will first determine
if the Karatsuba multiplier should be used. If not it will determine if the ``comba'' method can be used
and finally call the standard catch-all ``baseline'' method.
Throughout the rest of this manual several variants for various functions will be referenced to as
the ``comba'', ``baseline'', etc... method. Keep in mind you call one function to use any of the optimal
variants.
\subsection{Return Values} \subsection{Return Values}
All functions that return errors will return \textbf{MP\_OKAY} if the function was succesful. It will return All functions that return errors will return \textbf{MP\_OKAY} if the function was succesful. It will return
\textbf{MP\_MEM} if it ran out of heap memory or \textbf{MP\_VAL} if one of the arguements is out of range. \textbf{MP\_MEM} if it ran out of heap memory or \textbf{MP\_VAL} if one of the arguements is out of range.
@ -326,10 +342,53 @@ int mp_montgomery_setup(mp_int *a, mp_digit *mp);
/* computes xR^-1 == x (mod N) via Montgomery Reduction */ /* computes xR^-1 == x (mod N) via Montgomery Reduction */
int mp_montgomery_reduce(mp_int *a, mp_int *m, mp_digit mp); int mp_montgomery_reduce(mp_int *a, mp_int *m, mp_digit mp);
/* returns 1 if a is a valid DR modulus */
int mp_dr_is_modulus(mp_int *a);
/* sets the value of "d" required for mp_dr_reduce */
void mp_dr_setup(mp_int *a, mp_digit *d);
/* reduces a modulo b using the Diminished Radix method */
int mp_dr_reduce(mp_int *a, mp_int *b, mp_digit mp);
/* d = a^b (mod c) */ /* d = a^b (mod c) */
int mp_exptmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d); int mp_exptmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
\end{verbatim} \end{verbatim}
\subsection{Primality Routines}
\begin{verbatim}
/* ---> Primes <--- */
/* table of first 256 primes */
extern const mp_digit __prime_tab[];
/* result=1 if a is divisible by one of the first 256 primes */
int mp_prime_is_divisible(mp_int *a, int *result);
/* performs one Fermat test of "a" using base "b".
* Sets result to 0 if composite or 1 if probable prime
*/
int mp_prime_fermat(mp_int *a, mp_int *b, int *result);
/* performs one Miller-Rabin test of "a" using base "b".
* Sets result to 0 if composite or 1 if probable prime
*/
int mp_prime_miller_rabin(mp_int *a, mp_int *b, int *result);
/* 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
* of error no more than (1/4)^t.
*
* Sets result to 1 if probably prime, 0 otherwise
*/
int mp_prime_is_prime(mp_int *a, int t, int *result);
/* finds the next prime after the number "a" using "t" trials
* of Miller-Rabin.
*/
int mp_prime_next_prime(mp_int *a, int t);
\end{verbatim}
\subsection{Radix Conversions} \subsection{Radix Conversions}
To read or store integers in other formats there are the following functions. To read or store integers in other formats there are the following functions.
@ -533,23 +592,131 @@ $n$ is prime then $\left ( {a \over n} \right )$ is equal to $1$ if $a$ is a qua
it is not. it is not.
\subsubsection{mp\_exptmod(mp\_int *a, mp\_int *b, mp\_int *c, mp\_int *d)} \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 Computes $d \equiv 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 + 2^{k-1}$ multiplications. The value of $k$ is exponent it performs $\alpha$ squarings and at most $\lfloor \alpha/k \rfloor + 2^{k-1}$ multiplications. The value of $k$ is
chosen to minimize the number of multiplications required for a given value of $\alpha$. Barrett or Montgomery chosen to minimize the number of multiplications required for a given value of $\alpha$. Barrett, Montgomery or
reductions are used to reduce the squared or multiplied temporary results modulo $c$. Dimminished-Radix reductions are used to reduce the squared or multiplied temporary results modulo $c$.
\subsection{Fast Modular Reductions} \subsection{Fast Modular Reductions}
A modular reduction of $a \mbox{ (mod }b\mbox{)}$ means to divide $a$ by $b$ and obtain the remainder.
Typically modular reductions are popular in public key cryptography algorithms such as RSA,
Diffie-Hellman and Elliptic Curve. Modular reductions are also a large portion of modular exponentiation
(e.g. $a^b \mbox{ (mod }c\mbox{)}$).
In a simplistic sense a normal integer division could be used to compute reduction. Division is by far
the most complicated of routines in terms of the work required. As a result it is desirable to avoid
division as much as possible. This is evident in quite a few fields in computing. For example, often in
signal analysis uses multiplication by the reciprocal to approximate divisions. Number theory is no
different.
In most cases for the reduction of $a$ modulo $b$ the integer $a$ will be limited to the range
$0 \le a \le b^2$ which led to the invention of specialized algorithms to do the work.
The first algorithm is the most generic and is called the Barrett reduction. When the input is of the
limited form (e.g. $0 \le a \le b^2$) Barrett reduction is numerically equivalent to a full integer
division with remainder. For a $n$-digit value $b$ the Barrett reduction requires approximately $2n^2$
multiplications.
The second algorithm is the Montgomery reduction. It is slightly different since the result is not
numerically equivalent to a standard integer division with remainder. Also this algorithm only works for
odd moduli. The final result can be converted easily back to the desired for which makes the reduction
technique useful for algorithms where only the final output is desired. For a $n$-digit value $b$ the
Montgomery reduction requires approximately $n^2 + n$ multiplications, about half as many as the
Barrett algorithm.
The third algorithm is the Diminished Radix ``DR'' reduction. It is a highly optimized reduction algorithm
suitable only for a limited set of moduli. For the specific moduli it is numerically equivalent to
integer division with remainder. For a $n$-digit value $b$ the DR reduction rquires exactly $n$
multiplications which is considerably faster than either of the two previous algorithms.
All three algorithms are automatically used in the modular exponentiation function mp\_exptmod() when
appropriate moduli are detected.
\begin{figure}[here]
\begin{small}
\begin{center}
\begin{tabular}{|c|c|l|}
\hline \textbf{Algorithm} & \textbf{Multiplications} & \textbf{Limitations} \\
Barrett Reduction & $2n^2$ & Any modulus. \\
Montgomery Reduction & $n^2 + n$ & Any odd modulus. \\
DR Reduction & $n$ & Moduli of the form $p = \beta^k - p'$.\\
\hline
\end{tabular}
\caption{Summary of reduction techniques.}
\end{center}
\end{small}
\end{figure}
\subsubsection{mp\_reduce(mp\_int *a, mp\_int *b, mp\_int *c)} \subsubsection{mp\_reduce(mp\_int *a, mp\_int *b, mp\_int *c)}
Computes a Barrett reduction in-place of $a$ modulo $b$ with respect to $c$. In essence it computes Computes a Barrett reduction in-place of $a$ modulo $b$ with respect to $c$. In essence it computes
$a \equiv a \mbox{ (mod }b\mbox{)}$ provided $0 \le a \le b^2$. The value of $c$ is precomputed with the $a \mbox{ (mod }b\mbox{)}$ provided $0 \le a \le b^2$. The value of $c$ is precomputed with the
function mp\_reduce\_setup(). The modulus $b$ must be larger than zero. function mp\_reduce\_setup(). The modulus $b$ must be larger than zero.
This reduction function is much faster than simply calling mp\_mod() (\textit{Which simply uses mp\_div() anyways}) and is
desirable where ever an appropriate reduction is desired.
The Barrett reduction function has been optimized to use partial multipliers which means compared to MPI it performs The Barrett reduction function has been optimized to use partial multipliers which means compared to MPI it performs
have the number of single precision multipliers (\textit{provided they have the same size digits}). The partial have the number of single precision multipliers (\textit{provided they have the same size digits}). The partial
multipliers (\textit{one of which is shared with mp\_mul}) both have baseline and comba variants. Barrett reduction multipliers (\textit{one of which is shared with mp\_mul}) both have baseline and comba variants. Barrett reduction
can reduce a number modulo a $n-$digit modulus with approximately $2n^2$ single precision multiplications. can reduce a number modulo a $n-$digit modulus with approximately $2n^2$ single precision multiplications.
Consider the following snippet (from a BBS generator) using the more traditional approach:
\begin{small}
\begin{verbatim}
mp_int modulus, n;
unsigned char buf[128];
int ix, err;
/* ... init code ..., e.g. init modulus and n */
/* now output 128 bytes */
for (ix = 0; ix < 128; ix++) {
if ((err = mp_sqrmod(&n, &modulus, &n)) != MP_OKAY) {
printf("Err: %d\n", err);
exit(EXIT_FAILURE);
}
buf[ix] = n->dp[0] & 255;
}
\end{verbatim}
\end{small}
And now consider the same function using Barrett reductions:
\begin{small}
\begin{verbatim}
mp_int modulus, n, mp;
unsigned char buf[128];
int ix, err;
/* ... init code ... e.g. modulus and n */
/* now setup mp which is the Barrett param */
if ((err = mp_reduce_setup(&mp, &modulus)) != MP_OKAY) {
printf("Err: %d\n", err);
exit(EXIT_FAILURE);
}
/* now output 128 bytes */
for (ix = 0; ix < 128; ix++) {
/* square n */
if ((err = mp_sqr(&n, &n)) != MP_OKAY) {
printf("Err: %d\n", err);
exit(EXIT_FAILURE);
}
/* now reduce the square modulo modulus */
if ((err = mp_reduce(&n, &modulus, &mp)) != MP_OKAY) {
printf("Err: %d\n", err);
exit(EXIT_FAILURE);
}
buf[ix] = n->dp[0] & 255;
}
\end{verbatim}
\end{small}
Both routines will produce the same output provided the same initial values of $modulus$ and $n$. The Barrett
method seems like more work but the optimization stems from the use of the Barrett reduction instead of the normal
integer division.
\subsubsection{mp\_montgomery\_reduce(mp\_int *a, mp\_int *m, mp\_digit mp)} \subsubsection{mp\_montgomery\_reduce(mp\_int *a, mp\_int *m, mp\_digit mp)}
Computes a Montgomery reduction in-place of $a$ modulo $b$ with respect to $mp$. If $b$ is some $n-$digit modulus then Computes a Montgomery reduction in-place of $a$ modulo $b$ with respect to $mp$. If $b$ is some $n-$digit modulus then
$R = \beta^{n+1}$. The result of this function is $aR^{-1} \mbox{ (mod }b\mbox{)}$ provided that $0 \le a \le b^2$. $R = \beta^{n+1}$. The result of this function is $aR^{-1} \mbox{ (mod }b\mbox{)}$ provided that $0 \le a \le b^2$.
@ -578,6 +745,94 @@ two long divisions would be required to setup $\hat x$ and a multiplication foll
A very useful observation is that multiplying by $R = \beta^n$ amounts to performing a left shift by $n$ positions which A very useful observation is that multiplying by $R = \beta^n$ amounts to performing a left shift by $n$ positions which
requires no single precision multiplications. requires no single precision multiplications.
\subsubsection{mp\_dr\_reduce(mp\_int *a, mp\_int *b, mp\_digit mp)}
Computes the Diminished-Radix reduction of $a$ in place modulo $b$ with respect to $mp$. $a$ must be in the range
$0 \le a \le b^2$ and $mp$ must be precomputed with the function mp\_dr\_setup().
This reduction technique performs the reduction with $n$ multiplications and is much faster than either of the previous
reduction methods. Essentially it is very much like the Montgomery reduction except it is particularly optimized for
specific types of moduli. The moduli must be of the form $p = \beta^k - p'$ where $0 \le p' < \beta$ for $k \ge 2$.
This algorithm is suitable for several applications such as Diffie-Hellman public key cryptsystems where the prime $p$ is
of this form.
In appendix A several ``safe'' primes of various sizes are provided. These primes are DR moduli and of the form
$p = 2q + 1$ where both $p$ and $q$ are prime. A trivial observation is that $g = 4$ will be a generator for all of them
since the order of the multiplicative sub-group is at most $2q$. Since $2^2 \ne 1$ that means $4^q \equiv 2^{2q} \equiv 1$
and that $g = 4$ is a generator of order $q$.
These moduli can be used to construct a Diffie-Hellman public key cryptosystem. Since the moduli are of the
DR form the modular exponentiation steps will be efficient.
\subsection{Primality Testing and Generation}
\subsubsection{mp\_prime\_is\_divisible(mp\_int *a, int *result)}
Determines if $a$ is divisible by any of the first 256 primes. Sets $result$ to $1$ if true or $0$
otherwise. Also will set $result$ to $1$ if $a$ is equal to one of the first 256 primes.
\subsubsection{mp\_prime\_fermat(mp\_int *a, mp\_int *b, int *result)}
Determines if $b$ is a witness to the compositeness of $a$ using the Fermat test. Essentially this
computes $b^a \mbox{ (mod }a\mbox{)}$ and compares it to $b$. If they match $result$ is set
to $1$ otherwise it is set to $0$. If $a$ is prime and $1 < b < a$ then this function will set
$result$ to $1$ with a probability of one. If $a$ is composite then this function will set
$result$ to $1$ with a probability of no more than $1 \over 2$.
If this function is repeated $t$ times with different bases $b$ then the probability of a false positive
is no more than $2^{-t}$.
\subsubsection{mp\_prime\_miller\_rabin(mp\_int *a, mp\_int *b, int *result)}
Determines if $b$ is a witness to the compositeness of $a$ using the Miller-Rabin test. This test
works much (\textit{on an abstract level}) the same as the Fermat test except is more robust. The
set of pseudo-primes to any given base for the Miller-Rabin test is a proper subset of the pseudo-primes
for the Fermat test.
If $a$ is prime and $1 < b < a$ then this function will always set $result$ to $1$. If $a$ is composite
the trivial bound of error is $1 \over 4$. However, according to HAC\footnote{Handbook of Applied
Cryptography, Chapter 4, Section 4, pp. 147, Fact 4.48.} the following bounds are
known. For a test of $t$ trials on a $k$-bit number the probability $P_{k,t}$ of error is given as
follows.
\begin{enumerate}
\item $P_{k,1} < k^24^{2 - \sqrt{k}}$ for $k \ge 2$
\item $P_{k,t} < k^{3/2}2^tt^{-1/2}4^{2-\sqrt{tk}}$ for $(t = 2, k \ge 88)$ or $(3 \le t \le k/9, k \ge 21)$.
\item $P_{k,t} < {7 \over 20}k2^{-5t} + {1 \over 7}k^{15/4}2^{-k/2-2t} + 12k2^{-k/4-3t}$ for $k/9 \le t \le k/4, k \ge 21$.
\item $P_{k,t} < {1 \over 7}k^{15/4}2^{-k/2 - 2t}$ for $t \ge k/4, k \ge 21$.
\end{enumerate}
For instance, $P_{1024,1}$ which indicates the probability of failure of one test with a 1024-bit candidate
is no more than $2^{-40}$. However, ideally at least a couple of trials should be used. In LibTomCrypt
for instance eight tests are used. In this case $P_{1024,8}$ falls under the second rule which leads
to a probability of failure of no more than $2^{-155.52}$.
\begin{figure}[here]
\begin{small}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline \textbf{Size (k)} & \textbf{$t = 3$} & \textbf{$t = 4$} & \textbf{$t = 5$} & \textbf{$t = 6$} & \textbf{$t = 7$} & \textbf{$t = 8$}\\
\hline 512 & -58 & -70 & -79 & -88 & -96 & -104 \\
\hline 768 & -75 & -89 & -101 & -112 & -122 & -131\\
\hline 1024 & -89 & -106 & -120 & -133 & -144 & -155 \\
\hline 1280 & -102 & -120 & -136 & -151 & -164 & -176 \\
\hline 1536 & -113 & -133 & -151 & -167 & -181 & -195 \\
\hline 1792 & -124 & -146 & -165 & -182 & -198 & -212 \\
\hline 2048 & -134 & -157 & -178 & -196 & -213 & -228\\
\hline
\end{tabular}
\end{center}
\end{small}
\caption{Probability of error for a given random candidate of $k$ bits with $t$ trials. Denoted as
log$_2(P_{k,t})$. }
\end{figure}
\subsubsection{mp\_prime\_is\_prime(mp\_int *a, int t, int *result)}
This function determines if $a$ is probably prime by first performing trial division by the first 256
primes and then $t$ rounds of Miller-Rabin using the first $t$ primes as bases. If $a$ is prime this
function will always set $result$ to $1$. If $a$ is composite then it will almost always set $result$
to $0$. The probability of error is given in figure two.
\subsubsection{mp\_prime\_next\_prime(mp\_int *a, int t)}
This function will find the next prime \textbf{after} $a$ by using trial division and $t$ trials of
Miller-Rabin.
\section{Timing Analysis} \section{Timing Analysis}
\subsection{Digit Size} \subsection{Digit Size}
@ -662,8 +917,12 @@ MPI uses a binary square-multiply method for exponentiation. For the same expon
perform 8 squarings and 5 multiplications. There is a precomputation phase for the method LibTomMath uses but it perform 8 squarings and 5 multiplications. 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 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 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 and 512 multiplications.
technique on top of the savings the $k$-ary method provides.
Randomly the most probable event is that every $2k^2$ bits another multiplication is saved via the
sliding-window technique on top of the savings the $k$-ary method provides. This stems from the fact that each window
has a probability of $2^{-1}$ of being delayed by one bit. In reality the savings can be much more when the exponent
has an abundance of zero bits.
Both LibTomMath and MPI use Barrett reduction instead of division to reduce the numbers modulo the modulus given. 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 However, LibTomMath can take advantage of the fact that the multiplications required within the Barrett reduction
@ -671,12 +930,103 @@ do not have to give full precision. As a result the reduction step is much fast
code will automatically determine at run-time (e.g. when its called) whether the faster multiplier can be used. The 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 comba baseline). faster multipliers have also been optimized into the two variants (baseline and comba baseline).
LibTomMath also has a variant of the exptmod function that uses Montgomery reductions instead of Barrett reductions LibTomMath also has a variant of the exptmod function that uses Montgomery or Diminished-Radix reductions instead of
which is faster. The code will automatically detect when the Montgomery version can be used (\textit{Requires the Barrett reductions which are faster. The code will automatically detect when the Montgomery version can be used
modulus to be odd and below the MONTGOMERY\_EXPT\_CUTOFF size}). The Montgomery routine is essentially a copy of the (\textit{Requires the modulus to be odd and below the MONTGOMERY\_EXPT\_CUTOFF size}). The Montgomery routine is
Barrett exponentiation routine except it uses Montgomery reduction. essentially a copy of the Barrett exponentiation routine except it uses Montgomery reduction.
As a result of all these changes exponentiation in LibTomMath is much faster than compared to MPI. On most ALU-strong As a result of all these changes exponentiation in LibTomMath is much faster than compared to MPI. On most ALU-strong
processors (AMD Athlon for instance) exponentiation in LibTomMath is often more then ten times faster than MPI. processors (AMD Athlon for instance) exponentiation in LibTomMath is often more then ten times faster than MPI.
\newpage
\section*{Appendix A -- DR Safe Prime Moduli}
These are safe primes suitable for the DR reduction techniques.
\begin{small}
\begin{verbatim}
224-bit prime:
p == 26959946667150639794667015087019630673637144422540572481103341844143
532-bit prime:
p == 14059105607947488696282932836518693308967803494693489478439861164411
99243959839959474700214407465892859350284572975279726002583142341968
6528151609940203368691747
784-bit prime:
p == 10174582569701926077392351975587856746131528201775982910760891436407
52752352543956225804474009941755789631639189671820136396606697711084
75957692810857098847138903161308502419410142185759152435680068435915
159402496058513611411688900243039
1036-bit prime:
p == 73633510803960459580592340614718453088992337057476877219196961242207
30400993319449915739231125812675425079864519532271929704028930638504
85730703075899286013451337291468249027691733891486704001513279827771
74018362916106519487472796251714810077522836342108369176406547759082
3919364012917984605619526140821798437127
1540-bit prime:
p == 38564998830736521417281865696453025806593491967131023221754800625044
11826546885121070536038571753679461518026049420807660579867166071933
31995138078062523944232834134301060035963325132466829039948295286901
98205120921557533726473585751382193953592127439965050261476810842071
57368450587885458870662348457392592590350574754547108886771218500413
52012892734056144158994382765356263460989042410208779740029161680999
51885406379295536200413493190419727789712076165162175783
2072-bit prime:
p == 54218939133169617266167044061918053674999416641599333415160174539219
34845902966009796023786766248081296137779934662422030250545736925626
89251250471628358318743978285860720148446448885701001277560572526947
61939255157449083928645845499448866574499182283776991809511712954641
41244487770339412235658314203908468644295047744779491537946899487476
80362212954278693335653935890352619041936727463717926744868338358149
56836864340303776864961677852601361049369618605589931826833943267154
13281957242613296066998310166663594408748431030206661065682224010477
20269951530296879490444224546654729111504346660859907296364097126834
834235287147
\end{verbatim}
\newpage
\begin{verbatim}
3080-bit prime:
p == 14872591348147092640920326485259710388958656451489011805853404549855
24155135260217788758027400478312256339496385275012465661575576202252
06314569873207988029466422057976484876770407676185319721656326266004
66027039730507982182461708359620055985616697068444694474354610925422
65792444947706769615695252256130901271870341005768912974433684521436
21126335809752272646208391793909176002665892575707673348417320292714
14414925737999142402226287954056239531091315945236233530448983394814
94120112723445689647986475279242446083151413667587008191682564376412
34796414611389856588668313940700594138366932599747507691048808666325
63356891811579575714450674901879395531659037735542902605310091218790
44170766615232300936675369451260747671432073394867530820527479172464
10644245072764022650374658634027981631882139521072626829153564850619
07146160831634031899433344310568760382865303657571873671474460048559
12033137386225053275419626102417236133948503
4116-bit prime:
p == 10951211157166778028568112903923951285881685924091094949001780089679
55253005183831872715423151551999734857184538199864469605657805519106
71752965504405483319768745978263629725521974299473675154181526972794
07518606702687749033402960400061140139713092570283328496790968248002
50742691718610670812374272414086863715763724622797509437062518082383
05605014462496277630214789052124947706021514827516368830127584715531
60422794055576326393660668474428614221648326558746558242215778499288
63023018366835675399949740429332468186340518172487073360822220449055
34058256846156864525995487330361695377639385317484513208112197632746
27403549307444874296172025850155107442985301015477068215901887335158
80733527449780963163909830077616357506845523215289297624086914545378
51108253422962011656326016849452390656670941816601111275452976618355
45793212249409511773940884655967126200762400673705890369240247283750
76210477267488679008016579588696191194060127319035195370137160936882
40224439969917201783514453748848639690614421772002899286394128821718
53539149915834004216827510006035966557909908155251261543943446413363
97793791497068253936771017031980867706707490224041075826337383538651
82549367950377193483609465580277633166426163174014828176348776585274
6577808019633679
\end{verbatim}
\end{small}
\end{document} \end{document}

1
bn_fast_mp_invmod.c
View File

@ -80,7 +80,6 @@ fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c)
} }
mp_set (&D, 1); mp_set (&D, 1);
top: top:
/* 4. while u is even do */ /* 4. while u is even do */
while (mp_iseven (&u) == 1) { while (mp_iseven (&u) == 1) {

4
bn_mp_div.c
View File

@ -106,7 +106,7 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
/* step 3. for i from n down to (t + 1) */ /* step 3. for i from n down to (t + 1) */
for (i = n; i >= (t + 1); i--) { for (i = n; i >= (t + 1); i--) {
if (i > x.alloc) if (i > x.used)
continue; continue;
/* step 3.1 if xi == yt then set q{i-t-1} to b-1, otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */ /* step 3.1 if xi == yt then set q{i-t-1} to b-1, otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
@ -175,6 +175,7 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
/* now q is the quotient and x is the remainder [which we have to normalize] */ /* now q is the quotient and x is the remainder [which we have to normalize] */
/* get sign before writing to c */ /* get sign before writing to c */
x.sign = a->sign; x.sign = a->sign;
if (c != NULL) { if (c != NULL) {
mp_clamp (&q); mp_clamp (&q);
mp_exch (&q, c); mp_exch (&q, c);
@ -183,7 +184,6 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
if (d != NULL) { if (d != NULL) {
mp_div_2d (&x, norm, &x, NULL); mp_div_2d (&x, norm, &x, NULL);
mp_clamp (&x);
mp_exch (&x, d); mp_exch (&x, d);
} }

2
bn_mp_div_2d.c
View File

@ -52,7 +52,7 @@ mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d)
/* shift by as many digits in the bit count */ /* shift by as many digits in the bit count */
if (b >= DIGIT_BIT) { if (b >= DIGIT_BIT) {
mp_rshd (c, b / DIGIT_BIT); mp_rshd (c, b / DIGIT_BIT);
} }
/* shift any bit count < DIGIT_BIT */ /* shift any bit count < DIGIT_BIT */

1
bn_mp_div_d.c
View File

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

150
bn_mp_dr_reduce.c Normal file
View File

@ -0,0 +1,150 @@
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
* LibTomMath is library that provides for multiple-precision
* integer arithmetic as well as number theoretic functionality.
*
* The library is 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>
/* reduce "a" in place modulo "b" using the Diminished Radix algorithm.
*
* Based on algorithm from the paper
*
* "Generating Efficient Primes for Discrete Log Cryptosystems"
* Chae Hoon Lim, Pil Loong Lee,
* POSTECH Information Research Laboratories
*
* The modulus must be of a special format [see manual]
*/
int
mp_dr_reduce (mp_int * a, mp_int * b, mp_digit mp)
{
int err, i, j, k;
mp_word r;
mp_digit mu, *tmpj, *tmpi;
/* k = digits in modulus */
k = b->used;
/* ensure that "a" has at least 2k digits */
if (a->alloc < k + k) {
if ((err = mp_grow (a, k + k)) != MP_OKAY) {
return err;
}
}
/* alias for a->dp[i] */
tmpi = a->dp + k + k - 1;
/* for (i = 2k - 1; i >= k; i = i - 1)
*
* This is the main loop of the reduction. Note that at the end
* the words above position k are not zeroed as expected. The end
* result is that the digits from 0 to k-1 are the residue. So
* we have to clear those afterwards.
*/
for (i = k + k - 1; i >= k; i = i - 1) {
/* x[i - 1 : i - k] += x[i]*mp */
/* x[i] * mp */
r = ((mp_word) *tmpi--) * ((mp_word) mp);
/* now add r to x[i-1:i-k]
*
* First add it to the first digit x[i-k] then form the carry
* then enter the main loop
*/
j = i - k;
/* alias for a->dp[j] */
tmpj = a->dp + j;
/* add digit */
*tmpj += (mp_digit)(r & MP_MASK);
/* this is the carry */
mu = (r >> ((mp_word) DIGIT_BIT)) + (*tmpj >> DIGIT_BIT);
/* clear carry from a->dp[j] */
*tmpj++ &= MP_MASK;
/* now add rest of the digits
*
* Note this is basically a simple single digit addition to
* a larger multiple digit number. This is optimized somewhat
* because the propagation of carries is not likely to move
* more than a few digits.
*
*/
for (++j; mu != 0 && j <= (i - 1); ++j) {
*tmpj += mu;
mu = *tmpj >> DIGIT_BIT;
*tmpj++ &= MP_MASK;
}
/* if final carry */
if (mu != 0) {
/* add mp to this to correct */
j = i - k;
tmpj = a->dp + j;
*tmpj += mp;
mu = *tmpj >> DIGIT_BIT;
*tmpj++ &= MP_MASK;
/* now handle carries */
for (++j; mu != 0 && j <= (i - 1); j++) {
*tmpj += mu;
mu = *tmpj >> DIGIT_BIT;
*tmpj++ &= MP_MASK;
}
}
}
/* zero words above k */
tmpi = a->dp + k;
for (i = k; i < a->used; i++) {
*tmpi++ = 0;
}
/* clamp, sub and return */
mp_clamp (a);
if (mp_cmp_mag (a, b) != MP_LT) {
return s_mp_sub (a, b, a);
}
return MP_OKAY;
}
/* determines if a number is a valid DR modulus */
int mp_dr_is_modulus(mp_int *a)
{
int ix;
/* must be at least two digits */
if (a->used < 2) {
return 0;
}
for (ix = 1; ix < a->used; ix++) {
if (a->dp[ix] != MP_MASK) {
return 0;
}
}
return 1;
}
/* determines the setup value */
void mp_dr_setup(mp_int *a, mp_digit *d)
{
*d = (1 << DIGIT_BIT) - a->dp[0];
}

7
bn_mp_exptmod.c
View File

@ -24,9 +24,12 @@ static int f_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y);
int int
mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y) mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
{ {
int dr;
dr = mp_dr_is_modulus(P);
/* if the modulus is odd use the fast method */ /* if the modulus is odd use the fast method */
if (mp_isodd (P) == 1 && P->used > 4 && P->used < MONTGOMERY_EXPT_CUTOFF) { if (((mp_isodd (P) == 1 && P->used < MONTGOMERY_EXPT_CUTOFF) || dr == 1) && P->used > 4) {
return mp_exptmod_fast (G, X, P, Y); return mp_exptmod_fast (G, X, P, Y, dr);
} else { } else {
return f_mp_exptmod (G, X, P, Y); return f_mp_exptmod (G, X, P, Y);
} }

61
bn_mp_exptmod_fast.c
View File

@ -22,11 +22,13 @@
* Uses Montgomery reduction * Uses Montgomery reduction
*/ */
int int
mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y) mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
{ {
mp_int M[256], res; mp_int M[256], res;
mp_digit buf, mp; mp_digit buf, mp;
int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize; int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
int (*redux)(mp_int*,mp_int*,mp_digit);
/* find window size */ /* find window size */
x = mp_count_bits (X); x = mp_count_bits (X);
@ -56,9 +58,16 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
} }
} }
/* now setup montgomery */ if (redmode == 0) {
if ((err = mp_montgomery_setup (P, &mp)) != MP_OKAY) { /* now setup montgomery */
goto __M; if ((err = mp_montgomery_setup (P, &mp)) != MP_OKAY) {
goto __M;
}
redux = mp_montgomery_reduce;
} else {
/* setup DR reduction */
mp_dr_setup(P, &mp);
redux = mp_dr_reduce;
} }
/* setup result */ /* setup result */
@ -73,15 +82,23 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
* The first half of the table is not computed though accept for M[0] and M[1] * The first half of the table is not computed though accept for M[0] and M[1]
*/ */
/* now we need R mod m */ if (redmode == 0) {
if ((err = mp_montgomery_calc_normalization (&res, P)) != MP_OKAY) { /* now we need R mod m */
goto __RES; if ((err = mp_montgomery_calc_normalization (&res, P)) != MP_OKAY) {
goto __RES;
}
/* now set M[1] to G * R mod m */
if ((err = mp_mulmod (G, &res, P, &M[1])) != MP_OKAY) {
goto __RES;
}
} else {
mp_set(&res, 1);
if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) {
goto __RES;
}
} }
/* now set M[1] to G * R mod m */
if ((err = mp_mulmod (G, &res, P, &M[1])) != MP_OKAY) {
goto __RES;
}
/* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */ /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) { if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
goto __RES; goto __RES;
@ -91,7 +108,7 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
if ((err = mp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)])) != MP_OKAY) { if ((err = mp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)])) != MP_OKAY) {
goto __RES; goto __RES;
} }
if ((err = mp_montgomery_reduce (&M[1 << (winsize - 1)], P, mp)) != MP_OKAY) { if ((err = redux (&M[1 << (winsize - 1)], P, mp)) != MP_OKAY) {
goto __RES; goto __RES;
} }
} }
@ -101,7 +118,7 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) { if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
goto __RES; goto __RES;
} }
if ((err = mp_montgomery_reduce (&M[x], P, mp)) != MP_OKAY) { if ((err = redux (&M[x], P, mp)) != MP_OKAY) {
goto __RES; goto __RES;
} }
} }
@ -141,7 +158,7 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
if ((err = mp_sqr (&res, &res)) != MP_OKAY) { if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
goto __RES; goto __RES;
} }
if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) { if ((err = redux (&res, P, mp)) != MP_OKAY) {
goto __RES; goto __RES;
} }
continue; continue;
@ -158,7 +175,7 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
if ((err = mp_sqr (&res, &res)) != MP_OKAY) { if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
goto __RES; goto __RES;
} }
if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) { if ((err = redux (&res, P, mp)) != MP_OKAY) {
goto __RES; goto __RES;
} }
} }
@ -167,7 +184,7 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) { if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
goto __RES; goto __RES;
} }
if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) { if ((err = redux (&res, P, mp)) != MP_OKAY) {
goto __RES; goto __RES;
} }
@ -184,7 +201,7 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
if ((err = mp_sqr (&res, &res)) != MP_OKAY) { if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
goto __RES; goto __RES;
} }
if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) { if ((err = redux (&res, P, mp)) != MP_OKAY) {
goto __RES; goto __RES;
} }
@ -194,16 +211,18 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) { if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
goto __RES; goto __RES;
} }
if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) { if ((err = redux (&res, P, mp)) != MP_OKAY) {
goto __RES; goto __RES;
} }
} }
} }
} }
/* fixup result */ if (redmode == 0) {
if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) { /* fixup result */
goto __RES; if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) {
goto __RES;
}
} }
mp_exch (&res, Y); mp_exch (&res, Y);

2
bn_mp_grow.c
View File

@ -24,7 +24,7 @@ mp_grow (mp_int * a, int size)
if (a->alloc < size) { if (a->alloc < size) {
size += (MP_PREC * 2) - (size & (MP_PREC - 1)); /* ensure there are always at least MP_PREC digits extra on top */ size += (MP_PREC * 2) - (size & (MP_PREC - 1)); /* ensure there are always at least MP_PREC digits extra on top */
a->dp = realloc (a->dp, sizeof (mp_digit) * size); a->dp = OPT_CAST realloc (a->dp, sizeof (mp_digit) * size);
if (a->dp == NULL) { if (a->dp == NULL) {
return MP_MEM; return MP_MEM;
} }

2
bn_mp_init.c
View File

@ -20,7 +20,7 @@ mp_init (mp_int * a)
{ {
/* allocate ram required and clear it */ /* allocate ram required and clear it */
a->dp = calloc (sizeof (mp_digit), MP_PREC); a->dp = OPT_CAST calloc (sizeof (mp_digit), MP_PREC);
if (a->dp == NULL) { if (a->dp == NULL) {
return MP_MEM; return MP_MEM;
} }

2
bn_mp_init_size.c
View File

@ -21,7 +21,7 @@ mp_init_size (mp_int * a, int size)
/* pad up so there are at least 16 zero digits */ /* pad up so there are at least 16 zero digits */
size += (MP_PREC * 2) - (size & (MP_PREC - 1)); /* ensure there are always at least 16 digits extra on top */ size += (MP_PREC * 2) - (size & (MP_PREC - 1)); /* ensure there are always at least 16 digits extra on top */
a->dp = calloc (sizeof (mp_digit), size); a->dp = OPT_CAST calloc (sizeof (mp_digit), size);
if (a->dp == NULL) { if (a->dp == NULL) {
return MP_MEM; return MP_MEM;
} }

6
bn_mp_mul_2d.c
View File

@ -33,9 +33,9 @@ mp_mul_2d (mp_int * a, int b, mp_int * c)
/* shift by as many digits in the bit count */ /* shift by as many digits in the bit count */
if (b >= DIGIT_BIT) { if (b >= DIGIT_BIT) {
if ((res = mp_lshd (c, b / DIGIT_BIT)) != MP_OKAY) { if ((res = mp_lshd (c, b / DIGIT_BIT)) != MP_OKAY) {
return res; return res;
} }
} }
c->used = c->alloc; c->used = c->alloc;

52
bn_mp_prime_fermat.c Normal file
View File

@ -0,0 +1,52 @@
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
* LibTomMath is library that provides for multiple-precision
* integer arithmetic as well as number theoretic functionality.
*
* The library is 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>
/* performs one Fermat test.
*
* If "a" were prime then b^a == b (mod a) since the order of
* the multiplicative sub-group would be phi(a) = a-1. That means
* it would be the same as b^(a mod (a-1)) == b^1 == b (mod a).
*
* Sets result to 1 if the congruence holds, or zero otherwise.
*/
int
mp_prime_fermat (mp_int * a, mp_int * b, int *result)
{
mp_int t;
int err;
/* default to fail */
*result = 0;
/* init t */
if ((err = mp_init (&t)) != MP_OKAY) {
return err;
}
/* compute t = b^a mod a */
if ((err = mp_exptmod (b, a, a, &t)) != MP_OKAY) {
goto __T;
}
/* is it equal to b? */
if (mp_cmp (&t, b) == MP_EQ) {
*result = 1;
}
err = MP_OKAY;
__T:mp_clear (&t);
return err;
}

50
bn_mp_prime_is_divisible.c Normal file
View File

@ -0,0 +1,50 @@
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
* LibTomMath is library that provides for multiple-precision
* integer arithmetic as well as number theoretic functionality.
*
* The library is 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>
/* determines if an integers is divisible by one of the first 256 primes or not
*
* sets result to 0 if not, 1 if yes
*/
int
mp_prime_is_divisible (mp_int * a, int *result)
{
int err, ix;
mp_digit res;
/* default to not */
*result = 0;
for (ix = 0; ix < 256; ix++) {
/* is it equal to the prime? */
if (mp_cmp_d (a, __prime_tab[ix]) == MP_EQ) {
*result = 1;
return MP_OKAY;
}
/* what is a mod __prime_tab[ix] */
if ((err = mp_mod_d (a, __prime_tab[ix], &res)) != MP_OKAY) {
return err;
}
/* is the residue zero? */
if (res == 0) {
*result = 1;
return MP_OKAY;
}
}
return MP_OKAY;
}

68
bn_mp_prime_is_prime.c Normal file
View File

@ -0,0 +1,68 @@
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
* LibTomMath is library that provides for multiple-precision
* integer arithmetic as well as number theoretic functionality.
*
* The library is 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>
/* performs a variable number of rounds of Miller-Rabin
*
* Probability of error after t rounds is no more than
* (1/4)^t when 1 <= t <= 256
*
* Sets result to 1 if probably prime, 0 otherwise
*/
int
mp_prime_is_prime (mp_int * a, int t, int *result)
{
mp_int b;
int ix, err, res;
/* default to no */
*result = 0;
/* valid value of t? */
if (t < 1 || t > 256) {
return MP_VAL;
}
/* first perform trial division */
if ((err = mp_prime_is_divisible (a, &res)) != MP_OKAY) {
return err;
}
if (res == 1) {
return MP_OKAY;
}
/* now perform the miller-rabin rounds */
if ((err = mp_init (&b)) != MP_OKAY) {
return err;
}
for (ix = 0; ix < t; ix++) {
/* set the prime */
mp_set (&b, __prime_tab[ix]);
if ((err = mp_prime_miller_rabin (a, &b, &res)) != MP_OKAY) {
goto __B;
}
if (res == 0) {
goto __B;
}
}
/* passed the test */
*result = 1;
__B:mp_clear (&b);
return err;
}

90
bn_mp_prime_miller_rabin.c Normal file
View File

@ -0,0 +1,90 @@
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
* LibTomMath is library that provides for multiple-precision
* integer arithmetic as well as number theoretic functionality.
*
* The library is 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>
/* Miller-Rabin test of "a" to the base of "b" as described in
* HAC pp. 139 Algorithm 4.24
*
* Sets result to 0 if definitely composite or 1 if probably prime.
* 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)
{
mp_int n1, y, r;
int s, j, err;
/* default */
*result = 0;
/* get n1 = a - 1 */
if ((err = mp_init_copy (&n1, a)) != MP_OKAY) {
return err;
}
if ((err = mp_sub_d (&n1, 1, &n1)) != MP_OKAY) {
goto __N1;
}
/* set 2^s * r = n1 */
if ((err = mp_init_copy (&r, &n1)) != MP_OKAY) {
goto __N1;
}
s = 0;
while (mp_iseven (&r) == 1) {
++s;
if ((err = mp_div_2 (&r, &r)) != MP_OKAY) {
goto __R;
}
}
/* compute y = b^r mod a */
if ((err = mp_init (&y)) != MP_OKAY) {
goto __R;
}
if ((err = mp_exptmod (b, &r, a, &y)) != MP_OKAY) {
goto __Y;
}
/* if y != 1 and y != n1 do */
if (mp_cmp_d (&y, 1) != MP_EQ && mp_cmp (&y, &n1) != MP_EQ) {
j = 1;
/* while j <= s-1 and y != n1 */
while ((j <= (s - 1)) && mp_cmp (&y, &n1) != MP_EQ) {
if ((err = mp_sqrmod (&y, a, &y)) != MP_OKAY) {
goto __Y;
}
/* if y == 1 then composite */
if (mp_cmp_d (&y, 1) == MP_EQ) {
goto __Y;
}
++j;
}
/* if y != n1 then composite */
if (mp_cmp (&y, &n1) != MP_EQ) {
goto __Y;
}
}
/* probably prime now */
*result = 1;
__Y:mp_clear (&y);
__R:mp_clear (&r);
__N1:mp_clear (&n1);
return err;
}

54
bn_mp_prime_next_prime.c Normal file
View File

@ -0,0 +1,54 @@
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
* LibTomMath is library that provides for multiple-precision
* integer arithmetic as well as number theoretic functionality.
*
* The library is 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>
/* finds the next prime after the number "a" using "t" trials
* of Miller-Rabin.
*/
int mp_prime_next_prime(mp_int *a, int t)
{
int err, res;
if (mp_iseven(a) == 1) {
/* force odd */
if ((err = mp_add_d(a, 1, a)) != MP_OKAY) {
return err;
}
} else {
/* force to next number */
if ((err = mp_add_d(a, 2, a)) != MP_OKAY) {
return err;
}
}
for (;;) {
/* is this prime? */
if ((err = mp_prime_is_prime(a, t, &res)) != MP_OKAY) {
return err;
}
if (res == 1) {
break;
}
/* add two, next candidate */
if ((err = mp_add_d(a, 2, a)) != MP_OKAY) {
return err;
}
}
return MP_OKAY;
}

2
bn_mp_rshd.c
View File

@ -50,7 +50,7 @@ mp_rshd (mp_int * a, int b)
b-2 | b-1 | b0 | b1 | b2 | ... | bb | ----> b-2 | b-1 | b0 | b1 | b2 | ... | bb | ---->
/\ | ----> /\ | ---->
\-------------------/ ----> \-------------------/ ---->
*/ */
for (x = 0; x < (a->used - b); x++) { for (x = 0; x < (a->used - b); x++) {
*tmpa++ = *tmpaa++; *tmpa++ = *tmpaa++;
} }

2
bn_mp_shrink.c
View File

@ -19,7 +19,7 @@ int
mp_shrink (mp_int * a) mp_shrink (mp_int * a)
{ {
if (a->alloc != a->used) { if (a->alloc != a->used) {
if ((a->dp = realloc (a->dp, sizeof (mp_digit) * a->used)) == NULL) { if ((a->dp = OPT_CAST realloc (a->dp, sizeof (mp_digit) * a->used)) == NULL) {
return MP_MEM; return MP_MEM;
} }
a->alloc = a->used; a->alloc = a->used;

52
bn_prime_tab.c Normal file
View File

@ -0,0 +1,52 @@
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
* LibTomMath is library that provides for multiple-precision
* integer arithmetic as well as number theoretic functionality.
*
* The library is 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>
const mp_digit __prime_tab[] = {
0x0002, 0x0003, 0x0005, 0x0007, 0x000B, 0x000D, 0x0011, 0x0013,
0x0017, 0x001D, 0x001F, 0x0025, 0x0029, 0x002B, 0x002F, 0x0035,
0x003B, 0x003D, 0x0043, 0x0047, 0x0049, 0x004F, 0x0053, 0x0059,
0x0061, 0x0065, 0x0067, 0x006B, 0x006D, 0x0071, 0x007F, 0x0083,
0x0089, 0x008B, 0x0095, 0x0097, 0x009D, 0x00A3, 0x00A7, 0x00AD,
0x00B3, 0x00B5, 0x00BF, 0x00C1, 0x00C5, 0x00C7, 0x00D3, 0x00DF,
0x00E3, 0x00E5, 0x00E9, 0x00EF, 0x00F1, 0x00FB, 0x0101, 0x0107,
0x010D, 0x010F, 0x0115, 0x0119, 0x011B, 0x0125, 0x0133, 0x0137,
0x0139, 0x013D, 0x014B, 0x0151, 0x015B, 0x015D, 0x0161, 0x0167,
0x016F, 0x0175, 0x017B, 0x017F, 0x0185, 0x018D, 0x0191, 0x0199,
0x01A3, 0x01A5, 0x01AF, 0x01B1, 0x01B7, 0x01BB, 0x01C1, 0x01C9,
0x01CD, 0x01CF, 0x01D3, 0x01DF, 0x01E7, 0x01EB, 0x01F3, 0x01F7,
0x01FD, 0x0209, 0x020B, 0x021D, 0x0223, 0x022D, 0x0233, 0x0239,
0x023B, 0x0241, 0x024B, 0x0251, 0x0257, 0x0259, 0x025F, 0x0265,
0x0269, 0x026B, 0x0277, 0x0281, 0x0283, 0x0287, 0x028D, 0x0293,
0x0295, 0x02A1, 0x02A5, 0x02AB, 0x02B3, 0x02BD, 0x02C5, 0x02CF,
0x02D7, 0x02DD, 0x02E3, 0x02E7, 0x02EF, 0x02F5, 0x02F9, 0x0301,
0x0305, 0x0313, 0x031D, 0x0329, 0x032B, 0x0335, 0x0337, 0x033B,
0x033D, 0x0347, 0x0355, 0x0359, 0x035B, 0x035F, 0x036D, 0x0371,
0x0373, 0x0377, 0x038B, 0x038F, 0x0397, 0x03A1, 0x03A9, 0x03AD,
0x03B3, 0x03B9, 0x03C7, 0x03CB, 0x03D1, 0x03D7, 0x03DF, 0x03E5,
0x03F1, 0x03F5, 0x03FB, 0x03FD, 0x0407, 0x0409, 0x040F, 0x0419,
0x041B, 0x0425, 0x0427, 0x042D, 0x043F, 0x0443, 0x0445, 0x0449,
0x044F, 0x0455, 0x045D, 0x0463, 0x0469, 0x047F, 0x0481, 0x048B,
0x0493, 0x049D, 0x04A3, 0x04A9, 0x04B1, 0x04BD, 0x04C1, 0x04C7,
0x04CD, 0x04CF, 0x04D5, 0x04E1, 0x04EB, 0x04FD, 0x04FF, 0x0503,
0x0509, 0x050B, 0x0511, 0x0515, 0x0517, 0x051B, 0x0527, 0x0529,
0x052F, 0x0551, 0x0557, 0x055D, 0x0565, 0x0577, 0x0581, 0x058F,
0x0593, 0x0595, 0x0599, 0x059F, 0x05A7, 0x05AB, 0x05AD, 0x05B3,
0x05BF, 0x05C9, 0x05CB, 0x05CF, 0x05D1, 0x05D5, 0x05DB, 0x05E7,
0x05F3, 0x05FB, 0x0607, 0x060D, 0x0611, 0x0617, 0x061F, 0x0623,
0x062B, 0x062F, 0x063D, 0x0641, 0x0647, 0x0649, 0x064D, 0x0653
};

2
bn_radix.c
View File

@ -93,7 +93,7 @@ mp_toradix (mp_int * a, char *str, int radix)
*str++ = s_rmap[d]; *str++ = s_rmap[d];
++digs; ++digs;
} }
bn_reverse ((unsigned char *) _s, digs); bn_reverse ((unsigned char *)_s, digs);
*str++ = '\0'; *str++ = '\0';
mp_clear (&t); mp_clear (&t);
return MP_OKAY; return MP_OKAY;

2
bncore.c
View File

@ -18,5 +18,3 @@
int KARATSUBA_MUL_CUTOFF = 73, /* Min. number of digits before Karatsuba multiplication is used. */ int KARATSUBA_MUL_CUTOFF = 73, /* Min. number of digits before Karatsuba multiplication is used. */
KARATSUBA_SQR_CUTOFF = 121, /* Min. number of digits before Karatsuba squaring is used. */ KARATSUBA_SQR_CUTOFF = 121, /* Min. number of digits before Karatsuba squaring is used. */
MONTGOMERY_EXPT_CUTOFF = 128; /* max. number of digits that montgomery reductions will help for */ MONTGOMERY_EXPT_CUTOFF = 128; /* max. number of digits that montgomery reductions will help for */

11
changes.txt
View File

@ -1,3 +1,14 @@
Mar 22nd, 2003
v0.15 -- Added series of prime testing routines to lib
-- Fixed up etc/tune.c
-- Added DR reduction algorithm
-- Beefed up the manual more.
-- Fixed up demo/demo.c so it doesn't have so many warnings and it does the full series of
tests
-- Added "pre-gen" directory which will hold a "gen.pl"'ed copy of the entire lib [done at
zipup time so its always the latest]
-- Added conditional casts for C++ users [boo!]
Mar 15th, 2003 Mar 15th, 2003
v0.14 -- Tons of manual updates v0.14 -- Tons of manual updates
-- cleaned up the directory -- cleaned up the directory

52
demo/demo.c
View File

@ -89,7 +89,7 @@ int main(void)
unsigned long expt_n, add_n, sub_n, mul_n, div_n, sqr_n, mul2d_n, div2d_n, gcd_n, lcm_n, inv_n, unsigned long expt_n, add_n, sub_n, mul_n, div_n, sqr_n, mul2d_n, div2d_n, gcd_n, lcm_n, inv_n,
div2_n, mul2_n; div2_n, mul2_n;
unsigned rr; unsigned rr;
int cnt; int cnt, ix;
#ifdef TIMER #ifdef TIMER
int n; int n;
@ -104,9 +104,42 @@ int main(void)
mp_init(&e); mp_init(&e);
mp_init(&f); mp_init(&f);
/* test the DR reduction */
#if 0
srand(time(NULL));
for (cnt = 2; cnt < 32; cnt++) {
printf("%d digit modulus\n", cnt);
mp_grow(&a, cnt);
mp_zero(&a);
for (ix = 1; ix < cnt; ix++) {
a.dp[ix] = MP_MASK;
}
a.used = cnt;
mp_prime_next_prime(&a, 3);
mp_rand(&b, cnt - 1);
mp_copy(&b, &c);
rr = 0;
do {
if (!(rr & 127)) { printf("%9lu\r", rr); fflush(stdout); }
mp_sqr(&b, &b); mp_add_d(&b, 1, &b);
mp_copy(&b, &c);
mp_mod(&b, &a, &b);
mp_dr_reduce(&c, &a, (1<<DIGIT_BIT)-a.dp[0]);
if (mp_cmp(&b, &c) != MP_EQ) {
printf("Failed on trial %lu\n", rr); exit(-1);
}
} while (++rr < 1000000);
printf("Passed DR test for %d digits\n", cnt);
}
#endif
#ifdef TIMER #ifdef TIMER
printf("CLOCKS_PER_SEC == %lu\n", CLOCKS_PER_SEC); printf("CLOCKS_PER_SEC == %lu\n", CLOCKS_PER_SEC);
goto expttime;
log = fopen("add.log", "w"); log = fopen("add.log", "w");
for (cnt = 4; cnt <= 128; cnt += 4) { for (cnt = 4; cnt <= 128; cnt += 4) {
@ -136,7 +169,6 @@ goto expttime;
} }
fclose(log); fclose(log);
multtime:
log = fopen("sqr.log", "w"); log = fopen("sqr.log", "w");
for (cnt = 4; cnt <= 128; cnt += 4) { for (cnt = 4; cnt <= 128; cnt += 4) {
@ -165,9 +197,18 @@ multtime:
} }
fclose(log); fclose(log);
expttime:
{ {
char *primes[] = { char *primes[] = {
/* DR moduli */
"14059105607947488696282932836518693308967803494693489478439861164411992439598399594747002144074658928593502845729752797260025831423419686528151609940203368612079",
"101745825697019260773923519755878567461315282017759829107608914364075275235254395622580447400994175578963163918967182013639660669771108475957692810857098847138903161308502419410142185759152435680068435915159402496058513611411688900243039",
"736335108039604595805923406147184530889923370574768772191969612422073040099331944991573923112581267542507986451953227192970402893063850485730703075899286013451337291468249027691733891486704001513279827771740183629161065194874727962517148100775228363421083691764065477590823919364012917984605619526140821797602431",
"38564998830736521417281865696453025806593491967131023221754800625044118265468851210705360385717536794615180260494208076605798671660719333199513807806252394423283413430106003596332513246682903994829528690198205120921557533726473585751382193953592127439965050261476810842071573684505878854588706623484573925925903505747545471088867712185004135201289273405614415899438276535626346098904241020877974002916168099951885406379295536200413493190419727789712076165162175783",
"542189391331696172661670440619180536749994166415993334151601745392193484590296600979602378676624808129613777993466242203025054573692562689251250471628358318743978285860720148446448885701001277560572526947619392551574490839286458454994488665744991822837769918095117129546414124448777033941223565831420390846864429504774477949153794689948747680362212954278693335653935890352619041936727463717926744868338358149568368643403037768649616778526013610493696186055899318268339432671541328195724261329606699831016666359440874843103020666106568222401047720269951530296879490444224546654729111504346660859907296364097126834834235287147",
"1487259134814709264092032648525971038895865645148901180585340454985524155135260217788758027400478312256339496385275012465661575576202252063145698732079880294664220579764848767704076761853197216563262660046602703973050798218246170835962005598561669706844469447435461092542265792444947706769615695252256130901271870341005768912974433684521436211263358097522726462083917939091760026658925757076733484173202927141441492573799914240222628795405623953109131594523623353044898339481494120112723445689647986475279242446083151413667587008191682564376412347964146113898565886683139407005941383669325997475076910488086663256335689181157957571445067490187939553165903773554290260531009121879044170766615232300936675369451260747671432073394867530820527479172464106442450727640226503746586340279816318821395210726268291535648506190714616083163403189943334431056876038286530365757187367147446004855912033137386225053275419626102417236133948503",
"1095121115716677802856811290392395128588168592409109494900178008967955253005183831872715423151551999734857184538199864469605657805519106717529655044054833197687459782636297255219742994736751541815269727940751860670268774903340296040006114013971309257028332849679096824800250742691718610670812374272414086863715763724622797509437062518082383056050144624962776302147890521249477060215148275163688301275847155316042279405557632639366066847442861422164832655874655824221577849928863023018366835675399949740429332468186340518172487073360822220449055340582568461568645259954873303616953776393853174845132081121976327462740354930744487429617202585015510744298530101547706821590188733515880733527449780963163909830077616357506845523215289297624086914545378511082534229620116563260168494523906566709418166011112754529766183554579321224940951177394088465596712620076240067370589036924024728375076210477267488679008016579588696191194060127319035195370137160936882402244399699172017835144537488486396906144217720028992863941288217185353914991583400421682751000603596655790990815525126154394344641336397793791497068253936771017031980867706707490224041075826337383538651825493679503771934836094655802776331664261631740148281763487765852746577808019633679",
/* generic unrestricted moduli */
"17933601194860113372237070562165128350027320072176844226673287945873370751245439587792371960615073855669274087805055507977323024886880985062002853331424203", "17933601194860113372237070562165128350027320072176844226673287945873370751245439587792371960615073855669274087805055507977323024886880985062002853331424203",
"2893527720709661239493896562339544088620375736490408468011883030469939904368086092336458298221245707898933583190713188177399401852627749210994595974791782790253946539043962213027074922559572312141181787434278708783207966459019479487", "2893527720709661239493896562339544088620375736490408468011883030469939904368086092336458298221245707898933583190713188177399401852627749210994595974791782790253946539043962213027074922559572312141181787434278708783207966459019479487",
"347743159439876626079252796797422223177535447388206607607181663903045907591201940478223621722118173270898487582987137708656414344685816179420855160986340457973820182883508387588163122354089264395604796675278966117567294812714812796820596564876450716066283126720010859041484786529056457896367683122960411136319", "347743159439876626079252796797422223177535447388206607607181663903045907591201940478223621722118173270898487582987137708656414344685816179420855160986340457973820182883508387588163122354089264395604796675278966117567294812714812796820596564876450716066283126720010859041484786529056457896367683122960411136319",
@ -208,7 +249,7 @@ expttime:
} }
} }
fclose(log); fclose(log);
invtime:
log = fopen("invmod.log", "w"); log = fopen("invmod.log", "w");
for (cnt = 4; cnt <= 128; cnt += 4) { for (cnt = 4; cnt <= 128; cnt += 4) {
mp_rand(&a, cnt); mp_rand(&a, cnt);
@ -241,7 +282,6 @@ invtime:
div2_n = mul2_n = inv_n = expt_n = lcm_n = gcd_n = add_n = div2_n = mul2_n = inv_n = expt_n = lcm_n = gcd_n = add_n =
sub_n = mul_n = div_n = sqr_n = mul2d_n = div2d_n = cnt = 0; sub_n = mul_n = div_n = sqr_n = mul2d_n = div2d_n = cnt = 0;
for (;;) { for (;;) {
if (!(++cnt & 15)) sleep(3);
/* randomly clear and re-init one variable, this has the affect of triming the alloc space */ /* randomly clear and re-init one variable, this has the affect of triming the alloc space */
switch (abs(rand()) % 7) { switch (abs(rand()) % 7) {

53
etc/drprime.c Normal file
View File

@ -0,0 +1,53 @@
/* Makes safe primes of a DR nature */
#include <tommath.h>
const int sizes[] = { 8, 19, 28, 37, 55, 74, 110, 147 };
int main(void)
{
int res, x, y;
char buf[4096];
FILE *out;
mp_int a, b;
mp_init(&a);
mp_init(&b);
out = fopen("drprimes.txt", "w");
for (x = 0; x < (int)(sizeof(sizes)/sizeof(sizes[0])); x++) {
printf("Seeking a %d-bit safe prime\n", sizes[x] * DIGIT_BIT);
mp_grow(&a, sizes[x]);
mp_zero(&a);
for (y = 1; y < sizes[x]; y++) {
a.dp[y] = MP_MASK;
}
/* make a DR modulus */
a.dp[0] = 1;
a.used = sizes[x];
/* now loop */
do {
fflush(stdout);
mp_prime_next_prime(&a, 3);
printf(".");
mp_sub_d(&a, 1, &b);
mp_div_2(&b, &b);
mp_prime_is_prime(&b, 3, &res);
} while (res == 0);
if (mp_dr_is_modulus(&a) != 1) {
printf("Error not DR modulus\n");
} else {
mp_toradix(&a, buf, 10);
printf("\n\np == %s\n\n", buf);
fprintf(out, "%d-bit prime:\np == %s\n\n", mp_count_bits(&a), buf); fflush(out);
}
}
fclose(out);
mp_clear(&a);
mp_clear(&b);
return 0;
}

23
etc/drprimes.1 Normal file
View File

@ -0,0 +1,23 @@
224-bit prime:
p == 26959946667150639794667015087019630673637144422540572481103341844143
532-bit prime:
p == 14059105607947488696282932836518693308967803494693489478439861164411992439598399594747002144074658928593502845729752797260025831423419686528151609940203368691747
784-bit prime:
p == 101745825697019260773923519755878567461315282017759829107608914364075275235254395622580447400994175578963163918967182013639660669771108475957692810857098847138903161308502419410142185759152435680068435915159402496058513611411688900243039
1036-bit prime:
p == 736335108039604595805923406147184530889923370574768772191969612422073040099331944991573923112581267542507986451953227192970402893063850485730703075899286013451337291468249027691733891486704001513279827771740183629161065194874727962517148100775228363421083691764065477590823919364012917984605619526140821798437127
1540-bit prime:
p == 38564998830736521417281865696453025806593491967131023221754800625044118265468851210705360385717536794615180260494208076605798671660719333199513807806252394423283413430106003596332513246682903994829528690198205120921557533726473585751382193953592127439965050261476810842071573684505878854588706623484573925925903505747545471088867712185004135201289273405614415899438276535626346098904241020877974002916168099951885406379295536200413493190419727789712076165162175783
2072-bit prime:
p == 542189391331696172661670440619180536749994166415993334151601745392193484590296600979602378676624808129613777993466242203025054573692562689251250471628358318743978285860720148446448885701001277560572526947619392551574490839286458454994488665744991822837769918095117129546414124448777033941223565831420390846864429504774477949153794689948747680362212954278693335653935890352619041936727463717926744868338358149568368643403037768649616778526013610493696186055899318268339432671541328195724261329606699831016666359440874843103020666106568222401047720269951530296879490444224546654729111504346660859907296364097126834834235287147
3080-bit prime:
p == 1487259134814709264092032648525971038895865645148901180585340454985524155135260217788758027400478312256339496385275012465661575576202252063145698732079880294664220579764848767704076761853197216563262660046602703973050798218246170835962005598561669706844469447435461092542265792444947706769615695252256130901271870341005768912974433684521436211263358097522726462083917939091760026658925757076733484173202927141441492573799914240222628795405623953109131594523623353044898339481494120112723445689647986475279242446083151413667587008191682564376412347964146113898565886683139407005941383669325997475076910488086663256335689181157957571445067490187939553165903773554290260531009121879044170766615232300936675369451260747671432073394867530820527479172464106442450727640226503746586340279816318821395210726268291535648506190714616083163403189943334431056876038286530365757187367147446004855912033137386225053275419626102417236133948503
4116-bit prime:
p == 1095121115716677802856811290392395128588168592409109494900178008967955253005183831872715423151551999734857184538199864469605657805519106717529655044054833197687459782636297255219742994736751541815269727940751860670268774903340296040006114013971309257028332849679096824800250742691718610670812374272414086863715763724622797509437062518082383056050144624962776302147890521249477060215148275163688301275847155316042279405557632639366066847442861422164832655874655824221577849928863023018366835675399949740429332468186340518172487073360822220449055340582568461568645259954873303616953776393853174845132081121976327462740354930744487429617202585015510744298530101547706821590188733515880733527449780963163909830077616357506845523215289297624086914545378511082534229620116563260168494523906566709418166011112754529766183554579321224940951177394088465596712620076240067370589036924024728375076210477267488679008016579588696191194060127319035195370137160936882402244399699172017835144537488486396906144217720028992863941288217185353914991583400421682751000603596655790990815525126154394344641336397793791497068253936771017031980867706707490224041075826337383538651825493679503771934836094655802776331664261631740148281763487765852746577808019633679

5
etc/makefile
View File

@ -16,5 +16,8 @@ tune: tune.o
mersenne: mersenne.o mersenne: mersenne.o
$(CC) mersenne.o $(LIBNAME) -o mersenne $(CC) mersenne.o $(LIBNAME) -o mersenne
drprime: drprime.o
$(CC) drprime.o $(LIBNAME) -o drprime
clean: clean:
rm -f *.log *.o *.obj *.exe pprime tune mersenne rm -f *.log *.o *.obj *.exe pprime tune mersenne drprime

3
etc/makefile.msvc
View File

@ -12,3 +12,6 @@ mersenne: mersenne.obj
tune: tune.obj tune: tune.obj
cl tune.obj ../tommath.lib cl tune.obj ../tommath.lib
drprime: drprime.obj
cl drprime.obj ../tommath.lib

12
etc/tune.c
View File

@ -17,7 +17,7 @@ time_mult (void)
mp_init (&c); mp_init (&c);
t1 = clock (); t1 = clock ();
for (x = 4; x <= 128; x += 4) { for (x = 4; x <= 144; x += 4) {
mp_rand (&a, x); mp_rand (&a, x);
mp_rand (&b, x); mp_rand (&b, x);
for (y = 0; y < 10000; y++) { for (y = 0; y < 10000; y++) {
@ -41,7 +41,7 @@ time_sqr (void)
mp_init (&b); mp_init (&b);
t1 = clock (); t1 = clock ();
for (x = 4; x <= 128; x += 4) { for (x = 4; x <= 144; x += 4) {
mp_rand (&a, x); mp_rand (&a, x);
for (y = 0; y < 10000; y++) { for (y = 0; y < 10000; y++) {
mp_sqr (&a, &b); mp_sqr (&a, &b);
@ -65,7 +65,7 @@ time_expt (void)
mp_init (&d); mp_init (&d);
t1 = clock (); t1 = clock ();
for (x = 4; x <= 128; x += 4) { for (x = 4; x <= 144; x += 4) {
mp_rand (&a, x); mp_rand (&a, x);
mp_rand (&b, x); mp_rand (&b, x);
mp_rand (&c, x); mp_rand (&c, x);
@ -96,7 +96,7 @@ main (void)
/* tune multiplication first */ /* tune multiplication first */
log = fopen ("mult.log", "w"); log = fopen ("mult.log", "w");
best = CLOCKS_PER_SEC * 1000; best = CLOCKS_PER_SEC * 1000;
for (KARATSUBA_MUL_CUTOFF = 8; KARATSUBA_MUL_CUTOFF <= 128; KARATSUBA_MUL_CUTOFF++) { for (KARATSUBA_MUL_CUTOFF = 8; KARATSUBA_MUL_CUTOFF <= 144; KARATSUBA_MUL_CUTOFF++) {
ti = time_mult (); ti = time_mult ();
printf ("%4d : %9lu\r", KARATSUBA_MUL_CUTOFF, ti); printf ("%4d : %9lu\r", KARATSUBA_MUL_CUTOFF, ti);
fprintf (log, "%d, %lu\n", KARATSUBA_MUL_CUTOFF, ti); fprintf (log, "%d, %lu\n", KARATSUBA_MUL_CUTOFF, ti);
@ -112,7 +112,7 @@ main (void)
/* tune squaring */ /* tune squaring */
log = fopen ("sqr.log", "w"); log = fopen ("sqr.log", "w");
best = CLOCKS_PER_SEC * 1000; best = CLOCKS_PER_SEC * 1000;
for (KARATSUBA_SQR_CUTOFF = 8; KARATSUBA_SQR_CUTOFF <= 128; KARATSUBA_SQR_CUTOFF++) { for (KARATSUBA_SQR_CUTOFF = 8; KARATSUBA_SQR_CUTOFF <= 144; KARATSUBA_SQR_CUTOFF++) {
ti = time_sqr (); ti = time_sqr ();
printf ("%4d : %9lu\r", KARATSUBA_SQR_CUTOFF, ti); printf ("%4d : %9lu\r", KARATSUBA_SQR_CUTOFF, ti);
fprintf (log, "%d, %lu\n", KARATSUBA_SQR_CUTOFF, ti); fprintf (log, "%d, %lu\n", KARATSUBA_SQR_CUTOFF, ti);
@ -131,7 +131,7 @@ main (void)
log = fopen ("expt.log", "w"); log = fopen ("expt.log", "w");
best = CLOCKS_PER_SEC * 1000; best = CLOCKS_PER_SEC * 1000;
for (MONTGOMERY_EXPT_CUTOFF = 8; MONTGOMERY_EXPT_CUTOFF <= 192; MONTGOMERY_EXPT_CUTOFF++) { for (MONTGOMERY_EXPT_CUTOFF = 8; MONTGOMERY_EXPT_CUTOFF <= 144; MONTGOMERY_EXPT_CUTOFF++) {
ti = time_expt (); ti = time_expt ();
printf ("%4d : %9lu\r", MONTGOMERY_EXPT_CUTOFF, ti); printf ("%4d : %9lu\r", MONTGOMERY_EXPT_CUTOFF, ti);
fflush (stdout); fflush (stdout);

7
makefile
View File

@ -1,6 +1,6 @@
CFLAGS += -I./ -Wall -W -Wshadow -O3 -fomit-frame-pointer -funroll-loops CFLAGS += -I./ -Wall -W -Wshadow -O3 -fomit-frame-pointer -funroll-loops
VERSION=0.14 VERSION=0.15
default: libtommath.a default: libtommath.a
@ -30,7 +30,9 @@ bn_mp_reduce.o bn_mp_montgomery_setup.o bn_fast_mp_montgomery_reduce.o bn_mp_mon
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_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_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_radix.o \ bn_mp_to_signed_bin.o bn_mp_unsigned_bin_size.o bn_mp_signed_bin_size.o bn_radix.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_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
libtommath.a: $(OBJECTS) libtommath.a: $(OBJECTS)
$(AR) $(ARFLAGS) libtommath.a $(OBJECTS) $(AR) $(ARFLAGS) libtommath.a $(OBJECTS)
@ -65,6 +67,7 @@ clean:
cd etc ; make clean cd etc ; make clean
zipup: clean docs zipup: clean docs
perl gen.pl ; mv mpi.c pre_gen/ ; \
cd .. ; rm -rf ltm* libtommath-$(VERSION) ; mkdir libtommath-$(VERSION) ; \ cd .. ; rm -rf ltm* libtommath-$(VERSION) ; mkdir libtommath-$(VERSION) ; \
cp -R ./libtommath/* ./libtommath-$(VERSION)/ ; tar -c libtommath-$(VERSION)/* > ltm-$(VERSION).tar ; \ cp -R ./libtommath/* ./libtommath-$(VERSION)/ ; tar -c libtommath-$(VERSION)/* > ltm-$(VERSION).tar ; \
bzip2 -9vv ltm-$(VERSION).tar ; zip -9 -r ltm-$(VERSION).zip libtommath-$(VERSION)/* bzip2 -9vv ltm-$(VERSION).tar ; zip -9 -r ltm-$(VERSION).zip libtommath-$(VERSION)/*

5
makefile.msvc
View File

@ -20,7 +20,10 @@ bn_mp_reduce.obj bn_mp_montgomery_setup.obj bn_fast_mp_montgomery_reduce.obj bn_
bn_mp_exptmod_fast.obj bn_mp_exptmod.obj bn_mp_2expt.obj bn_mp_n_root.obj bn_mp_jacobi.obj bn_reverse.obj \ bn_mp_exptmod_fast.obj bn_mp_exptmod.obj bn_mp_2expt.obj bn_mp_n_root.obj bn_mp_jacobi.obj bn_reverse.obj \
bn_mp_count_bits.obj bn_mp_read_unsigned_bin.obj bn_mp_read_signed_bin.obj bn_mp_to_unsigned_bin.obj \ bn_mp_count_bits.obj bn_mp_read_unsigned_bin.obj bn_mp_read_signed_bin.obj bn_mp_to_unsigned_bin.obj \
bn_mp_to_signed_bin.obj bn_mp_unsigned_bin_size.obj bn_mp_signed_bin_size.obj bn_radix.obj \ bn_mp_to_signed_bin.obj bn_mp_unsigned_bin_size.obj bn_mp_signed_bin_size.obj bn_radix.obj \
bn_mp_xor.obj bn_mp_and.obj bn_mp_or.obj bn_mp_rand.obj bn_mp_montgomery_calc_normalization.obj bn_mp_xor.obj bn_mp_and.obj bn_mp_or.obj bn_mp_rand.obj bn_mp_montgomery_calc_normalization.obj \
bn_mp_prime_is_divisible.obj bn_prime_tab.obj bn_mp_prime_fermat.obj bn_mp_prime_miller_rabin.obj \
bn_mp_prime_is_prime.obj bn_mp_prime_next_prime.obj bn_mp_dr_reduce.obj
library: $(OBJECTS) library: $(OBJECTS)
lib /out:tommath.lib $(OBJECTS) lib /out:tommath.lib $(OBJECTS)

6
mtest/mtest.c
View File

@ -41,7 +41,7 @@ void rand_num(mp_int *a)
unsigned char buf[512]; unsigned char buf[512];
top: top:
size = 1 + ((fgetc(rng)*fgetc(rng)) % 512); size = 1 + ((fgetc(rng)*fgetc(rng)) % 96);
buf[0] = (fgetc(rng)&1)?1:0; buf[0] = (fgetc(rng)&1)?1:0;
fread(buf+1, 1, size, rng); fread(buf+1, 1, size, rng);
for (n = 0; n < size; n++) { for (n = 0; n < size; n++) {
@ -57,7 +57,7 @@ void rand_num2(mp_int *a)
unsigned char buf[512]; unsigned char buf[512];
top: top:
size = 1 + ((fgetc(rng)*fgetc(rng)) % 512); size = 1 + ((fgetc(rng)*fgetc(rng)) % 96);
buf[0] = (fgetc(rng)&1)?1:0; buf[0] = (fgetc(rng)&1)?1:0;
fread(buf+1, 1, size, rng); fread(buf+1, 1, size, rng);
for (n = 0; n < size; n++) { for (n = 0; n < size; n++) {
@ -73,8 +73,6 @@ int main(void)
mp_int a, b, c, d, e; mp_int a, b, c, d, e;
char buf[4096]; char buf[4096];
static int tests[] = { 11, 12 };
mp_init(&a); mp_init(&a);
mp_init(&b); mp_init(&b);
mp_init(&c); mp_init(&c);

5993
pre_gen/mpi.c Normal file
View File

File diff suppressed because it is too large Load Diff

56
tommath.h
View File

@ -28,8 +28,16 @@
#ifdef __cplusplus #ifdef __cplusplus
extern "C" { extern "C" {
#endif
/* C++ compilers don't like assigning void * to mp_digit * */
#define OPT_CAST (mp_digit *)
#else
/* C on the other hand dosen't care */
#define OPT_CAST
#endif
/* some default configurations. /* some default configurations.
* *
@ -202,7 +210,6 @@ int mp_cmp_mag(mp_int *a, mp_int *b);
/* c = a + b */ /* c = a + b */
int mp_add(mp_int *a, mp_int *b, mp_int *c); int mp_add(mp_int *a, mp_int *b, mp_int *c);
/* c = a - b */ /* c = a - b */
int mp_sub(mp_int *a, mp_int *b, mp_int *c); int mp_sub(mp_int *a, mp_int *b, mp_int *c);
@ -297,9 +304,52 @@ int mp_montgomery_calc_normalization(mp_int *a, mp_int *b);
/* computes xR^-1 == x (mod N) via Montgomery Reduction */ /* computes xR^-1 == x (mod N) via Montgomery Reduction */
int mp_montgomery_reduce(mp_int *a, mp_int *m, mp_digit mp); int mp_montgomery_reduce(mp_int *a, mp_int *m, mp_digit mp);
/* returns 1 if a is a valid DR modulus */
int mp_dr_is_modulus(mp_int *a);
/* sets the value of "d" required for mp_dr_reduce */
void mp_dr_setup(mp_int *a, mp_digit *d);
/* reduces a modulo b using the Diminished Radix method */
int mp_dr_reduce(mp_int *a, mp_int *b, mp_digit mp);
/* d = a^b (mod c) */ /* d = a^b (mod c) */
int mp_exptmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d); int mp_exptmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
/* ---> Primes <--- */
#define PRIME_SIZE 256 /* number of primes */
/* table of first 256 primes */
extern const mp_digit __prime_tab[];
/* result=1 if a is divisible by one of the first 256 primes */
int mp_prime_is_divisible(mp_int *a, int *result);
/* performs one Fermat test of "a" using base "b".
* Sets result to 0 if composite or 1 if probable prime
*/
int mp_prime_fermat(mp_int *a, mp_int *b, int *result);
/* performs one Miller-Rabin test of "a" using base "b".
* Sets result to 0 if composite or 1 if probable prime
*/
int mp_prime_miller_rabin(mp_int *a, mp_int *b, int *result);
/* 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
* of error no more than (1/4)^t.
*
* Sets result to 1 if probably prime, 0 otherwise
*/
int mp_prime_is_prime(mp_int *a, int t, int *result);
/* finds the next prime after the number "a" using "t" trials
* of Miller-Rabin.
*/
int mp_prime_next_prime(mp_int *a, int t);
/* ---> radix conversion <--- */ /* ---> radix conversion <--- */
int mp_count_bits(mp_int *a); int mp_count_bits(mp_int *a);
@ -341,7 +391,7 @@ int mp_karatsuba_mul(mp_int *a, mp_int *b, mp_int *c);
int mp_karatsuba_sqr(mp_int *a, mp_int *b); int mp_karatsuba_sqr(mp_int *a, mp_int *b);
int fast_mp_invmod(mp_int *a, mp_int *b, mp_int *c); int fast_mp_invmod(mp_int *a, mp_int *b, mp_int *c);
int fast_mp_montgomery_reduce(mp_int *a, mp_int *m, mp_digit mp); int fast_mp_montgomery_reduce(mp_int *a, mp_int *m, mp_digit mp);
int mp_exptmod_fast(mp_int *G, mp_int *X, mp_int *P, mp_int *Y); int mp_exptmod_fast(mp_int *G, mp_int *X, mp_int *P, mp_int *Y, int mode);
void bn_reverse(unsigned char *s, int len); void bn_reverse(unsigned char *s, int len);
#ifdef __cplusplus #ifdef __cplusplus