mpi.c

支持SSL v2/v3, TLS, PKCS #5, PKCS #7, PKCS #11, PKCS #12, S/MIME, X.509v3证书等安全协议或标准的开发库编译用到NSPR

      MP_CHECKOK( mp_sub(&v, &u, &v) );
      MP_CHECKOK( mp_sub(&C, &A, &C) );
      MP_CHECKOK( mp_sub(&D, &B, &D) );
    }
  } while (mp_cmp_z(&u) != 0);

  /* copy results to output */
  if(x)
    MP_CHECKOK( mp_copy(&C, x) );

  if(y)
    MP_CHECKOK( mp_copy(&D, y) );
      
  if(g)
    MP_CHECKOK( mp_mul(&gx, &v, g) );

 CLEANUP:
  while(last >= 0)
    mp_clear(clean[last--]);

  return res;

} /* end mp_xgcd() */

/* }}} */

mp_size mp_trailing_zeros(const mp_int *mp)
{
  mp_digit d;
  mp_size  n = 0;
  int      ix;

  if (!mp || !MP_DIGITS(mp) || !mp_cmp_z(mp))
    return n;

  for (ix = 0; !(d = MP_DIGIT(mp,ix)) && (ix < MP_USED(mp)); ++ix)
    n += MP_DIGIT_BIT;
  if (!d)
    return 0;	/* shouldn't happen, but ... */
#if (MP_DIGIT_MAX > MP_32BIT_MAX) 
  if (!(d & 0xffffffffU)) {
    d >>= 32;
    n  += 32;
  }
#endif
  if (!(d & 0xffffU)) {
    d >>= 16;
    n  += 16;
  }
  if (!(d & 0xffU)) {
    d >>= 8;
    n  += 8;
  }
  if (!(d & 0xfU)) {
    d >>= 4;
    n  += 4;
  }
  if (!(d & 0x3U)) {
    d >>= 2;
    n  += 2;
  }
  if (!(d & 0x1U)) {
    d >>= 1;
    n  += 1;
  }
#if MP_ARGCHK == 2
  assert(0 != (d & 1));
#endif
  return n;
}

/* Given a and prime p, computes c and k such that a*c == 2**k (mod p).
** Returns k (positive) or error (negative).
** This technique from the paper "Fast Modular Reciprocals" (unpublished)
** by Richard Schroeppel (a.k.a. Captain Nemo).
*/
mp_err s_mp_almost_inverse(const mp_int *a, const mp_int *p, mp_int *c)
{
  mp_err res;
  mp_err k    = 0;
  mp_int d, f, g;

  ARGCHK(a && p && c, MP_BADARG);

  MP_DIGITS(&d) = 0;
  MP_DIGITS(&f) = 0;
  MP_DIGITS(&g) = 0;
  MP_CHECKOK( mp_init(&d) );
  MP_CHECKOK( mp_init_copy(&f, a) );	/* f = a */
  MP_CHECKOK( mp_init_copy(&g, p) );	/* g = p */

  mp_set(c, 1);
  mp_zero(&d);

  if (mp_cmp_z(&f) == 0) {
    res = MP_UNDEF;
  } else 
  for (;;) {
    int diff_sign;
    while (mp_iseven(&f)) {
      mp_size n = mp_trailing_zeros(&f);
      if (!n) {
	res = MP_UNDEF;
	goto CLEANUP;
      }
      s_mp_div_2d(&f, n);
      MP_CHECKOK( s_mp_mul_2d(&d, n) );
      k += n;
    }
    if (mp_cmp_d(&f, 1) == MP_EQ) {	/* f == 1 */
      res = k;
      break;
    }
    diff_sign = mp_cmp(&f, &g);
    if (diff_sign < 0) {		/* f < g */
      s_mp_exch(&f, &g);
      s_mp_exch(c, &d);
    } else if (diff_sign == 0) {		/* f == g */
      res = MP_UNDEF;		/* a and p are not relatively prime */
      break;
    }
    if ((MP_DIGIT(&f,0) % 4) == (MP_DIGIT(&g,0) % 4)) {
      MP_CHECKOK( mp_sub(&f, &g, &f) );	/* f = f - g */
      MP_CHECKOK( mp_sub(c,  &d,  c) );	/* c = c - d */
    } else {
      MP_CHECKOK( mp_add(&f, &g, &f) );	/* f = f + g */
      MP_CHECKOK( mp_add(c,  &d,  c) );	/* c = c + d */
    }
  }
  if (res >= 0) {
    while (MP_SIGN(c) != MP_ZPOS) {
      MP_CHECKOK( mp_add(c, p, c) );
    }
    res = k;
  }

CLEANUP:
  mp_clear(&d);
  mp_clear(&f);
  mp_clear(&g);
  return res;
}

