AuroraMiddleware/libtommath
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

format with astyle (step 6)

Browse Source
This commit is contained in:
Francois Perrad 2017-08-30 20:23:46 +02:00
parent e2cd147a46
commit 2344bcea3a
20 changed files with 1742 additions and 1742 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

210
bn_fast_mp_invmod.c
View File

@ -23,124 +23,124 @@
*/
int fast_mp_invmod(mp_int *a, mp_int *b, mp_int *c)
{
mp_int x, y, u, v, B, D;
int res, neg;
mp_int x, y, u, v, B, D;
int res, neg;
/* 2. [modified] b must be odd */
if (mp_iseven(b) == MP_YES) {
return MP_VAL;
}
/* 2. [modified] b must be odd */
if (mp_iseven(b) == MP_YES) {
return MP_VAL;
}
/* init all our temps */
if ((res = mp_init_multi(&x, &y, &u, &v, &B, &D, NULL)) != MP_OKAY) {
return res;
}
/* init all our temps */
if ((res = mp_init_multi(&x, &y, &u, &v, &B, &D, NULL)) != MP_OKAY) {
return res;
}
/* x == modulus, y == value to invert */
if ((res = mp_copy(b, &x)) != MP_OKAY) {
goto LBL_ERR;
}
/* x == modulus, y == value to invert */
if ((res = mp_copy(b, &x)) != MP_OKAY) {
goto LBL_ERR;
}
/* we need y = |a| */
if ((res = mp_mod(a, b, &y)) != MP_OKAY) {
goto LBL_ERR;
}
/* we need y = |a| */
if ((res = mp_mod(a, b, &y)) != MP_OKAY) {
goto LBL_ERR;
}
/* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
if ((res = mp_copy(&x, &u)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_copy(&y, &v)) != MP_OKAY) {
goto LBL_ERR;
}
mp_set(&D, 1);
/* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
if ((res = mp_copy(&x, &u)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_copy(&y, &v)) != MP_OKAY) {
goto LBL_ERR;
}
mp_set(&D, 1);
top:
/* 4. while u is even do */
while (mp_iseven(&u) == MP_YES) {
/* 4.1 u = u/2 */
if ((res = mp_div_2(&u, &u)) != MP_OKAY) {
goto LBL_ERR;
}
/* 4.2 if B is odd then */
if (mp_isodd(&B) == MP_YES) {
if ((res = mp_sub(&B, &x, &B)) != MP_OKAY) {
goto LBL_ERR;
/* 4. while u is even do */
while (mp_iseven(&u) == MP_YES) {
/* 4.1 u = u/2 */
if ((res = mp_div_2(&u, &u)) != MP_OKAY) {
goto LBL_ERR;
}
}
/* B = B/2 */
if ((res = mp_div_2(&B, &B)) != MP_OKAY) {
goto LBL_ERR;
}
}
/* 5. while v is even do */
while (mp_iseven(&v) == MP_YES) {
/* 5.1 v = v/2 */
if ((res = mp_div_2(&v, &v)) != MP_OKAY) {
goto LBL_ERR;
}
/* 5.2 if D is odd then */
if (mp_isodd(&D) == MP_YES) {
/* D = (D-x)/2 */
if ((res = mp_sub(&D, &x, &D)) != MP_OKAY) {
goto LBL_ERR;
/* 4.2 if B is odd then */
if (mp_isodd(&B) == MP_YES) {
if ((res = mp_sub(&B, &x, &B)) != MP_OKAY) {
goto LBL_ERR;
}
}
}
/* D = D/2 */
if ((res = mp_div_2(&D, &D)) != MP_OKAY) {
/* B = B/2 */
if ((res = mp_div_2(&B, &B)) != MP_OKAY) {
goto LBL_ERR;
}
}
/* 5. while v is even do */
while (mp_iseven(&v) == MP_YES) {
/* 5.1 v = v/2 */
if ((res = mp_div_2(&v, &v)) != MP_OKAY) {
goto LBL_ERR;
}
/* 5.2 if D is odd then */
if (mp_isodd(&D) == MP_YES) {
/* D = (D-x)/2 */
if ((res = mp_sub(&D, &x, &D)) != MP_OKAY) {
goto LBL_ERR;
}
}
/* D = D/2 */
if ((res = mp_div_2(&D, &D)) != MP_OKAY) {
goto LBL_ERR;
}
}
/* 6. if u >= v then */
if (mp_cmp(&u, &v) != MP_LT) {
/* u = u - v, B = B - D */
if ((res = mp_sub(&u, &v, &u)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_sub(&B, &D, &B)) != MP_OKAY) {
goto LBL_ERR;
}
} else {
/* v - v - u, D = D - B */
if ((res = mp_sub(&v, &u, &v)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_sub(&D, &B, &D)) != MP_OKAY) {
goto LBL_ERR;
}
}
/* if not zero goto step 4 */
if (mp_iszero(&u) == MP_NO) {
goto top;
}
/* now a = C, b = D, gcd == g*v */
/* if v != 1 then there is no inverse */
if (mp_cmp_d(&v, 1) != MP_EQ) {
res = MP_VAL;
goto LBL_ERR;
}
}
}
/* 6. if u >= v then */
if (mp_cmp(&u, &v) != MP_LT) {
/* u = u - v, B = B - D */
if ((res = mp_sub(&u, &v, &u)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_sub(&B, &D, &B)) != MP_OKAY) {
goto LBL_ERR;
}
} else {
/* v - v - u, D = D - B */
if ((res = mp_sub(&v, &u, &v)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_sub(&D, &B, &D)) != MP_OKAY) {
goto LBL_ERR;
}
}
/* if not zero goto step 4 */
if (mp_iszero(&u) == MP_NO) {
goto top;
}
/* now a = C, b = D, gcd == g*v */
/* if v != 1 then there is no inverse */
if (mp_cmp_d(&v, 1) != MP_EQ) {
res = MP_VAL;
goto LBL_ERR;
}
/* b is now the inverse */
neg = a->sign;
while (D.sign == MP_NEG) {
if ((res = mp_add(&D, b, &D)) != MP_OKAY) {
goto LBL_ERR;
}
}
mp_exch(&D, c);
c->sign = neg;
res = MP_OKAY;
/* b is now the inverse */
neg = a->sign;
while (D.sign == MP_NEG) {
if ((res = mp_add(&D, b, &D)) != MP_OKAY) {
goto LBL_ERR;
}
}
mp_exch(&D, c);
c->sign = neg;
res = MP_OKAY;
LBL_ERR:
mp_clear_multi(&x, &y, &u, &v, &B, &D, NULL);
return res;
mp_clear_multi(&x, &y, &u, &v, &B, &D, NULL);
return res;
}
#endif

246
bn_fast_mp_montgomery_reduce.c
View File

@ -25,145 +25,145 @@
*/
int fast_mp_montgomery_reduce(mp_int *x, mp_int *n, mp_digit rho)
{
int ix, res, olduse;
mp_word W[MP_WARRAY];
int ix, res, olduse;
mp_word W[MP_WARRAY];
/* get old used count */
olduse = x->used;
/* get old used count */
olduse = x->used;
/* grow a as required */
if (x->alloc < (n->used + 1)) {
if ((res = mp_grow(x, n->used + 1)) != MP_OKAY) {
return res;
}
}
/* first we have to get the digits of the input into
* an array of double precision words W[...]
*/
{
mp_word *_W;
mp_digit *tmpx;
/* alias for the W[] array */
_W = W;
/* alias for the digits of x*/
tmpx = x->dp;
/* copy the digits of a into W[0..a->used-1] */
for (ix = 0; ix < x->used; ix++) {
*_W++ = *tmpx++;
}
/* zero the high words of W[a->used..m->used*2] */
for (; ix < ((n->used * 2) + 1); ix++) {
*_W++ = 0;
}
}
/* now we proceed to zero successive digits
* from the least significant upwards
*/
for (ix = 0; ix < n->used; ix++) {
/* mu = ai * m' mod b
*
* We avoid a double precision multiplication (which isn't required)
* by casting the value down to a mp_digit. Note this requires
* that W[ix-1] have the carry cleared (see after the inner loop)
*/
mp_digit mu;
mu = (mp_digit)(((W[ix] & MP_MASK) * rho) & MP_MASK);
/* a = a + mu * m * b**i
*
* This is computed in place and on the fly. The multiplication
* by b**i is handled by offseting which columns the results
* are added to.
*
* Note the comba method normally doesn't handle carries in the
* inner loop In this case we fix the carry from the previous
* column since the Montgomery reduction requires digits of the
* result (so far) [see above] to work. This is
* handled by fixing up one carry after the inner loop. The
* carry fixups are done in order so after these loops the
* first m->used words of W[] have the carries fixed
*/
{
int iy;
mp_digit *tmpn;
mp_word *_W;
/* alias for the digits of the modulus */
tmpn = n->dp;
/* Alias for the columns set by an offset of ix */
_W = W + ix;
/* inner loop */
for (iy = 0; iy < n->used; iy++) {
*_W++ += ((mp_word)mu) * ((mp_word)*tmpn++);
/* grow a as required */
if (x->alloc < (n->used + 1)) {
if ((res = mp_grow(x, n->used + 1)) != MP_OKAY) {
return res;
}
}
}
/* now fix carry for next digit, W[ix+1] */
W[ix + 1] += W[ix] >> ((mp_word) DIGIT_BIT);
}
/* first we have to get the digits of the input into
* an array of double precision words W[...]
*/
{
mp_word *_W;
mp_digit *tmpx;
/* now we have to propagate the carries and
* shift the words downward [all those least
* significant digits we zeroed].
*/
{
mp_digit *tmpx;
mp_word *_W, *_W1;
/* alias for the W[] array */
_W = W;
/* nox fix rest of carries */
/* alias for the digits of x*/
tmpx = x->dp;
/* alias for current word */
_W1 = W + ix;
/* copy the digits of a into W[0..a->used-1] */
for (ix = 0; ix < x->used; ix++) {
*_W++ = *tmpx++;
}
/* alias for next word, where the carry goes */
_W = W + ++ix;
/* zero the high words of W[a->used..m->used*2] */
for (; ix < ((n->used * 2) + 1); ix++) {
*_W++ = 0;
}
}
for (; ix <= ((n->used * 2) + 1); ix++) {
*_W++ += *_W1++ >> ((mp_word) DIGIT_BIT);
}
/* now we proceed to zero successive digits
* from the least significant upwards
*/
for (ix = 0; ix < n->used; ix++) {
/* mu = ai * m' mod b
*
* We avoid a double precision multiplication (which isn't required)
* by casting the value down to a mp_digit. Note this requires
* that W[ix-1] have the carry cleared (see after the inner loop)
*/
mp_digit mu;
mu = (mp_digit)(((W[ix] & MP_MASK) * rho) & MP_MASK);
/* copy out, A = A/b**n
*
* The result is A/b**n but instead of converting from an
* array of mp_word to mp_digit than calling mp_rshd
* we just copy them in the right order
*/
/* a = a + mu * m * b**i
*
* This is computed in place and on the fly. The multiplication
* by b**i is handled by offseting which columns the results
* are added to.
*
* Note the comba method normally doesn't handle carries in the
* inner loop In this case we fix the carry from the previous
* column since the Montgomery reduction requires digits of the
* result (so far) [see above] to work. This is
* handled by fixing up one carry after the inner loop. The
* carry fixups are done in order so after these loops the
* first m->used words of W[] have the carries fixed
*/
{
int iy;
mp_digit *tmpn;
mp_word *_W;
/* alias for destination word */
tmpx = x->dp;
/* alias for the digits of the modulus */
tmpn = n->dp;
/* alias for shifted double precision result */
_W = W + n->used;
/* Alias for the columns set by an offset of ix */
_W = W + ix;
for (ix = 0; ix < (n->used + 1); ix++) {
*tmpx++ = (mp_digit)(*_W++ & ((mp_word) MP_MASK));
}
/* inner loop */
for (iy = 0; iy < n->used; iy++) {
*_W++ += ((mp_word)mu) * ((mp_word)*tmpn++);
}
}
/* zero oldused digits, if the input a was larger than
* m->used+1 we'll have to clear the digits
*/
for (; ix < olduse; ix++) {
*tmpx++ = 0;
}
}
/* now fix carry for next digit, W[ix+1] */
W[ix + 1] += W[ix] >> ((mp_word) DIGIT_BIT);
}
/* set the max used and clamp */
x->used = n->used + 1;
mp_clamp(x);
/* now we have to propagate the carries and
* shift the words downward [all those least
* significant digits we zeroed].
*/
{
mp_digit *tmpx;
mp_word *_W, *_W1;
/* if A >= m then A = A - m */
if (mp_cmp_mag(x, n) != MP_LT) {
return s_mp_sub(x, n, x);
}
return MP_OKAY;
/* nox fix rest of carries */
/* alias for current word */
_W1 = W + ix;
/* alias for next word, where the carry goes */
_W = W + ++ix;
for (; ix <= ((n->used * 2) + 1); ix++) {
*_W++ += *_W1++ >> ((mp_word) DIGIT_BIT);
}
/* copy out, A = A/b**n
*
* The result is A/b**n but instead of converting from an
* array of mp_word to mp_digit than calling mp_rshd
* we just copy them in the right order
*/
/* alias for destination word */
tmpx = x->dp;
/* alias for shifted double precision result */
_W = W + n->used;
for (ix = 0; ix < (n->used + 1); ix++) {
*tmpx++ = (mp_digit)(*_W++ & ((mp_word) MP_MASK));
}
/* zero oldused digits, if the input a was larger than
* m->used+1 we'll have to clear the digits
*/
for (; ix < olduse; ix++) {
*tmpx++ = 0;
}
}
/* set the max used and clamp */
x->used = n->used + 1;
mp_clamp(x);
/* if A >= m then A = A - m */
if (mp_cmp_mag(x, n) != MP_LT) {
return s_mp_sub(x, n, x);
}
return MP_OKAY;
}
#endif

46
bn_mp_add.c
View File

@ -18,32 +18,32 @@
/* high level addition (handles signs) */
int mp_add(mp_int *a, mp_int *b, mp_int *c)
{
int sa, sb, res;
int sa, sb, res;
/* get sign of both inputs */
sa = a->sign;
sb = b->sign;
/* get sign of both inputs */
sa = a->sign;
sb = b->sign;
/* handle two cases, not four */
if (sa == sb) {
/* both positive or both negative */
/* add their magnitudes, copy the sign */
c->sign = sa;
res = s_mp_add(a, b, c);
} else {
/* one positive, the other negative */
/* subtract the one with the greater magnitude from */
/* the one of the lesser magnitude. The result gets */
/* the sign of the one with the greater magnitude. */
if (mp_cmp_mag(a, b) == MP_LT) {
c->sign = sb;
res = s_mp_sub(b, a, c);
} else {
/* handle two cases, not four */
if (sa == sb) {
/* both positive or both negative */
/* add their magnitudes, copy the sign */
c->sign = sa;
res = s_mp_sub(a, b, c);
}
}
return res;
res = s_mp_add(a, b, c);
} else {
/* one positive, the other negative */
/* subtract the one with the greater magnitude from */
/* the one of the lesser magnitude. The result gets */
/* the sign of the one with the greater magnitude. */
if (mp_cmp_mag(a, b) == MP_LT) {
c->sign = sb;
res = s_mp_sub(b, a, c);
} else {
c->sign = sa;
res = s_mp_sub(a, b, c);
}
}
return res;
}
#endif

430
bn_mp_div.c
View File

@ -23,63 +23,63 @@ int mp_div(mp_int *a, mp_int *b, mp_int *c, mp_int *d)
mp_int ta, tb, tq, q;
int res, n, n2;
/* is divisor zero ? */
if (mp_iszero(b) == MP_YES) {
return MP_VAL;
}
/* is divisor zero ? */
if (mp_iszero(b) == MP_YES) {
return MP_VAL;
}
/* if a < b then q=0, r = a */
if (mp_cmp_mag(a, b) == MP_LT) {
if (d != NULL) {
res = mp_copy(a, d);
} else {
res = MP_OKAY;
}
if (c != NULL) {
mp_zero(c);
}
return res;
}
/* if a < b then q=0, r = a */
if (mp_cmp_mag(a, b) == MP_LT) {
if (d != NULL) {
res = mp_copy(a, d);
} else {
res = MP_OKAY;
}
if (c != NULL) {
mp_zero(c);
}
return res;
}
/* init our temps */
if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL)) != MP_OKAY) {
return res;
}
/* init our temps */
if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL)) != MP_OKAY) {
return res;
}
mp_set(&tq, 1);
n = mp_count_bits(a) - mp_count_bits(b);
if (((res = mp_abs(a, &ta)) != MP_OKAY) ||
((res = mp_abs(b, &tb)) != MP_OKAY) ||
((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) ||
((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) {
mp_set(&tq, 1);
n = mp_count_bits(a) - mp_count_bits(b);
if (((res = mp_abs(a, &ta)) != MP_OKAY) ||
((res = mp_abs(b, &tb)) != MP_OKAY) ||
((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) ||
((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) {
goto LBL_ERR;
}
}
while (n-- >= 0) {
if (mp_cmp(&tb, &ta) != MP_GT) {
if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) ||
((res = mp_add(&q, &tq, &q)) != MP_OKAY)) {
goto LBL_ERR;
}
}
if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) ||
((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) {
goto LBL_ERR;
}
}
while (n-- >= 0) {
if (mp_cmp(&tb, &ta) != MP_GT) {
if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) ||
((res = mp_add(&q, &tq, &q)) != MP_OKAY)) {
goto LBL_ERR;
}
}
if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) ||
((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) {
goto LBL_ERR;
}
}
/* now q == quotient and ta == remainder */
n = a->sign;
n2 = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
if (c != NULL) {
mp_exch(c, &q);
c->sign = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2;
}
if (d != NULL) {
mp_exch(d, &ta);
d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n;
}
/* now q == quotient and ta == remainder */
n = a->sign;
n2 = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
if (c != NULL) {
mp_exch(c, &q);
c->sign = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2;
}
if (d != NULL) {
mp_exch(d, &ta);
d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n;
}
LBL_ERR:
mp_clear_multi(&ta, &tb, &tq, &q, NULL);
return res;
@ -100,195 +100,195 @@ LBL_ERR:
* The overall algorithm is as described as
* 14.20 from HAC but fixed to treat these cases.
*/
int mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
int mp_div(mp_int *a, mp_int *b, mp_int *c, mp_int *d)
{
mp_int q, x, y, t1, t2;
int res, n, t, i, norm, neg;
mp_int q, x, y, t1, t2;
int res, n, t, i, norm, neg;
/* is divisor zero ? */
if (mp_iszero(b) == MP_YES) {
return MP_VAL;
}
/* is divisor zero ? */
if (mp_iszero(b) == MP_YES) {
return MP_VAL;
}
/* if a < b then q=0, r = a */
if (mp_cmp_mag(a, b) == MP_LT) {
if (d != NULL) {
res = mp_copy(a, d);
} else {
res = MP_OKAY;
}
if (c != NULL) {
mp_zero(c);
}
return res;
}
if ((res = mp_init_size(&q, a->used + 2)) != MP_OKAY) {
return res;
}
q.used = a->used + 2;
if ((res = mp_init(&t1)) != MP_OKAY) {
goto LBL_Q;
}
if ((res = mp_init(&t2)) != MP_OKAY) {
goto LBL_T1;
}
if ((res = mp_init_copy(&x, a)) != MP_OKAY) {
goto LBL_T2;
}
if ((res = mp_init_copy(&y, b)) != MP_OKAY) {
goto LBL_X;
}
/* fix the sign */
neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
x.sign = y.sign = MP_ZPOS;
/* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
norm = mp_count_bits(&y) % DIGIT_BIT;
if (norm < (int)(DIGIT_BIT-1)) {
norm = (DIGIT_BIT-1) - norm;
if ((res = mp_mul_2d(&x, norm, &x)) != MP_OKAY) {
goto LBL_Y;
}
if ((res = mp_mul_2d(&y, norm, &y)) != MP_OKAY) {
goto LBL_Y;
}
} else {
norm = 0;
}
/* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
n = x.used - 1;
t = y.used - 1;
/* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
if ((res = mp_lshd(&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */
goto LBL_Y;
}
while (mp_cmp(&x, &y) != MP_LT) {
++(q.dp[n - t]);
if ((res = mp_sub(&x, &y, &x)) != MP_OKAY) {
goto LBL_Y;
}
}
/* reset y by shifting it back down */
mp_rshd(&y, n - t);
/* step 3. for i from n down to (t + 1) */
for (i = n; i >= (t + 1); i--) {
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 */
if (x.dp[i] == y.dp[t]) {
q.dp[(i - t) - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1);
} else {
mp_word tmp;
tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT);
tmp |= ((mp_word) x.dp[i - 1]);
tmp /= ((mp_word) y.dp[t]);
if (tmp > (mp_word) MP_MASK) {
tmp = MP_MASK;
/* if a < b then q=0, r = a */
if (mp_cmp_mag(a, b) == MP_LT) {
if (d != NULL) {
res = mp_copy(a, d);
} else {
res = MP_OKAY;
}
q.dp[(i - t) - 1] = (mp_digit)(tmp & (mp_word)(MP_MASK));
}
if (c != NULL) {
mp_zero(c);
}
return res;
}
/* while (q{i-t-1} * (yt * b + y{t-1})) >
xi * b**2 + xi-1 * b + xi-2
if ((res = mp_init_size(&q, a->used + 2)) != MP_OKAY) {
return res;
}
q.used = a->used + 2;
do q{i-t-1} -= 1;
*/
q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] + 1) & MP_MASK;
do {
q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1) & MP_MASK;
if ((res = mp_init(&t1)) != MP_OKAY) {
goto LBL_Q;
}
/* find left hand */
mp_zero(&t1);
t1.dp[0] = ((t - 1) < 0) ? 0 : y.dp[t - 1];
t1.dp[1] = y.dp[t];
t1.used = 2;
if ((res = mp_mul_d(&t1, q.dp[(i - t) - 1], &t1)) != MP_OKAY) {
goto LBL_Y;
if ((res = mp_init(&t2)) != MP_OKAY) {
goto LBL_T1;
}
if ((res = mp_init_copy(&x, a)) != MP_OKAY) {
goto LBL_T2;
}
if ((res = mp_init_copy(&y, b)) != MP_OKAY) {
goto LBL_X;
}
/* fix the sign */
neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
x.sign = y.sign = MP_ZPOS;
/* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
norm = mp_count_bits(&y) % DIGIT_BIT;
if (norm < (int)(DIGIT_BIT-1)) {
norm = (DIGIT_BIT-1) - norm;
if ((res = mp_mul_2d(&x, norm, &x)) != MP_OKAY) {
goto LBL_Y;
}
if ((res = mp_mul_2d(&y, norm, &y)) != MP_OKAY) {
goto LBL_Y;
}
} else {
norm = 0;
}
/* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
n = x.used - 1;
t = y.used - 1;
/* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
if ((res = mp_lshd(&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */
goto LBL_Y;
}
while (mp_cmp(&x, &y) != MP_LT) {
++(q.dp[n - t]);
if ((res = mp_sub(&x, &y, &x)) != MP_OKAY) {
goto LBL_Y;
}
}
/* reset y by shifting it back down */
mp_rshd(&y, n - t);
/* step 3. for i from n down to (t + 1) */
for (i = n; i >= (t + 1); i--) {
if (i > x.used) {
continue;
}
/* find right hand */
t2.dp[0] = ((i - 2) < 0) ? 0 : x.dp[i - 2];
t2.dp[1] = ((i - 1) < 0) ? 0 : x.dp[i - 1];
t2.dp[2] = x.dp[i];
t2.used = 3;
} while (mp_cmp_mag(&t1, &t2) == MP_GT);
/* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
if ((res = mp_mul_d(&y, q.dp[(i - t) - 1], &t1)) != MP_OKAY) {
goto LBL_Y;
}
if ((res = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY) {
goto LBL_Y;
}
if ((res = mp_sub(&x, &t1, &x)) != MP_OKAY) {
goto LBL_Y;
}
/* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
if (x.sign == MP_NEG) {
if ((res = mp_copy(&y, &t1)) != MP_OKAY) {
goto LBL_Y;
/* 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 */
if (x.dp[i] == y.dp[t]) {
q.dp[(i - t) - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1);
} else {
mp_word tmp;
tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT);
tmp |= ((mp_word) x.dp[i - 1]);
tmp /= ((mp_word) y.dp[t]);
if (tmp > (mp_word) MP_MASK) {
tmp = MP_MASK;
}
q.dp[(i - t) - 1] = (mp_digit)(tmp & (mp_word)(MP_MASK));
}
/* while (q{i-t-1} * (yt * b + y{t-1})) >
xi * b**2 + xi-1 * b + xi-2
do q{i-t-1} -= 1;
*/
q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] + 1) & MP_MASK;
do {
q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1) & MP_MASK;
/* find left hand */
mp_zero(&t1);
t1.dp[0] = ((t - 1) < 0) ? 0 : y.dp[t - 1];
t1.dp[1] = y.dp[t];
t1.used = 2;
if ((res = mp_mul_d(&t1, q.dp[(i - t) - 1], &t1)) != MP_OKAY) {
goto LBL_Y;
}
/* find right hand */
t2.dp[0] = ((i - 2) < 0) ? 0 : x.dp[i - 2];
t2.dp[1] = ((i - 1) < 0) ? 0 : x.dp[i - 1];
t2.dp[2] = x.dp[i];
t2.used = 3;
} while (mp_cmp_mag(&t1, &t2) == MP_GT);
/* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
if ((res = mp_mul_d(&y, q.dp[(i - t) - 1], &t1)) != MP_OKAY) {
goto LBL_Y;
}
if ((res = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY) {
goto LBL_Y;
}
if ((res = mp_add(&x, &t1, &x)) != MP_OKAY) {
goto LBL_Y;
goto LBL_Y;
}
q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1UL) & MP_MASK;
}
}
if ((res = mp_sub(&x, &t1, &x)) != MP_OKAY) {
goto LBL_Y;
}
/* now q is the quotient and x is the remainder
* [which we have to normalize]
*/
/* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
if (x.sign == MP_NEG) {
if ((res = mp_copy(&y, &t1)) != MP_OKAY) {
goto LBL_Y;
}
if ((res = mp_lshd(&t1, (i - t) - 1)) != MP_OKAY) {
goto LBL_Y;
}
if ((res = mp_add(&x, &t1, &x)) != MP_OKAY) {
goto LBL_Y;
}
/* get sign before writing to c */
x.sign = (x.used == 0) ? MP_ZPOS : a->sign;
q.dp[(i - t) - 1] = (q.dp[(i - t) - 1] - 1UL) & MP_MASK;
}
}
if (c != NULL) {
mp_clamp(&q);
mp_exch(&q, c);
c->sign = neg;
}
/* now q is the quotient and x is the remainder
* [which we have to normalize]
*/
if (d != NULL) {
if ((res = mp_div_2d(&x, norm, &x, NULL)) != MP_OKAY) {
goto LBL_Y;
}
mp_exch(&x, d);
}
/* get sign before writing to c */
x.sign = (x.used == 0) ? MP_ZPOS : a->sign;
res = MP_OKAY;
if (c != NULL) {
mp_clamp(&q);
mp_exch(&q, c);
c->sign = neg;
}
if (d != NULL) {
if ((res = mp_div_2d(&x, norm, &x, NULL)) != MP_OKAY) {
goto LBL_Y;
}
mp_exch(&x, d);
}
res = MP_OKAY;
LBL_Y:
mp_clear(&y);
mp_clear(&y);
LBL_X:
mp_clear(&x);
mp_clear(&x);
LBL_T2:
mp_clear(&t2);
mp_clear(&t2);
LBL_T1:
mp_clear(&t1);
mp_clear(&t1);
LBL_Q:
mp_clear(&q);
return res;
mp_clear(&q);
return res;
}
#endif

