mpi.c

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

      goto CLEANUP;
  }

  s_mp_exch(&s, c);

CLEANUP:
  mp_clear(&x);
X:
  mp_clear(&s);

  return res;

} /* end mp_expt_d() */

/* }}} */

/* }}} */

/*------------------------------------------------------------------------*/
/* {{{ Full arithmetic */

/* {{{ mp_abs(a, b) */

/*
  mp_abs(a, b)

  Compute b = |a|.  'a' and 'b' may be identical.
 */

mp_err mp_abs(const mp_int *a, mp_int *b)
{
  mp_err   res;

  ARGCHK(a != NULL && b != NULL, MP_BADARG);

  if((res = mp_copy(a, b)) != MP_OKAY)
    return res;

  SIGN(b) = ZPOS;

  return MP_OKAY;

} /* end mp_abs() */

/* }}} */

/* {{{ mp_neg(a, b) */

/*
  mp_neg(a, b)

  Compute b = -a.  'a' and 'b' may be identical.
 */

mp_err mp_neg(const mp_int *a, mp_int *b)
{
  mp_err   res;

  ARGCHK(a != NULL && b != NULL, MP_BADARG);

  if((res = mp_copy(a, b)) != MP_OKAY)
    return res;

  if(s_mp_cmp_d(b, 0) == MP_EQ) 
    SIGN(b) = ZPOS;
  else 
    SIGN(b) = (SIGN(b) == NEG) ? ZPOS : NEG;

  return MP_OKAY;

} /* end mp_neg() */

/* }}} */

/* {{{ mp_add(a, b, c) */

/*
  mp_add(a, b, c)

  Compute c = a + b.  All parameters may be identical.
 */

mp_err mp_add(const mp_int *a, const mp_int *b, mp_int *c)
{
  mp_err  res;

  ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);

  if(SIGN(a) == SIGN(b)) { /* same sign:  add values, keep sign */
    MP_CHECKOK( s_mp_add_3arg(a, b, c) );
  } else if(s_mp_cmp(a, b) >= 0) {  /* different sign: |a| >= |b|   */
    MP_CHECKOK( s_mp_sub_3arg(a, b, c) );
  } else {                          /* different sign: |a|  < |b|   */
    MP_CHECKOK( s_mp_sub_3arg(b, a, c) );
  }

  if (s_mp_cmp_d(c, 0) == MP_EQ)
    SIGN(c) = ZPOS;

CLEANUP:
  return res;

} /* end mp_add() */

/* }}} */

/* {{{ mp_sub(a, b, c) */

/*
  mp_sub(a, b, c)

  Compute c = a - b.  All parameters may be identical.
 */

mp_err mp_sub(const mp_int *a, const mp_int *b, mp_int *c)
{
  mp_err  res;
  int     magDiff;

  ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);

  if (a == b) {
    mp_zero(c);
    return MP_OKAY;
  }

  if (MP_SIGN(a) != MP_SIGN(b)) {
    MP_CHECKOK( s_mp_add_3arg(a, b, c) );
  } else if (!(magDiff = s_mp_cmp(a, b))) {
    mp_zero(c);
    res = MP_OKAY;
  } else if (magDiff > 0) {
    MP_CHECKOK( s_mp_sub_3arg(a, b, c) );
  } else {
    MP_CHECKOK( s_mp_sub_3arg(b, a, c) );
    MP_SIGN(c) = !MP_SIGN(a);
  }

  if (s_mp_cmp_d(c, 0) == MP_EQ)
    MP_SIGN(c) = MP_ZPOS;

CLEANUP:
  return res;

} /* end mp_sub() */

/* }}} */

/* {{{ mp_mul(a, b, c) */

/*
  mp_mul(a, b, c)

  Compute c = a * b.  All parameters may be identical.
 */