/* Compute T = (P ** -1) mod MP_RADIX.  Also works for 16-bit mp_digits.
** This technique from the paper "Fast Modular Reciprocals" (unpublished)
** by Richard Schroeppel (a.k.a. Captain Nemo).
*/
mp_digit  s_mp_invmod_radix(mp_digit P)
{
  mp_digit T = P;
  T *= 2 - (P * T);
  T *= 2 - (P * T);
  T *= 2 - (P * T);
  T *= 2 - (P * T);
#if MP_DIGIT_MAX > MP_32BIT_MAX
  T *= 2 - (P * T);
  T *= 2 - (P * T);
#endif
  return T;
}

/* Given c, k, and prime p, where a*c == 2**k (mod p), 
** Compute x = (a ** -1) mod p.  This is similar to Montgomery reduction.
** This technique from the paper "Fast Modular Reciprocals" (unpublished)
** by Richard Schroeppel (a.k.a. Captain Nemo).
*/
mp_err  s_mp_fixup_reciprocal(const mp_int *c, const mp_int *p, int k, mp_int *x)
{
  int      k_orig = k;
  mp_digit r;
  mp_size  ix;
  mp_err   res;

  if (mp_cmp_z(c) < 0) {		/* c < 0 */
    MP_CHECKOK( mp_add(c, p, x) );	/* x = c + p */
  } else {
    MP_CHECKOK( mp_copy(c, x) );	/* x = c */
  }

  /* make sure x is large enough */
  ix = MP_HOWMANY(k, MP_DIGIT_BIT) + MP_USED(p) + 1;
  ix = MP_MAX(ix, MP_USED(x));
  MP_CHECKOK( s_mp_pad(x, ix) );

  r = 0 - s_mp_invmod_radix(MP_DIGIT(p,0));

  for (ix = 0; k > 0; ix++) {
    int      j = MP_MIN(k, MP_DIGIT_BIT);
    mp_digit v = r * MP_DIGIT(x, ix);
    if (j < MP_DIGIT_BIT) {
      v &= ((mp_digit)1 << j) - 1;	/* v = v mod (2 ** j) */
    }
    s_mp_mul_d_add_offset(p, v, x, ix); /* x += p * v * (RADIX ** ix) */
    k -= j;
  }
  s_mp_clamp(x);
  s_mp_div_2d(x, k_orig);
  res = MP_OKAY;

CLEANUP:
  return res;
}

/* compute mod inverse using Schroeppel's method, only if m is odd */
mp_err s_mp_invmod_odd_m(const mp_int *a, const mp_int *m, mp_int *c)
{
  int k;
  mp_err  res;
  mp_int  x;

  ARGCHK(a && m && c, MP_BADARG);

  if(mp_cmp_z(a) == 0 || mp_cmp_z(m) == 0)
    return MP_RANGE;
  if (mp_iseven(m))
    return MP_UNDEF;

  MP_DIGITS(&x) = 0;

  if (a == c) {
    if ((res = mp_init_copy(&x, a)) != MP_OKAY)
      return res;
    if (a == m) 
      m = &x;
    a = &x;
  } else if (m == c) {
    if ((res = mp_init_copy(&x, m)) != MP_OKAY)
      return res;
    m = &x;
  } else {
    MP_DIGITS(&x) = 0;
  }

  MP_CHECKOK( s_mp_almost_inverse(a, m, c) );
  k = res;
  MP_CHECKOK( s_mp_fixup_reciprocal(c, m, k, c) );
CLEANUP:
  mp_clear(&x);
  return res;
}