100
bn_mp_div_2d.c
View File

@ -18,66 +18,66 @@
/* shift right by a certain bit count (store quotient in c, optional remainder in d) */
int mp_div_2d(mp_int *a, int b, mp_int *c, mp_int *d)
{
mp_digit D, r, rr;
int x, res;
mp_digit D, r, rr;
int x, res;
/* if the shift count is <= 0 then we do no work */
if (b <= 0) {
res = mp_copy(a, c);
if (d != NULL) {
mp_zero(d);
}
return res;
}
/* copy */
if ((res = mp_copy(a, c)) != MP_OKAY) {
return res;
}
/* 'a' should not be used after here - it might be the same as d */
/* get the remainder */
if (d != NULL) {
if ((res = mp_mod_2d(a, b, d)) != MP_OKAY) {
/* if the shift count is <= 0 then we do no work */
if (b <= 0) {
res = mp_copy(a, c);
if (d != NULL) {
mp_zero(d);
}
return res;
}
}
}
/* shift by as many digits in the bit count */
if (b >= (int)DIGIT_BIT) {
mp_rshd(c, b / DIGIT_BIT);
}
/* copy */
if ((res = mp_copy(a, c)) != MP_OKAY) {
return res;
}
/* 'a' should not be used after here - it might be the same as d */
/* shift any bit count < DIGIT_BIT */
D = (mp_digit)(b % DIGIT_BIT);
if (D != 0) {
mp_digit *tmpc, mask, shift;
/* get the remainder */
if (d != NULL) {
if ((res = mp_mod_2d(a, b, d)) != MP_OKAY) {
return res;
}
}
/* mask */
mask = (((mp_digit)1) << D) - 1;
/* shift by as many digits in the bit count */
if (b >= (int)DIGIT_BIT) {
mp_rshd(c, b / DIGIT_BIT);
}
/* shift for lsb */
shift = DIGIT_BIT - D;
/* shift any bit count < DIGIT_BIT */
D = (mp_digit)(b % DIGIT_BIT);
if (D != 0) {
mp_digit *tmpc, mask, shift;
/* alias */
tmpc = c->dp + (c->used - 1);
/* mask */
mask = (((mp_digit)1) << D) - 1;
/* carry */
r = 0;
for (x = c->used - 1; x >= 0; x--) {
/* get the lower bits of this word in a temp */
rr = *tmpc & mask;
/* shift for lsb */
shift = DIGIT_BIT - D;
/* shift the current word and mix in the carry bits from the previous word */
*tmpc = (*tmpc >> D) | (r << shift);
--tmpc;
/* alias */
tmpc = c->dp + (c->used - 1);
/* set the carry to the carry bits of the current word found above */
r = rr;
}
}
mp_clamp(c);
return MP_OKAY;
/* carry */
r = 0;
for (x = c->used - 1; x >= 0; x--) {
/* get the lower bits of this word in a temp */
rr = *tmpc & mask;
/* shift the current word and mix in the carry bits from the previous word */
*tmpc = (*tmpc >> D) | (r << shift);
--tmpc;
/* set the carry to the carry bits of the current word found above */
r = rr;
}
}
mp_clamp(c);
return MP_OKAY;
}
#endif

