imath.c

samba最新软件

{
  mp_size  i, j;
  mp_word  w;

  for(i = 0; i < size_a; ++i, dc += 2, ++da) {
    mp_digit  *dct = dc, *dat = da;

    if(*da == 0)
      continue;

    /* Take care of the first digit, no rollover */
    w = (mp_word)*dat * (mp_word)*dat + (mp_word)*dct;
    *dct = LOWER_HALF(w);
    w = UPPER_HALF(w);
    ++dat; ++dct;

    for(j = i + 1; j < size_a; ++j, ++dat, ++dct) {
      mp_word  t = (mp_word)*da * (mp_word)*dat;
      mp_word  u = w + (mp_word)*dct, ov = 0;

      /* Check if doubling t will overflow a word */
      if(HIGH_BIT_SET(t))
	ov = 1;

      w = t + t;

      /* Check if adding u to w will overflow a word */
      if(ADD_WILL_OVERFLOW(w, u))
	ov = 1;

      w += u;

      *dct = LOWER_HALF(w);
      w = UPPER_HALF(w);
      if(ov) {
	w += MP_DIGIT_MAX; /* MP_RADIX */
	++w;
      }
    }

    w = w + *dct;
    *dct = (mp_digit)w; 
    while((w = UPPER_HALF(w)) != 0) {
      ++dct; w = w + *dct;
      *dct = LOWER_HALF(w);
    }

    assert(w == 0);
  }
}

/* }}} */

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

static void      s_dadd(mp_int a, mp_digit b)
{
  mp_word   w = 0;
  mp_digit *da = MP_DIGITS(a);
  mp_size   ua = MP_USED(a);

  w = (mp_word)*da + b;
  *da++ = LOWER_HALF(w);
  w = UPPER_HALF(w);

  for(ua -= 1; ua > 0; --ua, ++da) {
    w = (mp_word)*da + w;

    *da = LOWER_HALF(w);
    w = UPPER_HALF(w);
  }

  if(w) {
    *da = (mp_digit)w;
    MP_USED(a) += 1;
  }
}

/* }}} */

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

static void      s_dmul(mp_int a, mp_digit b)
{
  mp_word   w = 0;
  mp_digit *da = MP_DIGITS(a);
  mp_size   ua = MP_USED(a);

  while(ua > 0) {
    w = (mp_word)*da * b + w;
    *da++ = LOWER_HALF(w);
    w = UPPER_HALF(w);
    --ua;
  }

  if(w) {
    *da = (mp_digit)w;
    MP_USED(a) += 1;
  }
}

/* }}} */

/* {{{ s_dbmul(da, b, dc, size_a) */

static void      s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, mp_size size_a)
{
  mp_word  w = 0;

  while(size_a > 0) {
    w = (mp_word)*da++ * (mp_word)b + w;

    *dc++ = LOWER_HALF(w);
    w = UPPER_HALF(w);
    --size_a;
  }

  if(w)
    *dc = LOWER_HALF(w);
}

/* }}} */

/* {{{ s_ddiv(da, d, dc, size_a) */

static mp_digit  s_ddiv(mp_int a, mp_digit b)
{
  mp_word   w = 0, qdigit;
  mp_size   ua = MP_USED(a);
  mp_digit *da = MP_DIGITS(a) + ua - 1;
  
  for(/* */; ua > 0; --ua, --da) {
    w = (w << MP_DIGIT_BIT) | *da;

    if(w >= b) {
      qdigit = w / b;
      w = w % b;
    } 
    else {
      qdigit = 0;
    }
      
    *da = (mp_digit)qdigit;
  }

  CLAMP(a);
  return (mp_digit)w;
}

/* }}} */

/* {{{ s_qdiv(z, p2) */