/* Known good algorithm for computing modular inverse.  But slow. */
mp_err mp_invmod_xgcd(const mp_int *a, const mp_int *m, mp_int *c)
{
  mp_int  g, x;
  mp_err  res;

  ARGCHK(a && m && c, MP_BADARG);

  if(mp_cmp_z(a) == 0 || mp_cmp_z(m) == 0)
    return MP_RANGE;

  MP_DIGITS(&g) = 0;
  MP_DIGITS(&x) = 0;
  MP_CHECKOK( mp_init(&x) );
  MP_CHECKOK( mp_init(&g) );

  MP_CHECKOK( mp_xgcd(a, m, &g, &x, NULL) );

  if (mp_cmp_d(&g, 1) != MP_EQ) {
    res = MP_UNDEF;
    goto CLEANUP;
  }

  res = mp_mod(&x, m, c);
  SIGN(c) = SIGN(a);

CLEANUP:
  mp_clear(&x);
  mp_clear(&g);

  return res;
}

/* modular inverse where modulus is 2**k. */
/* c = a**-1 mod 2**k */
mp_err s_mp_invmod_2d(const mp_int *a, mp_size k, mp_int *c)
{
  mp_err res;
  mp_size ix = k + 4;
  mp_int t0, t1, val, tmp, two2k;

  static const mp_digit d2 = 2;
  static const mp_int two = { MP_ZPOS, 1, 1, (mp_digit *)&d2 };

  if (mp_iseven(a))
    return MP_UNDEF;
  if (k <= MP_DIGIT_BIT) {
    mp_digit i = s_mp_invmod_radix(MP_DIGIT(a,0));
    if (k < MP_DIGIT_BIT)
      i &= ((mp_digit)1 << k) - (mp_digit)1;
    mp_set(c, i);
    return MP_OKAY;
  }
  MP_DIGITS(&t0) = 0;
  MP_DIGITS(&t1) = 0;
  MP_DIGITS(&val) = 0;
  MP_DIGITS(&tmp) = 0;
  MP_DIGITS(&two2k) = 0;
  MP_CHECKOK( mp_init_copy(&val, a) );
  s_mp_mod_2d(&val, k);
  MP_CHECKOK( mp_init_copy(&t0, &val) );
  MP_CHECKOK( mp_init_copy(&t1, &t0)  );
  MP_CHECKOK( mp_init(&tmp) );
  MP_CHECKOK( mp_init(&two2k) );
  MP_CHECKOK( s_mp_2expt(&two2k, k) );
  do {
    MP_CHECKOK( mp_mul(&val, &t1, &tmp)  );
    MP_CHECKOK( mp_sub(&two, &tmp, &tmp) );
    MP_CHECKOK( mp_mul(&t1, &tmp, &t1)   );
    s_mp_mod_2d(&t1, k);
    while (MP_SIGN(&t1) != MP_ZPOS) {
      MP_CHECKOK( mp_add(&t1, &two2k, &t1) );
    }
    if (mp_cmp(&t1, &t0) == MP_EQ) 
      break;
    MP_CHECKOK( mp_copy(&t1, &t0) );
  } while (--ix > 0);
  if (!ix) {
    res = MP_UNDEF;
  } else {
    mp_exch(c, &t1);
  }

CLEANUP:
  mp_clear(&t0);
  mp_clear(&t1);
  mp_clear(&val);
  mp_clear(&tmp);
  mp_clear(&two2k);
  return res;
}

