# moritz/libtommath forked from libtom/libtommath

1 parent 82f4858 commit b1756f2f9872a2981a6fea650d40fe5eedb561e7 Tom St Denis committed with sjaeckel Mar 22, 2003
BIN bn.pdf
Binary file not shown.
376 bn.tex
 @@ -1,7 +1,7 @@ \documentclass{article} \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} \maketitle \newpage @@ -100,6 +100,22 @@ \subsection{Calling Functions} mp_div_2(&x, &y); /* y = x / 2 */ \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} 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. @@ -326,10 +342,53 @@ \subsection{Modular Arithmetic} /* computes xR^-1 == x (mod N) via Montgomery Reduction */ 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) */ int mp_exptmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d); \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} To read or store integers in other formats there are the following functions. @@ -533,23 +592,131 @@ \subsubsection{mp\_jacobi(mp\_int *a, mp\_int *n, int *c)} it is not. \subsubsection{mp\_exptmod(mp\_int *a, mp\_int *b, mp\_int *c, mp\_int *d)} -Computes $d = a^b \mbox{ (mod }c\mbox{)}$ using a sliding window $k$-ary exponentiation algorithm. For an $\alpha$-bit +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 -chosen to minimize the number of multiplications required for a given value of $\alpha$. Barrett or Montgomery -reductions are used to reduce the squared or multiplied temporary results modulo $c$. +chosen to minimize the number of multiplications required for a given value of $\alpha$. Barrett, Montgomery or +Dimminished-Radix reductions are used to reduce the squared or multiplied temporary results modulo $c$. \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)} 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. +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 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 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)} 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$. @@ -576,7 +743,95 @@ \subsubsection{mp\_montgomery\_reduce(mp\_int *a, mp\_int *m, mp\_digit mp)} two long divisions would be required to setup $\hat x$ and a multiplication followed by reduction for each variable. 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} @@ -662,21 +917,116 @@ \subsection{Exponentiation Algorithms} 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 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. +and 512 multiplications. + +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. However, LibTomMath can take advantage of the fact that the multiplications required within the Barrett reduction 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 comba baseline). -LibTomMath also has a variant of the exptmod function that uses Montgomery reductions instead of Barrett reductions -which is faster. The code will automatically detect when the Montgomery version can be used (\textit{Requires the -modulus to be odd and below the MONTGOMERY\_EXPT\_CUTOFF size}). The Montgomery routine is essentially a copy of the -Barrett exponentiation routine except it uses Montgomery reduction. +LibTomMath also has a variant of the exptmod function that uses Montgomery or Diminished-Radix reductions instead of +Barrett reductions which are faster. The code will automatically detect when the Montgomery version can be used +(\textit{Requires the modulus to be odd and below the MONTGOMERY\_EXPT\_CUTOFF size}). The Montgomery routine is +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 -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}
1 bn_fast_mp_invmod.c
 @@ -80,7 +80,6 @@ fast_mp_invmod (mp_int * a, mp_int * b, mp_int * c) } mp_set (&D, 1); - top: /* 4. while u is even do */ while (mp_iseven (&u) == 1) {
6 bn_mp_div.c
 @@ -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) */ for (i = n; i >= (t + 1); i--) { - if (i > x.alloc) + if (i > x.used) 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 */ @@ -171,10 +171,11 @@ mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d) q.dp[i - t - 1] = (q.dp[i - t - 1] - 1UL) & MP_MASK; } } - + /* now q is the quotient and x is the remainder [which we have to normalize] */ /* get sign before writing to c */ x.sign = a->sign; + if (c != NULL) { mp_clamp (&q); 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) { mp_div_2d (&x, norm, &x, NULL); - mp_clamp (&x); mp_exch (&x, d); }