static void     s_qdiv(mp_int z, mp_size p2)
{
  mp_size ndig = p2 / MP_DIGIT_BIT, nbits = p2 % MP_DIGIT_BIT;
  mp_size uz = MP_USED(z);

  if(ndig) {
    mp_size  mark;
    mp_digit *to, *from;

    if(ndig >= uz) {
      mp_int_zero(z);
      return;
    }

    to = MP_DIGITS(z); from = to + ndig;

    for(mark = ndig; mark < uz; ++mark) 
      *to++ = *from++;

    MP_USED(z) = uz - ndig;
  }

  if(nbits) {
    mp_digit d = 0, *dz, save;
    mp_size  up = MP_DIGIT_BIT - nbits;

    uz = MP_USED(z);
    dz = MP_DIGITS(z) + uz - 1;

    for(/* */; uz > 0; --uz, --dz) {
      save = *dz;

      *dz = (*dz >> nbits) | (d << up);
      d = save;
    }

    CLAMP(z);
  }

  if(MP_USED(z) == 1 && z->digits[0] == 0)
    MP_SIGN(z) = MP_ZPOS;
}

/* }}} */

/* {{{ s_qmod(z, p2) */

static void     s_qmod(mp_int z, mp_size p2)
{
  mp_size   start = p2 / MP_DIGIT_BIT + 1, rest = p2 % MP_DIGIT_BIT;
  mp_size   uz = MP_USED(z);
  mp_digit  mask = (1 << rest) - 1;

  if(start <= uz) {
    MP_USED(z) = start;
    z->digits[start - 1] &= mask;
    CLAMP(z);
  }
}

/* }}} */

/* {{{ s_qmul(z, p2) */

static int      s_qmul(mp_int z, mp_size p2)
{
  mp_size   uz, need, rest, extra, i;
  mp_digit *from, *to, d;

  if(p2 == 0)
    return 1;

  uz = MP_USED(z); 
  need = p2 / MP_DIGIT_BIT; rest = p2 % MP_DIGIT_BIT;

  /* Figure out if we need an extra digit at the top end; this occurs
     if the topmost `rest' bits of the high-order digit of z are not
     zero, meaning they will be shifted off the end if not preserved */
  extra = 0;
  if(rest != 0) {
    mp_digit *dz = MP_DIGITS(z) + uz - 1;

    if((*dz >> (MP_DIGIT_BIT - rest)) != 0)
      extra = 1;
  }

  if(!s_pad(z, uz + need + extra))
    return 0;

  /* If we need to shift by whole digits, do that in one pass, then
     to back and shift by partial digits.
   */
  if(need > 0) {
    from = MP_DIGITS(z) + uz - 1;
    to = from + need;

    for(i = 0; i < uz; ++i)
      *to-- = *from--;

    ZERO(MP_DIGITS(z), need);
    uz += need;
  }

  if(rest) {
    d = 0;
    for(i = need, from = MP_DIGITS(z) + need; i < uz; ++i, ++from) {
      mp_digit save = *from;
      
      *from = (*from << rest) | (d >> (MP_DIGIT_BIT - rest));
      d = save;
    }

    d >>= (MP_DIGIT_BIT - rest);
    if(d != 0) {
      *from = d;
      uz += extra;
    }
  }

  MP_USED(z) = uz;
  CLAMP(z);

  return 1;
}

/* }}} */

/* {{{ s_qsub(z, p2) */

/* Compute z = 2^p2 - |z|; requires that 2^p2 >= |z|
   The sign of the result is always zero/positive.
 */
static int       s_qsub(mp_int z, mp_size p2)
{
  mp_digit hi = (1 << (p2 % MP_DIGIT_BIT)), *zp;
  mp_size  tdig = (p2 / MP_DIGIT_BIT), pos;
  mp_word  w = 0;

  if(!s_pad(z, tdig + 1))
    return 0;

  for(pos = 0, zp = MP_DIGITS(z); pos < tdig; ++pos, ++zp) {
    w = ((mp_word) MP_DIGIT_MAX + 1) - w - (mp_word)*zp;

    *zp = LOWER_HALF(w);
    w = UPPER_HALF(w) ? 0 : 1;
  }

  w = ((mp_word) MP_DIGIT_MAX + 1 + hi) - w - (mp_word)*zp;
  *zp = LOWER_HALF(w);

  assert(UPPER_HALF(w) != 0); /* no borrow out should be possible */
  
  MP_SIGN(z) = MP_ZPOS;
  CLAMP(z);

  return 1;
}