mp_err   mp_mul(const mp_int *a, const mp_int *b, mp_int * c)
{
  mp_digit *pb;
  mp_int   tmp;
  mp_err   res;
  mp_size  ib;
  mp_size  useda, usedb;

  ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);

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

  if (MP_USED(a) < MP_USED(b)) {
    const mp_int *xch = b;	/* switch a and b, to do fewer outer loops */
    b = a;
    a = xch;
  }

  MP_USED(c) = 1; MP_DIGIT(c, 0) = 0;
  if((res = s_mp_pad(c, USED(a) + USED(b))) != MP_OKAY)
    goto CLEANUP;

  pb = MP_DIGITS(b);
  s_mpv_mul_d(MP_DIGITS(a), MP_USED(a), *pb++, MP_DIGITS(c));

  /* Outer loop:  Digits of b */
  useda = MP_USED(a);
  usedb = MP_USED(b);
  for (ib = 1; ib < usedb; ib++) {
    mp_digit b_i    = *pb++;

    /* Inner product:  Digits of a */
    if (b_i)
      s_mpv_mul_d_add(MP_DIGITS(a), useda, b_i, MP_DIGITS(c) + ib);
    else
      MP_DIGIT(c, ib + useda) = b_i;
  }

  s_mp_clamp(c);

  if(SIGN(a) == SIGN(b) || s_mp_cmp_d(c, 0) == MP_EQ)
    SIGN(c) = ZPOS;
  else
    SIGN(c) = NEG;

CLEANUP:
  mp_clear(&tmp);
  return res;
} /* end mp_mul() */

/* }}} */

/* {{{ mp_sqr(a, sqr) */

#if MP_SQUARE
/*
  Computes the square of a.  This can be done more
  efficiently than a general multiplication, because many of the
  computation steps are redundant when squaring.  The inner product
  step is a bit more complicated, but we save a fair number of
  iterations of the multiplication loop.
 */

/* sqr = a^2;   Caller provides both a and tmp; */
mp_err   mp_sqr(const mp_int *a, mp_int *sqr)
{
  mp_digit *pa;
  mp_digit d;
  mp_err   res;
  mp_size  ix;
  mp_int   tmp;
  int      count;

  ARGCHK(a != NULL && sqr != NULL, MP_BADARG);

  if (a == sqr) {
    if((res = mp_init_copy(&tmp, a)) != MP_OKAY)
      return res;
    a = &tmp;
  } else {
    DIGITS(&tmp) = 0;
    res = MP_OKAY;
  }

  ix = 2 * MP_USED(a);
  if (ix > MP_ALLOC(sqr)) {
    MP_USED(sqr) = 1; 
    MP_CHECKOK( s_mp_grow(sqr, ix) );
  } 
  MP_USED(sqr) = ix;
  MP_DIGIT(sqr, 0) = 0;

  pa = MP_DIGITS(a);
  count = MP_USED(a) - 1;
  if (count > 0) {
    d = *pa++;
    s_mpv_mul_d(pa, count, d, MP_DIGITS(sqr) + 1);
    for (ix = 3; --count > 0; ix += 2) {
      d = *pa++;
      s_mpv_mul_d_add(pa, count, d, MP_DIGITS(sqr) + ix);
    } /* for(ix ...) */
    MP_DIGIT(sqr, MP_USED(sqr)-1) = 0; /* above loop stopped short of this. */

    /* now sqr *= 2 */
    s_mp_mul_2(sqr);
  } else {
    MP_DIGIT(sqr, 1) = 0;
  }

  /* now add the squares of the digits of a to sqr. */
  s_mpv_sqr_add_prop(MP_DIGITS(a), MP_USED(a), MP_DIGITS(sqr));

  SIGN(sqr) = ZPOS;
  s_mp_clamp(sqr);

CLEANUP:
  mp_clear(&tmp);
  return res;

} /* end mp_sqr() */
#endif

/* }}} */

/* {{{ mp_div(a, b, q, r) */

/*
  mp_div(a, b, q, r)

  Compute q = a / b and r = a mod b.  Input parameters may be re-used
  as output parameters.  If q or r is NULL, that portion of the
  computation will be discarded (although it will still be computed)
 */