mp_err s_mp_invmod_even_m(const mp_int *a, const mp_int *m, mp_int *c)
{
  mp_err res;
  mp_size k;
  mp_int oddFactor, evenFactor;	/* factors of the modulus */
  mp_int oddPart, evenPart;	/* parts to combine via CRT. */
  mp_int C2, tmp1, tmp2;

  static const mp_digit d1 = 1;
  static const mp_int one = { MP_ZPOS, 1, 1, (mp_digit *)&d1 };

  if ((res = s_mp_ispow2(m)) >= 0) {
    k = res;
    return s_mp_invmod_2d(a, k, c);
  }
  MP_DIGITS(&oddFactor) = 0;
  MP_DIGITS(&evenFactor) = 0;
  MP_DIGITS(&oddPart) = 0;
  MP_DIGITS(&evenPart) = 0;
  MP_DIGITS(&C2)     = 0;
  MP_DIGITS(&tmp1)   = 0;
  MP_DIGITS(&tmp2)   = 0;

  MP_CHECKOK( mp_init_copy(&oddFactor, m) );    /* oddFactor = m */
  MP_CHECKOK( mp_init(&evenFactor) );
  MP_CHECKOK( mp_init(&oddPart) );
  MP_CHECKOK( mp_init(&evenPart) );
  MP_CHECKOK( mp_init(&C2)     );
  MP_CHECKOK( mp_init(&tmp1)   );
  MP_CHECKOK( mp_init(&tmp2)   );

  k = mp_trailing_zeros(m);
  s_mp_div_2d(&oddFactor, k);
  MP_CHECKOK( s_mp_2expt(&evenFactor, k) );

  /* compute a**-1 mod oddFactor. */
  MP_CHECKOK( s_mp_invmod_odd_m(a, &oddFactor, &oddPart) );
  /* compute a**-1 mod evenFactor, where evenFactor == 2**k. */
  MP_CHECKOK( s_mp_invmod_2d(   a,       k,    &evenPart) );

  /* Use Chinese Remainer theorem to compute a**-1 mod m. */
  /* let m1 = oddFactor,  v1 = oddPart, 
   * let m2 = evenFactor, v2 = evenPart.
   */

  /* Compute C2 = m1**-1 mod m2. */
  MP_CHECKOK( s_mp_invmod_2d(&oddFactor, k,    &C2) );

  /* compute u = (v2 - v1)*C2 mod m2 */
  MP_CHECKOK( mp_sub(&evenPart, &oddPart,   &tmp1) );
  MP_CHECKOK( mp_mul(&tmp1,     &C2,        &tmp2) );
  s_mp_mod_2d(&tmp2, k);
  while (MP_SIGN(&tmp2) != MP_ZPOS) {
    MP_CHECKOK( mp_add(&tmp2, &evenFactor, &tmp2) );
  }

  /* compute answer = v1 + u*m1 */
  MP_CHECKOK( mp_mul(&tmp2,     &oddFactor, c) );
  MP_CHECKOK( mp_add(&oddPart,  c,          c) );
  /* not sure this is necessary, but it's low cost if not. */
  MP_CHECKOK( mp_mod(c,         m,          c) );

CLEANUP:
  mp_clear(&oddFactor);
  mp_clear(&evenFactor);
  mp_clear(&oddPart);
  mp_clear(&evenPart);
  mp_clear(&C2);
  mp_clear(&tmp1);
  mp_clear(&tmp2);
  return res;
}


/* {{{ mp_invmod(a, m, c) */

/*
  mp_invmod(a, m, c)

  Compute c = a^-1 (mod m), if there is an inverse for a (mod m).
  This is equivalent to the question of whether (a, m) = 1.  If not,
  MP_UNDEF is returned, and there is no inverse.
 */

mp_err mp_invmod(const mp_int *a, const mp_int *m, mp_int *c)
{

  ARGCHK(a && m && c, MP_BADARG);

  if(mp_cmp_z(a) == 0 || mp_cmp_z(m) == 0)
    return MP_RANGE;

  if (mp_isodd(m)) {
    return s_mp_invmod_odd_m(a, m, c);
  }
  if (mp_iseven(a))
    return MP_UNDEF;	/* not invertable */

  return s_mp_invmod_even_m(a, m, c);

} /* end mp_invmod() */

/* }}} */
#endif /* if MP_NUMTH */

/* }}} */

/*------------------------------------------------------------------------*/
/* {{{ mp_print(mp, ofp) */

#if MP_IOFUNC
/*
  mp_print(mp, ofp)

  Print a textual representation of the given mp_int on the output
  stream 'ofp'.  Output is generated using the internal radix.
 */

void   mp_print(mp_int *mp, FILE *ofp)
{
  int   ix;

  if(mp == NULL || ofp == NULL)
    return;

  fputc((SIGN(mp) == NEG) ? '-' : '+', ofp);

  for(ix = USED(mp) - 1; ix >= 0; ix--) {
    fprintf(ofp, DIGIT_FMT, DIGIT(mp, ix));
  }

} /* end mp_print() */

#endif /* if MP_IOFUNC */

/* }}} */

/*------------------------------------------------------------------------*/
/* {{{ More I/O Functions */

/* {{{ mp_read_raw(mp, str, len) */

/* 
   mp_read_raw(mp, str, len)

   Read in a raw value (base 256) into the given mp_int
 */