/* }}} */

/* {{{ s_dp2k(z) */

static int      s_dp2k(mp_int z)
{
  int       k = 0;
  mp_digit *dp = MP_DIGITS(z), d;

  if(MP_USED(z) == 1 && *dp == 0)
    return 1;

  while(*dp == 0) {
    k += MP_DIGIT_BIT;
    ++dp;
  }
  
  d = *dp;
  while((d & 1) == 0) {
    d >>= 1;
    ++k;
  }

  return k;
}

/* }}} */

/* {{{ s_isp2(z) */

static int       s_isp2(mp_int z)
{
  mp_size uz = MP_USED(z), k = 0;
  mp_digit *dz = MP_DIGITS(z), d;

  while(uz > 1) {
    if(*dz++ != 0)
      return -1;
    k += MP_DIGIT_BIT;
    --uz;
  }

  d = *dz;
  while(d > 1) {
    if(d & 1)
      return -1;
    ++k; d >>= 1;
  }

  return (int) k;
}

/* }}} */

/* {{{ s_2expt(z, k) */

static int       s_2expt(mp_int z, int k)
{
  mp_size  ndig, rest;
  mp_digit *dz;

  ndig = (k + MP_DIGIT_BIT) / MP_DIGIT_BIT;
  rest = k % MP_DIGIT_BIT;

  if(!s_pad(z, ndig))
    return 0;

  dz = MP_DIGITS(z);
  ZERO(dz, ndig);
  *(dz + ndig - 1) = (1 << rest);
  MP_USED(z) = ndig;

  return 1;
}

/* }}} */

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

static int      s_norm(mp_int a, mp_int b)
{
  mp_digit d = b->digits[MP_USED(b) - 1];
  int      k = 0;

  while(d < (mp_digit) (1 << (MP_DIGIT_BIT - 1))) { /* d < (MP_RADIX / 2) */
    d <<= 1;
    ++k;
  }

  /* These multiplications can't fail */
  if(k != 0) {
    (void) s_qmul(a, (mp_size) k);
    (void) s_qmul(b, (mp_size) k);
  }

  return k;
}

/* }}} */

/* {{{ s_brmu(z, m) */

static mp_result s_brmu(mp_int z, mp_int m)
{
  mp_size um = MP_USED(m) * 2;

  if(!s_pad(z, um))
    return MP_MEMORY;

  s_2expt(z, MP_DIGIT_BIT * um);
  return mp_int_div(z, m, z, NULL);
}

/* }}} */

/* {{{ s_reduce(x, m, mu, q1, q2) */

static int       s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2)
{
  mp_size   um = MP_USED(m), umb_p1, umb_m1;

  umb_p1 = (um + 1) * MP_DIGIT_BIT;
  umb_m1 = (um - 1) * MP_DIGIT_BIT;

  if(mp_int_copy(x, q1) != MP_OK)
    return 0;

  /* Compute q2 = floor((floor(x / b^(k-1)) * mu) / b^(k+1)) */
  s_qdiv(q1, umb_m1);
  UMUL(q1, mu, q2);
  s_qdiv(q2, umb_p1);

  /* Set x = x mod b^(k+1) */
  s_qmod(x, umb_p1);

  /* Now, q is a guess for the quotient a / m.
     Compute x - q * m mod b^(k+1), replacing x.  This may be off
     by a factor of 2m, but no more than that.
   */
  UMUL(q2, m, q1);
  s_qmod(q1, umb_p1);
  (void) mp_int_sub(x, q1, x); /* can't fail */

  /* The result may be < 0; if it is, add b^(k+1) to pin it in the
     proper range. */
  if((CMPZ(x) < 0) && !s_qsub(x, umb_p1))
    return 0;

  /* If x > m, we need to back it off until it is in range.
     This will be required at most twice.  */
  if(mp_int_compare(x, m) >= 0) {
    (void) mp_int_sub(x, m, x);
    if(mp_int_compare(x, m) >= 0)
      (void) mp_int_sub(x, m, x);
  }

  /* At this point, x has been properly reduced. */
  return 1;
}

/* }}} */

