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