1189251Ssam/* 2189251Ssam * Minimal code for RSA support from LibTomMath 0.41 3189251Ssam * http://libtom.org/ 4189251Ssam * http://libtom.org/files/ltm-0.41.tar.bz2 5189251Ssam * This library was released in public domain by Tom St Denis. 6189251Ssam * 7189251Ssam * The combination in this file may not use all of the optimized algorithms 8189251Ssam * from LibTomMath and may be considerable slower than the LibTomMath with its 9189251Ssam * default settings. The main purpose of having this version here is to make it 10189251Ssam * easier to build bignum.c wrapper without having to install and build an 11189251Ssam * external library. 12189251Ssam * 13189251Ssam * If CONFIG_INTERNAL_LIBTOMMATH is defined, bignum.c includes this 14189251Ssam * libtommath.c file instead of using the external LibTomMath library. 15189251Ssam */ 16189251Ssam 17189251Ssam#ifndef CHAR_BIT 18189251Ssam#define CHAR_BIT 8 19189251Ssam#endif 20189251Ssam 21189251Ssam#define BN_MP_INVMOD_C 22189251Ssam#define BN_S_MP_EXPTMOD_C /* Note: #undef in tommath_superclass.h; this would 23189251Ssam * require BN_MP_EXPTMOD_FAST_C instead */ 24189251Ssam#define BN_S_MP_MUL_DIGS_C 25189251Ssam#define BN_MP_INVMOD_SLOW_C 26189251Ssam#define BN_S_MP_SQR_C 27189251Ssam#define BN_S_MP_MUL_HIGH_DIGS_C /* Note: #undef in tommath_superclass.h; this 28189251Ssam * would require other than mp_reduce */ 29189251Ssam 30189251Ssam#ifdef LTM_FAST 31189251Ssam 32189251Ssam/* Use faster div at the cost of about 1 kB */ 33189251Ssam#define BN_MP_MUL_D_C 34189251Ssam 35189251Ssam/* Include faster exptmod (Montgomery) at the cost of about 2.5 kB in code */ 36189251Ssam#define BN_MP_EXPTMOD_FAST_C 37189251Ssam#define BN_MP_MONTGOMERY_SETUP_C 38189251Ssam#define BN_FAST_MP_MONTGOMERY_REDUCE_C 39189251Ssam#define BN_MP_MONTGOMERY_CALC_NORMALIZATION_C 40189251Ssam#define BN_MP_MUL_2_C 41189251Ssam 42189251Ssam/* Include faster sqr at the cost of about 0.5 kB in code */ 43189251Ssam#define BN_FAST_S_MP_SQR_C 44189251Ssam 45189251Ssam#else /* LTM_FAST */ 46189251Ssam 47189251Ssam#define BN_MP_DIV_SMALL 48189251Ssam#define BN_MP_INIT_MULTI_C 49189251Ssam#define BN_MP_CLEAR_MULTI_C 50189251Ssam#define BN_MP_ABS_C 51189251Ssam#endif /* LTM_FAST */ 52189251Ssam 53189251Ssam/* Current uses do not require support for negative exponent in exptmod, so we 54189251Ssam * can save about 1.5 kB in leaving out invmod. */ 55189251Ssam#define LTM_NO_NEG_EXP 56189251Ssam 57189251Ssam/* from tommath.h */ 58189251Ssam 59189251Ssam#ifndef MIN 60189251Ssam #define MIN(x,y) ((x)<(y)?(x):(y)) 61189251Ssam#endif 62189251Ssam 63189251Ssam#ifndef MAX 64189251Ssam #define MAX(x,y) ((x)>(y)?(x):(y)) 65189251Ssam#endif 66189251Ssam 67189251Ssam#define OPT_CAST(x) 68189251Ssam 69189251Ssamtypedef unsigned long mp_digit; 70189251Ssamtypedef u64 mp_word; 71189251Ssam 72189251Ssam#define DIGIT_BIT 28 73189251Ssam#define MP_28BIT 74189251Ssam 75189251Ssam 76189251Ssam#define XMALLOC os_malloc 77189251Ssam#define XFREE os_free 78189251Ssam#define XREALLOC os_realloc 79189251Ssam 80189251Ssam 81189251Ssam#define MP_MASK ((((mp_digit)1)<<((mp_digit)DIGIT_BIT))-((mp_digit)1)) 82189251Ssam 83189251Ssam#define MP_LT -1 /* less than */ 84189251Ssam#define MP_EQ 0 /* equal to */ 85189251Ssam#define MP_GT 1 /* greater than */ 86189251Ssam 87189251Ssam#define MP_ZPOS 0 /* positive integer */ 88189251Ssam#define MP_NEG 1 /* negative */ 89189251Ssam 90189251Ssam#define MP_OKAY 0 /* ok result */ 91189251Ssam#define MP_MEM -2 /* out of mem */ 92189251Ssam#define MP_VAL -3 /* invalid input */ 93189251Ssam 94189251Ssam#define MP_YES 1 /* yes response */ 95189251Ssam#define MP_NO 0 /* no response */ 96189251Ssam 97189251Ssamtypedef int mp_err; 98189251Ssam 99189251Ssam/* define this to use lower memory usage routines (exptmods mostly) */ 100189251Ssam#define MP_LOW_MEM 101189251Ssam 102189251Ssam/* default precision */ 103189251Ssam#ifndef MP_PREC 104189251Ssam #ifndef MP_LOW_MEM 105189251Ssam #define MP_PREC 32 /* default digits of precision */ 106189251Ssam #else 107189251Ssam #define MP_PREC 8 /* default digits of precision */ 108189251Ssam #endif 109189251Ssam#endif 110189251Ssam 111189251Ssam/* size of comba arrays, should be at least 2 * 2**(BITS_PER_WORD - BITS_PER_DIGIT*2) */ 112189251Ssam#define MP_WARRAY (1 << (sizeof(mp_word) * CHAR_BIT - 2 * DIGIT_BIT + 1)) 113189251Ssam 114189251Ssam/* the infamous mp_int structure */ 115189251Ssamtypedef struct { 116189251Ssam int used, alloc, sign; 117189251Ssam mp_digit *dp; 118189251Ssam} mp_int; 119189251Ssam 120189251Ssam 121189251Ssam/* ---> Basic Manipulations <--- */ 122189251Ssam#define mp_iszero(a) (((a)->used == 0) ? MP_YES : MP_NO) 123189251Ssam#define mp_iseven(a) (((a)->used > 0 && (((a)->dp[0] & 1) == 0)) ? MP_YES : MP_NO) 124189251Ssam#define mp_isodd(a) (((a)->used > 0 && (((a)->dp[0] & 1) == 1)) ? MP_YES : MP_NO) 125189251Ssam 126189251Ssam 127189251Ssam/* prototypes for copied functions */ 128189251Ssam#define s_mp_mul(a, b, c) s_mp_mul_digs(a, b, c, (a)->used + (b)->used + 1) 129189251Ssamstatic int s_mp_exptmod(mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode); 130189251Ssamstatic int s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs); 131189251Ssamstatic int s_mp_sqr(mp_int * a, mp_int * b); 132189251Ssamstatic int s_mp_mul_high_digs(mp_int * a, mp_int * b, mp_int * c, int digs); 133189251Ssam 134189251Ssamstatic int fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs); 135189251Ssam 136189251Ssam#ifdef BN_MP_INIT_MULTI_C 137189251Ssamstatic int mp_init_multi(mp_int *mp, ...); 138189251Ssam#endif 139189251Ssam#ifdef BN_MP_CLEAR_MULTI_C 140189251Ssamstatic void mp_clear_multi(mp_int *mp, ...); 141189251Ssam#endif 142189251Ssamstatic int mp_lshd(mp_int * a, int b); 143189251Ssamstatic void mp_set(mp_int * a, mp_digit b); 144189251Ssamstatic void mp_clamp(mp_int * a); 145189251Ssamstatic void mp_exch(mp_int * a, mp_int * b); 146189251Ssamstatic void mp_rshd(mp_int * a, int b); 147189251Ssamstatic void mp_zero(mp_int * a); 148189251Ssamstatic int mp_mod_2d(mp_int * a, int b, mp_int * c); 149189251Ssamstatic int mp_div_2d(mp_int * a, int b, mp_int * c, mp_int * d); 150189251Ssamstatic int mp_init_copy(mp_int * a, mp_int * b); 151189251Ssamstatic int mp_mul_2d(mp_int * a, int b, mp_int * c); 152189251Ssam#ifndef LTM_NO_NEG_EXP 153189251Ssamstatic int mp_div_2(mp_int * a, mp_int * b); 154189251Ssamstatic int mp_invmod(mp_int * a, mp_int * b, mp_int * c); 155189251Ssamstatic int mp_invmod_slow(mp_int * a, mp_int * b, mp_int * c); 156189251Ssam#endif /* LTM_NO_NEG_EXP */ 157189251Ssamstatic int mp_copy(mp_int * a, mp_int * b); 158189251Ssamstatic int mp_count_bits(mp_int * a); 159189251Ssamstatic int mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d); 160189251Ssamstatic int mp_mod(mp_int * a, mp_int * b, mp_int * c); 161189251Ssamstatic int mp_grow(mp_int * a, int size); 162189251Ssamstatic int mp_cmp_mag(mp_int * a, mp_int * b); 163189251Ssam#ifdef BN_MP_ABS_C 164189251Ssamstatic int mp_abs(mp_int * a, mp_int * b); 165189251Ssam#endif 166189251Ssamstatic int mp_sqr(mp_int * a, mp_int * b); 167189251Ssamstatic int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d); 168189251Ssamstatic int mp_reduce_2k_setup_l(mp_int *a, mp_int *d); 169189251Ssamstatic int mp_2expt(mp_int * a, int b); 170189251Ssamstatic int mp_reduce_setup(mp_int * a, mp_int * b); 171189251Ssamstatic int mp_reduce(mp_int * x, mp_int * m, mp_int * mu); 172189251Ssamstatic int mp_init_size(mp_int * a, int size); 173189251Ssam#ifdef BN_MP_EXPTMOD_FAST_C 174189251Ssamstatic int mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode); 175189251Ssam#endif /* BN_MP_EXPTMOD_FAST_C */ 176189251Ssam#ifdef BN_FAST_S_MP_SQR_C 177189251Ssamstatic int fast_s_mp_sqr (mp_int * a, mp_int * b); 178189251Ssam#endif /* BN_FAST_S_MP_SQR_C */ 179189251Ssam#ifdef BN_MP_MUL_D_C 180189251Ssamstatic int mp_mul_d (mp_int * a, mp_digit b, mp_int * c); 181189251Ssam#endif /* BN_MP_MUL_D_C */ 182189251Ssam 183189251Ssam 184189251Ssam 185189251Ssam/* functions from bn_<func name>.c */ 186189251Ssam 187189251Ssam 188189251Ssam/* reverse an array, used for radix code */ 189189251Ssamstatic void bn_reverse (unsigned char *s, int len) 190189251Ssam{ 191189251Ssam int ix, iy; 192189251Ssam unsigned char t; 193189251Ssam 194189251Ssam ix = 0; 195189251Ssam iy = len - 1; 196189251Ssam while (ix < iy) { 197189251Ssam t = s[ix]; 198189251Ssam s[ix] = s[iy]; 199189251Ssam s[iy] = t; 200189251Ssam ++ix; 201189251Ssam --iy; 202189251Ssam } 203189251Ssam} 204189251Ssam 205189251Ssam 206189251Ssam/* low level addition, based on HAC pp.594, Algorithm 14.7 */ 207189251Ssamstatic int s_mp_add (mp_int * a, mp_int * b, mp_int * c) 208189251Ssam{ 209189251Ssam mp_int *x; 210189251Ssam int olduse, res, min, max; 211189251Ssam 212189251Ssam /* find sizes, we let |a| <= |b| which means we have to sort 213189251Ssam * them. "x" will point to the input with the most digits 214189251Ssam */ 215189251Ssam if (a->used > b->used) { 216189251Ssam min = b->used; 217189251Ssam max = a->used; 218189251Ssam x = a; 219189251Ssam } else { 220189251Ssam min = a->used; 221189251Ssam max = b->used; 222189251Ssam x = b; 223189251Ssam } 224189251Ssam 225189251Ssam /* init result */ 226189251Ssam if (c->alloc < max + 1) { 227189251Ssam if ((res = mp_grow (c, max + 1)) != MP_OKAY) { 228189251Ssam return res; 229189251Ssam } 230189251Ssam } 231189251Ssam 232189251Ssam /* get old used digit count and set new one */ 233189251Ssam olduse = c->used; 234189251Ssam c->used = max + 1; 235189251Ssam 236189251Ssam { 237189251Ssam register mp_digit u, *tmpa, *tmpb, *tmpc; 238189251Ssam register int i; 239189251Ssam 240189251Ssam /* alias for digit pointers */ 241189251Ssam 242189251Ssam /* first input */ 243189251Ssam tmpa = a->dp; 244189251Ssam 245189251Ssam /* second input */ 246189251Ssam tmpb = b->dp; 247189251Ssam 248189251Ssam /* destination */ 249189251Ssam tmpc = c->dp; 250189251Ssam 251189251Ssam /* zero the carry */ 252189251Ssam u = 0; 253189251Ssam for (i = 0; i < min; i++) { 254189251Ssam /* Compute the sum at one digit, T[i] = A[i] + B[i] + U */ 255189251Ssam *tmpc = *tmpa++ + *tmpb++ + u; 256189251Ssam 257189251Ssam /* U = carry bit of T[i] */ 258189251Ssam u = *tmpc >> ((mp_digit)DIGIT_BIT); 259189251Ssam 260189251Ssam /* take away carry bit from T[i] */ 261189251Ssam *tmpc++ &= MP_MASK; 262189251Ssam } 263189251Ssam 264189251Ssam /* now copy higher words if any, that is in A+B 265189251Ssam * if A or B has more digits add those in 266189251Ssam */ 267189251Ssam if (min != max) { 268189251Ssam for (; i < max; i++) { 269189251Ssam /* T[i] = X[i] + U */ 270189251Ssam *tmpc = x->dp[i] + u; 271189251Ssam 272189251Ssam /* U = carry bit of T[i] */ 273189251Ssam u = *tmpc >> ((mp_digit)DIGIT_BIT); 274189251Ssam 275189251Ssam /* take away carry bit from T[i] */ 276189251Ssam *tmpc++ &= MP_MASK; 277189251Ssam } 278189251Ssam } 279189251Ssam 280189251Ssam /* add carry */ 281189251Ssam *tmpc++ = u; 282189251Ssam 283189251Ssam /* clear digits above oldused */ 284189251Ssam for (i = c->used; i < olduse; i++) { 285189251Ssam *tmpc++ = 0; 286189251Ssam } 287189251Ssam } 288189251Ssam 289189251Ssam mp_clamp (c); 290189251Ssam return MP_OKAY; 291189251Ssam} 292189251Ssam 293189251Ssam 294189251Ssam/* low level subtraction (assumes |a| > |b|), HAC pp.595 Algorithm 14.9 */ 295189251Ssamstatic int s_mp_sub (mp_int * a, mp_int * b, mp_int * c) 296189251Ssam{ 297189251Ssam int olduse, res, min, max; 298189251Ssam 299189251Ssam /* find sizes */ 300189251Ssam min = b->used; 301189251Ssam max = a->used; 302189251Ssam 303189251Ssam /* init result */ 304189251Ssam if (c->alloc < max) { 305189251Ssam if ((res = mp_grow (c, max)) != MP_OKAY) { 306189251Ssam return res; 307189251Ssam } 308189251Ssam } 309189251Ssam olduse = c->used; 310189251Ssam c->used = max; 311189251Ssam 312189251Ssam { 313189251Ssam register mp_digit u, *tmpa, *tmpb, *tmpc; 314189251Ssam register int i; 315189251Ssam 316189251Ssam /* alias for digit pointers */ 317189251Ssam tmpa = a->dp; 318189251Ssam tmpb = b->dp; 319189251Ssam tmpc = c->dp; 320189251Ssam 321189251Ssam /* set carry to zero */ 322189251Ssam u = 0; 323189251Ssam for (i = 0; i < min; i++) { 324189251Ssam /* T[i] = A[i] - B[i] - U */ 325189251Ssam *tmpc = *tmpa++ - *tmpb++ - u; 326189251Ssam 327189251Ssam /* U = carry bit of T[i] 328189251Ssam * Note this saves performing an AND operation since 329189251Ssam * if a carry does occur it will propagate all the way to the 330189251Ssam * MSB. As a result a single shift is enough to get the carry 331189251Ssam */ 332189251Ssam u = *tmpc >> ((mp_digit)(CHAR_BIT * sizeof (mp_digit) - 1)); 333189251Ssam 334189251Ssam /* Clear carry from T[i] */ 335189251Ssam *tmpc++ &= MP_MASK; 336189251Ssam } 337189251Ssam 338189251Ssam /* now copy higher words if any, e.g. if A has more digits than B */ 339189251Ssam for (; i < max; i++) { 340189251Ssam /* T[i] = A[i] - U */ 341189251Ssam *tmpc = *tmpa++ - u; 342189251Ssam 343189251Ssam /* U = carry bit of T[i] */ 344189251Ssam u = *tmpc >> ((mp_digit)(CHAR_BIT * sizeof (mp_digit) - 1)); 345189251Ssam 346189251Ssam /* Clear carry from T[i] */ 347189251Ssam *tmpc++ &= MP_MASK; 348189251Ssam } 349189251Ssam 350189251Ssam /* clear digits above used (since we may not have grown result above) */ 351189251Ssam for (i = c->used; i < olduse; i++) { 352189251Ssam *tmpc++ = 0; 353189251Ssam } 354189251Ssam } 355189251Ssam 356189251Ssam mp_clamp (c); 357189251Ssam return MP_OKAY; 358189251Ssam} 359189251Ssam 360189251Ssam 361189251Ssam/* init a new mp_int */ 362189251Ssamstatic int mp_init (mp_int * a) 363189251Ssam{ 364189251Ssam int i; 365189251Ssam 366189251Ssam /* allocate memory required and clear it */ 367189251Ssam a->dp = OPT_CAST(mp_digit) XMALLOC (sizeof (mp_digit) * MP_PREC); 368189251Ssam if (a->dp == NULL) { 369189251Ssam return MP_MEM; 370189251Ssam } 371189251Ssam 372189251Ssam /* set the digits to zero */ 373189251Ssam for (i = 0; i < MP_PREC; i++) { 374189251Ssam a->dp[i] = 0; 375189251Ssam } 376189251Ssam 377189251Ssam /* set the used to zero, allocated digits to the default precision 378189251Ssam * and sign to positive */ 379189251Ssam a->used = 0; 380189251Ssam a->alloc = MP_PREC; 381189251Ssam a->sign = MP_ZPOS; 382189251Ssam 383189251Ssam return MP_OKAY; 384189251Ssam} 385189251Ssam 386189251Ssam 387189251Ssam/* clear one (frees) */ 388189251Ssamstatic void mp_clear (mp_int * a) 389189251Ssam{ 390189251Ssam int i; 391189251Ssam 392189251Ssam /* only do anything if a hasn't been freed previously */ 393189251Ssam if (a->dp != NULL) { 394189251Ssam /* first zero the digits */ 395189251Ssam for (i = 0; i < a->used; i++) { 396189251Ssam a->dp[i] = 0; 397189251Ssam } 398189251Ssam 399189251Ssam /* free ram */ 400189251Ssam XFREE(a->dp); 401189251Ssam 402189251Ssam /* reset members to make debugging easier */ 403189251Ssam a->dp = NULL; 404189251Ssam a->alloc = a->used = 0; 405189251Ssam a->sign = MP_ZPOS; 406189251Ssam } 407189251Ssam} 408189251Ssam 409189251Ssam 410189251Ssam/* high level addition (handles signs) */ 411189251Ssamstatic int mp_add (mp_int * a, mp_int * b, mp_int * c) 412189251Ssam{ 413189251Ssam int sa, sb, res; 414189251Ssam 415189251Ssam /* get sign of both inputs */ 416189251Ssam sa = a->sign; 417189251Ssam sb = b->sign; 418189251Ssam 419189251Ssam /* handle two cases, not four */ 420189251Ssam if (sa == sb) { 421189251Ssam /* both positive or both negative */ 422189251Ssam /* add their magnitudes, copy the sign */ 423189251Ssam c->sign = sa; 424189251Ssam res = s_mp_add (a, b, c); 425189251Ssam } else { 426189251Ssam /* one positive, the other negative */ 427189251Ssam /* subtract the one with the greater magnitude from */ 428189251Ssam /* the one of the lesser magnitude. The result gets */ 429189251Ssam /* the sign of the one with the greater magnitude. */ 430189251Ssam if (mp_cmp_mag (a, b) == MP_LT) { 431189251Ssam c->sign = sb; 432189251Ssam res = s_mp_sub (b, a, c); 433189251Ssam } else { 434189251Ssam c->sign = sa; 435189251Ssam res = s_mp_sub (a, b, c); 436189251Ssam } 437189251Ssam } 438189251Ssam return res; 439189251Ssam} 440189251Ssam 441189251Ssam 442189251Ssam/* high level subtraction (handles signs) */ 443189251Ssamstatic int mp_sub (mp_int * a, mp_int * b, mp_int * c) 444189251Ssam{ 445189251Ssam int sa, sb, res; 446189251Ssam 447189251Ssam sa = a->sign; 448189251Ssam sb = b->sign; 449189251Ssam 450189251Ssam if (sa != sb) { 451189251Ssam /* subtract a negative from a positive, OR */ 452189251Ssam /* subtract a positive from a negative. */ 453189251Ssam /* In either case, ADD their magnitudes, */ 454189251Ssam /* and use the sign of the first number. */ 455189251Ssam c->sign = sa; 456189251Ssam res = s_mp_add (a, b, c); 457189251Ssam } else { 458189251Ssam /* subtract a positive from a positive, OR */ 459189251Ssam /* subtract a negative from a negative. */ 460189251Ssam /* First, take the difference between their */ 461189251Ssam /* magnitudes, then... */ 462189251Ssam if (mp_cmp_mag (a, b) != MP_LT) { 463189251Ssam /* Copy the sign from the first */ 464189251Ssam c->sign = sa; 465189251Ssam /* The first has a larger or equal magnitude */ 466189251Ssam res = s_mp_sub (a, b, c); 467189251Ssam } else { 468189251Ssam /* The result has the *opposite* sign from */ 469189251Ssam /* the first number. */ 470189251Ssam c->sign = (sa == MP_ZPOS) ? MP_NEG : MP_ZPOS; 471189251Ssam /* The second has a larger magnitude */ 472189251Ssam res = s_mp_sub (b, a, c); 473189251Ssam } 474189251Ssam } 475189251Ssam return res; 476189251Ssam} 477189251Ssam 478189251Ssam 479189251Ssam/* high level multiplication (handles sign) */ 480189251Ssamstatic int mp_mul (mp_int * a, mp_int * b, mp_int * c) 481189251Ssam{ 482189251Ssam int res, neg; 483189251Ssam neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG; 484189251Ssam 485189251Ssam /* use Toom-Cook? */ 486189251Ssam#ifdef BN_MP_TOOM_MUL_C 487189251Ssam if (MIN (a->used, b->used) >= TOOM_MUL_CUTOFF) { 488189251Ssam res = mp_toom_mul(a, b, c); 489189251Ssam } else 490189251Ssam#endif 491189251Ssam#ifdef BN_MP_KARATSUBA_MUL_C 492189251Ssam /* use Karatsuba? */ 493189251Ssam if (MIN (a->used, b->used) >= KARATSUBA_MUL_CUTOFF) { 494189251Ssam res = mp_karatsuba_mul (a, b, c); 495189251Ssam } else 496189251Ssam#endif 497189251Ssam { 498189251Ssam /* can we use the fast multiplier? 499189251Ssam * 500189251Ssam * The fast multiplier can be used if the output will 501189251Ssam * have less than MP_WARRAY digits and the number of 502189251Ssam * digits won't affect carry propagation 503189251Ssam */ 504189251Ssam#ifdef BN_FAST_S_MP_MUL_DIGS_C 505189251Ssam int digs = a->used + b->used + 1; 506189251Ssam 507189251Ssam if ((digs < MP_WARRAY) && 508189251Ssam MIN(a->used, b->used) <= 509189251Ssam (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) { 510189251Ssam res = fast_s_mp_mul_digs (a, b, c, digs); 511189251Ssam } else 512189251Ssam#endif 513189251Ssam#ifdef BN_S_MP_MUL_DIGS_C 514189251Ssam res = s_mp_mul (a, b, c); /* uses s_mp_mul_digs */ 515189251Ssam#else 516189251Ssam#error mp_mul could fail 517189251Ssam res = MP_VAL; 518189251Ssam#endif 519189251Ssam 520189251Ssam } 521189251Ssam c->sign = (c->used > 0) ? neg : MP_ZPOS; 522189251Ssam return res; 523189251Ssam} 524189251Ssam 525189251Ssam 526189251Ssam/* d = a * b (mod c) */ 527189251Ssamstatic int mp_mulmod (mp_int * a, mp_int * b, mp_int * c, mp_int * d) 528189251Ssam{ 529189251Ssam int res; 530189251Ssam mp_int t; 531189251Ssam 532189251Ssam if ((res = mp_init (&t)) != MP_OKAY) { 533189251Ssam return res; 534189251Ssam } 535189251Ssam 536189251Ssam if ((res = mp_mul (a, b, &t)) != MP_OKAY) { 537189251Ssam mp_clear (&t); 538189251Ssam return res; 539189251Ssam } 540189251Ssam res = mp_mod (&t, c, d); 541189251Ssam mp_clear (&t); 542189251Ssam return res; 543189251Ssam} 544189251Ssam 545189251Ssam 546189251Ssam/* c = a mod b, 0 <= c < b */ 547189251Ssamstatic int mp_mod (mp_int * a, mp_int * b, mp_int * c) 548189251Ssam{ 549189251Ssam mp_int t; 550189251Ssam int res; 551189251Ssam 552189251Ssam if ((res = mp_init (&t)) != MP_OKAY) { 553189251Ssam return res; 554189251Ssam } 555189251Ssam 556189251Ssam if ((res = mp_div (a, b, NULL, &t)) != MP_OKAY) { 557189251Ssam mp_clear (&t); 558189251Ssam return res; 559189251Ssam } 560189251Ssam 561189251Ssam if (t.sign != b->sign) { 562189251Ssam res = mp_add (b, &t, c); 563189251Ssam } else { 564189251Ssam res = MP_OKAY; 565189251Ssam mp_exch (&t, c); 566189251Ssam } 567189251Ssam 568189251Ssam mp_clear (&t); 569189251Ssam return res; 570189251Ssam} 571189251Ssam 572189251Ssam 573189251Ssam/* this is a shell function that calls either the normal or Montgomery 574189251Ssam * exptmod functions. Originally the call to the montgomery code was 575189251Ssam * embedded in the normal function but that wasted alot of stack space 576189251Ssam * for nothing (since 99% of the time the Montgomery code would be called) 577189251Ssam */ 578189251Ssamstatic int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y) 579189251Ssam{ 580189251Ssam int dr; 581189251Ssam 582189251Ssam /* modulus P must be positive */ 583189251Ssam if (P->sign == MP_NEG) { 584189251Ssam return MP_VAL; 585189251Ssam } 586189251Ssam 587189251Ssam /* if exponent X is negative we have to recurse */ 588189251Ssam if (X->sign == MP_NEG) { 589189251Ssam#ifdef LTM_NO_NEG_EXP 590189251Ssam return MP_VAL; 591189251Ssam#else /* LTM_NO_NEG_EXP */ 592189251Ssam#ifdef BN_MP_INVMOD_C 593189251Ssam mp_int tmpG, tmpX; 594189251Ssam int err; 595189251Ssam 596189251Ssam /* first compute 1/G mod P */ 597189251Ssam if ((err = mp_init(&tmpG)) != MP_OKAY) { 598189251Ssam return err; 599189251Ssam } 600189251Ssam if ((err = mp_invmod(G, P, &tmpG)) != MP_OKAY) { 601189251Ssam mp_clear(&tmpG); 602189251Ssam return err; 603189251Ssam } 604189251Ssam 605189251Ssam /* now get |X| */ 606189251Ssam if ((err = mp_init(&tmpX)) != MP_OKAY) { 607189251Ssam mp_clear(&tmpG); 608189251Ssam return err; 609189251Ssam } 610189251Ssam if ((err = mp_abs(X, &tmpX)) != MP_OKAY) { 611189251Ssam mp_clear_multi(&tmpG, &tmpX, NULL); 612189251Ssam return err; 613189251Ssam } 614189251Ssam 615189251Ssam /* and now compute (1/G)**|X| instead of G**X [X < 0] */ 616189251Ssam err = mp_exptmod(&tmpG, &tmpX, P, Y); 617189251Ssam mp_clear_multi(&tmpG, &tmpX, NULL); 618189251Ssam return err; 619189251Ssam#else 620189251Ssam#error mp_exptmod would always fail 621189251Ssam /* no invmod */ 622189251Ssam return MP_VAL; 623189251Ssam#endif 624189251Ssam#endif /* LTM_NO_NEG_EXP */ 625189251Ssam } 626189251Ssam 627189251Ssam/* modified diminished radix reduction */ 628189251Ssam#if defined(BN_MP_REDUCE_IS_2K_L_C) && defined(BN_MP_REDUCE_2K_L_C) && defined(BN_S_MP_EXPTMOD_C) 629189251Ssam if (mp_reduce_is_2k_l(P) == MP_YES) { 630189251Ssam return s_mp_exptmod(G, X, P, Y, 1); 631189251Ssam } 632189251Ssam#endif 633189251Ssam 634189251Ssam#ifdef BN_MP_DR_IS_MODULUS_C 635189251Ssam /* is it a DR modulus? */ 636189251Ssam dr = mp_dr_is_modulus(P); 637189251Ssam#else 638189251Ssam /* default to no */ 639189251Ssam dr = 0; 640189251Ssam#endif 641189251Ssam 642189251Ssam#ifdef BN_MP_REDUCE_IS_2K_C 643189251Ssam /* if not, is it a unrestricted DR modulus? */ 644189251Ssam if (dr == 0) { 645189251Ssam dr = mp_reduce_is_2k(P) << 1; 646189251Ssam } 647189251Ssam#endif 648189251Ssam 649189251Ssam /* if the modulus is odd or dr != 0 use the montgomery method */ 650189251Ssam#ifdef BN_MP_EXPTMOD_FAST_C 651189251Ssam if (mp_isodd (P) == 1 || dr != 0) { 652189251Ssam return mp_exptmod_fast (G, X, P, Y, dr); 653189251Ssam } else { 654189251Ssam#endif 655189251Ssam#ifdef BN_S_MP_EXPTMOD_C 656189251Ssam /* otherwise use the generic Barrett reduction technique */ 657189251Ssam return s_mp_exptmod (G, X, P, Y, 0); 658189251Ssam#else 659189251Ssam#error mp_exptmod could fail 660189251Ssam /* no exptmod for evens */ 661189251Ssam return MP_VAL; 662189251Ssam#endif 663189251Ssam#ifdef BN_MP_EXPTMOD_FAST_C 664189251Ssam } 665189251Ssam#endif 666189251Ssam} 667189251Ssam 668189251Ssam 669189251Ssam/* compare two ints (signed)*/ 670189251Ssamstatic int mp_cmp (mp_int * a, mp_int * b) 671189251Ssam{ 672189251Ssam /* compare based on sign */ 673189251Ssam if (a->sign != b->sign) { 674189251Ssam if (a->sign == MP_NEG) { 675189251Ssam return MP_LT; 676189251Ssam } else { 677189251Ssam return MP_GT; 678189251Ssam } 679189251Ssam } 680189251Ssam 681189251Ssam /* compare digits */ 682189251Ssam if (a->sign == MP_NEG) { 683189251Ssam /* if negative compare opposite direction */ 684189251Ssam return mp_cmp_mag(b, a); 685189251Ssam } else { 686189251Ssam return mp_cmp_mag(a, b); 687189251Ssam } 688189251Ssam} 689189251Ssam 690189251Ssam 691189251Ssam/* compare a digit */ 692189251Ssamstatic int mp_cmp_d(mp_int * a, mp_digit b) 693189251Ssam{ 694189251Ssam /* compare based on sign */ 695189251Ssam if (a->sign == MP_NEG) { 696189251Ssam return MP_LT; 697189251Ssam } 698189251Ssam 699189251Ssam /* compare based on magnitude */ 700189251Ssam if (a->used > 1) { 701189251Ssam return MP_GT; 702189251Ssam } 703189251Ssam 704189251Ssam /* compare the only digit of a to b */ 705189251Ssam if (a->dp[0] > b) { 706189251Ssam return MP_GT; 707189251Ssam } else if (a->dp[0] < b) { 708189251Ssam return MP_LT; 709189251Ssam } else { 710189251Ssam return MP_EQ; 711189251Ssam } 712189251Ssam} 713189251Ssam 714189251Ssam 715189251Ssam#ifndef LTM_NO_NEG_EXP 716189251Ssam/* hac 14.61, pp608 */ 717189251Ssamstatic int mp_invmod (mp_int * a, mp_int * b, mp_int * c) 718189251Ssam{ 719189251Ssam /* b cannot be negative */ 720189251Ssam if (b->sign == MP_NEG || mp_iszero(b) == 1) { 721189251Ssam return MP_VAL; 722189251Ssam } 723189251Ssam 724189251Ssam#ifdef BN_FAST_MP_INVMOD_C 725189251Ssam /* if the modulus is odd we can use a faster routine instead */ 726189251Ssam if (mp_isodd (b) == 1) { 727189251Ssam return fast_mp_invmod (a, b, c); 728189251Ssam } 729189251Ssam#endif 730189251Ssam 731189251Ssam#ifdef BN_MP_INVMOD_SLOW_C 732189251Ssam return mp_invmod_slow(a, b, c); 733189251Ssam#endif 734189251Ssam 735189251Ssam#ifndef BN_FAST_MP_INVMOD_C 736189251Ssam#ifndef BN_MP_INVMOD_SLOW_C 737189251Ssam#error mp_invmod would always fail 738189251Ssam#endif 739189251Ssam#endif 740189251Ssam return MP_VAL; 741189251Ssam} 742189251Ssam#endif /* LTM_NO_NEG_EXP */ 743189251Ssam 744189251Ssam 745189251Ssam/* get the size for an unsigned equivalent */ 746189251Ssamstatic int mp_unsigned_bin_size (mp_int * a) 747189251Ssam{ 748189251Ssam int size = mp_count_bits (a); 749189251Ssam return (size / 8 + ((size & 7) != 0 ? 1 : 0)); 750189251Ssam} 751189251Ssam 752189251Ssam 753189251Ssam#ifndef LTM_NO_NEG_EXP 754189251Ssam/* hac 14.61, pp608 */ 755189251Ssamstatic int mp_invmod_slow (mp_int * a, mp_int * b, mp_int * c) 756189251Ssam{ 757189251Ssam mp_int x, y, u, v, A, B, C, D; 758189251Ssam int res; 759189251Ssam 760189251Ssam /* b cannot be negative */ 761189251Ssam if (b->sign == MP_NEG || mp_iszero(b) == 1) { 762189251Ssam return MP_VAL; 763189251Ssam } 764189251Ssam 765189251Ssam /* init temps */ 766189251Ssam if ((res = mp_init_multi(&x, &y, &u, &v, 767189251Ssam &A, &B, &C, &D, NULL)) != MP_OKAY) { 768189251Ssam return res; 769189251Ssam } 770189251Ssam 771189251Ssam /* x = a, y = b */ 772189251Ssam if ((res = mp_mod(a, b, &x)) != MP_OKAY) { 773189251Ssam goto LBL_ERR; 774189251Ssam } 775189251Ssam if ((res = mp_copy (b, &y)) != MP_OKAY) { 776189251Ssam goto LBL_ERR; 777189251Ssam } 778189251Ssam 779189251Ssam /* 2. [modified] if x,y are both even then return an error! */ 780189251Ssam if (mp_iseven (&x) == 1 && mp_iseven (&y) == 1) { 781189251Ssam res = MP_VAL; 782189251Ssam goto LBL_ERR; 783189251Ssam } 784189251Ssam 785189251Ssam /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */ 786189251Ssam if ((res = mp_copy (&x, &u)) != MP_OKAY) { 787189251Ssam goto LBL_ERR; 788189251Ssam } 789189251Ssam if ((res = mp_copy (&y, &v)) != MP_OKAY) { 790189251Ssam goto LBL_ERR; 791189251Ssam } 792189251Ssam mp_set (&A, 1); 793189251Ssam mp_set (&D, 1); 794189251Ssam 795189251Ssamtop: 796189251Ssam /* 4. while u is even do */ 797189251Ssam while (mp_iseven (&u) == 1) { 798189251Ssam /* 4.1 u = u/2 */ 799189251Ssam if ((res = mp_div_2 (&u, &u)) != MP_OKAY) { 800189251Ssam goto LBL_ERR; 801189251Ssam } 802189251Ssam /* 4.2 if A or B is odd then */ 803189251Ssam if (mp_isodd (&A) == 1 || mp_isodd (&B) == 1) { 804189251Ssam /* A = (A+y)/2, B = (B-x)/2 */ 805189251Ssam if ((res = mp_add (&A, &y, &A)) != MP_OKAY) { 806189251Ssam goto LBL_ERR; 807189251Ssam } 808189251Ssam if ((res = mp_sub (&B, &x, &B)) != MP_OKAY) { 809189251Ssam goto LBL_ERR; 810189251Ssam } 811189251Ssam } 812189251Ssam /* A = A/2, B = B/2 */ 813189251Ssam if ((res = mp_div_2 (&A, &A)) != MP_OKAY) { 814189251Ssam goto LBL_ERR; 815189251Ssam } 816189251Ssam if ((res = mp_div_2 (&B, &B)) != MP_OKAY) { 817189251Ssam goto LBL_ERR; 818189251Ssam } 819189251Ssam } 820189251Ssam 821189251Ssam /* 5. while v is even do */ 822189251Ssam while (mp_iseven (&v) == 1) { 823189251Ssam /* 5.1 v = v/2 */ 824189251Ssam if ((res = mp_div_2 (&v, &v)) != MP_OKAY) { 825189251Ssam goto LBL_ERR; 826189251Ssam } 827189251Ssam /* 5.2 if C or D is odd then */ 828189251Ssam if (mp_isodd (&C) == 1 || mp_isodd (&D) == 1) { 829189251Ssam /* C = (C+y)/2, D = (D-x)/2 */ 830189251Ssam if ((res = mp_add (&C, &y, &C)) != MP_OKAY) { 831189251Ssam goto LBL_ERR; 832189251Ssam } 833189251Ssam if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) { 834189251Ssam goto LBL_ERR; 835189251Ssam } 836189251Ssam } 837189251Ssam /* C = C/2, D = D/2 */ 838189251Ssam if ((res = mp_div_2 (&C, &C)) != MP_OKAY) { 839189251Ssam goto LBL_ERR; 840189251Ssam } 841189251Ssam if ((res = mp_div_2 (&D, &D)) != MP_OKAY) { 842189251Ssam goto LBL_ERR; 843189251Ssam } 844189251Ssam } 845189251Ssam 846189251Ssam /* 6. if u >= v then */ 847189251Ssam if (mp_cmp (&u, &v) != MP_LT) { 848189251Ssam /* u = u - v, A = A - C, B = B - D */ 849189251Ssam if ((res = mp_sub (&u, &v, &u)) != MP_OKAY) { 850189251Ssam goto LBL_ERR; 851189251Ssam } 852189251Ssam 853189251Ssam if ((res = mp_sub (&A, &C, &A)) != MP_OKAY) { 854189251Ssam goto LBL_ERR; 855189251Ssam } 856189251Ssam 857189251Ssam if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) { 858189251Ssam goto LBL_ERR; 859189251Ssam } 860189251Ssam } else { 861189251Ssam /* v - v - u, C = C - A, D = D - B */ 862189251Ssam if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) { 863189251Ssam goto LBL_ERR; 864189251Ssam } 865189251Ssam 866189251Ssam if ((res = mp_sub (&C, &A, &C)) != MP_OKAY) { 867189251Ssam goto LBL_ERR; 868189251Ssam } 869189251Ssam 870189251Ssam if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) { 871189251Ssam goto LBL_ERR; 872189251Ssam } 873189251Ssam } 874189251Ssam 875189251Ssam /* if not zero goto step 4 */ 876189251Ssam if (mp_iszero (&u) == 0) 877189251Ssam goto top; 878189251Ssam 879189251Ssam /* now a = C, b = D, gcd == g*v */ 880189251Ssam 881189251Ssam /* if v != 1 then there is no inverse */ 882189251Ssam if (mp_cmp_d (&v, 1) != MP_EQ) { 883189251Ssam res = MP_VAL; 884189251Ssam goto LBL_ERR; 885189251Ssam } 886189251Ssam 887189251Ssam /* if its too low */ 888189251Ssam while (mp_cmp_d(&C, 0) == MP_LT) { 889189251Ssam if ((res = mp_add(&C, b, &C)) != MP_OKAY) { 890189251Ssam goto LBL_ERR; 891189251Ssam } 892189251Ssam } 893189251Ssam 894189251Ssam /* too big */ 895189251Ssam while (mp_cmp_mag(&C, b) != MP_LT) { 896189251Ssam if ((res = mp_sub(&C, b, &C)) != MP_OKAY) { 897189251Ssam goto LBL_ERR; 898189251Ssam } 899189251Ssam } 900189251Ssam 901189251Ssam /* C is now the inverse */ 902189251Ssam mp_exch (&C, c); 903189251Ssam res = MP_OKAY; 904189251SsamLBL_ERR:mp_clear_multi (&x, &y, &u, &v, &A, &B, &C, &D, NULL); 905189251Ssam return res; 906189251Ssam} 907189251Ssam#endif /* LTM_NO_NEG_EXP */ 908189251Ssam 909189251Ssam 910189251Ssam/* compare maginitude of two ints (unsigned) */ 911189251Ssamstatic int mp_cmp_mag (mp_int * a, mp_int * b) 912189251Ssam{ 913189251Ssam int n; 914189251Ssam mp_digit *tmpa, *tmpb; 915189251Ssam 916189251Ssam /* compare based on # of non-zero digits */ 917189251Ssam if (a->used > b->used) { 918189251Ssam return MP_GT; 919189251Ssam } 920189251Ssam 921189251Ssam if (a->used < b->used) { 922189251Ssam return MP_LT; 923189251Ssam } 924189251Ssam 925189251Ssam /* alias for a */ 926189251Ssam tmpa = a->dp + (a->used - 1); 927189251Ssam 928189251Ssam /* alias for b */ 929189251Ssam tmpb = b->dp + (a->used - 1); 930189251Ssam 931189251Ssam /* compare based on digits */ 932189251Ssam for (n = 0; n < a->used; ++n, --tmpa, --tmpb) { 933189251Ssam if (*tmpa > *tmpb) { 934189251Ssam return MP_GT; 935189251Ssam } 936189251Ssam 937189251Ssam if (*tmpa < *tmpb) { 938189251Ssam return MP_LT; 939189251Ssam } 940189251Ssam } 941189251Ssam return MP_EQ; 942189251Ssam} 943189251Ssam 944189251Ssam 945189251Ssam/* reads a unsigned char array, assumes the msb is stored first [big endian] */ 946189251Ssamstatic int mp_read_unsigned_bin (mp_int * a, const unsigned char *b, int c) 947189251Ssam{ 948189251Ssam int res; 949189251Ssam 950189251Ssam /* make sure there are at least two digits */ 951189251Ssam if (a->alloc < 2) { 952189251Ssam if ((res = mp_grow(a, 2)) != MP_OKAY) { 953189251Ssam return res; 954189251Ssam } 955189251Ssam } 956189251Ssam 957189251Ssam /* zero the int */ 958189251Ssam mp_zero (a); 959189251Ssam 960189251Ssam /* read the bytes in */ 961189251Ssam while (c-- > 0) { 962189251Ssam if ((res = mp_mul_2d (a, 8, a)) != MP_OKAY) { 963189251Ssam return res; 964189251Ssam } 965189251Ssam 966189251Ssam#ifndef MP_8BIT 967189251Ssam a->dp[0] |= *b++; 968189251Ssam a->used += 1; 969189251Ssam#else 970189251Ssam a->dp[0] = (*b & MP_MASK); 971189251Ssam a->dp[1] |= ((*b++ >> 7U) & 1); 972189251Ssam a->used += 2; 973189251Ssam#endif 974189251Ssam } 975189251Ssam mp_clamp (a); 976189251Ssam return MP_OKAY; 977189251Ssam} 978189251Ssam 979189251Ssam 980189251Ssam/* store in unsigned [big endian] format */ 981189251Ssamstatic int mp_to_unsigned_bin (mp_int * a, unsigned char *b) 982189251Ssam{ 983189251Ssam int x, res; 984189251Ssam mp_int t; 985189251Ssam 986189251Ssam if ((res = mp_init_copy (&t, a)) != MP_OKAY) { 987189251Ssam return res; 988189251Ssam } 989189251Ssam 990189251Ssam x = 0; 991189251Ssam while (mp_iszero (&t) == 0) { 992189251Ssam#ifndef MP_8BIT 993189251Ssam b[x++] = (unsigned char) (t.dp[0] & 255); 994189251Ssam#else 995189251Ssam b[x++] = (unsigned char) (t.dp[0] | ((t.dp[1] & 0x01) << 7)); 996189251Ssam#endif 997189251Ssam if ((res = mp_div_2d (&t, 8, &t, NULL)) != MP_OKAY) { 998189251Ssam mp_clear (&t); 999189251Ssam return res; 1000189251Ssam } 1001189251Ssam } 1002189251Ssam bn_reverse (b, x); 1003189251Ssam mp_clear (&t); 1004189251Ssam return MP_OKAY; 1005189251Ssam} 1006189251Ssam 1007189251Ssam 1008189251Ssam/* shift right by a certain bit count (store quotient in c, optional remainder in d) */ 1009189251Ssamstatic int mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d) 1010189251Ssam{ 1011189251Ssam mp_digit D, r, rr; 1012189251Ssam int x, res; 1013189251Ssam mp_int t; 1014189251Ssam 1015189251Ssam 1016189251Ssam /* if the shift count is <= 0 then we do no work */ 1017189251Ssam if (b <= 0) { 1018189251Ssam res = mp_copy (a, c); 1019189251Ssam if (d != NULL) { 1020189251Ssam mp_zero (d); 1021189251Ssam } 1022189251Ssam return res; 1023189251Ssam } 1024189251Ssam 1025189251Ssam if ((res = mp_init (&t)) != MP_OKAY) { 1026189251Ssam return res; 1027189251Ssam } 1028189251Ssam 1029189251Ssam /* get the remainder */ 1030189251Ssam if (d != NULL) { 1031189251Ssam if ((res = mp_mod_2d (a, b, &t)) != MP_OKAY) { 1032189251Ssam mp_clear (&t); 1033189251Ssam return res; 1034189251Ssam } 1035189251Ssam } 1036189251Ssam 1037189251Ssam /* copy */ 1038189251Ssam if ((res = mp_copy (a, c)) != MP_OKAY) { 1039189251Ssam mp_clear (&t); 1040189251Ssam return res; 1041189251Ssam } 1042189251Ssam 1043189251Ssam /* shift by as many digits in the bit count */ 1044189251Ssam if (b >= (int)DIGIT_BIT) { 1045189251Ssam mp_rshd (c, b / DIGIT_BIT); 1046189251Ssam } 1047189251Ssam 1048189251Ssam /* shift any bit count < DIGIT_BIT */ 1049189251Ssam D = (mp_digit) (b % DIGIT_BIT); 1050189251Ssam if (D != 0) { 1051189251Ssam register mp_digit *tmpc, mask, shift; 1052189251Ssam 1053189251Ssam /* mask */ 1054189251Ssam mask = (((mp_digit)1) << D) - 1; 1055189251Ssam 1056189251Ssam /* shift for lsb */ 1057189251Ssam shift = DIGIT_BIT - D; 1058189251Ssam 1059189251Ssam /* alias */ 1060189251Ssam tmpc = c->dp + (c->used - 1); 1061189251Ssam 1062189251Ssam /* carry */ 1063189251Ssam r = 0; 1064189251Ssam for (x = c->used - 1; x >= 0; x--) { 1065189251Ssam /* get the lower bits of this word in a temp */ 1066189251Ssam rr = *tmpc & mask; 1067189251Ssam 1068189251Ssam /* shift the current word and mix in the carry bits from the previous word */ 1069189251Ssam *tmpc = (*tmpc >> D) | (r << shift); 1070189251Ssam --tmpc; 1071189251Ssam 1072189251Ssam /* set the carry to the carry bits of the current word found above */ 1073189251Ssam r = rr; 1074189251Ssam } 1075189251Ssam } 1076189251Ssam mp_clamp (c); 1077189251Ssam if (d != NULL) { 1078189251Ssam mp_exch (&t, d); 1079189251Ssam } 1080189251Ssam mp_clear (&t); 1081189251Ssam return MP_OKAY; 1082189251Ssam} 1083189251Ssam 1084189251Ssam 1085189251Ssamstatic int mp_init_copy (mp_int * a, mp_int * b) 1086189251Ssam{ 1087189251Ssam int res; 1088189251Ssam 1089189251Ssam if ((res = mp_init (a)) != MP_OKAY) { 1090189251Ssam return res; 1091189251Ssam } 1092189251Ssam return mp_copy (b, a); 1093189251Ssam} 1094189251Ssam 1095189251Ssam 1096189251Ssam/* set to zero */ 1097189251Ssamstatic void mp_zero (mp_int * a) 1098189251Ssam{ 1099189251Ssam int n; 1100189251Ssam mp_digit *tmp; 1101189251Ssam 1102189251Ssam a->sign = MP_ZPOS; 1103189251Ssam a->used = 0; 1104189251Ssam 1105189251Ssam tmp = a->dp; 1106189251Ssam for (n = 0; n < a->alloc; n++) { 1107189251Ssam *tmp++ = 0; 1108189251Ssam } 1109189251Ssam} 1110189251Ssam 1111189251Ssam 1112189251Ssam/* copy, b = a */ 1113189251Ssamstatic int mp_copy (mp_int * a, mp_int * b) 1114189251Ssam{ 1115189251Ssam int res, n; 1116189251Ssam 1117189251Ssam /* if dst == src do nothing */ 1118189251Ssam if (a == b) { 1119189251Ssam return MP_OKAY; 1120189251Ssam } 1121189251Ssam 1122189251Ssam /* grow dest */ 1123189251Ssam if (b->alloc < a->used) { 1124189251Ssam if ((res = mp_grow (b, a->used)) != MP_OKAY) { 1125189251Ssam return res; 1126189251Ssam } 1127189251Ssam } 1128189251Ssam 1129189251Ssam /* zero b and copy the parameters over */ 1130189251Ssam { 1131189251Ssam register mp_digit *tmpa, *tmpb; 1132189251Ssam 1133189251Ssam /* pointer aliases */ 1134189251Ssam 1135189251Ssam /* source */ 1136189251Ssam tmpa = a->dp; 1137189251Ssam 1138189251Ssam /* destination */ 1139189251Ssam tmpb = b->dp; 1140189251Ssam 1141189251Ssam /* copy all the digits */ 1142189251Ssam for (n = 0; n < a->used; n++) { 1143189251Ssam *tmpb++ = *tmpa++; 1144189251Ssam } 1145189251Ssam 1146189251Ssam /* clear high digits */ 1147189251Ssam for (; n < b->used; n++) { 1148189251Ssam *tmpb++ = 0; 1149189251Ssam } 1150189251Ssam } 1151189251Ssam 1152189251Ssam /* copy used count and sign */ 1153189251Ssam b->used = a->used; 1154189251Ssam b->sign = a->sign; 1155189251Ssam return MP_OKAY; 1156189251Ssam} 1157189251Ssam 1158189251Ssam 1159189251Ssam/* shift right a certain amount of digits */ 1160189251Ssamstatic void mp_rshd (mp_int * a, int b) 1161189251Ssam{ 1162189251Ssam int x; 1163189251Ssam 1164189251Ssam /* if b <= 0 then ignore it */ 1165189251Ssam if (b <= 0) { 1166189251Ssam return; 1167189251Ssam } 1168189251Ssam 1169189251Ssam /* if b > used then simply zero it and return */ 1170189251Ssam if (a->used <= b) { 1171189251Ssam mp_zero (a); 1172189251Ssam return; 1173189251Ssam } 1174189251Ssam 1175189251Ssam { 1176189251Ssam register mp_digit *bottom, *top; 1177189251Ssam 1178189251Ssam /* shift the digits down */ 1179189251Ssam 1180189251Ssam /* bottom */ 1181189251Ssam bottom = a->dp; 1182189251Ssam 1183189251Ssam /* top [offset into digits] */ 1184189251Ssam top = a->dp + b; 1185189251Ssam 1186189251Ssam /* this is implemented as a sliding window where 1187189251Ssam * the window is b-digits long and digits from 1188189251Ssam * the top of the window are copied to the bottom 1189189251Ssam * 1190189251Ssam * e.g. 1191189251Ssam 1192189251Ssam b-2 | b-1 | b0 | b1 | b2 | ... | bb | ----> 1193189251Ssam /\ | ----> 1194189251Ssam \-------------------/ ----> 1195189251Ssam */ 1196189251Ssam for (x = 0; x < (a->used - b); x++) { 1197189251Ssam *bottom++ = *top++; 1198189251Ssam } 1199189251Ssam 1200189251Ssam /* zero the top digits */ 1201189251Ssam for (; x < a->used; x++) { 1202189251Ssam *bottom++ = 0; 1203189251Ssam } 1204189251Ssam } 1205189251Ssam 1206189251Ssam /* remove excess digits */ 1207189251Ssam a->used -= b; 1208189251Ssam} 1209189251Ssam 1210189251Ssam 1211189251Ssam/* swap the elements of two integers, for cases where you can't simply swap the 1212189251Ssam * mp_int pointers around 1213189251Ssam */ 1214189251Ssamstatic void mp_exch (mp_int * a, mp_int * b) 1215189251Ssam{ 1216189251Ssam mp_int t; 1217189251Ssam 1218189251Ssam t = *a; 1219189251Ssam *a = *b; 1220189251Ssam *b = t; 1221189251Ssam} 1222189251Ssam 1223189251Ssam 1224189251Ssam/* trim unused digits 1225189251Ssam * 1226189251Ssam * This is used to ensure that leading zero digits are 1227189251Ssam * trimed and the leading "used" digit will be non-zero 1228189251Ssam * Typically very fast. Also fixes the sign if there 1229189251Ssam * are no more leading digits 1230189251Ssam */ 1231189251Ssamstatic void mp_clamp (mp_int * a) 1232189251Ssam{ 1233189251Ssam /* decrease used while the most significant digit is 1234189251Ssam * zero. 1235189251Ssam */ 1236189251Ssam while (a->used > 0 && a->dp[a->used - 1] == 0) { 1237189251Ssam --(a->used); 1238189251Ssam } 1239189251Ssam 1240189251Ssam /* reset the sign flag if used == 0 */ 1241189251Ssam if (a->used == 0) { 1242189251Ssam a->sign = MP_ZPOS; 1243189251Ssam } 1244189251Ssam} 1245189251Ssam 1246189251Ssam 1247189251Ssam/* grow as required */ 1248189251Ssamstatic int mp_grow (mp_int * a, int size) 1249189251Ssam{ 1250189251Ssam int i; 1251189251Ssam mp_digit *tmp; 1252189251Ssam 1253189251Ssam /* if the alloc size is smaller alloc more ram */ 1254189251Ssam if (a->alloc < size) { 1255189251Ssam /* ensure there are always at least MP_PREC digits extra on top */ 1256189251Ssam size += (MP_PREC * 2) - (size % MP_PREC); 1257189251Ssam 1258189251Ssam /* reallocate the array a->dp 1259189251Ssam * 1260189251Ssam * We store the return in a temporary variable 1261189251Ssam * in case the operation failed we don't want 1262189251Ssam * to overwrite the dp member of a. 1263189251Ssam */ 1264189251Ssam tmp = OPT_CAST(mp_digit) XREALLOC (a->dp, sizeof (mp_digit) * size); 1265189251Ssam if (tmp == NULL) { 1266189251Ssam /* reallocation failed but "a" is still valid [can be freed] */ 1267189251Ssam return MP_MEM; 1268189251Ssam } 1269189251Ssam 1270189251Ssam /* reallocation succeeded so set a->dp */ 1271189251Ssam a->dp = tmp; 1272189251Ssam 1273189251Ssam /* zero excess digits */ 1274189251Ssam i = a->alloc; 1275189251Ssam a->alloc = size; 1276189251Ssam for (; i < a->alloc; i++) { 1277189251Ssam a->dp[i] = 0; 1278189251Ssam } 1279189251Ssam } 1280189251Ssam return MP_OKAY; 1281189251Ssam} 1282189251Ssam 1283189251Ssam 1284189251Ssam#ifdef BN_MP_ABS_C 1285189251Ssam/* b = |a| 1286189251Ssam * 1287189251Ssam * Simple function copies the input and fixes the sign to positive 1288189251Ssam */ 1289189251Ssamstatic int mp_abs (mp_int * a, mp_int * b) 1290189251Ssam{ 1291189251Ssam int res; 1292189251Ssam 1293189251Ssam /* copy a to b */ 1294189251Ssam if (a != b) { 1295189251Ssam if ((res = mp_copy (a, b)) != MP_OKAY) { 1296189251Ssam return res; 1297189251Ssam } 1298189251Ssam } 1299189251Ssam 1300189251Ssam /* force the sign of b to positive */ 1301189251Ssam b->sign = MP_ZPOS; 1302189251Ssam 1303189251Ssam return MP_OKAY; 1304189251Ssam} 1305189251Ssam#endif 1306189251Ssam 1307189251Ssam 1308189251Ssam/* set to a digit */ 1309189251Ssamstatic void mp_set (mp_int * a, mp_digit b) 1310189251Ssam{ 1311189251Ssam mp_zero (a); 1312189251Ssam a->dp[0] = b & MP_MASK; 1313189251Ssam a->used = (a->dp[0] != 0) ? 1 : 0; 1314189251Ssam} 1315189251Ssam 1316189251Ssam 1317189251Ssam#ifndef LTM_NO_NEG_EXP 1318189251Ssam/* b = a/2 */ 1319189251Ssamstatic int mp_div_2(mp_int * a, mp_int * b) 1320189251Ssam{ 1321189251Ssam int x, res, oldused; 1322189251Ssam 1323189251Ssam /* copy */ 1324189251Ssam if (b->alloc < a->used) { 1325189251Ssam if ((res = mp_grow (b, a->used)) != MP_OKAY) { 1326189251Ssam return res; 1327189251Ssam } 1328189251Ssam } 1329189251Ssam 1330189251Ssam oldused = b->used; 1331189251Ssam b->used = a->used; 1332189251Ssam { 1333189251Ssam register mp_digit r, rr, *tmpa, *tmpb; 1334189251Ssam 1335189251Ssam /* source alias */ 1336189251Ssam tmpa = a->dp + b->used - 1; 1337189251Ssam 1338189251Ssam /* dest alias */ 1339189251Ssam tmpb = b->dp + b->used - 1; 1340189251Ssam 1341189251Ssam /* carry */ 1342189251Ssam r = 0; 1343189251Ssam for (x = b->used - 1; x >= 0; x--) { 1344189251Ssam /* get the carry for the next iteration */ 1345189251Ssam rr = *tmpa & 1; 1346189251Ssam 1347189251Ssam /* shift the current digit, add in carry and store */ 1348189251Ssam *tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1)); 1349189251Ssam 1350189251Ssam /* forward carry to next iteration */ 1351189251Ssam r = rr; 1352189251Ssam } 1353189251Ssam 1354189251Ssam /* zero excess digits */ 1355189251Ssam tmpb = b->dp + b->used; 1356189251Ssam for (x = b->used; x < oldused; x++) { 1357189251Ssam *tmpb++ = 0; 1358189251Ssam } 1359189251Ssam } 1360189251Ssam b->sign = a->sign; 1361189251Ssam mp_clamp (b); 1362189251Ssam return MP_OKAY; 1363189251Ssam} 1364189251Ssam#endif /* LTM_NO_NEG_EXP */ 1365189251Ssam 1366189251Ssam 1367189251Ssam/* shift left by a certain bit count */ 1368189251Ssamstatic int mp_mul_2d (mp_int * a, int b, mp_int * c) 1369189251Ssam{ 1370189251Ssam mp_digit d; 1371189251Ssam int res; 1372189251Ssam 1373189251Ssam /* copy */ 1374189251Ssam if (a != c) { 1375189251Ssam if ((res = mp_copy (a, c)) != MP_OKAY) { 1376189251Ssam return res; 1377189251Ssam } 1378189251Ssam } 1379189251Ssam 1380189251Ssam if (c->alloc < (int)(c->used + b/DIGIT_BIT + 1)) { 1381189251Ssam if ((res = mp_grow (c, c->used + b / DIGIT_BIT + 1)) != MP_OKAY) { 1382189251Ssam return res; 1383189251Ssam } 1384189251Ssam } 1385189251Ssam 1386189251Ssam /* shift by as many digits in the bit count */ 1387189251Ssam if (b >= (int)DIGIT_BIT) { 1388189251Ssam if ((res = mp_lshd (c, b / DIGIT_BIT)) != MP_OKAY) { 1389189251Ssam return res; 1390189251Ssam } 1391189251Ssam } 1392189251Ssam 1393189251Ssam /* shift any bit count < DIGIT_BIT */ 1394189251Ssam d = (mp_digit) (b % DIGIT_BIT); 1395189251Ssam if (d != 0) { 1396189251Ssam register mp_digit *tmpc, shift, mask, r, rr; 1397189251Ssam register int x; 1398189251Ssam 1399189251Ssam /* bitmask for carries */ 1400189251Ssam mask = (((mp_digit)1) << d) - 1; 1401189251Ssam 1402189251Ssam /* shift for msbs */ 1403189251Ssam shift = DIGIT_BIT - d; 1404189251Ssam 1405189251Ssam /* alias */ 1406189251Ssam tmpc = c->dp; 1407189251Ssam 1408189251Ssam /* carry */ 1409189251Ssam r = 0; 1410189251Ssam for (x = 0; x < c->used; x++) { 1411189251Ssam /* get the higher bits of the current word */ 1412189251Ssam rr = (*tmpc >> shift) & mask; 1413189251Ssam 1414189251Ssam /* shift the current word and OR in the carry */ 1415189251Ssam *tmpc = ((*tmpc << d) | r) & MP_MASK; 1416189251Ssam ++tmpc; 1417189251Ssam 1418189251Ssam /* set the carry to the carry bits of the current word */ 1419189251Ssam r = rr; 1420189251Ssam } 1421189251Ssam 1422189251Ssam /* set final carry */ 1423189251Ssam if (r != 0) { 1424189251Ssam c->dp[(c->used)++] = r; 1425189251Ssam } 1426189251Ssam } 1427189251Ssam mp_clamp (c); 1428189251Ssam return MP_OKAY; 1429189251Ssam} 1430189251Ssam 1431189251Ssam 1432189251Ssam#ifdef BN_MP_INIT_MULTI_C 1433189251Ssamstatic int mp_init_multi(mp_int *mp, ...) 1434189251Ssam{ 1435189251Ssam mp_err res = MP_OKAY; /* Assume ok until proven otherwise */ 1436189251Ssam int n = 0; /* Number of ok inits */ 1437189251Ssam mp_int* cur_arg = mp; 1438189251Ssam va_list args; 1439189251Ssam 1440189251Ssam va_start(args, mp); /* init args to next argument from caller */ 1441189251Ssam while (cur_arg != NULL) { 1442189251Ssam if (mp_init(cur_arg) != MP_OKAY) { 1443189251Ssam /* Oops - error! Back-track and mp_clear what we already 1444189251Ssam succeeded in init-ing, then return error. 1445189251Ssam */ 1446189251Ssam va_list clean_args; 1447189251Ssam 1448189251Ssam /* end the current list */ 1449189251Ssam va_end(args); 1450189251Ssam 1451189251Ssam /* now start cleaning up */ 1452189251Ssam cur_arg = mp; 1453189251Ssam va_start(clean_args, mp); 1454189251Ssam while (n--) { 1455189251Ssam mp_clear(cur_arg); 1456189251Ssam cur_arg = va_arg(clean_args, mp_int*); 1457189251Ssam } 1458189251Ssam va_end(clean_args); 1459189251Ssam res = MP_MEM; 1460189251Ssam break; 1461189251Ssam } 1462189251Ssam n++; 1463189251Ssam cur_arg = va_arg(args, mp_int*); 1464189251Ssam } 1465189251Ssam va_end(args); 1466189251Ssam return res; /* Assumed ok, if error flagged above. */ 1467189251Ssam} 1468189251Ssam#endif 1469189251Ssam 1470189251Ssam 1471189251Ssam#ifdef BN_MP_CLEAR_MULTI_C 1472189251Ssamstatic void mp_clear_multi(mp_int *mp, ...) 1473189251Ssam{ 1474189251Ssam mp_int* next_mp = mp; 1475189251Ssam va_list args; 1476189251Ssam va_start(args, mp); 1477189251Ssam while (next_mp != NULL) { 1478189251Ssam mp_clear(next_mp); 1479189251Ssam next_mp = va_arg(args, mp_int*); 1480189251Ssam } 1481189251Ssam va_end(args); 1482189251Ssam} 1483189251Ssam#endif 1484189251Ssam 1485189251Ssam 1486189251Ssam/* shift left a certain amount of digits */ 1487189251Ssamstatic int mp_lshd (mp_int * a, int b) 1488189251Ssam{ 1489189251Ssam int x, res; 1490189251Ssam 1491189251Ssam /* if its less than zero return */ 1492189251Ssam if (b <= 0) { 1493189251Ssam return MP_OKAY; 1494189251Ssam } 1495189251Ssam 1496189251Ssam /* grow to fit the new digits */ 1497189251Ssam if (a->alloc < a->used + b) { 1498189251Ssam if ((res = mp_grow (a, a->used + b)) != MP_OKAY) { 1499189251Ssam return res; 1500189251Ssam } 1501189251Ssam } 1502189251Ssam 1503189251Ssam { 1504189251Ssam register mp_digit *top, *bottom; 1505189251Ssam 1506189251Ssam /* increment the used by the shift amount then copy upwards */ 1507189251Ssam a->used += b; 1508189251Ssam 1509189251Ssam /* top */ 1510189251Ssam top = a->dp + a->used - 1; 1511189251Ssam 1512189251Ssam /* base */ 1513189251Ssam bottom = a->dp + a->used - 1 - b; 1514189251Ssam 1515189251Ssam /* much like mp_rshd this is implemented using a sliding window 1516189251Ssam * except the window goes the otherway around. Copying from 1517189251Ssam * the bottom to the top. see bn_mp_rshd.c for more info. 1518189251Ssam */ 1519189251Ssam for (x = a->used - 1; x >= b; x--) { 1520189251Ssam *top-- = *bottom--; 1521189251Ssam } 1522189251Ssam 1523189251Ssam /* zero the lower digits */ 1524189251Ssam top = a->dp; 1525189251Ssam for (x = 0; x < b; x++) { 1526189251Ssam *top++ = 0; 1527189251Ssam } 1528189251Ssam } 1529189251Ssam return MP_OKAY; 1530189251Ssam} 1531189251Ssam 1532189251Ssam 1533189251Ssam/* returns the number of bits in an int */ 1534189251Ssamstatic int mp_count_bits (mp_int * a) 1535189251Ssam{ 1536189251Ssam int r; 1537189251Ssam mp_digit q; 1538189251Ssam 1539189251Ssam /* shortcut */ 1540189251Ssam if (a->used == 0) { 1541189251Ssam return 0; 1542189251Ssam } 1543189251Ssam 1544189251Ssam /* get number of digits and add that */ 1545189251Ssam r = (a->used - 1) * DIGIT_BIT; 1546189251Ssam 1547189251Ssam /* take the last digit and count the bits in it */ 1548189251Ssam q = a->dp[a->used - 1]; 1549189251Ssam while (q > ((mp_digit) 0)) { 1550189251Ssam ++r; 1551189251Ssam q >>= ((mp_digit) 1); 1552189251Ssam } 1553189251Ssam return r; 1554189251Ssam} 1555189251Ssam 1556189251Ssam 1557189251Ssam/* calc a value mod 2**b */ 1558189251Ssamstatic int mp_mod_2d (mp_int * a, int b, mp_int * c) 1559189251Ssam{ 1560189251Ssam int x, res; 1561189251Ssam 1562189251Ssam /* if b is <= 0 then zero the int */ 1563189251Ssam if (b <= 0) { 1564189251Ssam mp_zero (c); 1565189251Ssam return MP_OKAY; 1566189251Ssam } 1567189251Ssam 1568189251Ssam /* if the modulus is larger than the value than return */ 1569189251Ssam if (b >= (int) (a->used * DIGIT_BIT)) { 1570189251Ssam res = mp_copy (a, c); 1571189251Ssam return res; 1572189251Ssam } 1573189251Ssam 1574189251Ssam /* copy */ 1575189251Ssam if ((res = mp_copy (a, c)) != MP_OKAY) { 1576189251Ssam return res; 1577189251Ssam } 1578189251Ssam 1579189251Ssam /* zero digits above the last digit of the modulus */ 1580189251Ssam for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++) { 1581189251Ssam c->dp[x] = 0; 1582189251Ssam } 1583189251Ssam /* clear the digit that is not completely outside/inside the modulus */ 1584189251Ssam c->dp[b / DIGIT_BIT] &= 1585189251Ssam (mp_digit) ((((mp_digit) 1) << (((mp_digit) b) % DIGIT_BIT)) - ((mp_digit) 1)); 1586189251Ssam mp_clamp (c); 1587189251Ssam return MP_OKAY; 1588189251Ssam} 1589189251Ssam 1590189251Ssam 1591189251Ssam#ifdef BN_MP_DIV_SMALL 1592189251Ssam 1593189251Ssam/* slower bit-bang division... also smaller */ 1594189251Ssamstatic int mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d) 1595189251Ssam{ 1596189251Ssam mp_int ta, tb, tq, q; 1597189251Ssam int res, n, n2; 1598189251Ssam 1599189251Ssam /* is divisor zero ? */ 1600189251Ssam if (mp_iszero (b) == 1) { 1601189251Ssam return MP_VAL; 1602189251Ssam } 1603189251Ssam 1604189251Ssam /* if a < b then q=0, r = a */ 1605189251Ssam if (mp_cmp_mag (a, b) == MP_LT) { 1606189251Ssam if (d != NULL) { 1607189251Ssam res = mp_copy (a, d); 1608189251Ssam } else { 1609189251Ssam res = MP_OKAY; 1610189251Ssam } 1611189251Ssam if (c != NULL) { 1612189251Ssam mp_zero (c); 1613189251Ssam } 1614189251Ssam return res; 1615189251Ssam } 1616189251Ssam 1617189251Ssam /* init our temps */ 1618189251Ssam if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL) != MP_OKAY)) { 1619189251Ssam return res; 1620189251Ssam } 1621189251Ssam 1622189251Ssam 1623189251Ssam mp_set(&tq, 1); 1624189251Ssam n = mp_count_bits(a) - mp_count_bits(b); 1625189251Ssam if (((res = mp_abs(a, &ta)) != MP_OKAY) || 1626189251Ssam ((res = mp_abs(b, &tb)) != MP_OKAY) || 1627189251Ssam ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) || 1628189251Ssam ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) { 1629189251Ssam goto LBL_ERR; 1630189251Ssam } 1631189251Ssam 1632189251Ssam while (n-- >= 0) { 1633189251Ssam if (mp_cmp(&tb, &ta) != MP_GT) { 1634189251Ssam if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) || 1635189251Ssam ((res = mp_add(&q, &tq, &q)) != MP_OKAY)) { 1636189251Ssam goto LBL_ERR; 1637189251Ssam } 1638189251Ssam } 1639189251Ssam if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) || 1640189251Ssam ((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) { 1641189251Ssam goto LBL_ERR; 1642189251Ssam } 1643189251Ssam } 1644189251Ssam 1645189251Ssam /* now q == quotient and ta == remainder */ 1646189251Ssam n = a->sign; 1647189251Ssam n2 = (a->sign == b->sign ? MP_ZPOS : MP_NEG); 1648189251Ssam if (c != NULL) { 1649189251Ssam mp_exch(c, &q); 1650189251Ssam c->sign = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2; 1651189251Ssam } 1652189251Ssam if (d != NULL) { 1653189251Ssam mp_exch(d, &ta); 1654189251Ssam d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n; 1655189251Ssam } 1656189251SsamLBL_ERR: 1657189251Ssam mp_clear_multi(&ta, &tb, &tq, &q, NULL); 1658189251Ssam return res; 1659189251Ssam} 1660189251Ssam 1661189251Ssam#else 1662189251Ssam 1663189251Ssam/* integer signed division. 1664189251Ssam * c*b + d == a [e.g. a/b, c=quotient, d=remainder] 1665189251Ssam * HAC pp.598 Algorithm 14.20 1666189251Ssam * 1667189251Ssam * Note that the description in HAC is horribly 1668189251Ssam * incomplete. For example, it doesn't consider 1669189251Ssam * the case where digits are removed from 'x' in 1670189251Ssam * the inner loop. It also doesn't consider the 1671189251Ssam * case that y has fewer than three digits, etc.. 1672189251Ssam * 1673189251Ssam * The overall algorithm is as described as 1674189251Ssam * 14.20 from HAC but fixed to treat these cases. 1675189251Ssam*/ 1676189251Ssamstatic int mp_div (mp_int * a, mp_int * b, mp_int * c, mp_int * d) 1677189251Ssam{ 1678189251Ssam mp_int q, x, y, t1, t2; 1679189251Ssam int res, n, t, i, norm, neg; 1680189251Ssam 1681189251Ssam /* is divisor zero ? */ 1682189251Ssam if (mp_iszero (b) == 1) { 1683189251Ssam return MP_VAL; 1684189251Ssam } 1685189251Ssam 1686189251Ssam /* if a < b then q=0, r = a */ 1687189251Ssam if (mp_cmp_mag (a, b) == MP_LT) { 1688189251Ssam if (d != NULL) { 1689189251Ssam res = mp_copy (a, d); 1690189251Ssam } else { 1691189251Ssam res = MP_OKAY; 1692189251Ssam } 1693189251Ssam if (c != NULL) { 1694189251Ssam mp_zero (c); 1695189251Ssam } 1696189251Ssam return res; 1697189251Ssam } 1698189251Ssam 1699189251Ssam if ((res = mp_init_size (&q, a->used + 2)) != MP_OKAY) { 1700189251Ssam return res; 1701189251Ssam } 1702189251Ssam q.used = a->used + 2; 1703189251Ssam 1704189251Ssam if ((res = mp_init (&t1)) != MP_OKAY) { 1705189251Ssam goto LBL_Q; 1706189251Ssam } 1707189251Ssam 1708189251Ssam if ((res = mp_init (&t2)) != MP_OKAY) { 1709189251Ssam goto LBL_T1; 1710189251Ssam } 1711189251Ssam 1712189251Ssam if ((res = mp_init_copy (&x, a)) != MP_OKAY) { 1713189251Ssam goto LBL_T2; 1714189251Ssam } 1715189251Ssam 1716189251Ssam if ((res = mp_init_copy (&y, b)) != MP_OKAY) { 1717189251Ssam goto LBL_X; 1718189251Ssam } 1719189251Ssam 1720189251Ssam /* fix the sign */ 1721189251Ssam neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG; 1722189251Ssam x.sign = y.sign = MP_ZPOS; 1723189251Ssam 1724189251Ssam /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */ 1725189251Ssam norm = mp_count_bits(&y) % DIGIT_BIT; 1726189251Ssam if (norm < (int)(DIGIT_BIT-1)) { 1727189251Ssam norm = (DIGIT_BIT-1) - norm; 1728189251Ssam if ((res = mp_mul_2d (&x, norm, &x)) != MP_OKAY) { 1729189251Ssam goto LBL_Y; 1730189251Ssam } 1731189251Ssam if ((res = mp_mul_2d (&y, norm, &y)) != MP_OKAY) { 1732189251Ssam goto LBL_Y; 1733189251Ssam } 1734189251Ssam } else { 1735189251Ssam norm = 0; 1736189251Ssam } 1737189251Ssam 1738189251Ssam /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */ 1739189251Ssam n = x.used - 1; 1740189251Ssam t = y.used - 1; 1741189251Ssam 1742189251Ssam /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */ 1743189251Ssam if ((res = mp_lshd (&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */ 1744189251Ssam goto LBL_Y; 1745189251Ssam } 1746189251Ssam 1747189251Ssam while (mp_cmp (&x, &y) != MP_LT) { 1748189251Ssam ++(q.dp[n - t]); 1749189251Ssam if ((res = mp_sub (&x, &y, &x)) != MP_OKAY) { 1750189251Ssam goto LBL_Y; 1751189251Ssam } 1752189251Ssam } 1753189251Ssam 1754189251Ssam /* reset y by shifting it back down */ 1755189251Ssam mp_rshd (&y, n - t); 1756189251Ssam 1757189251Ssam /* step 3. for i from n down to (t + 1) */ 1758189251Ssam for (i = n; i >= (t + 1); i--) { 1759189251Ssam if (i > x.used) { 1760189251Ssam continue; 1761189251Ssam } 1762189251Ssam 1763189251Ssam /* step 3.1 if xi == yt then set q{i-t-1} to b-1, 1764189251Ssam * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */ 1765189251Ssam if (x.dp[i] == y.dp[t]) { 1766189251Ssam q.dp[i - t - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1); 1767189251Ssam } else { 1768189251Ssam mp_word tmp; 1769189251Ssam tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT); 1770189251Ssam tmp |= ((mp_word) x.dp[i - 1]); 1771189251Ssam tmp /= ((mp_word) y.dp[t]); 1772189251Ssam if (tmp > (mp_word) MP_MASK) 1773189251Ssam tmp = MP_MASK; 1774189251Ssam q.dp[i - t - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK)); 1775189251Ssam } 1776189251Ssam 1777189251Ssam /* while (q{i-t-1} * (yt * b + y{t-1})) > 1778189251Ssam xi * b**2 + xi-1 * b + xi-2 1779189251Ssam 1780189251Ssam do q{i-t-1} -= 1; 1781189251Ssam */ 1782189251Ssam q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK; 1783189251Ssam do { 1784189251Ssam q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK; 1785189251Ssam 1786189251Ssam /* find left hand */ 1787189251Ssam mp_zero (&t1); 1788189251Ssam t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1]; 1789189251Ssam t1.dp[1] = y.dp[t]; 1790189251Ssam t1.used = 2; 1791189251Ssam if ((res = mp_mul_d (&t1, q.dp[i - t - 1], &t1)) != MP_OKAY) { 1792189251Ssam goto LBL_Y; 1793189251Ssam } 1794189251Ssam 1795189251Ssam /* find right hand */ 1796189251Ssam t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2]; 1797189251Ssam t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1]; 1798189251Ssam t2.dp[2] = x.dp[i]; 1799189251Ssam t2.used = 3; 1800189251Ssam } while (mp_cmp_mag(&t1, &t2) == MP_GT); 1801189251Ssam 1802189251Ssam /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */ 1803189251Ssam if ((res = mp_mul_d (&y, q.dp[i - t - 1], &t1)) != MP_OKAY) { 1804189251Ssam goto LBL_Y; 1805189251Ssam } 1806189251Ssam 1807189251Ssam if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) { 1808189251Ssam goto LBL_Y; 1809189251Ssam } 1810189251Ssam 1811189251Ssam if ((res = mp_sub (&x, &t1, &x)) != MP_OKAY) { 1812189251Ssam goto LBL_Y; 1813189251Ssam } 1814189251Ssam 1815189251Ssam /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */ 1816189251Ssam if (x.sign == MP_NEG) { 1817189251Ssam if ((res = mp_copy (&y, &t1)) != MP_OKAY) { 1818189251Ssam goto LBL_Y; 1819189251Ssam } 1820189251Ssam if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) { 1821189251Ssam goto LBL_Y; 1822189251Ssam } 1823189251Ssam if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) { 1824189251Ssam goto LBL_Y; 1825189251Ssam } 1826189251Ssam 1827189251Ssam q.dp[i - t - 1] = (q.dp[i - t - 1] - 1UL) & MP_MASK; 1828189251Ssam } 1829189251Ssam } 1830189251Ssam 1831189251Ssam /* now q is the quotient and x is the remainder 1832189251Ssam * [which we have to normalize] 1833189251Ssam */ 1834189251Ssam 1835189251Ssam /* get sign before writing to c */ 1836189251Ssam x.sign = x.used == 0 ? MP_ZPOS : a->sign; 1837189251Ssam 1838189251Ssam if (c != NULL) { 1839189251Ssam mp_clamp (&q); 1840189251Ssam mp_exch (&q, c); 1841189251Ssam c->sign = neg; 1842189251Ssam } 1843189251Ssam 1844189251Ssam if (d != NULL) { 1845189251Ssam mp_div_2d (&x, norm, &x, NULL); 1846189251Ssam mp_exch (&x, d); 1847189251Ssam } 1848189251Ssam 1849189251Ssam res = MP_OKAY; 1850189251Ssam 1851189251SsamLBL_Y:mp_clear (&y); 1852189251SsamLBL_X:mp_clear (&x); 1853189251SsamLBL_T2:mp_clear (&t2); 1854189251SsamLBL_T1:mp_clear (&t1); 1855189251SsamLBL_Q:mp_clear (&q); 1856189251Ssam return res; 1857189251Ssam} 1858189251Ssam 1859189251Ssam#endif 1860189251Ssam 1861189251Ssam 1862189251Ssam#ifdef MP_LOW_MEM 1863189251Ssam #define TAB_SIZE 32 1864189251Ssam#else 1865189251Ssam #define TAB_SIZE 256 1866189251Ssam#endif 1867189251Ssam 1868189251Ssamstatic int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode) 1869189251Ssam{ 1870189251Ssam mp_int M[TAB_SIZE], res, mu; 1871189251Ssam mp_digit buf; 1872189251Ssam int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize; 1873189251Ssam int (*redux)(mp_int*,mp_int*,mp_int*); 1874189251Ssam 1875189251Ssam /* find window size */ 1876189251Ssam x = mp_count_bits (X); 1877189251Ssam if (x <= 7) { 1878189251Ssam winsize = 2; 1879189251Ssam } else if (x <= 36) { 1880189251Ssam winsize = 3; 1881189251Ssam } else if (x <= 140) { 1882189251Ssam winsize = 4; 1883189251Ssam } else if (x <= 450) { 1884189251Ssam winsize = 5; 1885189251Ssam } else if (x <= 1303) { 1886189251Ssam winsize = 6; 1887189251Ssam } else if (x <= 3529) { 1888189251Ssam winsize = 7; 1889189251Ssam } else { 1890189251Ssam winsize = 8; 1891189251Ssam } 1892189251Ssam 1893189251Ssam#ifdef MP_LOW_MEM 1894189251Ssam if (winsize > 5) { 1895189251Ssam winsize = 5; 1896189251Ssam } 1897189251Ssam#endif 1898189251Ssam 1899189251Ssam /* init M array */ 1900189251Ssam /* init first cell */ 1901189251Ssam if ((err = mp_init(&M[1])) != MP_OKAY) { 1902189251Ssam return err; 1903189251Ssam } 1904189251Ssam 1905189251Ssam /* now init the second half of the array */ 1906189251Ssam for (x = 1<<(winsize-1); x < (1 << winsize); x++) { 1907189251Ssam if ((err = mp_init(&M[x])) != MP_OKAY) { 1908189251Ssam for (y = 1<<(winsize-1); y < x; y++) { 1909189251Ssam mp_clear (&M[y]); 1910189251Ssam } 1911189251Ssam mp_clear(&M[1]); 1912189251Ssam return err; 1913189251Ssam } 1914189251Ssam } 1915189251Ssam 1916189251Ssam /* create mu, used for Barrett reduction */ 1917189251Ssam if ((err = mp_init (&mu)) != MP_OKAY) { 1918189251Ssam goto LBL_M; 1919189251Ssam } 1920189251Ssam 1921189251Ssam if (redmode == 0) { 1922189251Ssam if ((err = mp_reduce_setup (&mu, P)) != MP_OKAY) { 1923189251Ssam goto LBL_MU; 1924189251Ssam } 1925189251Ssam redux = mp_reduce; 1926189251Ssam } else { 1927189251Ssam if ((err = mp_reduce_2k_setup_l (P, &mu)) != MP_OKAY) { 1928189251Ssam goto LBL_MU; 1929189251Ssam } 1930189251Ssam redux = mp_reduce_2k_l; 1931189251Ssam } 1932189251Ssam 1933189251Ssam /* create M table 1934189251Ssam * 1935189251Ssam * The M table contains powers of the base, 1936189251Ssam * e.g. M[x] = G**x mod P 1937189251Ssam * 1938189251Ssam * The first half of the table is not 1939189251Ssam * computed though accept for M[0] and M[1] 1940189251Ssam */ 1941189251Ssam if ((err = mp_mod (G, P, &M[1])) != MP_OKAY) { 1942189251Ssam goto LBL_MU; 1943189251Ssam } 1944189251Ssam 1945189251Ssam /* compute the value at M[1<<(winsize-1)] by squaring 1946189251Ssam * M[1] (winsize-1) times 1947189251Ssam */ 1948189251Ssam if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) { 1949189251Ssam goto LBL_MU; 1950189251Ssam } 1951189251Ssam 1952189251Ssam for (x = 0; x < (winsize - 1); x++) { 1953189251Ssam /* square it */ 1954189251Ssam if ((err = mp_sqr (&M[1 << (winsize - 1)], 1955189251Ssam &M[1 << (winsize - 1)])) != MP_OKAY) { 1956189251Ssam goto LBL_MU; 1957189251Ssam } 1958189251Ssam 1959189251Ssam /* reduce modulo P */ 1960189251Ssam if ((err = redux (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) { 1961189251Ssam goto LBL_MU; 1962189251Ssam } 1963189251Ssam } 1964189251Ssam 1965189251Ssam /* create upper table, that is M[x] = M[x-1] * M[1] (mod P) 1966189251Ssam * for x = (2**(winsize - 1) + 1) to (2**winsize - 1) 1967189251Ssam */ 1968189251Ssam for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) { 1969189251Ssam if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) { 1970189251Ssam goto LBL_MU; 1971189251Ssam } 1972189251Ssam if ((err = redux (&M[x], P, &mu)) != MP_OKAY) { 1973189251Ssam goto LBL_MU; 1974189251Ssam } 1975189251Ssam } 1976189251Ssam 1977189251Ssam /* setup result */ 1978189251Ssam if ((err = mp_init (&res)) != MP_OKAY) { 1979189251Ssam goto LBL_MU; 1980189251Ssam } 1981189251Ssam mp_set (&res, 1); 1982189251Ssam 1983189251Ssam /* set initial mode and bit cnt */ 1984189251Ssam mode = 0; 1985189251Ssam bitcnt = 1; 1986189251Ssam buf = 0; 1987189251Ssam digidx = X->used - 1; 1988189251Ssam bitcpy = 0; 1989189251Ssam bitbuf = 0; 1990189251Ssam 1991189251Ssam for (;;) { 1992189251Ssam /* grab next digit as required */ 1993189251Ssam if (--bitcnt == 0) { 1994189251Ssam /* if digidx == -1 we are out of digits */ 1995189251Ssam if (digidx == -1) { 1996189251Ssam break; 1997189251Ssam } 1998189251Ssam /* read next digit and reset the bitcnt */ 1999189251Ssam buf = X->dp[digidx--]; 2000189251Ssam bitcnt = (int) DIGIT_BIT; 2001189251Ssam } 2002189251Ssam 2003189251Ssam /* grab the next msb from the exponent */ 2004189251Ssam y = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1; 2005189251Ssam buf <<= (mp_digit)1; 2006189251Ssam 2007189251Ssam /* if the bit is zero and mode == 0 then we ignore it 2008189251Ssam * These represent the leading zero bits before the first 1 bit 2009189251Ssam * in the exponent. Technically this opt is not required but it 2010189251Ssam * does lower the # of trivial squaring/reductions used 2011189251Ssam */ 2012189251Ssam if (mode == 0 && y == 0) { 2013189251Ssam continue; 2014189251Ssam } 2015189251Ssam 2016189251Ssam /* if the bit is zero and mode == 1 then we square */ 2017189251Ssam if (mode == 1 && y == 0) { 2018189251Ssam if ((err = mp_sqr (&res, &res)) != MP_OKAY) { 2019189251Ssam goto LBL_RES; 2020189251Ssam } 2021189251Ssam if ((err = redux (&res, P, &mu)) != MP_OKAY) { 2022189251Ssam goto LBL_RES; 2023189251Ssam } 2024189251Ssam continue; 2025189251Ssam } 2026189251Ssam 2027189251Ssam /* else we add it to the window */ 2028189251Ssam bitbuf |= (y << (winsize - ++bitcpy)); 2029189251Ssam mode = 2; 2030189251Ssam 2031189251Ssam if (bitcpy == winsize) { 2032189251Ssam /* ok window is filled so square as required and multiply */ 2033189251Ssam /* square first */ 2034189251Ssam for (x = 0; x < winsize; x++) { 2035189251Ssam if ((err = mp_sqr (&res, &res)) != MP_OKAY) { 2036189251Ssam goto LBL_RES; 2037189251Ssam } 2038189251Ssam if ((err = redux (&res, P, &mu)) != MP_OKAY) { 2039189251Ssam goto LBL_RES; 2040189251Ssam } 2041189251Ssam } 2042189251Ssam 2043189251Ssam /* then multiply */ 2044189251Ssam if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) { 2045189251Ssam goto LBL_RES; 2046189251Ssam } 2047189251Ssam if ((err = redux (&res, P, &mu)) != MP_OKAY) { 2048189251Ssam goto LBL_RES; 2049189251Ssam } 2050189251Ssam 2051189251Ssam /* empty window and reset */ 2052189251Ssam bitcpy = 0; 2053189251Ssam bitbuf = 0; 2054189251Ssam mode = 1; 2055189251Ssam } 2056189251Ssam } 2057189251Ssam 2058189251Ssam /* if bits remain then square/multiply */ 2059189251Ssam if (mode == 2 && bitcpy > 0) { 2060189251Ssam /* square then multiply if the bit is set */ 2061189251Ssam for (x = 0; x < bitcpy; x++) { 2062189251Ssam if ((err = mp_sqr (&res, &res)) != MP_OKAY) { 2063189251Ssam goto LBL_RES; 2064189251Ssam } 2065189251Ssam if ((err = redux (&res, P, &mu)) != MP_OKAY) { 2066189251Ssam goto LBL_RES; 2067189251Ssam } 2068189251Ssam 2069189251Ssam bitbuf <<= 1; 2070189251Ssam if ((bitbuf & (1 << winsize)) != 0) { 2071189251Ssam /* then multiply */ 2072189251Ssam if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) { 2073189251Ssam goto LBL_RES; 2074189251Ssam } 2075189251Ssam if ((err = redux (&res, P, &mu)) != MP_OKAY) { 2076189251Ssam goto LBL_RES; 2077189251Ssam } 2078189251Ssam } 2079189251Ssam } 2080189251Ssam } 2081189251Ssam 2082189251Ssam mp_exch (&res, Y); 2083189251Ssam err = MP_OKAY; 2084189251SsamLBL_RES:mp_clear (&res); 2085189251SsamLBL_MU:mp_clear (&mu); 2086189251SsamLBL_M: 2087189251Ssam mp_clear(&M[1]); 2088189251Ssam for (x = 1<<(winsize-1); x < (1 << winsize); x++) { 2089189251Ssam mp_clear (&M[x]); 2090189251Ssam } 2091189251Ssam return err; 2092189251Ssam} 2093189251Ssam 2094189251Ssam 2095189251Ssam/* computes b = a*a */ 2096189251Ssamstatic int mp_sqr (mp_int * a, mp_int * b) 2097189251Ssam{ 2098189251Ssam int res; 2099189251Ssam 2100189251Ssam#ifdef BN_MP_TOOM_SQR_C 2101189251Ssam /* use Toom-Cook? */ 2102189251Ssam if (a->used >= TOOM_SQR_CUTOFF) { 2103189251Ssam res = mp_toom_sqr(a, b); 2104189251Ssam /* Karatsuba? */ 2105189251Ssam } else 2106189251Ssam#endif 2107189251Ssam#ifdef BN_MP_KARATSUBA_SQR_C 2108189251Ssamif (a->used >= KARATSUBA_SQR_CUTOFF) { 2109189251Ssam res = mp_karatsuba_sqr (a, b); 2110189251Ssam } else 2111189251Ssam#endif 2112189251Ssam { 2113189251Ssam#ifdef BN_FAST_S_MP_SQR_C 2114189251Ssam /* can we use the fast comba multiplier? */ 2115189251Ssam if ((a->used * 2 + 1) < MP_WARRAY && 2116189251Ssam a->used < 2117189251Ssam (1 << (sizeof(mp_word) * CHAR_BIT - 2*DIGIT_BIT - 1))) { 2118189251Ssam res = fast_s_mp_sqr (a, b); 2119189251Ssam } else 2120189251Ssam#endif 2121189251Ssam#ifdef BN_S_MP_SQR_C 2122189251Ssam res = s_mp_sqr (a, b); 2123189251Ssam#else 2124189251Ssam#error mp_sqr could fail 2125189251Ssam res = MP_VAL; 2126189251Ssam#endif 2127189251Ssam } 2128189251Ssam b->sign = MP_ZPOS; 2129189251Ssam return res; 2130189251Ssam} 2131189251Ssam 2132189251Ssam 2133189251Ssam/* reduces a modulo n where n is of the form 2**p - d 2134189251Ssam This differs from reduce_2k since "d" can be larger 2135189251Ssam than a single digit. 2136189251Ssam*/ 2137189251Ssamstatic int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d) 2138189251Ssam{ 2139189251Ssam mp_int q; 2140189251Ssam int p, res; 2141189251Ssam 2142189251Ssam if ((res = mp_init(&q)) != MP_OKAY) { 2143189251Ssam return res; 2144189251Ssam } 2145189251Ssam 2146189251Ssam p = mp_count_bits(n); 2147189251Ssamtop: 2148189251Ssam /* q = a/2**p, a = a mod 2**p */ 2149189251Ssam if ((res = mp_div_2d(a, p, &q, a)) != MP_OKAY) { 2150189251Ssam goto ERR; 2151189251Ssam } 2152189251Ssam 2153189251Ssam /* q = q * d */ 2154189251Ssam if ((res = mp_mul(&q, d, &q)) != MP_OKAY) { 2155189251Ssam goto ERR; 2156189251Ssam } 2157189251Ssam 2158189251Ssam /* a = a + q */ 2159189251Ssam if ((res = s_mp_add(a, &q, a)) != MP_OKAY) { 2160189251Ssam goto ERR; 2161189251Ssam } 2162189251Ssam 2163189251Ssam if (mp_cmp_mag(a, n) != MP_LT) { 2164189251Ssam s_mp_sub(a, n, a); 2165189251Ssam goto top; 2166189251Ssam } 2167189251Ssam 2168189251SsamERR: 2169189251Ssam mp_clear(&q); 2170189251Ssam return res; 2171189251Ssam} 2172189251Ssam 2173189251Ssam 2174189251Ssam/* determines the setup value */ 2175189251Ssamstatic int mp_reduce_2k_setup_l(mp_int *a, mp_int *d) 2176189251Ssam{ 2177189251Ssam int res; 2178189251Ssam mp_int tmp; 2179189251Ssam 2180189251Ssam if ((res = mp_init(&tmp)) != MP_OKAY) { 2181189251Ssam return res; 2182189251Ssam } 2183189251Ssam 2184189251Ssam if ((res = mp_2expt(&tmp, mp_count_bits(a))) != MP_OKAY) { 2185189251Ssam goto ERR; 2186189251Ssam } 2187189251Ssam 2188189251Ssam if ((res = s_mp_sub(&tmp, a, d)) != MP_OKAY) { 2189189251Ssam goto ERR; 2190189251Ssam } 2191189251Ssam 2192189251SsamERR: 2193189251Ssam mp_clear(&tmp); 2194189251Ssam return res; 2195189251Ssam} 2196189251Ssam 2197189251Ssam 2198189251Ssam/* computes a = 2**b 2199189251Ssam * 2200189251Ssam * Simple algorithm which zeroes the int, grows it then just sets one bit 2201189251Ssam * as required. 2202189251Ssam */ 2203189251Ssamstatic int mp_2expt (mp_int * a, int b) 2204189251Ssam{ 2205189251Ssam int res; 2206189251Ssam 2207189251Ssam /* zero a as per default */ 2208189251Ssam mp_zero (a); 2209189251Ssam 2210189251Ssam /* grow a to accomodate the single bit */ 2211189251Ssam if ((res = mp_grow (a, b / DIGIT_BIT + 1)) != MP_OKAY) { 2212189251Ssam return res; 2213189251Ssam } 2214189251Ssam 2215189251Ssam /* set the used count of where the bit will go */ 2216189251Ssam a->used = b / DIGIT_BIT + 1; 2217189251Ssam 2218189251Ssam /* put the single bit in its place */ 2219189251Ssam a->dp[b / DIGIT_BIT] = ((mp_digit)1) << (b % DIGIT_BIT); 2220189251Ssam 2221189251Ssam return MP_OKAY; 2222189251Ssam} 2223189251Ssam 2224189251Ssam 2225189251Ssam/* pre-calculate the value required for Barrett reduction 2226189251Ssam * For a given modulus "b" it calulates the value required in "a" 2227189251Ssam */ 2228189251Ssamstatic int mp_reduce_setup (mp_int * a, mp_int * b) 2229189251Ssam{ 2230189251Ssam int res; 2231189251Ssam 2232189251Ssam if ((res = mp_2expt (a, b->used * 2 * DIGIT_BIT)) != MP_OKAY) { 2233189251Ssam return res; 2234189251Ssam } 2235189251Ssam return mp_div (a, b, a, NULL); 2236189251Ssam} 2237189251Ssam 2238189251Ssam 2239189251Ssam/* reduces x mod m, assumes 0 < x < m**2, mu is 2240189251Ssam * precomputed via mp_reduce_setup. 2241189251Ssam * From HAC pp.604 Algorithm 14.42 2242189251Ssam */ 2243189251Ssamstatic int mp_reduce (mp_int * x, mp_int * m, mp_int * mu) 2244189251Ssam{ 2245189251Ssam mp_int q; 2246189251Ssam int res, um = m->used; 2247189251Ssam 2248189251Ssam /* q = x */ 2249189251Ssam if ((res = mp_init_copy (&q, x)) != MP_OKAY) { 2250189251Ssam return res; 2251189251Ssam } 2252189251Ssam 2253189251Ssam /* q1 = x / b**(k-1) */ 2254189251Ssam mp_rshd (&q, um - 1); 2255189251Ssam 2256189251Ssam /* according to HAC this optimization is ok */ 2257189251Ssam if (((unsigned long) um) > (((mp_digit)1) << (DIGIT_BIT - 1))) { 2258189251Ssam if ((res = mp_mul (&q, mu, &q)) != MP_OKAY) { 2259189251Ssam goto CLEANUP; 2260189251Ssam } 2261189251Ssam } else { 2262189251Ssam#ifdef BN_S_MP_MUL_HIGH_DIGS_C 2263189251Ssam if ((res = s_mp_mul_high_digs (&q, mu, &q, um)) != MP_OKAY) { 2264189251Ssam goto CLEANUP; 2265189251Ssam } 2266189251Ssam#elif defined(BN_FAST_S_MP_MUL_HIGH_DIGS_C) 2267189251Ssam if ((res = fast_s_mp_mul_high_digs (&q, mu, &q, um)) != MP_OKAY) { 2268189251Ssam goto CLEANUP; 2269189251Ssam } 2270189251Ssam#else 2271189251Ssam { 2272189251Ssam#error mp_reduce would always fail 2273189251Ssam res = MP_VAL; 2274189251Ssam goto CLEANUP; 2275189251Ssam } 2276189251Ssam#endif 2277189251Ssam } 2278189251Ssam 2279189251Ssam /* q3 = q2 / b**(k+1) */ 2280189251Ssam mp_rshd (&q, um + 1); 2281189251Ssam 2282189251Ssam /* x = x mod b**(k+1), quick (no division) */ 2283189251Ssam if ((res = mp_mod_2d (x, DIGIT_BIT * (um + 1), x)) != MP_OKAY) { 2284189251Ssam goto CLEANUP; 2285189251Ssam } 2286189251Ssam 2287189251Ssam /* q = q * m mod b**(k+1), quick (no division) */ 2288189251Ssam if ((res = s_mp_mul_digs (&q, m, &q, um + 1)) != MP_OKAY) { 2289189251Ssam goto CLEANUP; 2290189251Ssam } 2291189251Ssam 2292189251Ssam /* x = x - q */ 2293189251Ssam if ((res = mp_sub (x, &q, x)) != MP_OKAY) { 2294189251Ssam goto CLEANUP; 2295189251Ssam } 2296189251Ssam 2297189251Ssam /* If x < 0, add b**(k+1) to it */ 2298189251Ssam if (mp_cmp_d (x, 0) == MP_LT) { 2299189251Ssam mp_set (&q, 1); 2300189251Ssam if ((res = mp_lshd (&q, um + 1)) != MP_OKAY) { 2301189251Ssam goto CLEANUP; 2302189251Ssam } 2303189251Ssam if ((res = mp_add (x, &q, x)) != MP_OKAY) { 2304189251Ssam goto CLEANUP; 2305189251Ssam } 2306189251Ssam } 2307189251Ssam 2308189251Ssam /* Back off if it's too big */ 2309189251Ssam while (mp_cmp (x, m) != MP_LT) { 2310189251Ssam if ((res = s_mp_sub (x, m, x)) != MP_OKAY) { 2311189251Ssam goto CLEANUP; 2312189251Ssam } 2313189251Ssam } 2314189251Ssam 2315189251SsamCLEANUP: 2316189251Ssam mp_clear (&q); 2317189251Ssam 2318189251Ssam return res; 2319189251Ssam} 2320189251Ssam 2321189251Ssam 2322189251Ssam/* multiplies |a| * |b| and only computes upto digs digits of result 2323189251Ssam * HAC pp. 595, Algorithm 14.12 Modified so you can control how 2324189251Ssam * many digits of output are created. 2325189251Ssam */ 2326189251Ssamstatic int s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs) 2327189251Ssam{ 2328189251Ssam mp_int t; 2329189251Ssam int res, pa, pb, ix, iy; 2330189251Ssam mp_digit u; 2331189251Ssam mp_word r; 2332189251Ssam mp_digit tmpx, *tmpt, *tmpy; 2333189251Ssam 2334189251Ssam /* can we use the fast multiplier? */ 2335189251Ssam if (((digs) < MP_WARRAY) && 2336189251Ssam MIN (a->used, b->used) < 2337189251Ssam (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) { 2338189251Ssam return fast_s_mp_mul_digs (a, b, c, digs); 2339189251Ssam } 2340189251Ssam 2341189251Ssam if ((res = mp_init_size (&t, digs)) != MP_OKAY) { 2342189251Ssam return res; 2343189251Ssam } 2344189251Ssam t.used = digs; 2345189251Ssam 2346189251Ssam /* compute the digits of the product directly */ 2347189251Ssam pa = a->used; 2348189251Ssam for (ix = 0; ix < pa; ix++) { 2349189251Ssam /* set the carry to zero */ 2350189251Ssam u = 0; 2351189251Ssam 2352189251Ssam /* limit ourselves to making digs digits of output */ 2353189251Ssam pb = MIN (b->used, digs - ix); 2354189251Ssam 2355189251Ssam /* setup some aliases */ 2356189251Ssam /* copy of the digit from a used within the nested loop */ 2357189251Ssam tmpx = a->dp[ix]; 2358189251Ssam 2359189251Ssam /* an alias for the destination shifted ix places */ 2360189251Ssam tmpt = t.dp + ix; 2361189251Ssam 2362189251Ssam /* an alias for the digits of b */ 2363189251Ssam tmpy = b->dp; 2364189251Ssam 2365189251Ssam /* compute the columns of the output and propagate the carry */ 2366189251Ssam for (iy = 0; iy < pb; iy++) { 2367189251Ssam /* compute the column as a mp_word */ 2368189251Ssam r = ((mp_word)*tmpt) + 2369189251Ssam ((mp_word)tmpx) * ((mp_word)*tmpy++) + 2370189251Ssam ((mp_word) u); 2371189251Ssam 2372189251Ssam /* the new column is the lower part of the result */ 2373189251Ssam *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK)); 2374189251Ssam 2375189251Ssam /* get the carry word from the result */ 2376189251Ssam u = (mp_digit) (r >> ((mp_word) DIGIT_BIT)); 2377189251Ssam } 2378189251Ssam /* set carry if it is placed below digs */ 2379189251Ssam if (ix + iy < digs) { 2380189251Ssam *tmpt = u; 2381189251Ssam } 2382189251Ssam } 2383189251Ssam 2384189251Ssam mp_clamp (&t); 2385189251Ssam mp_exch (&t, c); 2386189251Ssam 2387189251Ssam mp_clear (&t); 2388189251Ssam return MP_OKAY; 2389189251Ssam} 2390189251Ssam 2391189251Ssam 2392189251Ssam/* Fast (comba) multiplier 2393189251Ssam * 2394189251Ssam * This is the fast column-array [comba] multiplier. It is 2395189251Ssam * designed to compute the columns of the product first 2396189251Ssam * then handle the carries afterwards. This has the effect 2397189251Ssam * of making the nested loops that compute the columns very 2398189251Ssam * simple and schedulable on super-scalar processors. 2399189251Ssam * 2400189251Ssam * This has been modified to produce a variable number of 2401189251Ssam * digits of output so if say only a half-product is required 2402189251Ssam * you don't have to compute the upper half (a feature 2403189251Ssam * required for fast Barrett reduction). 2404189251Ssam * 2405189251Ssam * Based on Algorithm 14.12 on pp.595 of HAC. 2406189251Ssam * 2407189251Ssam */ 2408189251Ssamstatic int fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs) 2409189251Ssam{ 2410189251Ssam int olduse, res, pa, ix, iz; 2411189251Ssam mp_digit W[MP_WARRAY]; 2412189251Ssam register mp_word _W; 2413189251Ssam 2414189251Ssam /* grow the destination as required */ 2415189251Ssam if (c->alloc < digs) { 2416189251Ssam if ((res = mp_grow (c, digs)) != MP_OKAY) { 2417189251Ssam return res; 2418189251Ssam } 2419189251Ssam } 2420189251Ssam 2421189251Ssam /* number of output digits to produce */ 2422189251Ssam pa = MIN(digs, a->used + b->used); 2423189251Ssam 2424189251Ssam /* clear the carry */ 2425189251Ssam _W = 0; 2426189251Ssam for (ix = 0; ix < pa; ix++) { 2427189251Ssam int tx, ty; 2428189251Ssam int iy; 2429189251Ssam mp_digit *tmpx, *tmpy; 2430189251Ssam 2431189251Ssam /* get offsets into the two bignums */ 2432189251Ssam ty = MIN(b->used-1, ix); 2433189251Ssam tx = ix - ty; 2434189251Ssam 2435189251Ssam /* setup temp aliases */ 2436189251Ssam tmpx = a->dp + tx; 2437189251Ssam tmpy = b->dp + ty; 2438189251Ssam 2439189251Ssam /* this is the number of times the loop will iterrate, essentially 2440189251Ssam while (tx++ < a->used && ty-- >= 0) { ... } 2441189251Ssam */ 2442189251Ssam iy = MIN(a->used-tx, ty+1); 2443189251Ssam 2444189251Ssam /* execute loop */ 2445189251Ssam for (iz = 0; iz < iy; ++iz) { 2446189251Ssam _W += ((mp_word)*tmpx++)*((mp_word)*tmpy--); 2447189251Ssam 2448189251Ssam } 2449189251Ssam 2450189251Ssam /* store term */ 2451189251Ssam W[ix] = ((mp_digit)_W) & MP_MASK; 2452189251Ssam 2453189251Ssam /* make next carry */ 2454189251Ssam _W = _W >> ((mp_word)DIGIT_BIT); 2455189251Ssam } 2456189251Ssam 2457189251Ssam /* setup dest */ 2458189251Ssam olduse = c->used; 2459189251Ssam c->used = pa; 2460189251Ssam 2461189251Ssam { 2462189251Ssam register mp_digit *tmpc; 2463189251Ssam tmpc = c->dp; 2464189251Ssam for (ix = 0; ix < pa+1; ix++) { 2465189251Ssam /* now extract the previous digit [below the carry] */ 2466189251Ssam *tmpc++ = W[ix]; 2467189251Ssam } 2468189251Ssam 2469189251Ssam /* clear unused digits [that existed in the old copy of c] */ 2470189251Ssam for (; ix < olduse; ix++) { 2471189251Ssam *tmpc++ = 0; 2472189251Ssam } 2473189251Ssam } 2474189251Ssam mp_clamp (c); 2475189251Ssam return MP_OKAY; 2476189251Ssam} 2477189251Ssam 2478189251Ssam 2479189251Ssam/* init an mp_init for a given size */ 2480189251Ssamstatic int mp_init_size (mp_int * a, int size) 2481189251Ssam{ 2482189251Ssam int x; 2483189251Ssam 2484189251Ssam /* pad size so there are always extra digits */ 2485189251Ssam size += (MP_PREC * 2) - (size % MP_PREC); 2486189251Ssam 2487189251Ssam /* alloc mem */ 2488189251Ssam a->dp = OPT_CAST(mp_digit) XMALLOC (sizeof (mp_digit) * size); 2489189251Ssam if (a->dp == NULL) { 2490189251Ssam return MP_MEM; 2491189251Ssam } 2492189251Ssam 2493189251Ssam /* set the members */ 2494189251Ssam a->used = 0; 2495189251Ssam a->alloc = size; 2496189251Ssam a->sign = MP_ZPOS; 2497189251Ssam 2498189251Ssam /* zero the digits */ 2499189251Ssam for (x = 0; x < size; x++) { 2500189251Ssam a->dp[x] = 0; 2501189251Ssam } 2502189251Ssam 2503189251Ssam return MP_OKAY; 2504189251Ssam} 2505189251Ssam 2506189251Ssam 2507189251Ssam/* low level squaring, b = a*a, HAC pp.596-597, Algorithm 14.16 */ 2508189251Ssamstatic int s_mp_sqr (mp_int * a, mp_int * b) 2509189251Ssam{ 2510189251Ssam mp_int t; 2511189251Ssam int res, ix, iy, pa; 2512189251Ssam mp_word r; 2513189251Ssam mp_digit u, tmpx, *tmpt; 2514189251Ssam 2515189251Ssam pa = a->used; 2516189251Ssam if ((res = mp_init_size (&t, 2*pa + 1)) != MP_OKAY) { 2517189251Ssam return res; 2518189251Ssam } 2519189251Ssam 2520189251Ssam /* default used is maximum possible size */ 2521189251Ssam t.used = 2*pa + 1; 2522189251Ssam 2523189251Ssam for (ix = 0; ix < pa; ix++) { 2524189251Ssam /* first calculate the digit at 2*ix */ 2525189251Ssam /* calculate double precision result */ 2526189251Ssam r = ((mp_word) t.dp[2*ix]) + 2527189251Ssam ((mp_word)a->dp[ix])*((mp_word)a->dp[ix]); 2528189251Ssam 2529189251Ssam /* store lower part in result */ 2530189251Ssam t.dp[ix+ix] = (mp_digit) (r & ((mp_word) MP_MASK)); 2531189251Ssam 2532189251Ssam /* get the carry */ 2533189251Ssam u = (mp_digit)(r >> ((mp_word) DIGIT_BIT)); 2534189251Ssam 2535189251Ssam /* left hand side of A[ix] * A[iy] */ 2536189251Ssam tmpx = a->dp[ix]; 2537189251Ssam 2538189251Ssam /* alias for where to store the results */ 2539189251Ssam tmpt = t.dp + (2*ix + 1); 2540189251Ssam 2541189251Ssam for (iy = ix + 1; iy < pa; iy++) { 2542189251Ssam /* first calculate the product */ 2543189251Ssam r = ((mp_word)tmpx) * ((mp_word)a->dp[iy]); 2544189251Ssam 2545189251Ssam /* now calculate the double precision result, note we use 2546189251Ssam * addition instead of *2 since it's easier to optimize 2547189251Ssam */ 2548189251Ssam r = ((mp_word) *tmpt) + r + r + ((mp_word) u); 2549189251Ssam 2550189251Ssam /* store lower part */ 2551189251Ssam *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK)); 2552189251Ssam 2553189251Ssam /* get carry */ 2554189251Ssam u = (mp_digit)(r >> ((mp_word) DIGIT_BIT)); 2555189251Ssam } 2556189251Ssam /* propagate upwards */ 2557189251Ssam while (u != ((mp_digit) 0)) { 2558189251Ssam r = ((mp_word) *tmpt) + ((mp_word) u); 2559189251Ssam *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK)); 2560189251Ssam u = (mp_digit)(r >> ((mp_word) DIGIT_BIT)); 2561189251Ssam } 2562189251Ssam } 2563189251Ssam 2564189251Ssam mp_clamp (&t); 2565189251Ssam mp_exch (&t, b); 2566189251Ssam mp_clear (&t); 2567189251Ssam return MP_OKAY; 2568189251Ssam} 2569189251Ssam 2570189251Ssam 2571189251Ssam/* multiplies |a| * |b| and does not compute the lower digs digits 2572189251Ssam * [meant to get the higher part of the product] 2573189251Ssam */ 2574189251Ssamstatic int s_mp_mul_high_digs (mp_int * a, mp_int * b, mp_int * c, int digs) 2575189251Ssam{ 2576189251Ssam mp_int t; 2577189251Ssam int res, pa, pb, ix, iy; 2578189251Ssam mp_digit u; 2579189251Ssam mp_word r; 2580189251Ssam mp_digit tmpx, *tmpt, *tmpy; 2581189251Ssam 2582189251Ssam /* can we use the fast multiplier? */ 2583189251Ssam#ifdef BN_FAST_S_MP_MUL_HIGH_DIGS_C 2584189251Ssam if (((a->used + b->used + 1) < MP_WARRAY) 2585189251Ssam && MIN (a->used, b->used) < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) { 2586189251Ssam return fast_s_mp_mul_high_digs (a, b, c, digs); 2587189251Ssam } 2588189251Ssam#endif 2589189251Ssam 2590189251Ssam if ((res = mp_init_size (&t, a->used + b->used + 1)) != MP_OKAY) { 2591189251Ssam return res; 2592189251Ssam } 2593189251Ssam t.used = a->used + b->used + 1; 2594189251Ssam 2595189251Ssam pa = a->used; 2596189251Ssam pb = b->used; 2597189251Ssam for (ix = 0; ix < pa; ix++) { 2598189251Ssam /* clear the carry */ 2599189251Ssam u = 0; 2600189251Ssam 2601189251Ssam /* left hand side of A[ix] * B[iy] */ 2602189251Ssam tmpx = a->dp[ix]; 2603189251Ssam 2604189251Ssam /* alias to the address of where the digits will be stored */ 2605189251Ssam tmpt = &(t.dp[digs]); 2606189251Ssam 2607189251Ssam /* alias for where to read the right hand side from */ 2608189251Ssam tmpy = b->dp + (digs - ix); 2609189251Ssam 2610189251Ssam for (iy = digs - ix; iy < pb; iy++) { 2611189251Ssam /* calculate the double precision result */ 2612189251Ssam r = ((mp_word)*tmpt) + 2613189251Ssam ((mp_word)tmpx) * ((mp_word)*tmpy++) + 2614189251Ssam ((mp_word) u); 2615189251Ssam 2616189251Ssam /* get the lower part */ 2617189251Ssam *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK)); 2618189251Ssam 2619189251Ssam /* carry the carry */ 2620189251Ssam u = (mp_digit) (r >> ((mp_word) DIGIT_BIT)); 2621189251Ssam } 2622189251Ssam *tmpt = u; 2623189251Ssam } 2624189251Ssam mp_clamp (&t); 2625189251Ssam mp_exch (&t, c); 2626189251Ssam mp_clear (&t); 2627189251Ssam return MP_OKAY; 2628189251Ssam} 2629189251Ssam 2630189251Ssam 2631189251Ssam#ifdef BN_MP_MONTGOMERY_SETUP_C 2632189251Ssam/* setups the montgomery reduction stuff */ 2633189251Ssamstatic int 2634189251Ssammp_montgomery_setup (mp_int * n, mp_digit * rho) 2635189251Ssam{ 2636189251Ssam mp_digit x, b; 2637189251Ssam 2638189251Ssam/* fast inversion mod 2**k 2639189251Ssam * 2640189251Ssam * Based on the fact that 2641189251Ssam * 2642189251Ssam * XA = 1 (mod 2**n) => (X(2-XA)) A = 1 (mod 2**2n) 2643189251Ssam * => 2*X*A - X*X*A*A = 1 2644189251Ssam * => 2*(1) - (1) = 1 2645189251Ssam */ 2646189251Ssam b = n->dp[0]; 2647189251Ssam 2648189251Ssam if ((b & 1) == 0) { 2649189251Ssam return MP_VAL; 2650189251Ssam } 2651189251Ssam 2652189251Ssam x = (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2**4 */ 2653189251Ssam x *= 2 - b * x; /* here x*a==1 mod 2**8 */ 2654189251Ssam#if !defined(MP_8BIT) 2655189251Ssam x *= 2 - b * x; /* here x*a==1 mod 2**16 */ 2656189251Ssam#endif 2657189251Ssam#if defined(MP_64BIT) || !(defined(MP_8BIT) || defined(MP_16BIT)) 2658189251Ssam x *= 2 - b * x; /* here x*a==1 mod 2**32 */ 2659189251Ssam#endif 2660189251Ssam#ifdef MP_64BIT 2661189251Ssam x *= 2 - b * x; /* here x*a==1 mod 2**64 */ 2662189251Ssam#endif 2663189251Ssam 2664189251Ssam /* rho = -1/m mod b */ 2665189251Ssam *rho = (unsigned long)(((mp_word)1 << ((mp_word) DIGIT_BIT)) - x) & MP_MASK; 2666189251Ssam 2667189251Ssam return MP_OKAY; 2668189251Ssam} 2669189251Ssam#endif 2670189251Ssam 2671189251Ssam 2672189251Ssam#ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C 2673189251Ssam/* computes xR**-1 == x (mod N) via Montgomery Reduction 2674189251Ssam * 2675189251Ssam * This is an optimized implementation of montgomery_reduce 2676189251Ssam * which uses the comba method to quickly calculate the columns of the 2677189251Ssam * reduction. 2678189251Ssam * 2679189251Ssam * Based on Algorithm 14.32 on pp.601 of HAC. 2680189251Ssam*/ 2681189251Ssamint fast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho) 2682189251Ssam{ 2683189251Ssam int ix, res, olduse; 2684189251Ssam mp_word W[MP_WARRAY]; 2685189251Ssam 2686189251Ssam /* get old used count */ 2687189251Ssam olduse = x->used; 2688189251Ssam 2689189251Ssam /* grow a as required */ 2690189251Ssam if (x->alloc < n->used + 1) { 2691189251Ssam if ((res = mp_grow (x, n->used + 1)) != MP_OKAY) { 2692189251Ssam return res; 2693189251Ssam } 2694189251Ssam } 2695189251Ssam 2696189251Ssam /* first we have to get the digits of the input into 2697189251Ssam * an array of double precision words W[...] 2698189251Ssam */ 2699189251Ssam { 2700189251Ssam register mp_word *_W; 2701189251Ssam register mp_digit *tmpx; 2702189251Ssam 2703189251Ssam /* alias for the W[] array */ 2704189251Ssam _W = W; 2705189251Ssam 2706189251Ssam /* alias for the digits of x*/ 2707189251Ssam tmpx = x->dp; 2708189251Ssam 2709189251Ssam /* copy the digits of a into W[0..a->used-1] */ 2710189251Ssam for (ix = 0; ix < x->used; ix++) { 2711189251Ssam *_W++ = *tmpx++; 2712189251Ssam } 2713189251Ssam 2714189251Ssam /* zero the high words of W[a->used..m->used*2] */ 2715189251Ssam for (; ix < n->used * 2 + 1; ix++) { 2716189251Ssam *_W++ = 0; 2717189251Ssam } 2718189251Ssam } 2719189251Ssam 2720189251Ssam /* now we proceed to zero successive digits 2721189251Ssam * from the least significant upwards 2722189251Ssam */ 2723189251Ssam for (ix = 0; ix < n->used; ix++) { 2724189251Ssam /* mu = ai * m' mod b 2725189251Ssam * 2726189251Ssam * We avoid a double precision multiplication (which isn't required) 2727189251Ssam * by casting the value down to a mp_digit. Note this requires 2728189251Ssam * that W[ix-1] have the carry cleared (see after the inner loop) 2729189251Ssam */ 2730189251Ssam register mp_digit mu; 2731189251Ssam mu = (mp_digit) (((W[ix] & MP_MASK) * rho) & MP_MASK); 2732189251Ssam 2733189251Ssam /* a = a + mu * m * b**i 2734189251Ssam * 2735189251Ssam * This is computed in place and on the fly. The multiplication 2736189251Ssam * by b**i is handled by offseting which columns the results 2737189251Ssam * are added to. 2738189251Ssam * 2739189251Ssam * Note the comba method normally doesn't handle carries in the 2740189251Ssam * inner loop In this case we fix the carry from the previous 2741189251Ssam * column since the Montgomery reduction requires digits of the 2742189251Ssam * result (so far) [see above] to work. This is 2743189251Ssam * handled by fixing up one carry after the inner loop. The 2744189251Ssam * carry fixups are done in order so after these loops the 2745189251Ssam * first m->used words of W[] have the carries fixed 2746189251Ssam */ 2747189251Ssam { 2748189251Ssam register int iy; 2749189251Ssam register mp_digit *tmpn; 2750189251Ssam register mp_word *_W; 2751189251Ssam 2752189251Ssam /* alias for the digits of the modulus */ 2753189251Ssam tmpn = n->dp; 2754189251Ssam 2755189251Ssam /* Alias for the columns set by an offset of ix */ 2756189251Ssam _W = W + ix; 2757189251Ssam 2758189251Ssam /* inner loop */ 2759189251Ssam for (iy = 0; iy < n->used; iy++) { 2760189251Ssam *_W++ += ((mp_word)mu) * ((mp_word)*tmpn++); 2761189251Ssam } 2762189251Ssam } 2763189251Ssam 2764189251Ssam /* now fix carry for next digit, W[ix+1] */ 2765189251Ssam W[ix + 1] += W[ix] >> ((mp_word) DIGIT_BIT); 2766189251Ssam } 2767189251Ssam 2768189251Ssam /* now we have to propagate the carries and 2769189251Ssam * shift the words downward [all those least 2770189251Ssam * significant digits we zeroed]. 2771189251Ssam */ 2772189251Ssam { 2773189251Ssam register mp_digit *tmpx; 2774189251Ssam register mp_word *_W, *_W1; 2775189251Ssam 2776189251Ssam /* nox fix rest of carries */ 2777189251Ssam 2778189251Ssam /* alias for current word */ 2779189251Ssam _W1 = W + ix; 2780189251Ssam 2781189251Ssam /* alias for next word, where the carry goes */ 2782189251Ssam _W = W + ++ix; 2783189251Ssam 2784189251Ssam for (; ix <= n->used * 2 + 1; ix++) { 2785189251Ssam *_W++ += *_W1++ >> ((mp_word) DIGIT_BIT); 2786189251Ssam } 2787189251Ssam 2788189251Ssam /* copy out, A = A/b**n 2789189251Ssam * 2790189251Ssam * The result is A/b**n but instead of converting from an 2791189251Ssam * array of mp_word to mp_digit than calling mp_rshd 2792189251Ssam * we just copy them in the right order 2793189251Ssam */ 2794189251Ssam 2795189251Ssam /* alias for destination word */ 2796189251Ssam tmpx = x->dp; 2797189251Ssam 2798189251Ssam /* alias for shifted double precision result */ 2799189251Ssam _W = W + n->used; 2800189251Ssam 2801189251Ssam for (ix = 0; ix < n->used + 1; ix++) { 2802189251Ssam *tmpx++ = (mp_digit)(*_W++ & ((mp_word) MP_MASK)); 2803189251Ssam } 2804189251Ssam 2805189251Ssam /* zero oldused digits, if the input a was larger than 2806189251Ssam * m->used+1 we'll have to clear the digits 2807189251Ssam */ 2808189251Ssam for (; ix < olduse; ix++) { 2809189251Ssam *tmpx++ = 0; 2810189251Ssam } 2811189251Ssam } 2812189251Ssam 2813189251Ssam /* set the max used and clamp */ 2814189251Ssam x->used = n->used + 1; 2815189251Ssam mp_clamp (x); 2816189251Ssam 2817189251Ssam /* if A >= m then A = A - m */ 2818189251Ssam if (mp_cmp_mag (x, n) != MP_LT) { 2819189251Ssam return s_mp_sub (x, n, x); 2820189251Ssam } 2821189251Ssam return MP_OKAY; 2822189251Ssam} 2823189251Ssam#endif 2824189251Ssam 2825189251Ssam 2826189251Ssam#ifdef BN_MP_MUL_2_C 2827189251Ssam/* b = a*2 */ 2828189251Ssamstatic int mp_mul_2(mp_int * a, mp_int * b) 2829189251Ssam{ 2830189251Ssam int x, res, oldused; 2831189251Ssam 2832189251Ssam /* grow to accomodate result */ 2833189251Ssam if (b->alloc < a->used + 1) { 2834189251Ssam if ((res = mp_grow (b, a->used + 1)) != MP_OKAY) { 2835189251Ssam return res; 2836189251Ssam } 2837189251Ssam } 2838189251Ssam 2839189251Ssam oldused = b->used; 2840189251Ssam b->used = a->used; 2841189251Ssam 2842189251Ssam { 2843189251Ssam register mp_digit r, rr, *tmpa, *tmpb; 2844189251Ssam 2845189251Ssam /* alias for source */ 2846189251Ssam tmpa = a->dp; 2847189251Ssam 2848189251Ssam /* alias for dest */ 2849189251Ssam tmpb = b->dp; 2850189251Ssam 2851189251Ssam /* carry */ 2852189251Ssam r = 0; 2853189251Ssam for (x = 0; x < a->used; x++) { 2854189251Ssam 2855189251Ssam /* get what will be the *next* carry bit from the 2856189251Ssam * MSB of the current digit 2857189251Ssam */ 2858189251Ssam rr = *tmpa >> ((mp_digit)(DIGIT_BIT - 1)); 2859189251Ssam 2860189251Ssam /* now shift up this digit, add in the carry [from the previous] */ 2861189251Ssam *tmpb++ = ((*tmpa++ << ((mp_digit)1)) | r) & MP_MASK; 2862189251Ssam 2863189251Ssam /* copy the carry that would be from the source 2864189251Ssam * digit into the next iteration 2865189251Ssam */ 2866189251Ssam r = rr; 2867189251Ssam } 2868189251Ssam 2869189251Ssam /* new leading digit? */ 2870189251Ssam if (r != 0) { 2871189251Ssam /* add a MSB which is always 1 at this point */ 2872189251Ssam *tmpb = 1; 2873189251Ssam ++(b->used); 2874189251Ssam } 2875189251Ssam 2876189251Ssam /* now zero any excess digits on the destination 2877189251Ssam * that we didn't write to 2878189251Ssam */ 2879189251Ssam tmpb = b->dp + b->used; 2880189251Ssam for (x = b->used; x < oldused; x++) { 2881189251Ssam *tmpb++ = 0; 2882189251Ssam } 2883189251Ssam } 2884189251Ssam b->sign = a->sign; 2885189251Ssam return MP_OKAY; 2886189251Ssam} 2887189251Ssam#endif 2888189251Ssam 2889189251Ssam 2890189251Ssam#ifdef BN_MP_MONTGOMERY_CALC_NORMALIZATION_C 2891189251Ssam/* 2892189251Ssam * shifts with subtractions when the result is greater than b. 2893189251Ssam * 2894189251Ssam * The method is slightly modified to shift B unconditionally upto just under 2895189251Ssam * the leading bit of b. This saves alot of multiple precision shifting. 2896189251Ssam */ 2897189251Ssamstatic int mp_montgomery_calc_normalization (mp_int * a, mp_int * b) 2898189251Ssam{ 2899189251Ssam int x, bits, res; 2900189251Ssam 2901189251Ssam /* how many bits of last digit does b use */ 2902189251Ssam bits = mp_count_bits (b) % DIGIT_BIT; 2903189251Ssam 2904189251Ssam if (b->used > 1) { 2905189251Ssam if ((res = mp_2expt (a, (b->used - 1) * DIGIT_BIT + bits - 1)) != MP_OKAY) { 2906189251Ssam return res; 2907189251Ssam } 2908189251Ssam } else { 2909189251Ssam mp_set(a, 1); 2910189251Ssam bits = 1; 2911189251Ssam } 2912189251Ssam 2913189251Ssam 2914189251Ssam /* now compute C = A * B mod b */ 2915189251Ssam for (x = bits - 1; x < (int)DIGIT_BIT; x++) { 2916189251Ssam if ((res = mp_mul_2 (a, a)) != MP_OKAY) { 2917189251Ssam return res; 2918189251Ssam } 2919189251Ssam if (mp_cmp_mag (a, b) != MP_LT) { 2920189251Ssam if ((res = s_mp_sub (a, b, a)) != MP_OKAY) { 2921189251Ssam return res; 2922189251Ssam } 2923189251Ssam } 2924189251Ssam } 2925189251Ssam 2926189251Ssam return MP_OKAY; 2927189251Ssam} 2928189251Ssam#endif 2929189251Ssam 2930189251Ssam 2931189251Ssam#ifdef BN_MP_EXPTMOD_FAST_C 2932189251Ssam/* computes Y == G**X mod P, HAC pp.616, Algorithm 14.85 2933189251Ssam * 2934189251Ssam * Uses a left-to-right k-ary sliding window to compute the modular exponentiation. 2935189251Ssam * The value of k changes based on the size of the exponent. 2936189251Ssam * 2937189251Ssam * Uses Montgomery or Diminished Radix reduction [whichever appropriate] 2938189251Ssam */ 2939189251Ssam 2940189251Ssamstatic int mp_exptmod_fast (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode) 2941189251Ssam{ 2942189251Ssam mp_int M[TAB_SIZE], res; 2943189251Ssam mp_digit buf, mp; 2944189251Ssam int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize; 2945189251Ssam 2946189251Ssam /* use a pointer to the reduction algorithm. This allows us to use 2947189251Ssam * one of many reduction algorithms without modding the guts of 2948189251Ssam * the code with if statements everywhere. 2949189251Ssam */ 2950189251Ssam int (*redux)(mp_int*,mp_int*,mp_digit); 2951189251Ssam 2952189251Ssam /* find window size */ 2953189251Ssam x = mp_count_bits (X); 2954189251Ssam if (x <= 7) { 2955189251Ssam winsize = 2; 2956189251Ssam } else if (x <= 36) { 2957189251Ssam winsize = 3; 2958189251Ssam } else if (x <= 140) { 2959189251Ssam winsize = 4; 2960189251Ssam } else if (x <= 450) { 2961189251Ssam winsize = 5; 2962189251Ssam } else if (x <= 1303) { 2963189251Ssam winsize = 6; 2964189251Ssam } else if (x <= 3529) { 2965189251Ssam winsize = 7; 2966189251Ssam } else { 2967189251Ssam winsize = 8; 2968189251Ssam } 2969189251Ssam 2970189251Ssam#ifdef MP_LOW_MEM 2971189251Ssam if (winsize > 5) { 2972189251Ssam winsize = 5; 2973189251Ssam } 2974189251Ssam#endif 2975189251Ssam 2976189251Ssam /* init M array */ 2977189251Ssam /* init first cell */ 2978189251Ssam if ((err = mp_init(&M[1])) != MP_OKAY) { 2979189251Ssam return err; 2980189251Ssam } 2981189251Ssam 2982189251Ssam /* now init the second half of the array */ 2983189251Ssam for (x = 1<<(winsize-1); x < (1 << winsize); x++) { 2984189251Ssam if ((err = mp_init(&M[x])) != MP_OKAY) { 2985189251Ssam for (y = 1<<(winsize-1); y < x; y++) { 2986189251Ssam mp_clear (&M[y]); 2987189251Ssam } 2988189251Ssam mp_clear(&M[1]); 2989189251Ssam return err; 2990189251Ssam } 2991189251Ssam } 2992189251Ssam 2993189251Ssam /* determine and setup reduction code */ 2994189251Ssam if (redmode == 0) { 2995189251Ssam#ifdef BN_MP_MONTGOMERY_SETUP_C 2996189251Ssam /* now setup montgomery */ 2997189251Ssam if ((err = mp_montgomery_setup (P, &mp)) != MP_OKAY) { 2998189251Ssam goto LBL_M; 2999189251Ssam } 3000189251Ssam#else 3001189251Ssam err = MP_VAL; 3002189251Ssam goto LBL_M; 3003189251Ssam#endif 3004189251Ssam 3005189251Ssam /* automatically pick the comba one if available (saves quite a few calls/ifs) */ 3006189251Ssam#ifdef BN_FAST_MP_MONTGOMERY_REDUCE_C 3007189251Ssam if (((P->used * 2 + 1) < MP_WARRAY) && 3008189251Ssam P->used < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) { 3009189251Ssam redux = fast_mp_montgomery_reduce; 3010189251Ssam } else 3011189251Ssam#endif 3012189251Ssam { 3013189251Ssam#ifdef BN_MP_MONTGOMERY_REDUCE_C 3014189251Ssam /* use slower baseline Montgomery method */ 3015189251Ssam redux = mp_montgomery_reduce; 3016189251Ssam#else 3017189251Ssam err = MP_VAL; 3018189251Ssam goto LBL_M; 3019189251Ssam#endif 3020189251Ssam } 3021189251Ssam } else if (redmode == 1) { 3022189251Ssam#if defined(BN_MP_DR_SETUP_C) && defined(BN_MP_DR_REDUCE_C) 3023189251Ssam /* setup DR reduction for moduli of the form B**k - b */ 3024189251Ssam mp_dr_setup(P, &mp); 3025189251Ssam redux = mp_dr_reduce; 3026189251Ssam#else 3027189251Ssam err = MP_VAL; 3028189251Ssam goto LBL_M; 3029189251Ssam#endif 3030189251Ssam } else { 3031189251Ssam#if defined(BN_MP_REDUCE_2K_SETUP_C) && defined(BN_MP_REDUCE_2K_C) 3032189251Ssam /* setup DR reduction for moduli of the form 2**k - b */ 3033189251Ssam if ((err = mp_reduce_2k_setup(P, &mp)) != MP_OKAY) { 3034189251Ssam goto LBL_M; 3035189251Ssam } 3036189251Ssam redux = mp_reduce_2k; 3037189251Ssam#else 3038189251Ssam err = MP_VAL; 3039189251Ssam goto LBL_M; 3040189251Ssam#endif 3041189251Ssam } 3042189251Ssam 3043189251Ssam /* setup result */ 3044189251Ssam if ((err = mp_init (&res)) != MP_OKAY) { 3045189251Ssam goto LBL_M; 3046189251Ssam } 3047189251Ssam 3048189251Ssam /* create M table 3049189251Ssam * 3050189251Ssam 3051189251Ssam * 3052189251Ssam * The first half of the table is not computed though accept for M[0] and M[1] 3053189251Ssam */ 3054189251Ssam 3055189251Ssam if (redmode == 0) { 3056189251Ssam#ifdef BN_MP_MONTGOMERY_CALC_NORMALIZATION_C 3057189251Ssam /* now we need R mod m */ 3058189251Ssam if ((err = mp_montgomery_calc_normalization (&res, P)) != MP_OKAY) { 3059189251Ssam goto LBL_RES; 3060189251Ssam } 3061189251Ssam#else 3062189251Ssam err = MP_VAL; 3063189251Ssam goto LBL_RES; 3064189251Ssam#endif 3065189251Ssam 3066189251Ssam /* now set M[1] to G * R mod m */ 3067189251Ssam if ((err = mp_mulmod (G, &res, P, &M[1])) != MP_OKAY) { 3068189251Ssam goto LBL_RES; 3069189251Ssam } 3070189251Ssam } else { 3071189251Ssam mp_set(&res, 1); 3072189251Ssam if ((err = mp_mod(G, P, &M[1])) != MP_OKAY) { 3073189251Ssam goto LBL_RES; 3074189251Ssam } 3075189251Ssam } 3076189251Ssam 3077189251Ssam /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */ 3078189251Ssam if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) { 3079189251Ssam goto LBL_RES; 3080189251Ssam } 3081189251Ssam 3082189251Ssam for (x = 0; x < (winsize - 1); x++) { 3083189251Ssam if ((err = mp_sqr (&M[1 << (winsize - 1)], &M[1 << (winsize - 1)])) != MP_OKAY) { 3084189251Ssam goto LBL_RES; 3085189251Ssam } 3086189251Ssam if ((err = redux (&M[1 << (winsize - 1)], P, mp)) != MP_OKAY) { 3087189251Ssam goto LBL_RES; 3088189251Ssam } 3089189251Ssam } 3090189251Ssam 3091189251Ssam /* create upper table */ 3092189251Ssam for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) { 3093189251Ssam if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) { 3094189251Ssam goto LBL_RES; 3095189251Ssam } 3096189251Ssam if ((err = redux (&M[x], P, mp)) != MP_OKAY) { 3097189251Ssam goto LBL_RES; 3098189251Ssam } 3099189251Ssam } 3100189251Ssam 3101189251Ssam /* set initial mode and bit cnt */ 3102189251Ssam mode = 0; 3103189251Ssam bitcnt = 1; 3104189251Ssam buf = 0; 3105189251Ssam digidx = X->used - 1; 3106189251Ssam bitcpy = 0; 3107189251Ssam bitbuf = 0; 3108189251Ssam 3109189251Ssam for (;;) { 3110189251Ssam /* grab next digit as required */ 3111189251Ssam if (--bitcnt == 0) { 3112189251Ssam /* if digidx == -1 we are out of digits so break */ 3113189251Ssam if (digidx == -1) { 3114189251Ssam break; 3115189251Ssam } 3116189251Ssam /* read next digit and reset bitcnt */ 3117189251Ssam buf = X->dp[digidx--]; 3118189251Ssam bitcnt = (int)DIGIT_BIT; 3119189251Ssam } 3120189251Ssam 3121189251Ssam /* grab the next msb from the exponent */ 3122189251Ssam y = (mp_digit)(buf >> (DIGIT_BIT - 1)) & 1; 3123189251Ssam buf <<= (mp_digit)1; 3124189251Ssam 3125189251Ssam /* if the bit is zero and mode == 0 then we ignore it 3126189251Ssam * These represent the leading zero bits before the first 1 bit 3127189251Ssam * in the exponent. Technically this opt is not required but it 3128189251Ssam * does lower the # of trivial squaring/reductions used 3129189251Ssam */ 3130189251Ssam if (mode == 0 && y == 0) { 3131189251Ssam continue; 3132189251Ssam } 3133189251Ssam 3134189251Ssam /* if the bit is zero and mode == 1 then we square */ 3135189251Ssam if (mode == 1 && y == 0) { 3136189251Ssam if ((err = mp_sqr (&res, &res)) != MP_OKAY) { 3137189251Ssam goto LBL_RES; 3138189251Ssam } 3139189251Ssam if ((err = redux (&res, P, mp)) != MP_OKAY) { 3140189251Ssam goto LBL_RES; 3141189251Ssam } 3142189251Ssam continue; 3143189251Ssam } 3144189251Ssam 3145189251Ssam /* else we add it to the window */ 3146189251Ssam bitbuf |= (y << (winsize - ++bitcpy)); 3147189251Ssam mode = 2; 3148189251Ssam 3149189251Ssam if (bitcpy == winsize) { 3150189251Ssam /* ok window is filled so square as required and multiply */ 3151189251Ssam /* square first */ 3152189251Ssam for (x = 0; x < winsize; x++) { 3153189251Ssam if ((err = mp_sqr (&res, &res)) != MP_OKAY) { 3154189251Ssam goto LBL_RES; 3155189251Ssam } 3156189251Ssam if ((err = redux (&res, P, mp)) != MP_OKAY) { 3157189251Ssam goto LBL_RES; 3158189251Ssam } 3159189251Ssam } 3160189251Ssam 3161189251Ssam /* then multiply */ 3162189251Ssam if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) { 3163189251Ssam goto LBL_RES; 3164189251Ssam } 3165189251Ssam if ((err = redux (&res, P, mp)) != MP_OKAY) { 3166189251Ssam goto LBL_RES; 3167189251Ssam } 3168189251Ssam 3169189251Ssam /* empty window and reset */ 3170189251Ssam bitcpy = 0; 3171189251Ssam bitbuf = 0; 3172189251Ssam mode = 1; 3173189251Ssam } 3174189251Ssam } 3175189251Ssam 3176189251Ssam /* if bits remain then square/multiply */ 3177189251Ssam if (mode == 2 && bitcpy > 0) { 3178189251Ssam /* square then multiply if the bit is set */ 3179189251Ssam for (x = 0; x < bitcpy; x++) { 3180189251Ssam if ((err = mp_sqr (&res, &res)) != MP_OKAY) { 3181189251Ssam goto LBL_RES; 3182189251Ssam } 3183189251Ssam if ((err = redux (&res, P, mp)) != MP_OKAY) { 3184189251Ssam goto LBL_RES; 3185189251Ssam } 3186189251Ssam 3187189251Ssam /* get next bit of the window */ 3188189251Ssam bitbuf <<= 1; 3189189251Ssam if ((bitbuf & (1 << winsize)) != 0) { 3190189251Ssam /* then multiply */ 3191189251Ssam if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) { 3192189251Ssam goto LBL_RES; 3193189251Ssam } 3194189251Ssam if ((err = redux (&res, P, mp)) != MP_OKAY) { 3195189251Ssam goto LBL_RES; 3196189251Ssam } 3197189251Ssam } 3198189251Ssam } 3199189251Ssam } 3200189251Ssam 3201189251Ssam if (redmode == 0) { 3202189251Ssam /* fixup result if Montgomery reduction is used 3203189251Ssam * recall that any value in a Montgomery system is 3204189251Ssam * actually multiplied by R mod n. So we have 3205189251Ssam * to reduce one more time to cancel out the factor 3206189251Ssam * of R. 3207189251Ssam */ 3208189251Ssam if ((err = redux(&res, P, mp)) != MP_OKAY) { 3209189251Ssam goto LBL_RES; 3210189251Ssam } 3211189251Ssam } 3212189251Ssam 3213189251Ssam /* swap res with Y */ 3214189251Ssam mp_exch (&res, Y); 3215189251Ssam err = MP_OKAY; 3216189251SsamLBL_RES:mp_clear (&res); 3217189251SsamLBL_M: 3218189251Ssam mp_clear(&M[1]); 3219189251Ssam for (x = 1<<(winsize-1); x < (1 << winsize); x++) { 3220189251Ssam mp_clear (&M[x]); 3221189251Ssam } 3222189251Ssam return err; 3223189251Ssam} 3224189251Ssam#endif 3225189251Ssam 3226189251Ssam 3227189251Ssam#ifdef BN_FAST_S_MP_SQR_C 3228189251Ssam/* the jist of squaring... 3229189251Ssam * you do like mult except the offset of the tmpx [one that 3230189251Ssam * starts closer to zero] can't equal the offset of tmpy. 3231189251Ssam * So basically you set up iy like before then you min it with 3232189251Ssam * (ty-tx) so that it never happens. You double all those 3233189251Ssam * you add in the inner loop 3234189251Ssam 3235189251SsamAfter that loop you do the squares and add them in. 3236189251Ssam*/ 3237189251Ssam 3238189251Ssamstatic int fast_s_mp_sqr (mp_int * a, mp_int * b) 3239189251Ssam{ 3240189251Ssam int olduse, res, pa, ix, iz; 3241189251Ssam mp_digit W[MP_WARRAY], *tmpx; 3242189251Ssam mp_word W1; 3243189251Ssam 3244189251Ssam /* grow the destination as required */ 3245189251Ssam pa = a->used + a->used; 3246189251Ssam if (b->alloc < pa) { 3247189251Ssam if ((res = mp_grow (b, pa)) != MP_OKAY) { 3248189251Ssam return res; 3249189251Ssam } 3250189251Ssam } 3251189251Ssam 3252189251Ssam /* number of output digits to produce */ 3253189251Ssam W1 = 0; 3254189251Ssam for (ix = 0; ix < pa; ix++) { 3255189251Ssam int tx, ty, iy; 3256189251Ssam mp_word _W; 3257189251Ssam mp_digit *tmpy; 3258189251Ssam 3259189251Ssam /* clear counter */ 3260189251Ssam _W = 0; 3261189251Ssam 3262189251Ssam /* get offsets into the two bignums */ 3263189251Ssam ty = MIN(a->used-1, ix); 3264189251Ssam tx = ix - ty; 3265189251Ssam 3266189251Ssam /* setup temp aliases */ 3267189251Ssam tmpx = a->dp + tx; 3268189251Ssam tmpy = a->dp + ty; 3269189251Ssam 3270189251Ssam /* this is the number of times the loop will iterrate, essentially 3271189251Ssam while (tx++ < a->used && ty-- >= 0) { ... } 3272189251Ssam */ 3273189251Ssam iy = MIN(a->used-tx, ty+1); 3274189251Ssam 3275189251Ssam /* now for squaring tx can never equal ty 3276189251Ssam * we halve the distance since they approach at a rate of 2x 3277189251Ssam * and we have to round because odd cases need to be executed 3278189251Ssam */ 3279189251Ssam iy = MIN(iy, (ty-tx+1)>>1); 3280189251Ssam 3281189251Ssam /* execute loop */ 3282189251Ssam for (iz = 0; iz < iy; iz++) { 3283189251Ssam _W += ((mp_word)*tmpx++)*((mp_word)*tmpy--); 3284189251Ssam } 3285189251Ssam 3286189251Ssam /* double the inner product and add carry */ 3287189251Ssam _W = _W + _W + W1; 3288189251Ssam 3289189251Ssam /* even columns have the square term in them */ 3290189251Ssam if ((ix&1) == 0) { 3291189251Ssam _W += ((mp_word)a->dp[ix>>1])*((mp_word)a->dp[ix>>1]); 3292189251Ssam } 3293189251Ssam 3294189251Ssam /* store it */ 3295189251Ssam W[ix] = (mp_digit)(_W & MP_MASK); 3296189251Ssam 3297189251Ssam /* make next carry */ 3298189251Ssam W1 = _W >> ((mp_word)DIGIT_BIT); 3299189251Ssam } 3300189251Ssam 3301189251Ssam /* setup dest */ 3302189251Ssam olduse = b->used; 3303189251Ssam b->used = a->used+a->used; 3304189251Ssam 3305189251Ssam { 3306189251Ssam mp_digit *tmpb; 3307189251Ssam tmpb = b->dp; 3308189251Ssam for (ix = 0; ix < pa; ix++) { 3309189251Ssam *tmpb++ = W[ix] & MP_MASK; 3310189251Ssam } 3311189251Ssam 3312189251Ssam /* clear unused digits [that existed in the old copy of c] */ 3313189251Ssam for (; ix < olduse; ix++) { 3314189251Ssam *tmpb++ = 0; 3315189251Ssam } 3316189251Ssam } 3317189251Ssam mp_clamp (b); 3318189251Ssam return MP_OKAY; 3319189251Ssam} 3320189251Ssam#endif 3321189251Ssam 3322189251Ssam 3323189251Ssam#ifdef BN_MP_MUL_D_C 3324189251Ssam/* multiply by a digit */ 3325189251Ssamstatic int 3326189251Ssammp_mul_d (mp_int * a, mp_digit b, mp_int * c) 3327189251Ssam{ 3328189251Ssam mp_digit u, *tmpa, *tmpc; 3329189251Ssam mp_word r; 3330189251Ssam int ix, res, olduse; 3331189251Ssam 3332189251Ssam /* make sure c is big enough to hold a*b */ 3333189251Ssam if (c->alloc < a->used + 1) { 3334189251Ssam if ((res = mp_grow (c, a->used + 1)) != MP_OKAY) { 3335189251Ssam return res; 3336189251Ssam } 3337189251Ssam } 3338189251Ssam 3339189251Ssam /* get the original destinations used count */ 3340189251Ssam olduse = c->used; 3341189251Ssam 3342189251Ssam /* set the sign */ 3343189251Ssam c->sign = a->sign; 3344189251Ssam 3345189251Ssam /* alias for a->dp [source] */ 3346189251Ssam tmpa = a->dp; 3347189251Ssam 3348189251Ssam /* alias for c->dp [dest] */ 3349189251Ssam tmpc = c->dp; 3350189251Ssam 3351189251Ssam /* zero carry */ 3352189251Ssam u = 0; 3353189251Ssam 3354189251Ssam /* compute columns */ 3355189251Ssam for (ix = 0; ix < a->used; ix++) { 3356189251Ssam /* compute product and carry sum for this term */ 3357189251Ssam r = ((mp_word) u) + ((mp_word)*tmpa++) * ((mp_word)b); 3358189251Ssam 3359189251Ssam /* mask off higher bits to get a single digit */ 3360189251Ssam *tmpc++ = (mp_digit) (r & ((mp_word) MP_MASK)); 3361189251Ssam 3362189251Ssam /* send carry into next iteration */ 3363189251Ssam u = (mp_digit) (r >> ((mp_word) DIGIT_BIT)); 3364189251Ssam } 3365189251Ssam 3366189251Ssam /* store final carry [if any] and increment ix offset */ 3367189251Ssam *tmpc++ = u; 3368189251Ssam ++ix; 3369189251Ssam 3370189251Ssam /* now zero digits above the top */ 3371189251Ssam while (ix++ < olduse) { 3372189251Ssam *tmpc++ = 0; 3373189251Ssam } 3374189251Ssam 3375189251Ssam /* set used count */ 3376189251Ssam c->used = a->used + 1; 3377189251Ssam mp_clamp(c); 3378189251Ssam 3379189251Ssam return MP_OKAY; 3380189251Ssam} 3381189251Ssam#endif 3382