mp_err mp_div(const mp_int *a, const mp_int *b, mp_int *q, mp_int *r)
{
  mp_err   res;
  mp_int   *pQ, *pR;
  mp_int   qtmp, rtmp;
  int      cmp;

  ARGCHK(a != NULL && b != NULL, MP_BADARG);

  if(mp_cmp_z(b) == MP_EQ)
    return MP_RANGE;

  DIGITS(&qtmp) = 0;
  DIGITS(&rtmp) = 0;
  /* Set up some temporaries... */
  if (!q || q == a || q == b) {
    if((res = mp_init_copy(&qtmp, a)) != MP_OKAY)
      return res;
    pQ = &qtmp;
  } else {
    if((res = mp_copy(a, q)) != MP_OKAY)
      return res;
    pQ = q;
  }
  if (!r || r == a || r == b) {
    if((res = mp_init_copy(&rtmp, b)) != MP_OKAY)
      goto CLEANUP;
    pR = &rtmp;
  } else {
    if((res = mp_copy(b, r)) != MP_OKAY)
      goto CLEANUP;
    pR = r;
  }

  /*
    If a <= b, we can compute the solution without division;
    otherwise, we actually do the work required.
   */
  if((cmp = s_mp_cmp(pQ, pR)) < 0) {
    s_mp_exch(pQ, pR);
    mp_zero(pQ);
  } else if(cmp == 0) {
    mp_set(pQ, 1);
    mp_zero(pR);
  } else {
    if((res = s_mp_div(pQ, pR)) != MP_OKAY)
      goto CLEANUP;
  }

  /* Compute the signs for the output  */
  SIGN(pR) = SIGN(a); /* Sr = Sa              */
  if(SIGN(a) == SIGN(b))
    SIGN(pQ) = ZPOS;  /* Sq = ZPOS if Sa = Sb */
  else
    SIGN(pQ) = NEG;   /* Sq = NEG if Sa != Sb */

  if(s_mp_cmp_d(pQ, 0) == MP_EQ)
    SIGN(pQ) = ZPOS;
  if(s_mp_cmp_d(pR, 0) == MP_EQ)
    SIGN(pR) = ZPOS;

  /* Copy output, if it is needed      */
  if(q && q != pQ) 
    s_mp_exch(pQ, q);

  if(r && r != pR) 
    s_mp_exch(pR, r);

CLEANUP:
  mp_clear(&rtmp);
  mp_clear(&qtmp);

  return res;

} /* end mp_div() */

/* }}} */

/* {{{ mp_div_2d(a, d, q, r) */

mp_err mp_div_2d(const mp_int *a, mp_digit d, mp_int *q, mp_int *r)
{
  mp_err  res;

  ARGCHK(a != NULL, MP_BADARG);

  if(q) {
    if((res = mp_copy(a, q)) != MP_OKAY)
      return res;
  }
  if(r) {
    if((res = mp_copy(a, r)) != MP_OKAY)
      return res;
  }
  if(q) {
    s_mp_div_2d(q, d);
  }
  if(r) {
    s_mp_mod_2d(r, d);
  }

  return MP_OKAY;

} /* end mp_div_2d() */

/* }}} */

/* {{{ mp_expt(a, b, c) */

/*
  mp_expt(a, b, c)

  Compute c = a ** b, that is, raise a to the b power.  Uses a
  standard iterative square-and-multiply technique.
 */

mp_err mp_expt(mp_int *a, mp_int *b, mp_int *c)
{
  mp_int   s, x;
  mp_err   res;
  mp_digit d;
  int      dig, bit;

  ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);

  if(mp_cmp_z(b) < 0)
    return MP_RANGE;

  if((res = mp_init(&s)) != MP_OKAY)
    return res;

  mp_set(&s, 1);

  if((res = mp_init_copy(&x, a)) != MP_OKAY)
    goto X;

  /* Loop over low-order digits in ascending order */
  for(dig = 0; dig < (USED(b) - 1); dig++) {
    d = DIGIT(b, dig);

    /* Loop over bits of each non-maximal digit */
    for(bit = 0; bit < DIGIT_BIT; bit++) {
      if(d & 1) {
	if((res = s_mp_mul(&s, &x)) != MP_OKAY) 
	  goto CLEANUP;
      }

      d >>= 1;
      
      if((res = s_mp_sqr(&x)) != MP_OKAY)
	goto CLEANUP;
    }
  }

  /* Consider now the last digit... */
  d = DIGIT(b, dig);

  while(d) {
    if(d & 1) {
      if((res = s_mp_mul(&s, &x)) != MP_OKAY)
	goto CLEANUP;
    }

    d >>= 1;

    if((res = s_mp_sqr(&x)) != MP_OKAY)
      goto CLEANUP;
  }
  
  if(mp_iseven(b))
    SIGN(&s) = SIGN(a);

  res = mp_copy(&s, c);

CLEANUP:
  mp_clear(&x);
X:
  mp_clear(&s);

  return res;

} /* end mp_expt() */

/* }}} */

