/*
* Use is subject to license terms.
*
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*
* Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
* or visit www.oracle.com if you need additional information or have any
* questions.
*/
/* *********************************************************************
*
* The Original Code is the MPI Arbitrary Precision Integer Arithmetic library.
*
* The Initial Developer of the Original Code is
* Michael J. Fromberger.
* Portions created by the Initial Developer are Copyright (C) 1998
* the Initial Developer. All Rights Reserved.
*
* Contributor(s):
* Netscape Communications Corporation
* Douglas Stebila <douglas@stebila.ca> of Sun Laboratories.
*
*********************************************************************** */
/* Arbitrary precision integer arithmetic library */
#include "mpi-priv.h"
#if defined(OSF1)
#include <c_asm.h>
#endif
#if MP_LOGTAB
/*
A table of the logs of 2 for various bases (the 0 and 1 entries of
this table are meaningless and should not be referenced).
This table is used to compute output lengths for the mp_toradix()
function. Since a number n in radix r takes up about log_r(n)
digits, we estimate the output size by taking the least integer
greater than log_r(n), where:
log_r(n) = log_2(n) * log_r(2)
This table, therefore, is a table of log_r(2) for 2 <= r <= 36,
which are the output bases supported.
*/
#include "logtab.h"
#endif
/* {{{ Constant strings */
/* Constant strings returned by mp_strerror() */
static const char *mp_err_string[] = {
"unknown result code", /* say what? */
"boolean true", /* MP_OKAY, MP_YES */
"boolean false", /* MP_NO */
"out of memory", /* MP_MEM */
"argument out of range", /* MP_RANGE */
"invalid input parameter", /* MP_BADARG */
"result is undefined" /* MP_UNDEF */
};
/* Value to digit maps for radix conversion */
/* s_dmap_1 - standard digits and letters */
static const char *s_dmap_1 =
"0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";
/* }}} */
unsigned long mp_allocs;
unsigned long mp_frees;
unsigned long mp_copies;
/* {{{ Default precision manipulation */
/* Default precision for newly created mp_int's */
{
return s_mp_defprec;
} /* end mp_get_prec() */
{
if(prec == 0)
else
s_mp_defprec = prec;
} /* end mp_set_prec() */
/* }}} */
/*------------------------------------------------------------------------*/
/* {{{ mp_init(mp, kmflag) */
/*
mp_init(mp, kmflag)
Initialize a new zero-valued mp_int. Returns MP_OKAY if successful,
MP_MEM if memory could not be allocated for the structure.
*/
{
} /* end mp_init() */
/* }}} */
/* {{{ mp_init_size(mp, prec, kmflag) */
/*
mp_init_size(mp, prec, kmflag)
Initialize a new zero-valued mp_int with at least the given
precision; returns MP_OKAY if successful, or MP_MEM if memory could
not be allocated for the structure.
*/
{
return MP_MEM;
return MP_OKAY;
} /* end mp_init_size() */
/* }}} */
/* {{{ mp_init_copy(mp, from) */
/*
mp_init_copy(mp, from)
Initialize mp as an exact copy of from. Returns MP_OKAY if
successful, MP_MEM if memory could not be allocated for the new
structure.
*/
{
return MP_OKAY;
return MP_MEM;
#ifndef _WIN32
#endif /* _WIN32 */
return MP_OKAY;
} /* end mp_init_copy() */
/* }}} */
/* {{{ mp_copy(from, to) */
/*
mp_copy(from, to)
Copies the mp_int 'from' to the mp_int 'to'. It is presumed that
'to' has already been initialized (if not, use mp_init_copy()
instead). If 'from' and 'to' are identical, nothing happens.
*/
{
return MP_OKAY;
++mp_copies;
{ /* copy */
/*
If the allocated buffer in 'to' already has enough space to hold
all the used digits of 'from', we'll re-use it to avoid hitting
the memory allocater more than necessary; otherwise, we'd have
to grow anyway, so we just allocate a hunk and make the copy as
usual
*/
} else {
return MP_MEM;
#if MP_CRYPTO
#endif
}
}
/* Copy the precision and sign from the original */
} /* end copy */
return MP_OKAY;
} /* end mp_copy() */
/* }}} */
/* {{{ mp_exch(mp1, mp2) */
/*
mp_exch(mp1, mp2)
Exchange mp1 and mp2 without allocating any intermediate memory
(well, unless you count the stack space needed for this call and the
locals it creates...). This cannot fail.
*/
{
#if MP_ARGCHK == 2
#else
return;
#endif
} /* end mp_exch() */
/* }}} */
/* {{{ mp_clear(mp) */
/*
mp_clear(mp)
Release the storage used by an mp_int, and void its fields so that
if someone calls mp_clear() again for the same int later, we won't
get tollchocked.
*/
{
return;
#if MP_CRYPTO
#endif
}
} /* end mp_clear() */
/* }}} */
/* {{{ mp_zero(mp) */
/*
mp_zero(mp)
Set mp to zero. Does not change the allocated size of the structure,
and therefore cannot fail (except on a bad argument, which we ignore)
*/
{
return;
} /* end mp_zero() */
/* }}} */
/* {{{ mp_set(mp, d) */
{
return;
} /* end mp_set() */
/* }}} */
/* {{{ mp_set_int(mp, z) */
{
int ix;
unsigned long v = labs(z);
if(z == 0)
return MP_OKAY; /* shortcut for zero */
if (sizeof v <= sizeof(mp_digit)) {
} else {
return res;
return res;
}
}
if(z < 0)
return MP_OKAY;
} /* end mp_set_int() */
/* }}} */
/* {{{ mp_set_ulong(mp, z) */
{
int ix;
if(z == 0)
return MP_OKAY; /* shortcut for zero */
if (sizeof z <= sizeof(mp_digit)) {
} else {
return res;
return res;
}
}
return MP_OKAY;
} /* end mp_set_ulong() */
/* }}} */
/*------------------------------------------------------------------------*/
/* {{{ Digit arithmetic */
/* {{{ mp_add_d(a, d, b) */
/*
mp_add_d(a, d, b)
Compute the sum b = a + d, for a single digit d. Respects the sign of
its primary addend (single digits are unsigned anyway).
*/
{
return res;
goto CLEANUP;
} else if(s_mp_cmp_d(&tmp, d) >= 0) {
goto CLEANUP;
} else {
}
if(s_mp_cmp_d(&tmp, 0) == 0)
return res;
} /* end mp_add_d() */
/* }}} */
/* {{{ mp_sub_d(a, d, b) */
/*
mp_sub_d(a, d, b)
Compute the difference b = a - d, for a single digit d. Respects the
sign of its subtrahend (single digits are unsigned anyway).
*/
{
return res;
goto CLEANUP;
} else if(s_mp_cmp_d(&tmp, d) >= 0) {
goto CLEANUP;
} else {
}
if(s_mp_cmp_d(&tmp, 0) == 0)
return res;
} /* end mp_sub_d() */
/* }}} */
/* {{{ mp_mul_d(a, d, b) */
/*
mp_mul_d(a, d, b)
Compute the product b = a * d, for a single digit d. Respects the sign
of its multiplicand (single digits are unsigned anyway)
*/
{
if(d == 0) {
mp_zero(b);
return MP_OKAY;
}
return res;
res = s_mp_mul_d(b, d);
return res;
} /* end mp_mul_d() */
/* }}} */
/* {{{ mp_mul_2(a, c) */
{
return res;
return s_mp_mul_2(c);
} /* end mp_mul_2() */
/* }}} */
/* {{{ mp_div_d(a, d, q, r) */
/*
mp_div_d(a, d, q, r)
Compute the quotient q = a / d and remainder r = a mod d, for a
single digit d. Respects the sign of its divisor (single digits are
unsigned anyway).
*/
{
int pow;
if(d == 0)
return MP_RANGE;
/* Shortcut for powers of two ... */
if((pow = s_mp_ispow2d(d)) >= 0) {
if(q) {
mp_copy(a, q);
s_mp_div_2d(q, pow);
}
if(r)
*r = rem;
return MP_OKAY;
}
return res;
if(s_mp_cmp_d(&qp, 0) == 0)
if(r)
*r = rem;
if(q)
return res;
} /* end mp_div_d() */
/* }}} */
/* {{{ mp_div_2(a, c) */
/*
mp_div_2(a, c)
Compute c = a / 2, disregarding the remainder.
*/
{
return res;
s_mp_div_2(c);
return MP_OKAY;
} /* end mp_div_2() */
/* }}} */
/* {{{ mp_expt_d(a, d, b) */
{
mp_int s, x;
return res;
goto X;
DIGIT(&s, 0) = 1;
while(d != 0) {
if(d & 1) {
goto CLEANUP;
}
d /= 2;
goto CLEANUP;
}
s_mp_exch(&s, c);
mp_clear(&x);
X:
mp_clear(&s);
return res;
} /* end mp_expt_d() */
/* }}} */
/* }}} */
/*------------------------------------------------------------------------*/
/* {{{ Full arithmetic */
/* {{{ mp_abs(a, b) */
/*
mp_abs(a, b)
Compute b = |a|. 'a' and 'b' may be identical.
*/
{
return res;
return MP_OKAY;
} /* end mp_abs() */
/* }}} */
/* {{{ mp_neg(a, b) */
/*
mp_neg(a, b)
Compute b = -a. 'a' and 'b' may be identical.
*/
{
return res;
if(s_mp_cmp_d(b, 0) == MP_EQ)
else
return MP_OKAY;
} /* end mp_neg() */
/* }}} */
/* {{{ mp_add(a, b, c) */
/*
mp_add(a, b, c)
Compute c = a + b. All parameters may be identical.
*/
{
MP_CHECKOK( s_mp_add_3arg(a, b, c) );
} else if(s_mp_cmp(a, b) >= 0) { /* different sign: |a| >= |b| */
MP_CHECKOK( s_mp_sub_3arg(a, b, c) );
} else { /* different sign: |a| < |b| */
MP_CHECKOK( s_mp_sub_3arg(b, a, c) );
}
if (s_mp_cmp_d(c, 0) == MP_EQ)
return res;
} /* end mp_add() */
/* }}} */
/* {{{ mp_sub(a, b, c) */
/*
mp_sub(a, b, c)
Compute c = a - b. All parameters may be identical.
*/
{
int magDiff;
if (a == b) {
mp_zero(c);
return MP_OKAY;
}
MP_CHECKOK( s_mp_add_3arg(a, b, c) );
mp_zero(c);
} else if (magDiff > 0) {
MP_CHECKOK( s_mp_sub_3arg(a, b, c) );
} else {
MP_CHECKOK( s_mp_sub_3arg(b, a, c) );
}
if (s_mp_cmp_d(c, 0) == MP_EQ)
return res;
} /* end mp_sub() */
/* }}} */
/* {{{ mp_mul(a, b, c) */
/*
mp_mul(a, b, c)
Compute c = a * b. All parameters may be identical.
*/
{
if (a == c) {
return res;
if (a == b)
b = &tmp;
a = &tmp;
} else if (b == c) {
return res;
b = &tmp;
} else {
}
b = a;
a = xch;
}
goto CLEANUP;
#ifdef NSS_USE_COMBA
if (MP_USED(a) == 4) {
s_mp_mul_comba_4(a, b, c);
goto CLEANUP;
}
if (MP_USED(a) == 8) {
s_mp_mul_comba_8(a, b, c);
goto CLEANUP;
}
if (MP_USED(a) == 16) {
s_mp_mul_comba_16(a, b, c);
goto CLEANUP;
}
if (MP_USED(a) == 32) {
s_mp_mul_comba_32(a, b, c);
goto CLEANUP;
}
}
#endif
/* Outer loop: Digits of b */
/* Inner product: Digits of a */
if (b_i)
else
}
s_mp_clamp(c);
else
return res;
} /* end mp_mul() */
/* }}} */
/* {{{ mp_sqr(a, sqr) */
#if MP_SQUARE
/*
Computes the square of a. This can be done more
efficiently than a general multiplication, because many of the
computation steps are redundant when squaring. The inner product
step is a bit more complicated, but we save a fair number of
iterations of the multiplication loop.
*/
/* sqr = a^2; Caller provides both a and tmp; */
{
mp_digit d;
int count;
if (a == sqr) {
return res;
a = &tmp;
} else {
}
}
#ifdef NSS_USE_COMBA
if (IS_POWER_OF_2(MP_USED(a))) {
if (MP_USED(a) == 4) {
s_mp_sqr_comba_4(a, sqr);
goto CLEANUP;
}
if (MP_USED(a) == 8) {
s_mp_sqr_comba_8(a, sqr);
goto CLEANUP;
}
if (MP_USED(a) == 16) {
s_mp_sqr_comba_16(a, sqr);
goto CLEANUP;
}
if (MP_USED(a) == 32) {
s_mp_sqr_comba_32(a, sqr);
goto CLEANUP;
}
}
#endif
if (count > 0) {
d = *pa++;
d = *pa++;
} /* for(ix ...) */
/* now sqr *= 2 */
} else {
}
/* now add the squares of the digits of a to sqr. */
return res;
} /* end mp_sqr() */
#endif
/* }}} */
/* {{{ mp_div(a, b, q, r) */
/*
mp_div(a, b, q, r)
Compute q = a / b and r = a mod b. Input parameters may be re-used
as output parameters. If q or r is NULL, that portion of the
computation will be discarded (although it will still be computed)
*/
{
int cmp;
return MP_RANGE;
/* Set up some temporaries... */
if (!r || r == a || r == b) {
} else {
MP_CHECKOK( mp_copy(a, r) );
pR = r;
}
if (!q || q == a || q == b) {
} else {
pQ = q;
}
/*
If |a| <= |b|, we can compute the solution without division;
otherwise, we actually do the work required.
*/
if (cmp) {
/* r was set to a above. */
} else {
}
} else {
}
/* Compute the signs for the output */
/* Sq = ZPOS if Sa == Sb */ /* Sq = NEG if Sa != Sb */
/* Copy output, if it is needed */
if(q && q != pQ)
if(r && r != pR)
return res;
} /* end mp_div() */
/* }}} */
/* {{{ mp_div_2d(a, d, q, r) */
{
if(q) {
return res;
}
if(r) {
return res;
}
if(q) {
s_mp_div_2d(q, d);
}
if(r) {
s_mp_mod_2d(r, d);
}
return MP_OKAY;
} /* end mp_div_2d() */
/* }}} */
/* {{{ mp_expt(a, b, c) */
/*
mp_expt(a, b, c)
Compute c = a ** b, that is, raise a to the b power. Uses a
standard iterative square-and-multiply technique.
*/
{
mp_int s, x;
mp_digit d;
if(mp_cmp_z(b) < 0)
return MP_RANGE;
return res;
mp_set(&s, 1);
goto X;
/* Loop over low-order digits in ascending order */
/* Loop over bits of each non-maximal digit */
if(d & 1) {
goto CLEANUP;
}
d >>= 1;
goto CLEANUP;
}
}
/* Consider now the last digit... */
while(d) {
if(d & 1) {
goto CLEANUP;
}
d >>= 1;
goto CLEANUP;
}
if(mp_iseven(b))
mp_clear(&x);
X:
mp_clear(&s);
return res;
} /* end mp_expt() */
/* }}} */
/* {{{ mp_2expt(a, k) */
/* Compute a = 2^k */
{
return s_mp_2expt(a, k);
} /* end mp_2expt() */
/* }}} */
/* {{{ mp_mod(a, m, c) */
/*
mp_mod(a, m, c)
Compute c = a (mod m). Result will always be 0 <= c < m.
*/
{
int mag;
return MP_RANGE;
/*
If |a| > m, we need to divide to get the remainder and take the
absolute value.
If |a| < m, we don't need to do any division, just copy and adjust
the sign (if a is negative).
If |a| == m, we can simply set the result to zero.
This order is intended to minimize the average path length of the
comparison chain on common workloads -- the most frequent cases are
that |a| != m, so we do those first.
*/
return res;
return res;
}
} else if(mag < 0) {
return res;
if(mp_cmp_z(a) < 0) {
return res;
}
} else {
mp_zero(c);
}
return MP_OKAY;
} /* end mp_mod() */
/* }}} */
/* {{{ mp_mod_d(a, d, c) */
/*
mp_mod_d(a, d, c)
Compute c = a (mod d). Result will always be 0 <= c < d
*/
{
if(s_mp_cmp_d(a, d) > 0) {
return res;
} else {
else
}
if(c)
*c = rem;
return MP_OKAY;
} /* end mp_mod_d() */
/* }}} */
/* {{{ mp_sqrt(a, b) */
/*
mp_sqrt(a, b)
Compute the integer square root of a, and store the result in b.
Uses an integer-arithmetic version of Newton's iterative linear
approximation technique to determine this value; the result has the
following two properties:
b^2 <= a
(b+1)^2 >= a
It is a range error to pass a negative value.
*/
{
mp_int x, t;
/* Cannot take square root of a negative value */
return MP_RANGE;
/* Special cases for zero and one, trivial */
if(mp_cmp_d(a, 1) <= 0)
return mp_copy(a, b);
/* Initialize the temporaries we'll use below */
return res;
/* Compute an initial guess for the iteration as a itself */
goto X;
if (used > 1) {
}
for(;;) {
/* t = (x * x) - a */
mp_copy(&x, &t); /* can't fail, t is big enough for original x */
goto CLEANUP;
/* t = t / 2x */
s_mp_mul_2(&x);
goto CLEANUP;
s_mp_div_2(&x);
/* Terminate the loop, if the quotient is zero */
break;
/* x = x - t */
goto CLEANUP;
}
/* Copy result to output parameter */
mp_sub_d(&x, 1, &x);
s_mp_exch(&x, b);
mp_clear(&x);
X:
mp_clear(&t);
return res;
} /* end mp_sqrt() */
/* }}} */
/* }}} */
/*------------------------------------------------------------------------*/
/* {{{ Modular arithmetic */
#if MP_MODARITH
/* {{{ mp_addmod(a, b, m, c) */
/*
mp_addmod(a, b, m, c)
Compute c = (a + b) mod m
*/
{
return res;
return res;
return MP_OKAY;
}
/* }}} */
/* {{{ mp_submod(a, b, m, c) */
/*
mp_submod(a, b, m, c)
Compute c = (a - b) mod m
*/
{
return res;
return res;
return MP_OKAY;
}
/* }}} */
/* {{{ mp_mulmod(a, b, m, c) */
/*
mp_mulmod(a, b, m, c)
Compute c = (a * b) mod m
*/
{
return res;
return res;
return MP_OKAY;
}
/* }}} */
/* {{{ mp_sqrmod(a, m, c) */
#if MP_SQUARE
{
return res;
return res;
return MP_OKAY;
} /* end mp_sqrmod() */
#endif
/* }}} */
/* {{{ s_mp_exptmod(a, b, m, c) */
/*
s_mp_exptmod(a, b, m, c)
Compute c = (a ** b) mod m. Uses a standard square-and-multiply
method with modular reductions at each step. (This is basically the
same code as mp_expt(), except for the addition of the reductions)
The modular reductions are done using Barrett's algorithm (see
s_mp_reduce() below for details)
*/
{
mp_digit d;
return MP_RANGE;
return res;
goto X;
goto MU;
mp_set(&s, 1);
/* mu = b^2k / m */
goto CLEANUP;
/* Loop over digits of b in ascending order, except highest order */
/* Loop over the bits of the lower-order digits */
if(d & 1) {
goto CLEANUP;
goto CLEANUP;
}
d >>= 1;
goto CLEANUP;
goto CLEANUP;
}
}
/* Now do the last digit... */
while(d) {
if(d & 1) {
goto CLEANUP;
goto CLEANUP;
}
d >>= 1;
goto CLEANUP;
goto CLEANUP;
}
s_mp_exch(&s, c);
MU:
mp_clear(&x);
X:
mp_clear(&s);
return res;
} /* end s_mp_exptmod() */
/* }}} */
/* {{{ mp_exptmod_d(a, d, m, c) */
{
mp_int s, x;
return res;
goto X;
mp_set(&s, 1);
while(d != 0) {
if(d & 1) {
goto CLEANUP;
}
d /= 2;
goto CLEANUP;
}
s_mp_exch(&s, c);
mp_clear(&x);
X:
mp_clear(&s);
return res;
} /* end mp_exptmod_d() */
/* }}} */
#endif /* if MP_MODARITH */
/* }}} */
/*------------------------------------------------------------------------*/
/* {{{ Comparison functions */
/* {{{ mp_cmp_z(a) */
/*
mp_cmp_z(a)
Compare a <=> 0. Returns <0 if a<0, 0 if a=0, >0 if a>0.
*/
{
return MP_LT;
return MP_EQ;
else
return MP_GT;
} /* end mp_cmp_z() */
/* }}} */
/* {{{ mp_cmp_d(a, d) */
/*
mp_cmp_d(a, d)
Compare a <=> d. Returns <0 if a<d, 0 if a=d, >0 if a>d
*/
{
return MP_LT;
return s_mp_cmp_d(a, d);
} /* end mp_cmp_d() */
/* }}} */
/* {{{ mp_cmp(a, b) */
{
int mag;
return MP_EQ;
return mag;
else
return -mag;
return MP_GT;
} else {
return MP_LT;
}
} /* end mp_cmp() */
/* }}} */
/* {{{ mp_cmp_mag(a, b) */
/*
mp_cmp_mag(a, b)
Compares |a| <=> |b|, and returns an appropriate comparison result
*/
{
return s_mp_cmp(a, b);
} /* end mp_cmp_mag() */
/* }}} */
/* {{{ mp_cmp_int(a, z, kmflag) */
/*
This just converts z to an mp_int, and uses the existing comparison
routines. This is sort of inefficient, but it's not clear to me how
frequently this wil get used anyway. For small positive constants,
you can always use mp_cmp_d(), and for zero, there is mp_cmp_z().
*/
{
int out;
return out;
} /* end mp_cmp_int() */
/* }}} */
/* {{{ mp_isodd(a) */
/*
mp_isodd(a)
Returns a true (non-zero) value if a is odd, false (zero) otherwise.
*/
{
return (int)(DIGIT(a, 0) & 1);
} /* end mp_isodd() */
/* }}} */
/* {{{ mp_iseven(a) */
{
return !mp_isodd(a);
} /* end mp_iseven() */
/* }}} */
/* }}} */
/*------------------------------------------------------------------------*/
/* {{{ Number theoretic functions */
#if MP_NUMTH
/* {{{ mp_gcd(a, b, c) */
/*
Like the old mp_gcd() function, except computes the GCD using the
binary algorithm due to Josef Stein in 1961 (via Knuth).
*/
{
mp_int u, v, t;
mp_size k = 0;
return MP_RANGE;
return mp_copy(b, c);
return mp_copy(a, c);
}
return res;
goto U;
goto V;
/* Divide out common factors of 2 until at least 1 of a, b is even */
s_mp_div_2(&u);
s_mp_div_2(&v);
++k;
}
/* Initialize t */
if(mp_isodd(&u)) {
goto CLEANUP;
/* t = -v */
else
} else {
goto CLEANUP;
}
for(;;) {
while(mp_iseven(&t)) {
s_mp_div_2(&t);
}
goto CLEANUP;
} else {
goto CLEANUP;
/* v = -t */
else
}
goto CLEANUP;
if(s_mp_cmp_d(&t, 0) == MP_EQ)
break;
}
s_mp_2expt(&v, k); /* v = 2^k */
mp_clear(&v);
V:
mp_clear(&u);
U:
mp_clear(&t);
return res;
} /* end mp_gcd() */
/* }}} */
/* {{{ mp_lcm(a, b, c) */
/* We compute the least common multiple using the rule:
ab = [a, b](a, b)
... by computing the product, and dividing out the gcd.
*/
{
/* Set up temporaries */
return res;
goto GCD;
goto CLEANUP;
goto CLEANUP;
GCD:
return res;
} /* end mp_lcm() */
/* }}} */
/* {{{ mp_xgcd(a, b, g, x, y) */
/*
mp_xgcd(a, b, g, x, y)
Compute g = (a, b) and values x and y satisfying Bezout's identity
(that is, ax + by = g). This uses the binary extended GCD algorithm
based on the Stein algorithm used for mp_gcd()
See algorithm 14.61 in Handbook of Applied Cryptogrpahy.
*/
{
if(mp_cmp_z(b) == 0)
return MP_RANGE;
/* Initialize all these variables we need */
/* Divide by two until at least one of them is odd */
s_mp_div_2d(&xc,n);
s_mp_div_2d(&yc,n);
}
/* Loop through binary GCD algorithm */
do {
while(mp_iseven(&u)) {
s_mp_div_2(&u);
s_mp_div_2(&A); s_mp_div_2(&B);
} else {
s_mp_div_2(&A);
s_mp_div_2(&B);
}
}
while(mp_iseven(&v)) {
s_mp_div_2(&v);
s_mp_div_2(&C); s_mp_div_2(&D);
} else {
s_mp_div_2(&C);
s_mp_div_2(&D);
}
}
if(mp_cmp(&u, &v) >= 0) {
MP_CHECKOK( mp_sub(&u, &v, &u) );
MP_CHECKOK( mp_sub(&A, &C, &A) );
MP_CHECKOK( mp_sub(&B, &D, &B) );
} else {
MP_CHECKOK( mp_sub(&v, &u, &v) );
MP_CHECKOK( mp_sub(&C, &A, &C) );
MP_CHECKOK( mp_sub(&D, &B, &D) );
}
} while (mp_cmp_z(&u) != 0);
/* copy results to output */
if(x)
MP_CHECKOK( mp_copy(&C, x) );
if(y)
MP_CHECKOK( mp_copy(&D, y) );
if(g)
while(last >= 0)
return res;
} /* end mp_xgcd() */
/* }}} */
{
mp_digit d;
mp_size n = 0;
unsigned int ix;
return n;
n += MP_DIGIT_BIT;
if (!d)
return 0; /* shouldn't happen, but ... */
#if !defined(MP_USE_UINT_DIGIT)
if (!(d & 0xffffffffU)) {
d >>= 32;
n += 32;
}
#endif
if (!(d & 0xffffU)) {
d >>= 16;
n += 16;
}
if (!(d & 0xffU)) {
d >>= 8;
n += 8;
}
if (!(d & 0xfU)) {
d >>= 4;
n += 4;
}
if (!(d & 0x3U)) {
d >>= 2;
n += 2;
}
if (!(d & 0x1U)) {
d >>= 1;
n += 1;
}
#if MP_ARGCHK == 2
assert(0 != (d & 1));
#endif
return n;
}
/* Given a and prime p, computes c and k such that a*c == 2**k (mod p).
** Returns k (positive) or error (negative).
** This technique from the paper "Fast Modular Reciprocals" (unpublished)
** by Richard Schroeppel (a.k.a. Captain Nemo).
*/
{
mp_err k = 0;
mp_int d, f, g;
MP_DIGITS(&d) = 0;
MP_DIGITS(&f) = 0;
MP_DIGITS(&g) = 0;
mp_set(c, 1);
mp_zero(&d);
if (mp_cmp_z(&f) == 0) {
} else
for (;;) {
int diff_sign;
while (mp_iseven(&f)) {
mp_size n = mp_trailing_zeros(&f);
if (!n) {
goto CLEANUP;
}
s_mp_div_2d(&f, n);
MP_CHECKOK( s_mp_mul_2d(&d, n) );
k += n;
}
res = k;
break;
}
if (diff_sign < 0) { /* f < g */
s_mp_exch(&f, &g);
s_mp_exch(c, &d);
} else if (diff_sign == 0) { /* f == g */
break;
}
} else {
}
}
if (res >= 0) {
MP_CHECKOK( mp_add(c, p, c) );
}
res = k;
}
mp_clear(&d);
mp_clear(&f);
mp_clear(&g);
return res;
}
/* Compute T = (P ** -1) mod MP_RADIX. Also works for 16-bit mp_digits.
** This technique from the paper "Fast Modular Reciprocals" (unpublished)
** by Richard Schroeppel (a.k.a. Captain Nemo).
*/
{
mp_digit T = P;
T *= 2 - (P * T);
T *= 2 - (P * T);
T *= 2 - (P * T);
T *= 2 - (P * T);
#if !defined(MP_USE_UINT_DIGIT)
T *= 2 - (P * T);
T *= 2 - (P * T);
#endif
return T;
}
/* Given c, k, and prime p, where a*c == 2**k (mod p),
** Compute x = (a ** -1) mod p. This is similar to Montgomery reduction.
** This technique from the paper "Fast Modular Reciprocals" (unpublished)
** by Richard Schroeppel (a.k.a. Captain Nemo).
*/
{
int k_orig = k;
mp_digit r;
if (mp_cmp_z(c) < 0) { /* c < 0 */
} else {
}
/* make sure x is large enough */
r = 0 - s_mp_invmod_radix(MP_DIGIT(p,0));
int j = MP_MIN(k, MP_DIGIT_BIT);
if (j < MP_DIGIT_BIT) {
}
k -= j;
}
s_mp_clamp(x);
s_mp_div_2d(x, k_orig);
return res;
}
/* compute mod inverse using Schroeppel's method, only if m is odd */
{
int k;
mp_int x;
return MP_RANGE;
if (mp_iseven(m))
return MP_UNDEF;
MP_DIGITS(&x) = 0;
if (a == c) {
return res;
if (a == m)
m = &x;
a = &x;
} else if (m == c) {
return res;
m = &x;
} else {
MP_DIGITS(&x) = 0;
}
MP_CHECKOK( s_mp_almost_inverse(a, m, c) );
k = res;
MP_CHECKOK( s_mp_fixup_reciprocal(c, m, k, c) );
mp_clear(&x);
return res;
}
/* Known good algorithm for computing modular inverse. But slow. */
{
mp_int g, x;
return MP_RANGE;
MP_DIGITS(&g) = 0;
MP_DIGITS(&x) = 0;
goto CLEANUP;
}
mp_clear(&x);
mp_clear(&g);
return res;
}
/* modular inverse where modulus is 2**k. */
/* c = a**-1 mod 2**k */
{
if (mp_iseven(a))
return MP_UNDEF;
if (k <= MP_DIGIT_BIT) {
if (k < MP_DIGIT_BIT)
mp_set(c, i);
return MP_OKAY;
}
s_mp_mod_2d(&val, k);
do {
s_mp_mod_2d(&t1, k);
}
break;
} while (--ix > 0);
if (!ix) {
} else {
}
return res;
}
{
mp_size k;
/*static const mp_digit d1 = 1; */
/*static const mp_int one = { MP_ZPOS, 1, 1, (mp_digit *)&d1 }; */
if ((res = s_mp_ispow2(m)) >= 0) {
k = res;
return s_mp_invmod_2d(a, k, c);
}
MP_DIGITS(&evenFactor) = 0;
k = mp_trailing_zeros(m);
s_mp_div_2d(&oddFactor, k);
/* compute a**-1 mod oddFactor. */
/* compute a**-1 mod evenFactor, where evenFactor == 2**k. */
/* Use Chinese Remainer theorem to compute a**-1 mod m. */
/* let m1 = oddFactor, v1 = oddPart,
* let m2 = evenFactor, v2 = evenPart.
*/
/* Compute C2 = m1**-1 mod m2. */
/* compute u = (v2 - v1)*C2 mod m2 */
s_mp_mod_2d(&tmp2, k);
}
/* compute answer = v1 + u*m1 */
/* not sure this is necessary, but it's low cost if not. */
MP_CHECKOK( mp_mod(c, m, c) );
return res;
}
/* {{{ mp_invmod(a, m, c) */
/*
mp_invmod(a, m, c)
Compute c = a^-1 (mod m), if there is an inverse for a (mod m).
This is equivalent to the question of whether (a, m) = 1. If not,
MP_UNDEF is returned, and there is no inverse.
*/
{
return MP_RANGE;
if (mp_isodd(m)) {
return s_mp_invmod_odd_m(a, m, c);
}
if (mp_iseven(a))
return MP_UNDEF; /* not invertable */
return s_mp_invmod_even_m(a, m, c);
} /* end mp_invmod() */
/* }}} */
#endif /* if MP_NUMTH */
/* }}} */
/*------------------------------------------------------------------------*/
/* {{{ mp_print(mp, ofp) */
#if MP_IOFUNC
/*
mp_print(mp, ofp)
Print a textual representation of the given mp_int on the output
stream 'ofp'. Output is generated using the internal radix.
*/
{
int ix;
return;
}
} /* end mp_print() */
#endif /* if MP_IOFUNC */
/* }}} */
/*------------------------------------------------------------------------*/
/* {{{ More I/O Functions */
/* {{{ mp_read_raw(mp, str, len) */
/*
mp_read_raw(mp, str, len)
Read in a raw value (base 256) into the given mp_int
*/
{
int ix;
/* Get sign from first byte */
if(ustr[0])
else
/* Read the rest of the digits */
return res;
return res;
}
return MP_OKAY;
} /* end mp_read_raw() */
/* }}} */
/* {{{ mp_raw_size(mp) */
{
} /* end mp_raw_size() */
/* }}} */
/* {{{ mp_toraw(mp, str) */
{
/* Iterate over each digit... */
/* Unpack digit bytes, high order first */
}
}
return MP_OKAY;
} /* end mp_toraw() */
/* }}} */
/* {{{ mp_read_radix(mp, str, radix) */
/*
mp_read_radix(mp, str, radix)
Read an integer from the given string, and set mp to the resulting
value. The input is presumed to be in base 10. Leading non-digit
characters are ignored, and the function reads until a non-digit
character or the end of the string.
*/
{
/* Skip leading non-digit characters until a digit or '-' or '+' */
++ix;
}
++ix;
++ix;
}
return res;
return res;
++ix;
}
else
return MP_OKAY;
} /* end mp_read_radix() */
{
int cx;
/* Skip leading non-digit characters until a digit or '-' or '+' */
cx != '-' &&
cx != '+') {
++str;
}
if (cx == '-') {
++str;
} else if (cx == '+') {
++str;
}
if (str[0] == '0') {
radix = 16;
str += 2;
} else {
radix = 8;
str++;
}
}
}
return res;
}
/* }}} */
/* {{{ mp_radix_size(mp, radix) */
{
int bits;
return 0;
} /* end mp_radix_size() */
/* }}} */
/* {{{ mp_toradix(mp, str, radix) */
{
str[0] = '0';
} else {
char ch;
return res;
/* Save sign for later, and take absolute value */
/* Generate output digits in reverse order */
return res;
}
/* Generate digits, use capital letters */
}
/* Add - sign if original value was negative */
/* Add trailing NUL to end the string */
/* Reverse the digits and sign indicator */
ix = 0;
++ix;
--pos;
}
}
return MP_OKAY;
} /* end mp_toradix() */
/* }}} */
/* {{{ mp_tovalue(ch, r) */
{
return s_mp_tovalue(ch, r);
} /* end mp_tovalue() */
/* }}} */
/* }}} */
/* {{{ mp_strerror(ec) */
/*
mp_strerror(ec)
Return a string describing the meaning of error code 'ec'. The
string returned is allocated in static memory, so the caller should
not attempt to modify or free the memory associated with this
string.
*/
{
/* Code values are negative, so the senses of these comparisons
are accurate */
return mp_err_string[0]; /* unknown error code */
} else {
}
} /* end mp_strerror() */
/* }}} */
/*========================================================================*/
/*------------------------------------------------------------------------*/
/* Static function definitions (internal use only) */
/* {{{ Memory management */
/* {{{ s_mp_grow(mp, min) */
/* Make sure there are at least 'min' digits allocated to mp */
{
/* Set min to next nearest default precision block size */
return MP_MEM;
#if MP_CRYPTO
#endif
}
return MP_OKAY;
} /* end s_mp_grow() */
/* }}} */
/* {{{ s_mp_pad(mp, min) */
/* Make sure the used size of mp is at least 'min', growing if needed */
{
/* Make sure there is room to increase precision */
return res;
} else {
}
/* Increase precision; should already be 0-filled */
}
return MP_OKAY;
} /* end s_mp_pad() */
/* }}} */
/* {{{ s_mp_setz(dp, count) */
#if MP_MACRO == 0
/* Set 'count' digits pointed to by dp to be zeroes */
{
#if MP_MEMSET == 0
int ix;
#else
#endif
} /* end s_mp_setz() */
#endif
/* }}} */
/* {{{ s_mp_copy(sp, dp, count) */
#if MP_MACRO == 0
/* Copy 'count' digits from sp to dp */
{
#if MP_MEMCPY == 0
int ix;
#else
#endif
} /* end s_mp_copy() */
#endif
/* }}} */
/* {{{ s_mp_alloc(nb, ni, kmflag) */
#if MP_MACRO == 0
/* Allocate ni records of nb bytes each, and return a pointer to that */
{
++mp_allocs;
#ifdef _KERNEL
return (mp);
#else
#endif
} /* end s_mp_alloc() */
#endif
/* }}} */
/* {{{ s_mp_free(ptr) */
#if MP_MACRO == 0
/* Free the memory pointed to by ptr */
{
if(ptr) {
++mp_frees;
#ifdef _KERNEL
#else
#endif
}
} /* end s_mp_free() */
#endif
/* }}} */
/* {{{ s_mp_clamp(mp) */
#if MP_MACRO == 0
/* Remove leading zeroes from the given value */
{
--used;
} /* end s_mp_clamp() */
#endif
/* }}} */
/* {{{ s_mp_exch(a, b) */
/* Exchange the data for a and b; (b, a) = (a, b) */
{
tmp = *a;
*a = *b;
*b = tmp;
} /* end s_mp_exch() */
/* }}} */
/* }}} */
/* {{{ Arithmetic helpers */
/* {{{ s_mp_lshd(mp, p) */
/*
Shift mp leftward by p digits, growing if needed, and zero-filling
the in-shifted digits at the right end. This is a convenient
alternative to multiplication by powers of the radix
The value of USED(mp) must already have been set to the value for
the shifted result.
*/
{
int ix;
if(p == 0)
return MP_OKAY;
return MP_OKAY;
return res;
/* Shift all the significant figures over as needed */
/* Fill the bottom digits with zeroes */
return MP_OKAY;
} /* end s_mp_lshd() */
/* }}} */
/* {{{ s_mp_mul_2d(mp, d) */
/*
Multiply the integer by 2^d, where d is a number of bits. This
amounts to a bitwise shift of the value.
*/
{
dshift = d / MP_DIGIT_BIT;
bshift = d % MP_DIGIT_BIT;
/* bits to be shifted out of the top word */
return res;
return res;
if (bshift) {
}
}
s_mp_clamp(mp);
return MP_OKAY;
} /* end s_mp_mul_2d() */
/* {{{ s_mp_rshd(mp, p) */
/*
Shift mp rightward by p digits. Maintains the invariant that
digits above the precision are all zero. Digits shifted off the
end are lost. Cannot fail.
*/
{
if(p == 0)
return;
/* Shortcut when all digits are to be shifted off */
return;
}
/* Shift all the significant figures over as needed */
/* Fill the top digits with zeroes */
while (p-- > 0)
*dst++ = 0;
#if 0
/* Strip off any leading zeroes */
s_mp_clamp(mp);
#endif
} /* end s_mp_rshd() */
/* }}} */
/* {{{ s_mp_div_2(mp) */
/* Divide by two -- take advantage of radix properties to do it fast */
{
} /* end s_mp_div_2() */
/* }}} */
/* {{{ s_mp_mul_2(mp) */
{
/* Shift digits leftward by 1 bit */
}
/* Deal with rollover from last digit */
if (kin) {
return res;
}
}
return MP_OKAY;
} /* end s_mp_mul_2() */
/* }}} */
/* {{{ s_mp_mod_2d(mp, d) */
/*
Remainder the integer by 2^d, where d is a number of bits. This
amounts to a bitwise AND of the value, and does not require the full
division code
*/
{
return;
/* Flush all the bits above 2^d in its digit */
/* Flush all digits above the one with 2^d in it */
s_mp_clamp(mp);
} /* end s_mp_mod_2d() */
/* }}} */
/* {{{ s_mp_div_2d(mp, d) */
/*
Divide the integer by 2^d, where d is a number of bits. This
amounts to a bitwise shift of the value, and does not require the
full division code (used in Barrett reduction, see below)
*/
{
int ix;
d %= DIGIT_BIT;
if (d) {
save = 0;
}
}
s_mp_clamp(mp);
} /* end s_mp_div_2d() */
/* }}} */
/* {{{ s_mp_norm(a, b, *d) */
/*
s_mp_norm(a, b, *d)
Normalize a and b for division, where b is the divisor. In order
that we might make good guesses for quotient digits, we want the
leading digit of b to be at least half the radix, which we
accomplish by multiplying a and b by a power of 2. The exponent
(shift count) is placed in *pd, so that the remainder can be shifted
back at the end of the division process.
*/
{
mp_digit d;
d = 0;
b_msd <<= 1;
++d;
}
if (d) {
MP_CHECKOK( s_mp_mul_2d(a, d) );
MP_CHECKOK( s_mp_mul_2d(b, d) );
}
*pd = d;
return res;
} /* end s_mp_norm() */
/* }}} */
/* }}} */
/* {{{ Primitive digit arithmetic */
/* {{{ s_mp_add_d(mp, d) */
/* Add d to |mp| in place */
{
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
mp_word w, k = 0;
k = CARRYOUT(w);
k = CARRYOUT(w);
++ix;
}
if(k != 0) {
return res;
}
return MP_OKAY;
#else
}
/* mp is growing */
}
return res;
#endif
} /* end s_mp_add_d() */
/* }}} */
/* {{{ s_mp_sub_d(mp, d) */
/* Subtract d from |mp| in place, assumes |mp| > d */
{
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
mp_word w, b = 0;
/* Compute initial subtraction */
b = CARRYOUT(w) ? 0 : 1;
/* Propagate borrows leftward */
b = CARRYOUT(w) ? 0 : 1;
++ix;
}
/* Remove leading zeroes */
s_mp_clamp(mp);
/* If we have a borrow out, it's a violation of the input invariant */
if(b)
return MP_RANGE;
else
return MP_OKAY;
#else
}
s_mp_clamp(mp);
#endif
} /* end s_mp_sub_d() */
/* }}} */
/* {{{ s_mp_mul_d(a, d) */
/* Compute a = a * d, single digit multiplication */
{
int pow;
if (!d) {
mp_zero(a);
return MP_OKAY;
}
if (d == 1)
return MP_OKAY;
if (0 <= (pow = s_mp_ispow2d(d))) {
}
s_mp_clamp(a);
return res;
} /* end s_mp_mul_d() */
/* }}} */
/* {{{ s_mp_div_d(mp, d, r) */
/*
s_mp_div_d(mp, d, r)
Compute the quotient mp = mp / d and remainder r = mp mod d, for a
single digit d. If r is null, the remainder will be discarded.
*/
{
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_DIV_WORD)
mp_word w = 0, q;
#else
mp_digit w, q;
#endif
int ix;
if(d == 0)
return MP_RANGE;
if (d == 1) {
if (r)
*r = 0;
return MP_OKAY;
}
/* could check for power of 2 here, but mp_div_d does that. */
q = n / d;
rem = n % d;
if (r)
*r = rem;
return MP_OKAY;
}
/* Make room for the quotient */
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_DIV_WORD)
if(w >= d) {
q = w / d;
w = w % d;
} else {
q = 0;
}
}
#else
{
mp_digit p;
#if !defined(MP_ASSEMBLY_DIV_2DX1D)
#endif
#if !defined(MP_ASSEMBLY_DIV_2DX1D)
if (norm)
d <<= norm;
#endif
p = 0;
if (p) {
MP_CHECKOK( s_mpv_div_2dx1d(p, w, d, &q, &w) );
} else if (w >= d) {
q = w / d;
w = w % d;
} else {
q = 0;
}
p = w;
}
#if !defined(MP_ASSEMBLY_DIV_2DX1D)
if (norm)
w >>= norm;
#endif
}
#endif
/* Deliver the remainder, if desired */
if(r)
*r = (mp_digit)w;
s_mp_clamp(");
return res;
} /* end s_mp_div_d() */
/* }}} */
/* }}} */
/* {{{ Primitive full arithmetic */
/* {{{ s_mp_add(a, b) */
/* Compute a = |a| + |b| */
{
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
mp_word w = 0;
#else
#endif
/* Make sure a has enough precision for the output value */
return res;
/*
Add up all digits up to the precision of b. If b had initially
the same precision as a, or greater, we took care of it by the
padding step above, so there is no problem. If b had initially
less precision, we'll have to make sure the carry out is duly
propagated upward among the higher-order digits of the sum.
*/
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
w = CARRYOUT(w);
#else
d = *pa;
d = (sum < d); /* detect overflow */
#endif
}
/* If we run out of 'b' digits before we're actually done, make
sure the carries get propagated upward...
*/
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
w = w + *pa;
w = CARRYOUT(w);
++ix;
}
#else
++ix;
}
#endif
/* If there's an overall carry out, increase precision and include
it. We could have done this initially, but why touch the memory
allocator unless we're sure we have to?
*/
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
if (w) {
return res;
}
#else
if (carry) {
return res;
}
#endif
return MP_OKAY;
} /* end s_mp_add() */
/* }}} */
/* Compute c = |a| + |b| */ /* magnitude addition */
{
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
mp_word w = 0;
#else
#endif
a = b;
b = xch;
}
/* Make sure a has enough precision for the output value */
return res;
/*
Add up all digits up to the precision of b. If b had initially
the same precision as a, or greater, we took care of it by the
exchange step above, so there is no problem. If b had initially
less precision, we'll have to make sure the carry out is duly
propagated upward among the higher-order digits of the sum.
*/
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
w = CARRYOUT(w);
#else
d = *pa++;
d = (sum < d); /* detect overflow */
#endif
}
/* If we run out of 'b' digits before we're actually done, make
sure the carries get propagated upward...
*/
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
w = w + *pa++;
w = CARRYOUT(w);
#else
#endif
}
/* If there's an overall carry out, increase precision and include
it. We could have done this initially, but why touch the memory
allocator unless we're sure we have to?
*/
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
if (w) {
return res;
++used;
}
#else
if (carry) {
return res;
++used;
}
#endif
return MP_OKAY;
}
/* {{{ s_mp_add_offset(a, b, offset) */
/* Compute a = |a| + ( |b| * (RADIX ** offset) ) */
{
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
mp_word w, k = 0;
#else
#endif
/* Make sure a has enough precision for the output value */
return res;
/*
Add up all digits up to the precision of b. If b had initially
the same precision as a, or greater, we took care of it by the
padding step above, so there is no problem. If b had initially
less precision, we'll have to make sure the carry out is duly
propagated upward among the higher-order digits of the sum.
*/
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
k = CARRYOUT(w);
#else
d = (sum < d);
#endif
}
/* If we run out of 'b' digits before we're actually done, make
sure the carries get propagated upward...
*/
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
k = CARRYOUT(w);
}
#else
}
#endif
/* If there's an overall carry out, increase precision and include
it. We could have done this initially, but why touch the memory
allocator unless we're sure we have to?
*/
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
if(k) {
return res;
}
#else
if (carry) {
return res;
}
#endif
s_mp_clamp(a);
return MP_OKAY;
} /* end s_mp_add_offset() */
/* }}} */
/* {{{ s_mp_sub(a, b) */
/* Compute a = |a| - |b|, assumes |a| >= |b| */
{
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
mp_sword w = 0;
#else
#endif
/*
Subtract and propagate borrow. Up to the precision of b, this
accounts for the digits of b; after that, we just make sure the
carries get to the right place. This saves having to pad b out to
the precision of a just to make the loops work right...
*/
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
w >>= MP_DIGIT_BIT;
#else
d = *pa;
d = (diff > d); /* detect borrow */
++d;
borrow = d;
#endif
}
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
w = w + *pa;
w >>= MP_DIGIT_BIT;
}
#else
d = *pa;
}
#endif
/* Clobber any leading zeroes we created */
s_mp_clamp(a);
/*
If there was a borrow out, then |b| > |a| in violation
of our input invariant. We've already done the work,
but we'll at least complain about it...
*/
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
#else
#endif
} /* end s_mp_sub() */
/* }}} */
/* Compute c = |a| - |b|, assumes |a| >= |b| */ /* magnitude subtract */
{
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
mp_sword w = 0;
#else
#endif
/* Make sure a has enough precision for the output value */
return res;
/*
Subtract and propagate borrow. Up to the precision of b, this
accounts for the digits of b; after that, we just make sure the
carries get to the right place. This saves having to pad b out to
the precision of a just to make the loops work right...
*/
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
w >>= MP_DIGIT_BIT;
#else
d = *pa++;
d = (diff > d);
++d;
borrow = d;
#endif
}
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
w = w + *pa++;
w >>= MP_DIGIT_BIT;
#else
d = *pa++;
#endif
}
/* Clobber any leading zeroes we created */
s_mp_clamp(c);
/*
If there was a borrow out, then |b| > |a| in violation
of our input invariant. We've already done the work,
but we'll at least complain about it...
*/
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
#else
#endif
}
/* {{{ s_mp_mul(a, b) */
/* Compute a = |a| * |b| */
{
return mp_mul(a, b, a);
} /* end s_mp_mul() */
/* }}} */
#if defined(MP_USE_UINT_DIGIT) && defined(MP_USE_LONG_LONG_MULTIPLY)
/* This trick works on Sparc V8 CPUs with the Workshop compilers. */
{ unsigned long long product = (unsigned long long)a * b; \
{ Plo = asm ("mulq %a0, %a1, %v0", a, b);\
Phi = asm ("umulh %a0, %a1, %v0", a, b); }
#else
Phi += MP_HALF_RADIX; \
a1b0 <<= MP_HALF_DIGIT_BIT; \
++Phi; \
}
#endif
#if !defined(MP_ASSEMBLY_MULTIPLY)
/* c = a * b */
{
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_MUL_WORD)
mp_digit d = 0;
/* Inner product: Digits of a */
while (a_len--) {
*c++ = ACCUM(w);
d = CARRYOUT(w);
}
*c = d;
#else
while (a_len--) {
++a1b1;
*c++ = a0b0;
}
*c = carry;
#endif
}
/* c += a * b */
mp_digit *c)
{
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_MUL_WORD)
mp_digit d = 0;
/* Inner product: Digits of a */
while (a_len--) {
*c++ = ACCUM(w);
d = CARRYOUT(w);
}
*c = d;
#else
while (a_len--) {
++a1b1;
++a1b1;
*c++ = a0b0;
}
*c = carry;
#endif
}
/* Presently, this is only used by the Montgomery arithmetic code. */
/* c += a * b */
{
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_MUL_WORD)
mp_digit d = 0;
/* Inner product: Digits of a */
while (a_len--) {
*c++ = ACCUM(w);
d = CARRYOUT(w);
}
while (d) {
*c++ = ACCUM(w);
d = CARRYOUT(w);
}
#else
while (a_len--) {
++a1b1;
++a1b1;
*c++ = a0b0;
}
while (carry) {
*c++ = carry;
}
#endif
}
#endif
#if defined(MP_USE_UINT_DIGIT) && defined(MP_USE_LONG_LONG_MULTIPLY)
/* This trick works on Sparc V8 CPUs with the Workshop compilers. */
{ unsigned long long square = (unsigned long long)a * a; \
{ Plo = asm ("mulq %a0, %a0, %v0", a);\
Phi = asm ("umulh %a0, %a0, %v0", a); }
#else
++Phi; \
}
#endif
#if !defined(MP_ASSEMBLY_SQUARE)
/* Add the squares of the digits of a to the digits of b. */
{
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_MUL_WORD)
mp_word w;
mp_digit d;
w = 0;
#define ADD_SQUARE(n) \
d = pa[n]; \
w = (w >> DIGIT_BIT)
ADD_SQUARE(0);
ADD_SQUARE(1);
ADD_SQUARE(2);
ADD_SQUARE(3);
pa += 4;
ps += 8;
}
if (ix) {
switch (ix) {
case 0: break;
}
}
while (w) {
w += *ps;
w = (w >> DIGIT_BIT);
}
#else
while (a_len--) {
/* here a1a1 and a0a0 constitute a_i ** 2 */
++a1a1;
/* now add to ps */
++a1a1;
}
while (carry) {
}
#endif
}
#endif
#if (defined(MP_NO_MP_WORD) || defined(MP_NO_DIV_WORD)) \
&& !defined(MP_ASSEMBLY_DIV_2DX1D)
/*
** Divide 64-bit (Nhi,Nlo) by 32-bit divisor, which must be normalized
** so its high bit is 1. This code is from NSPR.
*/
{
if (r1 < m) {
}
}
r1 -= m;
if (r0 < m) {
}
}
if (qp)
if (rp)
return MP_OKAY;
}
#endif
#if MP_SQUARE
/* {{{ s_mp_sqr(a) */
{
return res;
}
return res;
}
/* }}} */
#endif
/* {{{ s_mp_div(a, b) */
/*
s_mp_div(a, b)
Compute a = a / b and b = a mod b. Assumes b > a.
*/
{
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_DIV_WORD)
#else
#endif
mp_digit d;
int ix;
return MP_RANGE;
/* Shortcut if divisor is power of two */
return MP_OKAY;
}
DIGITS(&t) = 0;
/* A working temporary for division */
/* Normalize to optimize guessing */
/* Perform the division itself...woo! */
/* Find a partial substring of rem which is at least div */
/* If we didn't find one, we're finished dividing */
int i;
int unusedRem;
-- unusedRem;
#if MP_ARGCHK == 2
#endif
}
/* Compute a guess for the next quotient digit */
q_msd = 1;
#if !defined(MP_NO_MP_WORD) && !defined(MP_NO_DIV_WORD)
--q_msd;
#else
mp_digit r;
#endif
} else {
q_msd = 0;
}
#if MP_ARGCHK == 2
#endif
if (q_msd <= 0)
break;
/* See what that multiplies out to */
/*
If it's too big, back it off. We should not have to do this
more than once, or, in rare cases, twice. Knuth describes a
method by which this could be reduced to a maximum of once, but
I didn't implement that here.
* When using s_mpv_div_2dx1d, we may have to do this 3 times.
*/
--q_msd;
}
if (i < 0) {
goto CLEANUP;
}
/* At this point, q_msd should be the right next digit */
/*
Include the digit in the quotient. We allocated enough memory
for any quotient we could ever possibly get, so we should not
have to check for failures here
*/
}
/* Denormalize remainder */
if (d) {
s_mp_div_2d(rem, d);
}
mp_clear(&t);
return res;
} /* end s_mp_div() */
/* }}} */
/* {{{ s_mp_2expt(a, k) */
{
mp_zero(a);
return res;
return MP_OKAY;
} /* end s_mp_2expt() */
/* }}} */
/* {{{ s_mp_reduce(x, m, mu) */
/*
Compute Barrett reduction, x (mod m), given a precomputed value for
mu = b^2k / m, where b = RADIX and k = #digits(m). This should be
faster than straight division, when many reductions by the same
value of m are required (such as in modular exponentiation). This
can nearly halve the time required to do modular exponentiation,
as compared to using the full integer divide to reduce.
This algorithm was derived from the _Handbook of Applied
Cryptography_ by Menezes, Oorschot and VanStone, Ch. 14,
pp. 603-604.
*/
{
mp_int q;
return res;
/* x = x mod b^(k+1), quick (no division) */
/* q = q * m mod b^(k+1), quick (no division) */
s_mp_mul(&q, m);
/* x = x - q */
goto CLEANUP;
/* If x < 0, add b^(k+1) to it */
if(mp_cmp_z(x) < 0) {
mp_set(&q, 1);
goto CLEANUP;
goto CLEANUP;
}
/* Back off if it's too big */
while(mp_cmp(x, m) >= 0) {
break;
}
mp_clear(&q);
return res;
} /* end s_mp_reduce() */
/* }}} */
/* }}} */
/* {{{ Primitive comparisons */
/* {{{ s_mp_cmp(a, b) */
/* Compare |a| <=> |b|, return 0 if equal, <0 if a<b, >0 if a>b */
{
{
goto IS_GT;
goto IS_LT;
}
{
while (used_a >= 4) {
pa -= 4;
pb -= 4;
used_a -= 4;
CMP_AB(3);
CMP_AB(2);
CMP_AB(1);
CMP_AB(0);
}
/* do nothing */;
done:
goto IS_GT;
goto IS_LT;
}
return MP_EQ;
return MP_LT;
return MP_GT;
} /* end s_mp_cmp() */
/* }}} */
/* {{{ s_mp_cmp_d(a, d) */
/* Compare |a| <=> d, return 0 if equal, <0 if a<d, >0 if a>d */
{
if(USED(a) > 1)
return MP_GT;
if(DIGIT(a, 0) < d)
return MP_LT;
else if(DIGIT(a, 0) > d)
return MP_GT;
else
return MP_EQ;
} /* end s_mp_cmp_d() */
/* }}} */
/* {{{ s_mp_ispow2(v) */
/*
Returns -1 if the value is not a power of two; otherwise, it returns
k such that v = 2^k, i.e. lg(v).
*/
{
mp_digit d;
extra = s_mp_ispow2d(d);
return extra;
while (--ix >= 0) {
return -1; /* not a power of two */
extra += MP_DIGIT_BIT;
}
return extra;
} /* end s_mp_ispow2() */
/* }}} */
/* {{{ s_mp_ispow2d(d) */
{
if ((d != 0) && ((d & (d-1)) == 0)) { /* d is a power of 2 */
int pow = 0;
#if defined (MP_USE_UINT_DIGIT)
if (d & 0xffff0000U)
pow += 16;
if (d & 0xff00ff00U)
pow += 8;
if (d & 0xf0f0f0f0U)
pow += 4;
if (d & 0xccccccccU)
pow += 2;
if (d & 0xaaaaaaaaU)
pow += 1;
#elif defined(MP_USE_LONG_LONG_DIGIT)
if (d & 0xffffffff00000000ULL)
pow += 32;
if (d & 0xffff0000ffff0000ULL)
pow += 16;
if (d & 0xff00ff00ff00ff00ULL)
pow += 8;
if (d & 0xf0f0f0f0f0f0f0f0ULL)
pow += 4;
if (d & 0xccccccccccccccccULL)
pow += 2;
if (d & 0xaaaaaaaaaaaaaaaaULL)
pow += 1;
#elif defined(MP_USE_LONG_DIGIT)
if (d & 0xffffffff00000000UL)
pow += 32;
if (d & 0xffff0000ffff0000UL)
pow += 16;
if (d & 0xff00ff00ff00ff00UL)
pow += 8;
if (d & 0xf0f0f0f0f0f0f0f0UL)
pow += 4;
if (d & 0xccccccccccccccccUL)
pow += 2;
if (d & 0xaaaaaaaaaaaaaaaaUL)
pow += 1;
#else
#error "unknown type for mp_digit"
#endif
return pow;
}
return -1;
} /* end s_mp_ispow2d() */
/* }}} */
/* }}} */
/* {{{ Primitive I/O helpers */
/* {{{ s_mp_tovalue(ch, r) */
/*
Convert the given character to its digit value, in the given radix.
If the given character is not understood in the given radix, -1 is
returned. Otherwise the digit's numeric value is returned.
The results will be odd if you use a radix < 2 or > 62, you are
expected to know what you're up to.
*/
{
if(r > 36)
else
else if(xch == '+')
val = 62;
else if(xch == '/')
val = 63;
else
return -1;
return -1;
return val;
} /* end s_mp_tovalue() */
/* }}} */
/* {{{ s_mp_todigit(val, r, low) */
/*
Convert val to a radix-r digit, if possible. If val is out of range
for r, returns zero. Otherwise, returns an ASCII character denoting
the value in the given radix.
The results may be odd if you use a radix < 2 or > 64, you are
expected to know what you're doing.
*/
{
char ch;
if(val >= (unsigned int)r)
return 0;
if(r <= 36 && low)
return ch;
} /* end s_mp_todigit() */
/* }}} */
/* {{{ s_mp_outlen(bits, radix) */
/*
Return an estimate for how long a string is needed to hold a radix
r representation of a number with 'bits' significant bits, plus an
extra for a zero terminator (assuming C style strings here)
*/
{
} /* end s_mp_outlen() */
/* }}} */
/* }}} */
/* {{{ mp_read_unsigned_octets(mp, str, len) */
/* mp_read_unsigned_octets(mp, str, len)
Read in a raw value (base 256) into the given mp_int
No sign bit, number is positive. Leading zeros ignored.
*/
{
int count;
mp_digit d;
if (count) {
d = (d << 8) | *str++;
}
}
/* Read the rest of the digits */
d = (d << 8) | *str++;
}
if (!d)
continue;
} else {
return res;
}
}
return MP_OKAY;
} /* end mp_read_unsigned_octets() */
/* }}} */
/* {{{ mp_unsigned_octet_size(mp) */
int
{
int bytes;
int ix;
mp_digit d = 0;
/* subtract leading zeros. */
/* Iterate over each digit... */
if (d)
break;
bytes -= sizeof(d);
}
if (!bytes)
return 1;
/* Have MSD, check digit bytes, high order first */
if (x)
break;
--bytes;
}
return bytes;
} /* end mp_unsigned_octet_size() */
/* }}} */
/* {{{ mp_to_unsigned_octets(mp, str) */
/* output a buffer of big endian octets no longer than specified. */
{
unsigned int bytes;
/* Iterate over each digit... */
int jx;
/* Unpack digit bytes, high order first */
if (!pos && !x) /* suppress leading zeros */
continue;
}
}
if (!pos)
return pos;
} /* end mp_to_unsigned_octets() */
/* }}} */
/* {{{ mp_to_signed_octets(mp, str) */
/* output a buffer of big endian octets no longer than specified. */
{
unsigned int bytes;
/* Iterate over each digit... */
int jx;
/* Unpack digit bytes, high order first */
if (!pos) {
if (!x) /* suppress leading zeros */
continue;
if (x & 0x80) { /* add one leading zero to make output positive. */
return MP_BADARG;
}
}
}
}
if (!pos)
return pos;
} /* end mp_to_signed_octets() */
/* }}} */
/* {{{ mp_to_fixlen_octets(mp, str) */
/* output a buffer of big endian octets exactly as long as requested. */
{
unsigned int bytes;
/* place any needed leading zeros */
*str++ = 0;
}
/* Iterate over each digit... */
int jx;
/* Unpack digit bytes, high order first */
if (!pos && !x) /* suppress leading zeros */
continue;
}
}
if (!pos)
return MP_OKAY;
} /* end mp_to_fixlen_octets() */
/* }}} */
/*------------------------------------------------------------------------*/
/* HERE THERE BE DRAGONS */