4 bn_mp_div_2d.c
 @@ -52,8 +52,8 @@ mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d) /* shift by as many digits in the bit count */ if (b >= DIGIT_BIT) { - mp_rshd (c, b / DIGIT_BIT); - } + mp_rshd (c, b / DIGIT_BIT); + } /* shift any bit count < DIGIT_BIT */ D = (mp_digit) (b % DIGIT_BIT);
1 bn_mp_div_d.c
 @@ -21,7 +21,6 @@ mp_div_d (mp_int * a, mp_digit b, mp_int * c, mp_digit * d) mp_int t, t2; int res; - if ((res = mp_init (&t)) != MP_OKAY) { return res; }
150 bn_mp_dr_reduce.c
7 bn_mp_exptmod.c
 @@ -24,9 +24,12 @@ static int f_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y); int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y) { + int dr; + + dr = mp_dr_is_modulus(P); /* if the modulus is odd use the fast method */ - if (mp_isodd (P) == 1 && P->used > 4 && P->used < MONTGOMERY_EXPT_CUTOFF) { - return mp_exptmod_fast (G, X, P, Y); + if (((mp_isodd (P) == 1 && P->used < MONTGOMERY_EXPT_CUTOFF) || dr == 1) && P->used > 4) { + return mp_exptmod_fast (G, X, P, Y, dr); } else { return f_mp_exptmod (G, X, P, Y); }
65 bn_mp_exptmod_fast.c
 @@ -22,11 +22,13 @@ * Uses Montgomery reduction */ 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_digit buf, mp; int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize; + int (*redux)(mp_int*,mp_int*,mp_digit); + /* find window size */ x = mp_count_bits (X); @@ -55,10 +57,17 @@ mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y) return err; } } - - /* now setup montgomery */ - if ((err = mp_montgomery_setup (P, &mp)) != MP_OKAY) { - goto __M; + + if (redmode == 0) { + /* now setup montgomery */ + 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 */ @@ -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] */ - /* now we need R mod m */ - if ((err = mp_montgomery_calc_normalization (&res, P)) != MP_OKAY) { - goto __RES; - } + if (redmode == 0) { + /* now we need R mod m */ + 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; + /* 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; + } } + /* 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) { 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) { 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; } } @@ -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) { goto __RES; } - if ((err = mp_montgomery_reduce (&M[x], P, mp)) != MP_OKAY) { + if ((err = redux (&M[x], P, mp)) != MP_OKAY) { 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) { goto __RES; } - if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) { + if ((err = redux (&res, P, mp)) != MP_OKAY) { goto __RES; } 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) { goto __RES; } - if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) { + if ((err = redux (&res, P, mp)) != MP_OKAY) { 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) { goto __RES; } - if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) { + if ((err = redux (&res, P, mp)) != MP_OKAY) { 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) { goto __RES; } - if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) { + if ((err = redux (&res, P, mp)) != MP_OKAY) { goto __RES; } @@ -194,17 +211,19 @@ 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) { goto __RES; } - if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) { + if ((err = redux (&res, P, mp)) != MP_OKAY) { goto __RES; } } } } - /* fixup result */ - if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) { - goto __RES; - } + if (redmode == 0) { + /* fixup result */ + if ((err = mp_montgomery_reduce (&res, P, mp)) != MP_OKAY) { + goto __RES; + } + } mp_exch (&res, Y); err = MP_OKAY;
2 bn_mp_grow.c
 @@ -24,7 +24,7 @@ mp_grow (mp_int * a, int 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 */ - a->dp = realloc (a->dp, sizeof (mp_digit) * size); + a->dp = OPT_CAST realloc (a->dp, sizeof (mp_digit) * size); if (a->dp == NULL) { return MP_MEM; }
2 bn_mp_init.c
 @@ -20,7 +20,7 @@ mp_init (mp_int * a) { /* 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) { return MP_MEM; }
2 bn_mp_init_size.c
 @@ -21,7 +21,7 @@ mp_init_size (mp_int * a, int size) /* 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 */ - a->dp = calloc (sizeof (mp_digit), size); + a->dp = OPT_CAST calloc (sizeof (mp_digit), size); if (a->dp == NULL) { return MP_MEM; }
4 bn_mp_lshd.c
 @@ -36,10 +36,10 @@ mp_lshd (mp_int * a, int b) /* increment the used by the shift amount than copy upwards */ a->used += b; - + /* top */ tmpa = a->dp + a->used - 1; - + /* base */ tmpaa = a->dp + a->used - 1 - b;
8 bn_mp_mul_2d.c
 @@ -33,10 +33,10 @@ mp_mul_2d (mp_int * a, int b, mp_int * c) /* shift by as many digits in the bit count */ if (b >= DIGIT_BIT) { - if ((res = mp_lshd (c, b / DIGIT_BIT)) != MP_OKAY) { - return res; - } - } + if ((res = mp_lshd (c, b / DIGIT_BIT)) != MP_OKAY) { + return res; + } + } c->used = c->alloc; /* shift any bit count < DIGIT_BIT */
52 bn_mp_prime_fermat.c
 @@ -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 + +/* 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
 @@ -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 + +/* 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
 @@ -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 + +/* 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
 @@ -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 + +/* 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
 @@ -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 + +/* 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; +} +
8 bn_mp_rshd.c
 @@ -38,19 +38,19 @@ mp_rshd (mp_int * a, int b) /* base */ tmpa = a->dp; - + /* offset into digits */ tmpaa = a->dp + b; - + /* this is implemented as a sliding window where the window is b-digits long * and digits from the top of the window are copied to the bottom * * e.g. - + b-2 | b-1 | b0 | b1 | b2 | ... | bb | ----> /\ | ----> \-------------------/ ----> - */ + */ for (x = 0; x < (a->used - b); x++) { *tmpa++ = *tmpaa++; }
2 bn_mp_shrink.c
 @@ -19,7 +19,7 @@ int mp_shrink (mp_int * a) { 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; } a->alloc = a->used;
52 bn_prime_tab.c
 @@ -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 +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 +};
 @@ -93,7 +93,7 @@ mp_toradix (mp_int * a, char *str, int radix) *str++ = s_rmap[d]; ++digs; } - bn_reverse ((unsigned char *) _s, digs); + bn_reverse ((unsigned char *)_s, digs); *str++ = '\0'; mp_clear (&t); return MP_OKAY;
 @@ -55,13 +55,13 @@ s_mp_add (mp_int * a, mp_int * b, mp_int * c) register int i; /* alias for digit pointers */ - + /* first input */ tmpa = a->dp; - + /* second input */ tmpb = b->dp; - + /* destination */ tmpc = c->dp;