/* {{{ mp_2expt(a, k) */

/* Compute a = 2^k */

mp_err mp_2expt(mp_int *a, mp_digit k)
{
  ARGCHK(a != NULL, MP_BADARG);

  return s_mp_2expt(a, k);

} /* end mp_2expt() */

/* }}} */

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

/*
  mp_mod(a, m, c)

  Compute c = a (mod m).  Result will always be 0 <= c < m.
 */

mp_err mp_mod(const mp_int *a, const mp_int *m, mp_int *c)
{
  mp_err  res;
  int     mag;

  ARGCHK(a != NULL && m != NULL && c != NULL, MP_BADARG);

  if(SIGN(m) == NEG)
    return MP_RANGE;

  /*
     If |a| > m, we need to divide to get the remainder and take the
     absolute value.  

     If |a| < m, we don't need to do any division, just copy and adjust
     the sign (if a is negative).

     If |a| == m, we can simply set the result to zero.

     This order is intended to minimize the average path length of the
     comparison chain on common workloads -- the most frequent cases are
     that |a| != m, so we do those first.
   */
  if((mag = s_mp_cmp(a, m)) > 0) {
    if((res = mp_div(a, m, NULL, c)) != MP_OKAY)
      return res;
    
    if(SIGN(c) == NEG) {
      if((res = mp_add(c, m, c)) != MP_OKAY)
	return res;
    }

  } else if(mag < 0) {
    if((res = mp_copy(a, c)) != MP_OKAY)
      return res;

    if(mp_cmp_z(a) < 0) {
      if((res = mp_add(c, m, c)) != MP_OKAY)
	return res;

    }
    
  } else {
    mp_zero(c);

  }

  return MP_OKAY;

} /* end mp_mod() */

/* }}} */

/* {{{ mp_mod_d(a, d, c) */

/*
  mp_mod_d(a, d, c)

  Compute c = a (mod d).  Result will always be 0 <= c < d
 */
mp_err mp_mod_d(const mp_int *a, mp_digit d, mp_digit *c)
{
  mp_err   res;
  mp_digit rem;

  ARGCHK(a != NULL && c != NULL, MP_BADARG);

  if(s_mp_cmp_d(a, d) > 0) {
    if((res = mp_div_d(a, d, NULL, &rem)) != MP_OKAY)
      return res;

  } else {
    if(SIGN(a) == NEG)
      rem = d - DIGIT(a, 0);
    else
      rem = DIGIT(a, 0);
  }

  if(c)
    *c = rem;

  return MP_OKAY;

} /* end mp_mod_d() */

/* }}} */

/* {{{ mp_sqrt(a, b) */

/*
  mp_sqrt(a, b)

  Compute the integer square root of a, and store the result in b.
  Uses an integer-arithmetic version of Newton's iterative linear
  approximation technique to determine this value; the result has the
  following two properties:

     b^2 <= a
     (b+1)^2 >= a

  It is a range error to pass a negative value.
 */
mp_err mp_sqrt(const mp_int *a, mp_int *b)
{
  mp_int   x, t;
  mp_err   res;
  mp_size  used;

  ARGCHK(a != NULL && b != NULL, MP_BADARG);

  /* Cannot take square root of a negative value */
  if(SIGN(a) == NEG)
    return MP_RANGE;

  /* Special cases for zero and one, trivial     */
  if(mp_cmp_d(a, 1) <= 0)
    return mp_copy(a, b);
    
  /* Initialize the temporaries we'll use below  */
  if((res = mp_init_size(&t, USED(a))) != MP_OKAY)
    return res;

  /* Compute an initial guess for the iteration as a itself */
  if((res = mp_init_copy(&x, a)) != MP_OKAY)
    goto X;

  used = MP_USED(&x);
  if (used > 1) {
    s_mp_rshd(&x, used / 2);
  }

  for(;;) {
    /* t = (x * x) - a */
    mp_copy(&x, &t);      /* can't fail, t is big enough for original x */
    if((res = mp_sqr(&t, &t)) != MP_OKAY ||
       (res = mp_sub(&t, a, &t)) != MP_OKAY)
      goto CLEANUP;

    /* t = t / 2x       */
    s_mp_mul_2(&x);
    if((res = mp_div(&t, &x, &t, NULL)) != MP_OKAY)
      goto CLEANUP;
    s_mp_div_2(&x);