mp_err  mp_read_raw(mp_int *mp, char *str, int len)
{
  int            ix;
  mp_err         res;
  unsigned char *ustr = (unsigned char *)str;

  ARGCHK(mp != NULL && str != NULL && len > 0, MP_BADARG);

  mp_zero(mp);

  /* Get sign from first byte */
  if(ustr[0])
    SIGN(mp) = NEG;
  else
    SIGN(mp) = ZPOS;

  /* Read the rest of the digits */
  for(ix = 1; ix < len; ix++) {
#if DIGIT_MAX < 256
    if((res = s_mp_lshd(mp, 1)) != MP_OKAY)
      return res;
#else
    if((res = mp_mul_d(mp, 256, mp)) != MP_OKAY)
      return res;
#endif
    if((res = mp_add_d(mp, ustr[ix], mp)) != MP_OKAY)
      return res;
  }

  return MP_OKAY;

} /* end mp_read_raw() */

/* }}} */

/* {{{ mp_raw_size(mp) */

int    mp_raw_size(mp_int *mp)
{
  ARGCHK(mp != NULL, 0);

  return (USED(mp) * sizeof(mp_digit)) + 1;

} /* end mp_raw_size() */

/* }}} */

/* {{{ mp_toraw(mp, str) */

mp_err mp_toraw(mp_int *mp, char *str)
{
  int  ix, jx, pos = 1;

  ARGCHK(mp != NULL && str != NULL, MP_BADARG);

  str[0] = (char)SIGN(mp);

  /* Iterate over each digit... */
  for(ix = USED(mp) - 1; ix >= 0; ix--) {
    mp_digit  d = DIGIT(mp, ix);

    /* Unpack digit bytes, high order first */
    for(jx = sizeof(mp_digit) - 1; jx >= 0; jx--) {
      str[pos++] = (char)(d >> (jx * CHAR_BIT));
    }
  }

  return MP_OKAY;

} /* end mp_toraw() */

/* }}} */

/* {{{ mp_read_radix(mp, str, radix) */

/*
  mp_read_radix(mp, str, radix)

  Read an integer from the given string, and set mp to the resulting
  value.  The input is presumed to be in base 10.  Leading non-digit
  characters are ignored, and the function reads until a non-digit
  character or the end of the string.
 */

mp_err  mp_read_radix(mp_int *mp, const char *str, int radix)
{
  int     ix = 0, val = 0;
  mp_err  res;
  mp_sign sig = ZPOS;

  ARGCHK(mp != NULL && str != NULL && radix >= 2 && radix <= MAX_RADIX, 
	 MP_BADARG);

  mp_zero(mp);

  /* Skip leading non-digit characters until a digit or '-' or '+' */
  while(str[ix] && 
	(s_mp_tovalue(str[ix], radix) < 0) && 
	str[ix] != '-' &&
	str[ix] != '+') {
    ++ix;
  }

  if(str[ix] == '-') {
    sig = NEG;
    ++ix;
  } else if(str[ix] == '+') {
    sig = ZPOS; /* this is the default anyway... */
    ++ix;
  }

  while((val = s_mp_tovalue(str[ix], radix)) >= 0) {
    if((res = s_mp_mul_d(mp, radix)) != MP_OKAY)
      return res;
    if((res = s_mp_add_d(mp, val)) != MP_OKAY)
      return res;
    ++ix;
  }

  if(s_mp_cmp_d(mp, 0) == MP_EQ)
    SIGN(mp) = ZPOS;
  else
    SIGN(mp) = sig;

  return MP_OKAY;

} /* end mp_read_radix() */

/* }}} */

/* {{{ mp_radix_size(mp, radix) */

int    mp_radix_size(mp_int *mp, int radix)
{
  int  bits;

  if(!mp || radix < 2 || radix > MAX_RADIX)
    return 0;

  bits = USED(mp) * DIGIT_BIT - 1;
 
  return s_mp_outlen(bits, radix);

} /* end mp_radix_size() */

/* }}} */

/* {{{ mp_toradix(mp, str, radix) */

mp_err mp_toradix(mp_int *mp, char *str, int radix)
{
  int  ix, pos = 0;

  ARGCHK(mp != NULL && str != NULL, MP_BADARG);
  ARGCHK(radix > 1 && radix <= MAX_RADIX, MP_RANGE);

  if(mp_cmp_z(mp) == MP_EQ) {
    str[0] = '0';
    str[1] = '\0';
  } else {
    mp_err   res;
    mp_int   tmp;
    mp_sign  sgn;
    mp_digit rem, rdx = (mp_digit)radix;
    char     ch;

    if((res = mp_init_copy(&tmp, mp)) != MP_OKAY)
      return res;