2 bncore.c
 @@ -18,5 +18,3 @@ 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. */ MONTGOMERY_EXPT_CUTOFF = 128; /* max. number of digits that montgomery reductions will help for */ - -
11 changes.txt
 @@ -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 v0.14 -- Tons of manual updates -- cleaned up the directory
54 demo/demo.c
 @@ -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, div2_n, mul2_n; unsigned rr; - int cnt; + int cnt, ix; #ifdef TIMER int n; @@ -103,10 +103,43 @@ int main(void) mp_init(&d); mp_init(&e); 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<
53 etc/drprime.c
 @@ -0,0 +1,53 @@ +/* Makes safe primes of a DR nature */ +#include + +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
 @@ -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
 @@ -15,6 +15,9 @@ tune: tune.o mersenne: mersenne.o $(CC) mersenne.o$(LIBNAME) -o mersenne + +drprime: drprime.o + $(CC) drprime.o$(LIBNAME) -o drprime clean: - rm -f *.log *.o *.obj *.exe pprime tune mersenne + rm -f *.log *.o *.obj *.exe pprime tune mersenne drprime
5 etc/makefile.msvc
 @@ -11,4 +11,7 @@ mersenne: mersenne.obj cl mersenne.obj ../tommath.lib 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
 @@ -17,7 +17,7 @@ time_mult (void) mp_init (&c); t1 = clock (); - for (x = 4; x <= 128; x += 4) { + for (x = 4; x <= 144; x += 4) { mp_rand (&a, x); mp_rand (&b, x); for (y = 0; y < 10000; y++) { @@ -41,7 +41,7 @@ time_sqr (void) mp_init (&b); t1 = clock (); - for (x = 4; x <= 128; x += 4) { + for (x = 4; x <= 144; x += 4) { mp_rand (&a, x); for (y = 0; y < 10000; y++) { mp_sqr (&a, &b); @@ -65,7 +65,7 @@ time_expt (void) mp_init (&d); t1 = clock (); - for (x = 4; x <= 128; x += 4) { + for (x = 4; x <= 144; x += 4) { mp_rand (&a, x); mp_rand (&b, x); mp_rand (&c, x); @@ -96,7 +96,7 @@ main (void) /* tune multiplication first */ log = fopen ("mult.log", "w"); 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 (); printf ("%4d : %9lu\r", KARATSUBA_MUL_CUTOFF, ti); fprintf (log, "%d, %lu\n", KARATSUBA_MUL_CUTOFF, ti); @@ -112,7 +112,7 @@ main (void) /* tune squaring */ log = fopen ("sqr.log", "w"); 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 (); printf ("%4d : %9lu\r", KARATSUBA_SQR_CUTOFF, ti); fprintf (log, "%d, %lu\n", KARATSUBA_SQR_CUTOFF, ti); @@ -131,7 +131,7 @@ main (void) log = fopen ("expt.log", "w"); 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 (); printf ("%4d : %9lu\r", MONTGOMERY_EXPT_CUTOFF, ti); fflush (stdout);
7 makefile
 @@ -1,6 +1,6 @@ CFLAGS += -I./ -Wall -W -Wshadow -O3 -fomit-frame-pointer -funroll-loops -VERSION=0.14 +VERSION=0.15 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_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_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)$(AR) $(ARFLAGS) libtommath.a$(OBJECTS) @@ -65,6 +67,7 @@ clean: cd etc ; make clean zipup: clean docs + perl gen.pl ; mv mpi.c pre_gen/ ; \ cd .. ; rm -rf ltm* libtommath-$(VERSION) ; mkdir libtommath-$(VERSION) ; \ 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)/*
5 makefile.msvc
 @@ -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_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_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) lib /out:tommath.lib$(OBJECTS)
6 mtest/mtest.c
 @@ -41,7 +41,7 @@ void rand_num(mp_int *a) unsigned char buf[512]; top: - size = 1 + ((fgetc(rng)*fgetc(rng)) % 512); + size = 1 + ((fgetc(rng)*fgetc(rng)) % 96); buf[0] = (fgetc(rng)&1)?1:0; fread(buf+1, 1, size, rng); for (n = 0; n < size; n++) { @@ -57,7 +57,7 @@ void rand_num2(mp_int *a) unsigned char buf[512]; top: - size = 1 + ((fgetc(rng)*fgetc(rng)) % 512); + size = 1 + ((fgetc(rng)*fgetc(rng)) % 96); buf[0] = (fgetc(rng)&1)?1:0; fread(buf+1, 1, size, rng); for (n = 0; n < size; n++) { @@ -73,8 +73,6 @@ int main(void) mp_int a, b, c, d, e; char buf[4096]; - static int tests[] = { 11, 12 }; - mp_init(&a); mp_init(&b); mp_init(&c);
5,993 pre_gen/mpi.c
5,993 additions, 0 deletions not shown because the diff is too large. Please use a local Git client to view these changes.
56 tommath.h
 @@ -28,8 +28,16 @@ #ifdef __cplusplus 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. * @@ -202,7 +210,6 @@ 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); - /* c = a - b */ 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 */ 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) */ 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 <--- */ 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 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 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); #ifdef __cplusplus