96
bn_mp_expt_d_ex.c
View File

@ -18,62 +18,62 @@
/* calculate c = a**b using a square-multiply algorithm */
int mp_expt_d_ex(mp_int *a, mp_digit b, mp_int *c, int fast)
{
int res;
unsigned int x;
int res;
unsigned int x;
mp_int g;
mp_int g;
if ((res = mp_init_copy(&g, a)) != MP_OKAY) {
return res;
}
if ((res = mp_init_copy(&g, a)) != MP_OKAY) {
return res;
}
/* set initial result */
mp_set(c, 1);
/* set initial result */
mp_set(c, 1);
if (fast != 0) {
while (b > 0) {
/* if the bit is set multiply */
if ((b & 1) != 0) {
if ((res = mp_mul(c, &g, c)) != MP_OKAY) {
mp_clear(&g);
return res;
}
if (fast != 0) {
while (b > 0) {
/* if the bit is set multiply */
if ((b & 1) != 0) {
if ((res = mp_mul(c, &g, c)) != MP_OKAY) {
mp_clear(&g);
return res;
}
}
/* square */
if (b > 1) {
if ((res = mp_sqr(&g, &g)) != MP_OKAY) {
mp_clear(&g);
return res;
}
}
/* shift to next bit */
b >>= 1;
}
} else {
for (x = 0; x < DIGIT_BIT; x++) {
/* square */
if ((res = mp_sqr(c, c)) != MP_OKAY) {
mp_clear(&g);
return res;
}
/* square */
if (b > 1) {
if ((res = mp_sqr(&g, &g)) != MP_OKAY) {
mp_clear(&g);
return res;
}
/* if the bit is set multiply */
if ((b & (mp_digit)(((mp_digit)1) << (DIGIT_BIT - 1))) != 0) {
if ((res = mp_mul(c, &g, c)) != MP_OKAY) {
mp_clear(&g);
return res;
}
}
/* shift to next bit */
b <<= 1;
}
} /* if ... else */
/* shift to next bit */
b >>= 1;
}
} else {
for (x = 0; x < DIGIT_BIT; x++) {
/* square */
if ((res = mp_sqr(c, c)) != MP_OKAY) {
mp_clear(&g);
return res;
}
/* if the bit is set multiply */
if ((b & (mp_digit)(((mp_digit)1) << (DIGIT_BIT - 1))) != 0) {
if ((res = mp_mul(c, &g, c)) != MP_OKAY) {
mp_clear(&g);
return res;
}
}
/* shift to next bit */
b <<= 1;
}
} /* if ... else */
mp_clear(&g);
return MP_OKAY;
mp_clear(&g);
return MP_OKAY;
}
#endif

482
bn_mp_exptmod_fast.c
View File

@ -31,288 +31,288 @@
int mp_exptmod_fast(mp_int *G, mp_int *X, mp_int *P, mp_int *Y, int redmode)
{
mp_int M[TAB_SIZE], res;
mp_digit buf, mp;
int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
mp_int M[TAB_SIZE], res;
mp_digit buf, mp;
int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
/* use a pointer to the reduction algorithm. This allows us to use
* one of many reduction algorithms without modding the guts of
* the code with if statements everywhere.
*/
int (*redux)(mp_int*,mp_int*,mp_digit);
/* use a pointer to the reduction algorithm. This allows us to use
* one of many reduction algorithms without modding the guts of
* the code with if statements everywhere.
*/
int (*redux)(mp_int *,mp_int *,mp_digit);
/* find window size */
x = mp_count_bits(X);
if (x <= 7) {
winsize = 2;
} else if (x <= 36) {
winsize = 3;
} else if (x <= 140) {
winsize = 4;
} else if (x <= 450) {
winsize = 5;
} else if (x <= 1303) {
winsize = 6;
} else if (x <= 3529) {
winsize = 7;
} else {
winsize = 8;
}
/* find window size */
x = mp_count_bits(X);
if (x <= 7) {
winsize = 2;
} else if (x <= 36) {
winsize = 3;
} else if (x <= 140) {
winsize = 4;
} else if (x <= 450) {
winsize = 5;
} else if (x <= 1303) {
winsize = 6;
} else if (x <= 3529) {
winsize = 7;
} else {
winsize = 8;
}
#ifdef MP_LOW_MEM
if (winsize > 5) {
winsize = 5;
}
if (winsize > 5) {
winsize = 5;
}
#endif
/* init M array */
/* init first cell */
if ((err = mp_init_size(&M[1], P->alloc)) != MP_OKAY) {
return err;
}
/* now init the second half of the array */
for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
if ((err = mp_init_size(&M[x], P->alloc)) != MP_OKAY) {
for (y = 1<<(winsize-1); y < x; y++) {
mp_clear(&M[y]);
}
mp_clear(&M[1]);
/* init M array */
/* init first cell */
if ((err = mp_init_size(&M[1], P->alloc)) != MP_OKAY) {
return err;
}
}
}
/* determine and setup reduction code */
if (redmode == 0) {
/* now init the second half of the array */
for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
if ((err = mp_init_size(&M[x], P->alloc)) != MP_OKAY) {
for (y = 1<<(winsize-1); y < x; y++) {
mp_clear(&M[y]);
}
mp_clear(&M[1]);
return err;
}
}
/* determine and setup reduction code */
if (redmode == 0) {
#ifdef BN_MP_MONTGOMERY_SETUP_C
/* now setup montgomery */
if ((err = mp_montgomery_setup(P, &mp)) != MP_OKAY) {
goto LBL_M;
}
/* now setup montgomery */
if ((err = mp_montgomery_setup(P, &mp)) != MP_OKAY) {
goto LBL_M;
}
#else
err = MP_VAL;
goto LBL_M;
err = MP_VAL;
goto LBL_M;
#endif
/* automatically pick the comba one if available (saves quite a few calls/ifs) */
/* automatically pick the comba one if available (saves quite a few calls/ifs) */
#ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C
if ((((P->used * 2) + 1) < MP_WARRAY) &&
if ((((P->used * 2) + 1) < MP_WARRAY) &&
(P->used < (1 << ((CHAR_BIT * sizeof(mp_word)) - (2 * DIGIT_BIT))))) {
redux = fast_mp_montgomery_reduce;
} else
redux = fast_mp_montgomery_reduce;
} else
#endif
{
{
#ifdef BN_MP_MONTGOMERY_REDUCE_C
/* use slower baseline Montgomery method */
redux = mp_montgomery_reduce;
/* use slower baseline Montgomery method */
redux = mp_montgomery_reduce;
#else
err = MP_VAL;
goto LBL_M;
err = MP_VAL;
goto LBL_M;
#endif
}
} else if (redmode == 1) {
}
} else if (redmode == 1) {
#if defined(BN_MP_DR_SETUP_C) && defined(BN_MP_DR_REDUCE_C)
/* setup DR reduction for moduli of the form B**k - b */
mp_dr_setup(P, &mp);
redux = mp_dr_reduce;
/* setup DR reduction for moduli of the form B**k - b */
mp_dr_setup(P, &mp);
redux = mp_dr_reduce;
#else
err = MP_VAL;
goto LBL_M;
err = MP_VAL;
goto LBL_M;
#endif
} else {
} else {
#if defined(BN_MP_REDUCE_2K_SETUP_C) && defined(BN_MP_REDUCE_2K_C)
/* setup DR reduction for moduli of the form 2**k - b */
if ((err = mp_reduce_2k_setup(P, &mp)) != MP_OKAY) {
goto LBL_M;
}
redux = mp_reduce_2k;
/* setup DR reduction for moduli of the form 2**k - b */
if ((err = mp_reduce_2k_setup(P, &mp)) != MP_OKAY) {
goto LBL_M;
}
redux = mp_reduce_2k;
#else
err = MP_VAL;
goto LBL_M;
err = MP_VAL;
goto LBL_M;
#endif
}
}
/* setup result */
if ((err = mp_init_size(&res, P->alloc)) != MP_OKAY) {
goto LBL_M;
}
/* setup result */
if ((err = mp_init_size(&res, P->alloc)) != MP_OKAY) {
goto LBL_M;
}
/* create M table
*
/* create M table
*
*
* The first half of the table is not computed though accept for M[0] and M[1]
*/
*
* The first half of the table is not computed though accept for M[0] and M[1]
*/
if (redmode == 0) {
if (redmode == 0) {
#ifdef BN_MP_MONTGOMERY_CALC_NORMALIZATION_C
/* now we need R mod m */
if ((err = mp_montgomery_calc_normalization(&res, P)) != MP_OKAY) {
goto LBL_RES;
}
/* now we need R mod m */
if ((err = mp_montgomery_calc_normalization(&res, P)) != MP_OKAY) {
goto LBL_RES;
}
/* now set M[1] to G * R mod m */
if ((err = mp_mulmod(G, &res, P, &M[1])) != MP_OKAY) {
goto LBL_RES;
}
/* now set M[1] to G * R mod m */
if ((err = mp_mulmod(G, &res, P, &M[1])) != MP_OKAY) {
goto LBL_RES;
}
#else
err = MP_VAL;
goto LBL_RES;
err = MP_VAL;
goto LBL_RES;
#endif
} else {
mp_set(&res, 1);
if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) {
goto LBL_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 LBL_RES;
}
for (x = 0; x < (winsize - 1); x++) {
if ((err = mp_sqr(&M[1 << (winsize - 1)], &M[1 << (winsize - 1)])) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&M[1 << (winsize - 1)], P, mp)) != MP_OKAY) {
goto LBL_RES;
}
}
/* create upper table */
for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
if ((err = mp_mul(&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&M[x], P, mp)) != MP_OKAY) {
goto LBL_RES;
}
}
/* set initial mode and bit cnt */
mode = 0;
bitcnt = 1;
buf = 0;
digidx = X->used - 1;
bitcpy = 0;
bitbuf = 0;
for (;;) {
/* grab next digit as required */
if (--bitcnt == 0) {
/* if digidx == -1 we are out of digits so break */
if (digidx == -1) {
break;
} else {
mp_set(&res, 1);
if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) {
goto LBL_RES;
}
/* read next digit and reset bitcnt */
buf = X->dp[digidx--];
bitcnt = (int)DIGIT_BIT;
}
}
/* grab the next msb from the exponent */
y = (mp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
buf <<= (mp_digit)1;
/* 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 LBL_RES;
}
/* if the bit is zero and mode == 0 then we ignore it
* These represent the leading zero bits before the first 1 bit
* in the exponent. Technically this opt is not required but it
* does lower the # of trivial squaring/reductions used
*/
if ((mode == 0) && (y == 0)) {
continue;
}
/* if the bit is zero and mode == 1 then we square */
if ((mode == 1) && (y == 0)) {
if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
goto LBL_RES;
for (x = 0; x < (winsize - 1); x++) {
if ((err = mp_sqr(&M[1 << (winsize - 1)], &M[1 << (winsize - 1)])) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&M[1 << (winsize - 1)], P, mp)) != MP_OKAY) {
goto LBL_RES;
}
}
/* create upper table */
for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
if ((err = mp_mul(&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&M[x], P, mp)) != MP_OKAY) {
goto LBL_RES;
}
}
/* set initial mode and bit cnt */
mode = 0;
bitcnt = 1;
buf = 0;
digidx = X->used - 1;
bitcpy = 0;
bitbuf = 0;
for (;;) {
/* grab next digit as required */
if (--bitcnt == 0) {
/* if digidx == -1 we are out of digits so break */
if (digidx == -1) {
break;
}
/* read next digit and reset bitcnt */
buf = X->dp[digidx--];
bitcnt = (int)DIGIT_BIT;
}
/* grab the next msb from the exponent */
y = (mp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
buf <<= (mp_digit)1;
/* if the bit is zero and mode == 0 then we ignore it
* These represent the leading zero bits before the first 1 bit
* in the exponent. Technically this opt is not required but it
* does lower the # of trivial squaring/reductions used
*/
if ((mode == 0) && (y == 0)) {
continue;
}
/* if the bit is zero and mode == 1 then we square */
if ((mode == 1) && (y == 0)) {
if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, mp)) != MP_OKAY) {
goto LBL_RES;
}
continue;
}
/* else we add it to the window */
bitbuf |= (y << (winsize - ++bitcpy));
mode = 2;
if (bitcpy == winsize) {
/* ok window is filled so square as required and multiply */
/* square first */
for (x = 0; x < winsize; x++) {
if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, mp)) != MP_OKAY) {
goto LBL_RES;
}
}
/* then multiply */
if ((err = mp_mul(&res, &M[bitbuf], &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, mp)) != MP_OKAY) {
goto LBL_RES;
}
/* empty window and reset */
bitcpy = 0;
bitbuf = 0;
mode = 1;
}
}
/* if bits remain then square/multiply */
if ((mode == 2) && (bitcpy > 0)) {
/* square then multiply if the bit is set */
for (x = 0; x < bitcpy; x++) {
if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, mp)) != MP_OKAY) {
goto LBL_RES;
}
/* get next bit of the window */
bitbuf <<= 1;
if ((bitbuf & (1 << winsize)) != 0) {
/* then multiply */
if ((err = mp_mul(&res, &M[1], &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, mp)) != MP_OKAY) {
goto LBL_RES;
}
}
}
}
if (redmode == 0) {
/* fixup result if Montgomery reduction is used
* recall that any value in a Montgomery system is
* actually multiplied by R mod n. So we have
* to reduce one more time to cancel out the factor
* of R.
*/
if ((err = redux(&res, P, mp)) != MP_OKAY) {
goto LBL_RES;
goto LBL_RES;
}
continue;
}
}
/* else we add it to the window */
bitbuf |= (y << (winsize - ++bitcpy));
mode = 2;
if (bitcpy == winsize) {
/* ok window is filled so square as required and multiply */
/* square first */
for (x = 0; x < winsize; x++) {
if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, mp)) != MP_OKAY) {
goto LBL_RES;
}
}
/* then multiply */
if ((err = mp_mul(&res, &M[bitbuf], &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, mp)) != MP_OKAY) {
goto LBL_RES;
}
/* empty window and reset */
bitcpy = 0;
bitbuf = 0;
mode = 1;
}
}
/* if bits remain then square/multiply */
if ((mode == 2) && (bitcpy > 0)) {
/* square then multiply if the bit is set */
for (x = 0; x < bitcpy; x++) {
if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, mp)) != MP_OKAY) {
goto LBL_RES;
}
/* get next bit of the window */
bitbuf <<= 1;
if ((bitbuf & (1 << winsize)) != 0) {
/* then multiply */
if ((err = mp_mul(&res, &M[1], &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, mp)) != MP_OKAY) {
goto LBL_RES;
}
}
}
}
if (redmode == 0) {
/* fixup result if Montgomery reduction is used
* recall that any value in a Montgomery system is
* actually multiplied by R mod n. So we have
* to reduce one more time to cancel out the factor
* of R.
*/
if ((err = redux(&res, P, mp)) != MP_OKAY) {
goto LBL_RES;
}
}
/* swap res with Y */
mp_exch(&res, Y);
err = MP_OKAY;
/* swap res with Y */
mp_exch(&res, Y);
err = MP_OKAY;
LBL_RES:
mp_clear(&res);
mp_clear(&res);
LBL_M:
mp_clear(&M[1]);
for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
mp_clear(&M[x]);
}
return err;
mp_clear(&M[1]);
for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
mp_clear(&M[x]);
}
return err;
}
#endif