/* {{{ s_embar(a, b, m, mu, c) */

/* Perform modular exponentiation using Barrett's method, where mu is
   the reduction constant for m.  Assumes a < m, b > 0. */
static mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c)
{
  mp_digit  *db, *dbt, umu, d;
  mpz_t     temp[3]; 
  mp_result res;
  int       last = 0;

  umu = MP_USED(mu); db = MP_DIGITS(b); dbt = db + MP_USED(b) - 1;

  while(last < 3) {
    SETUP(mp_int_init_size(TEMP(last), 4 * umu), last);
    ZERO(MP_DIGITS(TEMP(last - 1)), MP_ALLOC(TEMP(last - 1)));
  }

  (void) mp_int_set_value(c, 1);

  /* Take care of low-order digits */
  while(db < dbt) {
    int      i;

    for(d = *db, i = MP_DIGIT_BIT; i > 0; --i, d >>= 1) {
      if(d & 1) {
	/* The use of a second temporary avoids allocation */
	UMUL(c, a, TEMP(0));
	if(!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) {
	  res = MP_MEMORY; goto CLEANUP;
	}
	mp_int_copy(TEMP(0), c);
      }


      USQR(a, TEMP(0));
      assert(MP_SIGN(TEMP(0)) == MP_ZPOS);
      if(!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) {
	res = MP_MEMORY; goto CLEANUP;
      }
      assert(MP_SIGN(TEMP(0)) == MP_ZPOS);
      mp_int_copy(TEMP(0), a);


    }

    ++db;
  }

  /* Take care of highest-order digit */
  d = *dbt;
  for(;;) {
    if(d & 1) {
      UMUL(c, a, TEMP(0));
      if(!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) {
	res = MP_MEMORY; goto CLEANUP;
      }
      mp_int_copy(TEMP(0), c);
    }
    
    d >>= 1;
    if(!d) break;

    USQR(a, TEMP(0));
    if(!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) {
      res = MP_MEMORY; goto CLEANUP;
    }
    (void) mp_int_copy(TEMP(0), a);
  }

 CLEANUP:
  while(--last >= 0)
    mp_int_clear(TEMP(last));
  
  return res;
}

/* }}} */

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

/* Precondition:  a >= b and b > 0
   Postcondition: a' = a / b, b' = a % b
 */
static mp_result s_udiv(mp_int a, mp_int b)
{
  mpz_t     q, r, t;
  mp_size   ua, ub, qpos = 0;
  mp_digit *da, btop;
  mp_result res = MP_OK;
  int       k, skip = 0;

  /* Force signs to positive */
  MP_SIGN(a) = MP_ZPOS;
  MP_SIGN(b) = MP_ZPOS;

  /* Normalize, per Knuth */
  k = s_norm(a, b);

  ua = MP_USED(a); ub = MP_USED(b); btop = b->digits[ub - 1];
  if((res = mp_int_init_size(&q, ua)) != MP_OK) return res;
  if((res = mp_int_init_size(&t, ua + 1)) != MP_OK) goto CLEANUP;

  da = MP_DIGITS(a);
  r.digits = da + ua - 1;  /* The contents of r are shared with a */
  r.used   = 1;
  r.sign   = MP_ZPOS;
  r.alloc  = MP_ALLOC(a);
  ZERO(t.digits, t.alloc);

  /* Solve for quotient digits, store in q.digits in reverse order */
  while(r.digits >= da) {
    assert(qpos <= q.alloc);

    if(s_ucmp(b, &r) > 0) {
      r.digits -= 1;
      r.used += 1;
      
      if(++skip > 1 && qpos > 0) 
	q.digits[qpos++] = 0;
      
      CLAMP(&r);
    }
    else {
      mp_word  pfx = r.digits[r.used - 1];
      mp_word  qdigit;
      
      if(r.used > 1 && pfx <= btop) {
	pfx <<= MP_DIGIT_BIT / 2;
	pfx <<= MP_DIGIT_BIT / 2;
	pfx |= r.digits[r.used - 2];
      }

      qdigit = pfx / btop;
      if(qdigit > MP_DIGIT_MAX