246
bn_mp_invmod_slow.c
View File

@ -18,156 +18,156 @@
/* hac 14.61, pp608 */
int mp_invmod_slow(mp_int *a, mp_int *b, mp_int *c)
{
mp_int x, y, u, v, A, B, C, D;
int res;
mp_int x, y, u, v, A, B, C, D;
int res;
/* b cannot be negative */
if ((b->sign == MP_NEG) || (mp_iszero(b) == MP_YES)) {
return MP_VAL;
}
/* b cannot be negative */
if ((b->sign == MP_NEG) || (mp_iszero(b) == MP_YES)) {
return MP_VAL;
}
/* init temps */
if ((res = mp_init_multi(&x, &y, &u, &v,
&A, &B, &C, &D, NULL)) != MP_OKAY) {
return res;
}
/* init temps */
if ((res = mp_init_multi(&x, &y, &u, &v,
&A, &B, &C, &D, NULL)) != MP_OKAY) {
return res;
}
/* x = a, y = b */
if ((res = mp_mod(a, b, &x)) != MP_OKAY) {
/* x = a, y = b */
if ((res = mp_mod(a, b, &x)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_copy(b, &y)) != MP_OKAY) {
goto LBL_ERR;
}
}
if ((res = mp_copy(b, &y)) != MP_OKAY) {
goto LBL_ERR;
}
/* 2. [modified] if x,y are both even then return an error! */
if ((mp_iseven(&x) == MP_YES) && (mp_iseven(&y) == MP_YES)) {
res = MP_VAL;
goto LBL_ERR;
}
/* 2. [modified] if x,y are both even then return an error! */
if ((mp_iseven(&x) == MP_YES) && (mp_iseven(&y) == MP_YES)) {
res = MP_VAL;
goto LBL_ERR;
}
/* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
if ((res = mp_copy(&x, &u)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_copy(&y, &v)) != MP_OKAY) {
goto LBL_ERR;
}
mp_set(&A, 1);
mp_set(&D, 1);
/* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
if ((res = mp_copy(&x, &u)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_copy(&y, &v)) != MP_OKAY) {
goto LBL_ERR;
}
mp_set(&A, 1);
mp_set(&D, 1);
top:
/* 4. while u is even do */
while (mp_iseven(&u) == MP_YES) {
/* 4.1 u = u/2 */
if ((res = mp_div_2(&u, &u)) != MP_OKAY) {
goto LBL_ERR;
}
/* 4.2 if A or B is odd then */
if ((mp_isodd(&A) == MP_YES) || (mp_isodd(&B) == MP_YES)) {
/* A = (A+y)/2, B = (B-x)/2 */
if ((res = mp_add(&A, &y, &A)) != MP_OKAY) {
/* 4. while u is even do */
while (mp_iseven(&u) == MP_YES) {
/* 4.1 u = u/2 */
if ((res = mp_div_2(&u, &u)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_sub(&B, &x, &B)) != MP_OKAY) {
/* 4.2 if A or B is odd then */
if ((mp_isodd(&A) == MP_YES) || (mp_isodd(&B) == MP_YES)) {
/* A = (A+y)/2, B = (B-x)/2 */
if ((res = mp_add(&A, &y, &A)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_sub(&B, &x, &B)) != MP_OKAY) {
goto LBL_ERR;
}
}
/* A = A/2, B = B/2 */
if ((res = mp_div_2(&A, &A)) != MP_OKAY) {
goto LBL_ERR;
}
}
/* A = A/2, B = B/2 */
if ((res = mp_div_2(&A, &A)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_div_2(&B, &B)) != MP_OKAY) {
goto LBL_ERR;
}
}
/* 5. while v is even do */
while (mp_iseven(&v) == MP_YES) {
/* 5.1 v = v/2 */
if ((res = mp_div_2(&v, &v)) != MP_OKAY) {
goto LBL_ERR;
}
/* 5.2 if C or D is odd then */
if ((mp_isodd(&C) == MP_YES) || (mp_isodd(&D) == MP_YES)) {
/* C = (C+y)/2, D = (D-x)/2 */
if ((res = mp_add(&C, &y, &C)) != MP_OKAY) {
if ((res = mp_div_2(&B, &B)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_sub(&D, &x, &D)) != MP_OKAY) {
}
/* 5. while v is even do */
while (mp_iseven(&v) == MP_YES) {
/* 5.1 v = v/2 */
if ((res = mp_div_2(&v, &v)) != MP_OKAY) {
goto LBL_ERR;
}
}
/* C = C/2, D = D/2 */
if ((res = mp_div_2(&C, &C)) != MP_OKAY) {
/* 5.2 if C or D is odd then */
if ((mp_isodd(&C) == MP_YES) || (mp_isodd(&D) == MP_YES)) {
/* C = (C+y)/2, D = (D-x)/2 */
if ((res = mp_add(&C, &y, &C)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_sub(&D, &x, &D)) != MP_OKAY) {
goto LBL_ERR;
}
}
/* C = C/2, D = D/2 */
if ((res = mp_div_2(&C, &C)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_div_2(&D, &D)) != MP_OKAY) {
goto LBL_ERR;
}
}
/* 6. if u >= v then */
if (mp_cmp(&u, &v) != MP_LT) {
/* u = u - v, A = A - C, B = B - D */
if ((res = mp_sub(&u, &v, &u)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_sub(&A, &C, &A)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_sub(&B, &D, &B)) != MP_OKAY) {
goto LBL_ERR;
}
} else {
/* v - v - u, C = C - A, D = D - B */
if ((res = mp_sub(&v, &u, &v)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_sub(&C, &A, &C)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_sub(&D, &B, &D)) != MP_OKAY) {
goto LBL_ERR;
}
}
/* if not zero goto step 4 */
if (mp_iszero(&u) == MP_NO)
goto top;
/* now a = C, b = D, gcd == g*v */
/* if v != 1 then there is no inverse */
if (mp_cmp_d(&v, 1) != MP_EQ) {
res = MP_VAL;
goto LBL_ERR;
}
if ((res = mp_div_2(&D, &D)) != MP_OKAY) {
goto LBL_ERR;
}
}
}
/* 6. if u >= v then */
if (mp_cmp(&u, &v) != MP_LT) {
/* u = u - v, A = A - C, B = B - D */
if ((res = mp_sub(&u, &v, &u)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_sub(&A, &C, &A)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_sub(&B, &D, &B)) != MP_OKAY) {
goto LBL_ERR;
}
} else {
/* v - v - u, C = C - A, D = D - B */
if ((res = mp_sub(&v, &u, &v)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_sub(&C, &A, &C)) != MP_OKAY) {
goto LBL_ERR;
}
if ((res = mp_sub(&D, &B, &D)) != MP_OKAY) {
goto LBL_ERR;
}
}
/* if not zero goto step 4 */
if (mp_iszero(&u) == MP_NO)
goto top;
/* now a = C, b = D, gcd == g*v */
/* if v != 1 then there is no inverse */
if (mp_cmp_d(&v, 1) != MP_EQ) {
res = MP_VAL;
goto LBL_ERR;
}
/* if its too low */
while (mp_cmp_d(&C, 0) == MP_LT) {
/* if its too low */
while (mp_cmp_d(&C, 0) == MP_LT) {
if ((res = mp_add(&C, b, &C)) != MP_OKAY) {
goto LBL_ERR;
}
}
}
/* too big */
while (mp_cmp_mag(&C, b) != MP_LT) {
/* too big */
while (mp_cmp_mag(&C, b) != MP_LT) {
if ((res = mp_sub(&C, b, &C)) != MP_OKAY) {
goto LBL_ERR;
}
}
}
/* C is now the inverse */
mp_exch(&C, c);
res = MP_OKAY;
/* C is now the inverse */
mp_exch(&C, c);
res = MP_OKAY;
LBL_ERR:
mp_clear_multi(&x, &y, &u, &v, &A, &B, &C, &D, NULL);
return res;
mp_clear_multi(&x, &y, &u, &v, &A, &B, &C, &D, NULL);
return res;
}
#endif

148
bn_mp_jacobi.c
View File

@ -22,95 +22,95 @@
*/
int mp_jacobi(mp_int *a, mp_int *n, int *c)
{
mp_int a1, p1;
int k, s, r, res;
mp_digit residue;
mp_int a1, p1;
int k, s, r, res;
mp_digit residue;
/* if a < 0 return MP_VAL */
if (mp_isneg(a) == MP_YES) {
return MP_VAL;
}
/* if a < 0 return MP_VAL */
if (mp_isneg(a) == MP_YES) {
return MP_VAL;
}
/* if n <= 0 return MP_VAL */
if (mp_cmp_d(n, 0) != MP_GT) {
return MP_VAL;
}
/* if n <= 0 return MP_VAL */
if (mp_cmp_d(n, 0) != MP_GT) {
return MP_VAL;
}
/* step 1. handle case of a == 0 */
if (mp_iszero(a) == MP_YES) {
/* special case of a == 0 and n == 1 */
if (mp_cmp_d(n, 1) == MP_EQ) {
*c = 1;
} else {
*c = 0;
}
return MP_OKAY;
}
/* step 1. handle case of a == 0 */
if (mp_iszero(a) == MP_YES) {
/* special case of a == 0 and n == 1 */
if (mp_cmp_d(n, 1) == MP_EQ) {
*c = 1;
} else {
*c = 0;
}
return MP_OKAY;
}
/* step 2. if a == 1, return 1 */
if (mp_cmp_d(a, 1) == MP_EQ) {
*c = 1;
return MP_OKAY;
}
/* step 2. if a == 1, return 1 */
if (mp_cmp_d(a, 1) == MP_EQ) {
*c = 1;
return MP_OKAY;
}
/* default */
s = 0;
/* default */
s = 0;
/* step 3. write a = a1 * 2**k */
if ((res = mp_init_copy(&a1, a)) != MP_OKAY) {
return res;
}
/* step 3. write a = a1 * 2**k */
if ((res = mp_init_copy(&a1, a)) != MP_OKAY) {
return res;
}
if ((res = mp_init(&p1)) != MP_OKAY) {
goto LBL_A1;
}
if ((res = mp_init(&p1)) != MP_OKAY) {
goto LBL_A1;
}
/* divide out larger power of two */
k = mp_cnt_lsb(&a1);
if ((res = mp_div_2d(&a1, k, &a1, NULL)) != MP_OKAY) {
goto LBL_P1;
}
/* divide out larger power of two */
k = mp_cnt_lsb(&a1);
if ((res = mp_div_2d(&a1, k, &a1, NULL)) != MP_OKAY) {
goto LBL_P1;
}
/* step 4. if e is even set s=1 */
if ((k & 1) == 0) {
s = 1;
} else {
/* else set s=1 if p = 1/7 (mod 8) or s=-1 if p = 3/5 (mod 8) */
residue = n->dp[0] & 7;
if ((residue == 1) || (residue == 7)) {
/* step 4. if e is even set s=1 */
if ((k & 1) == 0) {
s = 1;
} else if ((residue == 3) || (residue == 5)) {
s = -1;
}
}
} else {
/* else set s=1 if p = 1/7 (mod 8) or s=-1 if p = 3/5 (mod 8) */
residue = n->dp[0] & 7;
/* step 5. if p == 3 (mod 4) *and* a1 == 3 (mod 4) then s = -s */
if ( ((n->dp[0] & 3) == 3) && ((a1.dp[0] & 3) == 3)) {
s = -s;
}
if ((residue == 1) || (residue == 7)) {
s = 1;
} else if ((residue == 3) || (residue == 5)) {
s = -1;
}
}
/* if a1 == 1 we're done */
if (mp_cmp_d(&a1, 1) == MP_EQ) {
*c = s;
} else {
/* n1 = n mod a1 */
if ((res = mp_mod(n, &a1, &p1)) != MP_OKAY) {
goto LBL_P1;
}
if ((res = mp_jacobi(&p1, &a1, &r)) != MP_OKAY) {
goto LBL_P1;
}
*c = s * r;
}
/* step 5. if p == 3 (mod 4) *and* a1 == 3 (mod 4) then s = -s */
if (((n->dp[0] & 3) == 3) && ((a1.dp[0] & 3) == 3)) {
s = -s;
}
/* done */
res = MP_OKAY;
/* if a1 == 1 we're done */
if (mp_cmp_d(&a1, 1) == MP_EQ) {
*c = s;
} else {
/* n1 = n mod a1 */
if ((res = mp_mod(n, &a1, &p1)) != MP_OKAY) {
goto LBL_P1;
}
if ((res = mp_jacobi(&p1, &a1, &r)) != MP_OKAY) {
goto LBL_P1;
}
*c = s * r;
}
/* done */
res = MP_OKAY;
LBL_P1:
mp_clear(&p1);
mp_clear(&p1);
LBL_A1:
mp_clear(&a1);
return res;
mp_clear(&a1);
return res;
}
#endif

48
bn_mp_montgomery_calc_normalization.c
View File

@ -23,34 +23,34 @@
*/
int mp_montgomery_calc_normalization(mp_int *a, mp_int *b)
{
int x, bits, res;
int x, bits, res;
/* how many bits of last digit does b use */
bits = mp_count_bits(b) % DIGIT_BIT;
/* how many bits of last digit does b use */
bits = mp_count_bits(b) % DIGIT_BIT;
if (b->used > 1) {
if ((res = mp_2expt(a, ((b->used - 1) * DIGIT_BIT) + bits - 1)) != MP_OKAY) {
return res;
}
} else {
mp_set(a, 1);
bits = 1;
}
/* now compute C = A * B mod b */
for (x = bits - 1; x < (int)DIGIT_BIT; x++) {
if ((res = mp_mul_2(a, a)) != MP_OKAY) {
return res;
}
if (mp_cmp_mag(a, b) != MP_LT) {
if ((res = s_mp_sub(a, b, a)) != MP_OKAY) {
return res;
if (b->used > 1) {
if ((res = mp_2expt(a, ((b->used - 1) * DIGIT_BIT) + bits - 1)) != MP_OKAY) {
return res;
}
}
}
} else {
mp_set(a, 1);
bits = 1;
}
return MP_OKAY;
/* now compute C = A * B mod b */
for (x = bits - 1; x < (int)DIGIT_BIT; x++) {
if ((res = mp_mul_2(a, a)) != MP_OKAY) {
return res;
}
if (mp_cmp_mag(a, b) != MP_LT) {
if ((res = s_mp_sub(a, b, a)) != MP_OKAY) {
return res;
}
}
}
return MP_OKAY;
}
#endif

162
bn_mp_montgomery_reduce.c
View File

@ -18,97 +18,97 @@
/* computes xR**-1 == x (mod N) via Montgomery Reduction */
int mp_montgomery_reduce(mp_int *x, mp_int *n, mp_digit rho)
{
int ix, res, digs;
mp_digit mu;
int ix, res, digs;
mp_digit mu;
/* can the fast reduction [comba] method be used?
*
* Note that unlike in mul you're safely allowed *less*
* than the available columns [255 per default] since carries
* are fixed up in the inner loop.
*/
digs = (n->used * 2) + 1;
if ((digs < MP_WARRAY) &&
(n->used <
(1 << ((CHAR_BIT * sizeof(mp_word)) - (2 * DIGIT_BIT))))) {
return fast_mp_montgomery_reduce(x, n, rho);
}
/* can the fast reduction [comba] method be used?
*
* Note that unlike in mul you're safely allowed *less*
* than the available columns [255 per default] since carries
* are fixed up in the inner loop.
*/
digs = (n->used * 2) + 1;
if ((digs < MP_WARRAY) &&
(n->used <
(1 << ((CHAR_BIT * sizeof(mp_word)) - (2 * DIGIT_BIT))))) {
return fast_mp_montgomery_reduce(x, n, rho);
}
/* grow the input as required */
if (x->alloc < digs) {
if ((res = mp_grow(x, digs)) != MP_OKAY) {
return res;
}
}
x->used = digs;
for (ix = 0; ix < n->used; ix++) {
/* mu = ai * rho mod b
*
* The value of rho must be precalculated via
* montgomery_setup() such that
* it equals -1/n0 mod b this allows the
* following inner loop to reduce the
* input one digit at a time
*/
mu = (mp_digit)(((mp_word)x->dp[ix] * (mp_word)rho) & MP_MASK);
/* a = a + mu * m * b**i */
{
int iy;
mp_digit *tmpn, *tmpx, u;
mp_word r;
/* alias for digits of the modulus */
tmpn = n->dp;
/* alias for the digits of x [the input] */
tmpx = x->dp + ix;
/* set the carry to zero */
u = 0;
/* Multiply and add in place */
for (iy = 0; iy < n->used; iy++) {
/* compute product and sum */
r = ((mp_word)mu * (mp_word)*tmpn++) +
(mp_word) u + (mp_word) *tmpx;
/* get carry */
u = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
/* fix digit */
*tmpx++ = (mp_digit)(r & ((mp_word) MP_MASK));
/* grow the input as required */
if (x->alloc < digs) {
if ((res = mp_grow(x, digs)) != MP_OKAY) {
return res;
}
/* At this point the ix'th digit of x should be zero */
}
x->used = digs;
for (ix = 0; ix < n->used; ix++) {
/* mu = ai * rho mod b
*
* The value of rho must be precalculated via
* montgomery_setup() such that
* it equals -1/n0 mod b this allows the
* following inner loop to reduce the
* input one digit at a time
*/
mu = (mp_digit)(((mp_word)x->dp[ix] * (mp_word)rho) & MP_MASK);
/* a = a + mu * m * b**i */
{
int iy;
mp_digit *tmpn, *tmpx, u;
mp_word r;
/* alias for digits of the modulus */
tmpn = n->dp;
/* alias for the digits of x [the input] */
tmpx = x->dp + ix;
/* set the carry to zero */
u = 0;
/* Multiply and add in place */
for (iy = 0; iy < n->used; iy++) {
/* compute product and sum */
r = ((mp_word)mu * (mp_word)*tmpn++) +
(mp_word) u + (mp_word) *tmpx;
/* get carry */
u = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
/* fix digit */
*tmpx++ = (mp_digit)(r & ((mp_word) MP_MASK));
}
/* At this point the ix'th digit of x should be zero */
/* propagate carries upwards as required*/
while (u != 0) {
*tmpx += u;
u = *tmpx >> DIGIT_BIT;
*tmpx++ &= MP_MASK;
/* propagate carries upwards as required*/
while (u != 0) {
*tmpx += u;
u = *tmpx >> DIGIT_BIT;
*tmpx++ &= MP_MASK;
}
}
}
}
}
/* at this point the n.used'th least
* significant digits of x are all zero
* which means we can shift x to the
* right by n.used digits and the
* residue is unchanged.
*/
/* at this point the n.used'th least
* significant digits of x are all zero
* which means we can shift x to the
* right by n.used digits and the
* residue is unchanged.
*/
/* x = x/b**n.used */
mp_clamp(x);
mp_rshd(x, n->used);
/* x = x/b**n.used */
mp_clamp(x);
mp_rshd(x, n->used);
/* if x >= n then x = x - n */
if (mp_cmp_mag(x, n) != MP_LT) {
return s_mp_sub(x, n, x);
}
/* if x >= n then x = x - n */
if (mp_cmp_mag(x, n) != MP_LT) {
return s_mp_sub(x, n, x);
}
return MP_OKAY;
return MP_OKAY;
}
#endif

166
bn_mp_n_root_ex.c
View File

@ -27,106 +27,106 @@
*/
int mp_n_root_ex(mp_int *a, mp_digit b, mp_int *c, int fast)
{
mp_int t1, t2, t3;
int res, neg;
mp_int t1, t2, t3;
int res, neg;
/* input must be positive if b is even */
if (((b & 1) == 0) && (a->sign == MP_NEG)) {
return MP_VAL;
}
/* input must be positive if b is even */
if (((b & 1) == 0) && (a->sign == MP_NEG)) {
return MP_VAL;
}
if ((res = mp_init(&t1)) != MP_OKAY) {
return res;
}
if ((res = mp_init(&t1)) != MP_OKAY) {
return res;
}
if ((res = mp_init(&t2)) != MP_OKAY) {
goto LBL_T1;
}
if ((res = mp_init(&t2)) != MP_OKAY) {
goto LBL_T1;
}
if ((res = mp_init(&t3)) != MP_OKAY) {
goto LBL_T2;
}
if ((res = mp_init(&t3)) != MP_OKAY) {
goto LBL_T2;
}
/* if a is negative fudge the sign but keep track */
neg = a->sign;
a->sign = MP_ZPOS;
/* if a is negative fudge the sign but keep track */
neg = a->sign;
a->sign = MP_ZPOS;
/* t2 = 2 */
mp_set(&t2, 2);
/* t2 = 2 */
mp_set(&t2, 2);
do {
/* t1 = t2 */
if ((res = mp_copy(&t2, &t1)) != MP_OKAY) {
goto LBL_T3;
}
/* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */
/* t3 = t1**(b-1) */
if ((res = mp_expt_d_ex(&t1, b - 1, &t3, fast)) != MP_OKAY) {
goto LBL_T3;
}
/* numerator */
/* t2 = t1**b */
if ((res = mp_mul(&t3, &t1, &t2)) != MP_OKAY) {
goto LBL_T3;
}
/* t2 = t1**b - a */
if ((res = mp_sub(&t2, a, &t2)) != MP_OKAY) {
goto LBL_T3;
}
/* denominator */
/* t3 = t1**(b-1) * b */
if ((res = mp_mul_d(&t3, b, &t3)) != MP_OKAY) {
goto LBL_T3;
}
/* t3 = (t1**b - a)/(b * t1**(b-1)) */
if ((res = mp_div(&t2, &t3, &t3, NULL)) != MP_OKAY) {
goto LBL_T3;
}
if ((res = mp_sub(&t1, &t3, &t2)) != MP_OKAY) {
goto LBL_T3;
}
} while (mp_cmp(&t1, &t2) != MP_EQ);
/* result can be off by a few so check */
for (;;) {
if ((res = mp_expt_d_ex(&t1, b, &t2, fast)) != MP_OKAY) {
goto LBL_T3;
}
if (mp_cmp(&t2, a) == MP_GT) {
if ((res = mp_sub_d(&t1, 1, &t1)) != MP_OKAY) {
do {
/* t1 = t2 */
if ((res = mp_copy(&t2, &t1)) != MP_OKAY) {
goto LBL_T3;
}
} else {
break;
}
}
/* reset the sign of a first */
a->sign = neg;
/* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */
/* set the result */
mp_exch(&t1, c);
/* t3 = t1**(b-1) */
if ((res = mp_expt_d_ex(&t1, b - 1, &t3, fast)) != MP_OKAY) {
goto LBL_T3;
}
/* set the sign of the result */
c->sign = neg;
/* numerator */
/* t2 = t1**b */
if ((res = mp_mul(&t3, &t1, &t2)) != MP_OKAY) {
goto LBL_T3;
}
res = MP_OKAY;
/* t2 = t1**b - a */
if ((res = mp_sub(&t2, a, &t2)) != MP_OKAY) {
goto LBL_T3;
}
/* denominator */
/* t3 = t1**(b-1) * b */
if ((res = mp_mul_d(&t3, b, &t3)) != MP_OKAY) {
goto LBL_T3;
}
/* t3 = (t1**b - a)/(b * t1**(b-1)) */
if ((res = mp_div(&t2, &t3, &t3, NULL)) != MP_OKAY) {
goto LBL_T3;
}
if ((res = mp_sub(&t1, &t3, &t2)) != MP_OKAY) {
goto LBL_T3;
}
} while (mp_cmp(&t1, &t2) != MP_EQ);
/* result can be off by a few so check */
for (;;) {
if ((res = mp_expt_d_ex(&t1, b, &t2, fast)) != MP_OKAY) {
goto LBL_T3;
}
if (mp_cmp(&t2, a) == MP_GT) {
if ((res = mp_sub_d(&t1, 1, &t1)) != MP_OKAY) {
goto LBL_T3;
}
} else {
break;
}
}
/* reset the sign of a first */
a->sign = neg;
/* set the result */
mp_exch(&t1, c);
/* set the sign of the result */
c->sign = neg;
res = MP_OKAY;
LBL_T3:
mp_clear(&t3);
mp_clear(&t3);
LBL_T2:
mp_clear(&t2);
mp_clear(&t2);
LBL_T1:
mp_clear(&t1);
return res;
mp_clear(&t1);
return res;
}
#endif

124
bn_mp_prime_miller_rabin.c
View File

@ -24,80 +24,80 @@
*/
int mp_prime_miller_rabin(mp_int *a, mp_int *b, int *result)
{
mp_int n1, y, r;
int s, j, err;
mp_int n1, y, r;
int s, j, err;
/* default */
*result = MP_NO;
/* default */
*result = MP_NO;
/* ensure b > 1 */
if (mp_cmp_d(b, 1) != MP_GT) {
return MP_VAL;
}
/* ensure b > 1 */
if (mp_cmp_d(b, 1) != MP_GT) {
return MP_VAL;
}
/* 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 LBL_N1;
}
/* 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 LBL_N1;
}
/* set 2**s * r = n1 */
if ((err = mp_init_copy(&r, &n1)) != MP_OKAY) {
goto LBL_N1;
}
/* set 2**s * r = n1 */
if ((err = mp_init_copy(&r, &n1)) != MP_OKAY) {
goto LBL_N1;
}
/* count the number of least significant bits
* which are zero
*/
s = mp_cnt_lsb(&r);
/* count the number of least significant bits
* which are zero
*/
s = mp_cnt_lsb(&r);
/* now divide n - 1 by 2**s */
if ((err = mp_div_2d(&r, s, &r, NULL)) != MP_OKAY) {
goto LBL_R;
}
/* now divide n - 1 by 2**s */
if ((err = mp_div_2d(&r, s, &r, NULL)) != MP_OKAY) {
goto LBL_R;
}
/* compute y = b**r mod a */
if ((err = mp_init(&y)) != MP_OKAY) {
goto LBL_R;
}
if ((err = mp_exptmod(b, &r, a, &y)) != MP_OKAY) {
goto LBL_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 LBL_Y;
}
/* if y == 1 then composite */
if (mp_cmp_d(&y, 1) == MP_EQ) {
goto LBL_Y;
}
++j;
}
/* if y != n1 then composite */
if (mp_cmp(&y, &n1) != MP_EQ) {
/* compute y = b**r mod a */
if ((err = mp_init(&y)) != MP_OKAY) {
goto LBL_R;
}
if ((err = mp_exptmod(b, &r, a, &y)) != MP_OKAY) {
goto LBL_Y;
}
}
}
/* probably prime now */
*result = MP_YES;
/* 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 LBL_Y;
}
/* if y == 1 then composite */
if (mp_cmp_d(&y, 1) == MP_EQ) {
goto LBL_Y;
}
++j;
}
/* if y != n1 then composite */
if (mp_cmp(&y, &n1) != MP_EQ) {
goto LBL_Y;
}
}
/* probably prime now */
*result = MP_YES;
LBL_Y:
mp_clear(&y);
mp_clear(&y);
LBL_R:
mp_clear(&r);
mp_clear(&r);
LBL_N1:
mp_clear(&n1);
return err;
mp_clear(&n1);
return err;
}
#endif

72
bn_mp_rand.c
View File

@ -16,13 +16,13 @@
*/
#if MP_GEN_RANDOM_MAX == 0xffffffff
#define MP_GEN_RANDOM_SHIFT 32
#define MP_GEN_RANDOM_SHIFT 32
#elif MP_GEN_RANDOM_MAX == 32767
/* SHRT_MAX */
#define MP_GEN_RANDOM_SHIFT 15
/* SHRT_MAX */
#define MP_GEN_RANDOM_SHIFT 15
#elif MP_GEN_RANDOM_MAX == 2147483647
/* INT_MAX */
#define MP_GEN_RANDOM_SHIFT 31
/* INT_MAX */
#define MP_GEN_RANDOM_SHIFT 31
#elif !defined(MP_GEN_RANDOM_SHIFT)
#error Thou shalt define their own valid MP_GEN_RANDOM_SHIFT
#endif
@ -30,47 +30,47 @@
/* makes a pseudo-random int of a given size */
static mp_digit s_gen_random(void)
{
mp_digit d = 0, msk = 0;
do {
d <<= MP_GEN_RANDOM_SHIFT;
d |= ((mp_digit) MP_GEN_RANDOM());
msk <<= MP_GEN_RANDOM_SHIFT;
msk |= (MP_MASK & MP_GEN_RANDOM_MAX);
} while ((MP_MASK & msk) != MP_MASK);
d &= MP_MASK;
return d;
mp_digit d = 0, msk = 0;
do {
d <<= MP_GEN_RANDOM_SHIFT;
d |= ((mp_digit) MP_GEN_RANDOM());
msk <<= MP_GEN_RANDOM_SHIFT;
msk |= (MP_MASK & MP_GEN_RANDOM_MAX);
} while ((MP_MASK & msk) != MP_MASK);
d &= MP_MASK;
return d;
}
int mp_rand(mp_int *a, int digits)
{
int res;
mp_digit d;
int res;
mp_digit d;
mp_zero(a);
if (digits <= 0) {
return MP_OKAY;
}
mp_zero(a);
if (digits <= 0) {
return MP_OKAY;
}
/* first place a random non-zero digit */
do {
d = s_gen_random();
} while (d == 0);
/* first place a random non-zero digit */
do {
d = s_gen_random();
} while (d == 0);
if ((res = mp_add_d(a, d, a)) != MP_OKAY) {
return res;
}
while (--digits > 0) {
if ((res = mp_lshd(a, 1)) != MP_OKAY) {
if ((res = mp_add_d(a, d, a)) != MP_OKAY) {
return res;
}
}
if ((res = mp_add_d(a, s_gen_random(), a)) != MP_OKAY) {
return res;
}
}
while (--digits > 0) {
if ((res = mp_lshd(a, 1)) != MP_OKAY) {
return res;
}
return MP_OKAY;
if ((res = mp_add_d(a, s_gen_random(), a)) != MP_OKAY) {
return res;
}
}
return MP_OKAY;
}
#endif

108
bn_mp_read_radix.c
View File

@ -18,71 +18,71 @@
/* read a string [ASCII] in a given radix */
int mp_read_radix(mp_int *a, const char *str, int radix)
{
int y, res, neg;
char ch;
int y, res, neg;
char ch;
/* zero the digit bignum */
mp_zero(a);
/* zero the digit bignum */
mp_zero(a);
/* make sure the radix is ok */
if ((radix < 2) || (radix > 64)) {
return MP_VAL;
}
/* make sure the radix is ok */
if ((radix < 2) || (radix > 64)) {
return MP_VAL;
}
/* if the leading digit is a
* minus set the sign to negative.
*/
if (*str == '-') {
++str;
neg = MP_NEG;
} else {
neg = MP_ZPOS;
}
/* if the leading digit is a
* minus set the sign to negative.
*/
if (*str == '-') {
++str;
neg = MP_NEG;
} else {
neg = MP_ZPOS;
}
/* set the integer to the default of zero */
mp_zero(a);
/* set the integer to the default of zero */
mp_zero(a);
/* process each digit of the string */
while (*str != '\0') {
/* if the radix <= 36 the conversion is case insensitive
* this allows numbers like 1AB and 1ab to represent the same value
* [e.g. in hex]
*/
ch = (radix <= 36) ? (char)toupper((int)*str) : *str;
for (y = 0; y < 64; y++) {
if (ch == mp_s_rmap[y]) {
/* process each digit of the string */
while (*str != '\0') {
/* if the radix <= 36 the conversion is case insensitive
* this allows numbers like 1AB and 1ab to represent the same value
* [e.g. in hex]
*/
ch = (radix <= 36) ? (char)toupper((int)*str) : *str;
for (y = 0; y < 64; y++) {
if (ch == mp_s_rmap[y]) {
break;
}
}
/* if the char was found in the map
* and is less than the given radix add it
* to the number, otherwise exit the loop.
*/
if (y < radix) {
if ((res = mp_mul_d(a, (mp_digit)radix, a)) != MP_OKAY) {
return res;
}
if ((res = mp_add_d(a, (mp_digit)y, a)) != MP_OKAY) {
return res;
}
} else {
break;
}
}
++str;
}
/* if the char was found in the map
* and is less than the given radix add it
* to the number, otherwise exit the loop.
*/
if (y < radix) {
if ((res = mp_mul_d(a, (mp_digit)radix, a)) != MP_OKAY) {
return res;
}
if ((res = mp_add_d(a, (mp_digit)y, a)) != MP_OKAY) {
return res;
}
} else {
break;
}
++str;
}
/* if an illegal character was found, fail. */
if (!(*str == '\0' || *str == '\r' || *str == '\n')) {
/* if an illegal character was found, fail. */
if (!(*str == '\0' || *str == '\r' || *str == '\n')) {
mp_zero(a);
return MP_VAL;
}
}
/* set the sign only if a != 0 */
if (mp_iszero(a) != MP_YES) {
a->sign = neg;
}
return MP_OKAY;
/* set the sign only if a != 0 */
if (mp_iszero(a) != MP_YES) {
a->sign = neg;
}
return MP_OKAY;
}
#endif

112
bn_mp_reduce.c
View File

@ -21,77 +21,77 @@
*/
int mp_reduce(mp_int *x, mp_int *m, mp_int *mu)
{
mp_int q;
int res, um = m->used;
mp_int q;
int res, um = m->used;
/* q = x */
if ((res = mp_init_copy(&q, x)) != MP_OKAY) {
return res;
}
/* q = x */
if ((res = mp_init_copy(&q, x)) != MP_OKAY) {
return res;
}
/* q1 = x / b**(k-1) */
mp_rshd(&q, um - 1);
/* q1 = x / b**(k-1) */
mp_rshd(&q, um - 1);
/* according to HAC this optimization is ok */
if (((mp_digit) um) > (((mp_digit)1) << (DIGIT_BIT - 1))) {
if ((res = mp_mul(&q, mu, &q)) != MP_OKAY) {
goto CLEANUP;
}
} else {
/* according to HAC this optimization is ok */
if (((mp_digit) um) > (((mp_digit)1) << (DIGIT_BIT - 1))) {
if ((res = mp_mul(&q, mu, &q)) != MP_OKAY) {
goto CLEANUP;
}
} else {
#ifdef BN_S_MP_MUL_HIGH_DIGS_C
if ((res = s_mp_mul_high_digs(&q, mu, &q, um)) != MP_OKAY) {
goto CLEANUP;
}
if ((res = s_mp_mul_high_digs(&q, mu, &q, um)) != MP_OKAY) {
goto CLEANUP;
}
#elif defined(BN_FAST_S_MP_MUL_HIGH_DIGS_C)
if ((res = fast_s_mp_mul_high_digs(&q, mu, &q, um)) != MP_OKAY) {
goto CLEANUP;
}
if ((res = fast_s_mp_mul_high_digs(&q, mu, &q, um)) != MP_OKAY) {
goto CLEANUP;
}
#else
{
res = MP_VAL;
goto CLEANUP;
}
{
res = MP_VAL;
goto CLEANUP;
}
#endif
}
}
/* q3 = q2 / b**(k+1) */
mp_rshd(&q, um + 1);
/* q3 = q2 / b**(k+1) */
mp_rshd(&q, um + 1);
/* x = x mod b**(k+1), quick (no division) */
if ((res = mp_mod_2d(x, DIGIT_BIT * (um + 1), x)) != MP_OKAY) {
goto CLEANUP;
}
/* q = q * m mod b**(k+1), quick (no division) */
if ((res = s_mp_mul_digs(&q, m, &q, um + 1)) != MP_OKAY) {
goto CLEANUP;
}
/* x = x - q */
if ((res = mp_sub(x, &q, x)) != MP_OKAY) {
goto CLEANUP;
}
/* If x < 0, add b**(k+1) to it */
if (mp_cmp_d(x, 0) == MP_LT) {
mp_set(&q, 1);
if ((res = mp_lshd(&q, um + 1)) != MP_OKAY)
/* x = x mod b**(k+1), quick (no division) */
if ((res = mp_mod_2d(x, DIGIT_BIT * (um + 1), x)) != MP_OKAY) {
goto CLEANUP;
if ((res = mp_add(x, &q, x)) != MP_OKAY)
goto CLEANUP;
}
}
/* Back off if it's too big */
while (mp_cmp(x, m) != MP_LT) {
if ((res = s_mp_sub(x, m, x)) != MP_OKAY) {
/* q = q * m mod b**(k+1), quick (no division) */
if ((res = s_mp_mul_digs(&q, m, &q, um + 1)) != MP_OKAY) {
goto CLEANUP;
}
}
}
/* x = x - q */
if ((res = mp_sub(x, &q, x)) != MP_OKAY) {
goto CLEANUP;
}
/* If x < 0, add b**(k+1) to it */
if (mp_cmp_d(x, 0) == MP_LT) {
mp_set(&q, 1);
if ((res = mp_lshd(&q, um + 1)) != MP_OKAY)
goto CLEANUP;
if ((res = mp_add(x, &q, x)) != MP_OKAY)
goto CLEANUP;
}
/* Back off if it's too big */
while (mp_cmp(x, m) != MP_LT) {
if ((res = s_mp_sub(x, m, x)) != MP_OKAY) {
goto CLEANUP;
}
}
CLEANUP:
mp_clear(&q);
mp_clear(&q);
return res;
return res;
}
#endif

86
bn_mp_sqrt.c
View File

@ -18,62 +18,62 @@
/* this function is less generic than mp_n_root, simpler and faster */
int mp_sqrt(mp_int *arg, mp_int *ret)
{
int res;
mp_int t1, t2;
int res;
mp_int t1, t2;
/* must be positive */
if (arg->sign == MP_NEG) {
return MP_VAL;
}
/* must be positive */
if (arg->sign == MP_NEG) {
return MP_VAL;
}
/* easy out */
if (mp_iszero(arg) == MP_YES) {
mp_zero(ret);
return MP_OKAY;
}
/* easy out */
if (mp_iszero(arg) == MP_YES) {
mp_zero(ret);
return MP_OKAY;
}
if ((res = mp_init_copy(&t1, arg)) != MP_OKAY) {
return res;
}
if ((res = mp_init_copy(&t1, arg)) != MP_OKAY) {
return res;
}
if ((res = mp_init(&t2)) != MP_OKAY) {
goto E2;
}
if ((res = mp_init(&t2)) != MP_OKAY) {
goto E2;
}
/* First approx. (not very bad for large arg) */
mp_rshd(&t1, t1.used/2);
/* First approx. (not very bad for large arg) */
mp_rshd(&t1, t1.used/2);
/* t1 > 0 */
if ((res = mp_div(arg, &t1, &t2, NULL)) != MP_OKAY) {
goto E1;
}
if ((res = mp_add(&t1, &t2, &t1)) != MP_OKAY) {
goto E1;
}
if ((res = mp_div_2(&t1, &t1)) != MP_OKAY) {
goto E1;
}
/* And now t1 > sqrt(arg) */
do {
if ((res = mp_div(arg, &t1, &t2, NULL)) != MP_OKAY) {
/* t1 > 0 */
if ((res = mp_div(arg, &t1, &t2, NULL)) != MP_OKAY) {
goto E1;
}
if ((res = mp_add(&t1, &t2, &t1)) != MP_OKAY) {
}
if ((res = mp_add(&t1, &t2, &t1)) != MP_OKAY) {
goto E1;
}
if ((res = mp_div_2(&t1, &t1)) != MP_OKAY) {
}
if ((res = mp_div_2(&t1, &t1)) != MP_OKAY) {
goto E1;
}
/* t1 >= sqrt(arg) >= t2 at this point */
} while (mp_cmp_mag(&t1, &t2) == MP_GT);
}
/* And now t1 > sqrt(arg) */
do {
if ((res = mp_div(arg, &t1, &t2, NULL)) != MP_OKAY) {
goto E1;
}
if ((res = mp_add(&t1, &t2, &t1)) != MP_OKAY) {
goto E1;
}
if ((res = mp_div_2(&t1, &t1)) != MP_OKAY) {
goto E1;
}
/* t1 >= sqrt(arg) >= t2 at this point */
} while (mp_cmp_mag(&t1, &t2) == MP_GT);
mp_exch(&t1, ret);
mp_exch(&t1, ret);
E1:
mp_clear(&t2);
mp_clear(&t2);
E2:
mp_clear(&t1);
return res;
mp_clear(&t1);
return res;
}
#endif

56
bn_mp_sub.c
View File

@ -18,37 +18,37 @@
/* high level subtraction (handles signs) */
int mp_sub(mp_int *a, mp_int *b, mp_int *c)
{
int sa, sb, res;
int sa, sb, res;
sa = a->sign;
sb = b->sign;
sa = a->sign;
sb = b->sign;
if (sa != sb) {
/* subtract a negative from a positive, OR */
/* subtract a positive from a negative. */
/* In either case, ADD their magnitudes, */
/* and use the sign of the first number. */
c->sign = sa;
res = s_mp_add(a, b, c);
} else {
/* subtract a positive from a positive, OR */
/* subtract a negative from a negative. */
/* First, take the difference between their */
/* magnitudes, then... */
if (mp_cmp_mag(a, b) != MP_LT) {
/* Copy the sign from the first */
if (sa != sb) {
/* subtract a negative from a positive, OR */
/* subtract a positive from a negative. */
/* In either case, ADD their magnitudes, */
/* and use the sign of the first number. */
c->sign = sa;
/* The first has a larger or equal magnitude */
res = s_mp_sub(a, b, c);
} else {
/* The result has the *opposite* sign from */
/* the first number. */
c->sign = (sa == MP_ZPOS) ? MP_NEG : MP_ZPOS;
/* The second has a larger magnitude */
res = s_mp_sub(b, a, c);
}
}
return res;
res = s_mp_add(a, b, c);
} else {
/* subtract a positive from a positive, OR */
/* subtract a negative from a negative. */
/* First, take the difference between their */
/* magnitudes, then... */
if (mp_cmp_mag(a, b) != MP_LT) {
/* Copy the sign from the first */
c->sign = sa;
/* The first has a larger or equal magnitude */
res = s_mp_sub(a, b, c);
} else {
/* The result has the *opposite* sign from */
/* the first number. */
c->sign = (sa == MP_ZPOS) ? MP_NEG : MP_ZPOS;
/* The second has a larger magnitude */
res = s_mp_sub(b, a, c);
}
}
return res;
}
#endif

152
bn_s_mp_add.c
View File

@ -18,88 +18,88 @@
/* low level addition, based on HAC pp.594, Algorithm 14.7 */
int s_mp_add(mp_int *a, mp_int *b, mp_int *c)
{
mp_int *x;
int olduse, res, min, max;
mp_int *x;
int olduse, res, min, max;
/* find sizes, we let |a| <= |b| which means we have to sort
* them. "x" will point to the input with the most digits
*/
if (a->used > b->used) {
min = b->used;
max = a->used;
x = a;
} else {
min = a->used;
max = b->used;
x = b;
}
/* find sizes, we let |a| <= |b| which means we have to sort
* them. "x" will point to the input with the most digits
*/
if (a->used > b->used) {
min = b->used;
max = a->used;
x = a;
} else {
min = a->used;
max = b->used;
x = b;
}
/* init result */
if (c->alloc < (max + 1)) {
if ((res = mp_grow(c, max + 1)) != MP_OKAY) {
return res;
}
}
/* get old used digit count and set new one */
olduse = c->used;
c->used = max + 1;
{
mp_digit u, *tmpa, *tmpb, *tmpc;
int i;
/* alias for digit pointers */
/* first input */
tmpa = a->dp;
/* second input */
tmpb = b->dp;
/* destination */
tmpc = c->dp;
/* zero the carry */
u = 0;
for (i = 0; i < min; i++) {
/* Compute the sum at one digit, T[i] = A[i] + B[i] + U */
*tmpc = *tmpa++ + *tmpb++ + u;
/* U = carry bit of T[i] */
u = *tmpc >> ((mp_digit)DIGIT_BIT);
/* take away carry bit from T[i] */
*tmpc++ &= MP_MASK;
}
/* now copy higher words if any, that is in A+B
* if A or B has more digits add those in
*/
if (min != max) {
for (; i < max; i++) {
/* T[i] = X[i] + U */
*tmpc = x->dp[i] + u;
/* U = carry bit of T[i] */
u = *tmpc >> ((mp_digit)DIGIT_BIT);
/* take away carry bit from T[i] */
*tmpc++ &= MP_MASK;
/* init result */
if (c->alloc < (max + 1)) {
if ((res = mp_grow(c, max + 1)) != MP_OKAY) {
return res;
}
}
}
/* add carry */
*tmpc++ = u;
/* get old used digit count and set new one */
olduse = c->used;
c->used = max + 1;
/* clear digits above oldused */
for (i = c->used; i < olduse; i++) {
*tmpc++ = 0;
}
}
{
mp_digit u, *tmpa, *tmpb, *tmpc;
int i;
mp_clamp(c);
return MP_OKAY;
/* alias for digit pointers */
/* first input */
tmpa = a->dp;
/* second input */
tmpb = b->dp;
/* destination */
tmpc = c->dp;
/* zero the carry */
u = 0;
for (i = 0; i < min; i++) {
/* Compute the sum at one digit, T[i] = A[i] + B[i] + U */
*tmpc = *tmpa++ + *tmpb++ + u;
/* U = carry bit of T[i] */
u = *tmpc >> ((mp_digit)DIGIT_BIT);
/* take away carry bit from T[i] */
*tmpc++ &= MP_MASK;
}
/* now copy higher words if any, that is in A+B
* if A or B has more digits add those in
*/
if (min != max) {
for (; i < max; i++) {
/* T[i] = X[i] + U */
*tmpc = x->dp[i] + u;
/* U = carry bit of T[i] */
u = *tmpc >> ((mp_digit)DIGIT_BIT);
/* take away carry bit from T[i] */
*tmpc++ &= MP_MASK;
}
}
/* add carry */
*tmpc++ = u;
/* clear digits above oldused */
for (i = c->used; i < olduse; i++) {
*tmpc++ = 0;
}
}
mp_clamp(c);
return MP_OKAY;
}
#endif

394
bn_s_mp_exptmod.c
View File

@ -22,230 +22,230 @@
int s_mp_exptmod(mp_int *G, mp_int *X, mp_int *P, mp_int *Y, int redmode)
{
mp_int M[TAB_SIZE], res, mu;
mp_digit buf;
int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
int (*redux)(mp_int*,mp_int*,mp_int*);
mp_int M[TAB_SIZE], res, mu;
mp_digit buf;
int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
int (*redux)(mp_int *,mp_int *,mp_int *);
/* find window size */
x = mp_count_bits(X);
if (x <= 7) {
winsize = 2;
} else if (x <= 36) {
winsize = 3;
} else if (x <= 140) {
winsize = 4;
} else if (x <= 450) {
winsize = 5;
} else if (x <= 1303) {
winsize = 6;
} else if (x <= 3529) {
winsize = 7;
} else {
winsize = 8;
}
/* find window size */
x = mp_count_bits(X);
if (x <= 7) {
winsize = 2;
} else if (x <= 36) {
winsize = 3;
} else if (x <= 140) {
winsize = 4;
} else if (x <= 450) {
winsize = 5;
} else if (x <= 1303) {
winsize = 6;
} else if (x <= 3529) {
winsize = 7;
} else {
winsize = 8;
}
#ifdef MP_LOW_MEM
if (winsize > 5) {
winsize = 5;
}
if (winsize > 5) {
winsize = 5;
}
#endif
/* init M array */
/* init first cell */
if ((err = mp_init(&M[1])) != MP_OKAY) {
return err;
}
/* now init the second half of the array */
for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
if ((err = mp_init(&M[x])) != MP_OKAY) {
for (y = 1<<(winsize-1); y < x; y++) {
mp_clear(&M[y]);
}
mp_clear(&M[1]);
/* init M array */
/* init first cell */
if ((err = mp_init(&M[1])) != MP_OKAY) {
return err;
}
}
}
/* create mu, used for Barrett reduction */
if ((err = mp_init(&mu)) != MP_OKAY) {
goto LBL_M;
}
/* now init the second half of the array */
for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
if ((err = mp_init(&M[x])) != MP_OKAY) {
for (y = 1<<(winsize-1); y < x; y++) {
mp_clear(&M[y]);
}
mp_clear(&M[1]);
return err;
}
}
if (redmode == 0) {
if ((err = mp_reduce_setup(&mu, P)) != MP_OKAY) {
goto LBL_MU;
}
redux = mp_reduce;
} else {
if ((err = mp_reduce_2k_setup_l(P, &mu)) != MP_OKAY) {
goto LBL_MU;
}
redux = mp_reduce_2k_l;
}
/* create mu, used for Barrett reduction */
if ((err = mp_init(&mu)) != MP_OKAY) {
goto LBL_M;
}
/* create M table
*
* The M table contains powers of the base,
* e.g. M[x] = G**x mod P
*
* The first half of the table is not
* computed though accept for M[0] and M[1]
*/
if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) {
goto LBL_MU;
}
if (redmode == 0) {
if ((err = mp_reduce_setup(&mu, P)) != MP_OKAY) {
goto LBL_MU;
}
redux = mp_reduce;
} else {
if ((err = mp_reduce_2k_setup_l(P, &mu)) != MP_OKAY) {
goto LBL_MU;
}
redux = mp_reduce_2k_l;
}
/* 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 LBL_MU;
}
for (x = 0; x < (winsize - 1); x++) {
/* square it */
if ((err = mp_sqr(&M[1 << (winsize - 1)],
&M[1 << (winsize - 1)])) != MP_OKAY) {
/* create M table
*
* The M table contains powers of the base,
* e.g. M[x] = G**x mod P
*
* The first half of the table is not
* computed though accept for M[0] and M[1]
*/
if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) {
goto LBL_MU;
}
}
/* reduce modulo P */
if ((err = redux(&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
/* 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 LBL_MU;
}
}
}
/* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
* for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
*/
for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
if ((err = mp_mul(&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
for (x = 0; x < (winsize - 1); x++) {
/* square it */
if ((err = mp_sqr(&M[1 << (winsize - 1)],
&M[1 << (winsize - 1)])) != MP_OKAY) {
goto LBL_MU;
}
/* reduce modulo P */
if ((err = redux(&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
goto LBL_MU;
}
}
/* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
* for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
*/
for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
if ((err = mp_mul(&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
goto LBL_MU;
}
if ((err = redux(&M[x], P, &mu)) != MP_OKAY) {
goto LBL_MU;
}
}
/* setup result */
if ((err = mp_init(&res)) != MP_OKAY) {
goto LBL_MU;
}
if ((err = redux(&M[x], P, &mu)) != MP_OKAY) {
goto LBL_MU;
}
}
}
mp_set(&res, 1);
/* setup result */
if ((err = mp_init(&res)) != MP_OKAY) {
goto LBL_MU;
}
mp_set(&res, 1);
/* set initial mode and bit cnt */
mode = 0;
bitcnt = 1;
buf = 0;
digidx = X->used - 1;
bitcpy = 0;
bitbuf = 0;
/* set initial mode and bit cnt */
mode = 0;
bitcnt = 1;
buf = 0;
digidx = X->used - 1;
bitcpy = 0;
bitbuf = 0;
for (;;) {
/* grab next digit as required */
if (--bitcnt == 0) {
/* if digidx == -1 we are out of digits */
if (digidx == -1) {
break;
}
/* read next digit and reset the bitcnt */
buf = X->dp[digidx--];
bitcnt = (int)DIGIT_BIT;
}
/* grab the next msb from the exponent */
y = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1;
buf <<= (mp_digit)1;
/* if the bit is zero and mode == 0 then we ignore it
* These represent the leading zero bits before the first 1 bit
* in the exponent. Technically this opt is not required but it
* does lower the # of trivial squaring/reductions used
*/
if ((mode == 0) && (y == 0)) {
continue;
}
/* if the bit is zero and mode == 1 then we square */
if ((mode == 1) && (y == 0)) {
if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, &mu)) != MP_OKAY) {
goto LBL_RES;
}
continue;
}
/* else we add it to the window */
bitbuf |= (y << (winsize - ++bitcpy));
mode = 2;
if (bitcpy == winsize) {
/* ok window is filled so square as required and multiply */
/* square first */
for (x = 0; x < winsize; x++) {
if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, &mu)) != MP_OKAY) {
goto LBL_RES;
}
for (;;) {
/* grab next digit as required */
if (--bitcnt == 0) {
/* if digidx == -1 we are out of digits */
if (digidx == -1) {
break;
}
/* read next digit and reset the bitcnt */
buf = X->dp[digidx--];
bitcnt = (int)DIGIT_BIT;
}
/* then multiply */
if ((err = mp_mul(&res, &M[bitbuf], &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, &mu)) != MP_OKAY) {
goto LBL_RES;
/* grab the next msb from the exponent */
y = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1;
buf <<= (mp_digit)1;
/* if the bit is zero and mode == 0 then we ignore it
* These represent the leading zero bits before the first 1 bit
* in the exponent. Technically this opt is not required but it
* does lower the # of trivial squaring/reductions used
*/
if ((mode == 0) && (y == 0)) {
continue;
}
/* empty window and reset */
bitcpy = 0;
bitbuf = 0;
mode = 1;
}
}
/* if bits remain then square/multiply */
if ((mode == 2) && (bitcpy > 0)) {
/* square then multiply if the bit is set */
for (x = 0; x < bitcpy; x++) {
if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, &mu)) != MP_OKAY) {
goto LBL_RES;
/* if the bit is zero and mode == 1 then we square */
if ((mode == 1) && (y == 0)) {
if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, &mu)) != MP_OKAY) {
goto LBL_RES;
}
continue;
}
bitbuf <<= 1;
if ((bitbuf & (1 << winsize)) != 0) {
/* then multiply */
if ((err = mp_mul(&res, &M[1], &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, &mu)) != MP_OKAY) {
goto LBL_RES;
}
}
}
}
/* else we add it to the window */
bitbuf |= (y << (winsize - ++bitcpy));
mode = 2;
mp_exch(&res, Y);
err = MP_OKAY;
if (bitcpy == winsize) {
/* ok window is filled so square as required and multiply */
/* square first */
for (x = 0; x < winsize; x++) {
if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, &mu)) != MP_OKAY) {
goto LBL_RES;
}
}
/* then multiply */
if ((err = mp_mul(&res, &M[bitbuf], &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, &mu)) != MP_OKAY) {
goto LBL_RES;
}
/* empty window and reset */
bitcpy = 0;
bitbuf = 0;
mode = 1;
}
}
/* if bits remain then square/multiply */
if ((mode == 2) && (bitcpy > 0)) {
/* square then multiply if the bit is set */
for (x = 0; x < bitcpy; x++) {
if ((err = mp_sqr(&res, &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, &mu)) != MP_OKAY) {
goto LBL_RES;
}
bitbuf <<= 1;
if ((bitbuf & (1 << winsize)) != 0) {
/* then multiply */
if ((err = mp_mul(&res, &M[1], &res)) != MP_OKAY) {
goto LBL_RES;
}
if ((err = redux(&res, P, &mu)) != MP_OKAY) {
goto LBL_RES;
}
}
}
}
mp_exch(&res, Y);
err = MP_OKAY;
LBL_RES:
mp_clear(&res);
mp_clear(&res);
LBL_MU:
mp_clear(&mu);
mp_clear(&mu);
LBL_M:
mp_clear(&M[1]);
for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
mp_clear(&M[x]);
}
return err;
mp_clear(&M[1]);
for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
mp_clear(&M[x]);
}
return err;
}
#endif