imath.c

samba最新软件

/* {{{ mp_int_exptmod_evalue(a, value, m, c) */

mp_result mp_int_exptmod_evalue(mp_int a, int value, mp_int m, mp_int c)
{
  mpz_t    vtmp;
  mp_digit vbuf[MP_VALUE_DIGITS(value)];

  s_fake(&vtmp, value, vbuf);

  return mp_int_exptmod(a, &vtmp, m, c);
}

/* }}} */

/* {{{ mp_int_exptmod_bvalue(v, b, m, c) */

mp_result mp_int_exptmod_bvalue(int value, mp_int b,
				mp_int m, mp_int c)
{
  mpz_t    vtmp;
  mp_digit vbuf[MP_VALUE_DIGITS(value)];

  s_fake(&vtmp, value, vbuf);

  return mp_int_exptmod(&vtmp, b, m, c);
}

/* }}} */

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

mp_result mp_int_exptmod_known(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c)
{
  mp_result res;
  mp_size   um;
  mpz_t     temp[2];
  mp_int    s;
  int       last = 0;

  CHECK(a && b && m && c);

  /* Zero moduli and negative exponents are not considered. */
  if(CMPZ(m) == 0)
    return MP_UNDEF;
  if(CMPZ(b) < 0)
    return MP_RANGE;

  um = MP_USED(m);
  SETUP(mp_int_init_size(TEMP(0), 2 * um), last);

  if(c == b || c == m) {
    SETUP(mp_int_init_size(TEMP(1), 2 * um), last);
    s = TEMP(1);
  } 
  else {
    s = c;
  }
  
  if((res = mp_int_mod(a, m, TEMP(0))) != MP_OK) goto CLEANUP;

  if((res = s_embar(TEMP(0), b, m, mu, s)) != MP_OK)
    goto CLEANUP;

  res = mp_int_copy(s, c);

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

  return res;
}

/* }}} */

/* {{{ mp_int_redux_const(m, c) */

mp_result mp_int_redux_const(mp_int m, mp_int c)
{
  CHECK(m != NULL && c != NULL && m != c);

  return s_brmu(c, m);
}

/* }}} */

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

mp_result mp_int_invmod(mp_int a, mp_int m, mp_int c)
{
  mp_result res;
  mp_sign   sa;
  int       last = 0;
  mpz_t     temp[2];

  CHECK(a != NULL && m != NULL && c != NULL);

  if(CMPZ(a) == 0 || CMPZ(m) <= 0)
    return MP_RANGE;

  sa = MP_SIGN(a); /* need this for the result later */

  for(last = 0; last < 2; ++last)
    mp_int_init(TEMP(last));

  if((res = mp_int_egcd(a, m, TEMP(0), TEMP(1), NULL)) != MP_OK) 
    goto CLEANUP;

  if(mp_int_compare_value(TEMP(0), 1) != 0) {
    res = MP_UNDEF;
    goto CLEANUP;
  }

  /* It is first necessary to constrain the value to the proper range */
  if((res = mp_int_mod(TEMP(1), m, TEMP(1))) != MP_OK)
    goto CLEANUP;

  /* Now, if 'a' was originally negative, the value we have is
     actually the magnitude of the negative representative; to get the
     positive value we have to subtract from the modulus.  Otherwise,
     the value is okay as it stands.
   */
  if(sa == MP_NEG)
    res = mp_int_sub(m, TEMP(1), c);
  else
    res = mp_int_copy(TEMP(1), c);

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

  return res;
}

/* }}} */

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

/* Binary GCD algorithm due to Josef Stein, 1961 */
mp_result mp_int_gcd(mp_int a, mp_int b, mp_int c)
{ 
  int       ca, cb, k = 0;
  mpz_t     u, v, t;
  mp_result res;

  CHECK(a != NULL && b != NULL && c != NULL);

  ca = CMPZ(a);
  cb = CMPZ(b);
  if(ca == 0 && cb == 0)
    return MP_UNDEF;
  else if(ca == 0) 
    return mp_int_abs(b, c);
  else if(cb == 0) 
    return mp_int_abs(a, c);

  mp_int_init(&t);
  if((res = mp_int_init_copy(&u, a)) != MP_OK)
    goto U;
  if((res = mp_int_init_copy(&v, b)) != MP_OK)
    goto V;

  MP_SIGN(&u) = MP_ZPOS; MP_SIGN(&v) = MP_ZPOS;

  { /* Divide out common factors of 2 from u and v */
    int div2_u = s_dp2k(&u), div2_v = s_dp2k(&v);
   
    k = MIN(div2_u, div2_v);
    s_qdiv(&u, (mp_size) k);
    s_qdiv(&v, (mp_size) k);
  }
  
  if(mp_int_is_odd(&u)) {
    if((res = mp_int_neg(&v, &t)) != MP_OK)
      goto CLEANUP;
  } 
  else {
    if((res = mp_int_copy(&u, &t)) != MP_OK)
      goto CLEANUP;
  }

  for(;;) {
    s_qdiv(&t, s_dp2k(&t));

    if(CMPZ(&t) > 0) {
      if((res = mp_int_copy(&t, &u)) != MP_OK)
	goto CLEANUP;
    } 
    else {
      if((res = mp_int_neg(&t, &v)) != MP_OK)
	goto CLEANUP;
    }

    if((res = mp_int_sub(&u, &v, &t)) != MP_OK)
      goto CLEANUP;

    if(CMPZ(&t) == 0)
      break;
  } 

  if((res = mp_int_abs(&u, c)) != MP_OK)
    goto CLEANUP;
  if(!s_qmul(c, (mp_size) k))
    res = MP_MEMORY;
  
 CLEANUP:
  mp_int_clear(&v);
 V: mp_int_clear(&u);
 U: mp_int_clear(&t);

  return res;
}

/* }}} */

/* {{{ mp_int_egcd(a, b, c, x, y) */

/* This is the binary GCD algorithm again, but this time we keep track
   of the elementary matrix operations as we go, so we can get values
   x and y satisfying c = ax + by.
 */
mp_result mp_int_egcd(mp_int a, mp_int b, mp_int c, 
		      mp_int x, mp_int y)
{ 
  int       k, last = 0, ca, cb;
  mpz_t     temp[8];
  mp_result res;
  
  CHECK(a != NULL && b != NULL && c != NULL && 
	(x != NULL || y != NULL));

  ca = CMPZ(a);
  cb = CMPZ(b);
  if(ca == 0 && cb == 0)
    return MP_UNDEF;
  else if(ca == 0) {
    if((res = mp_int_abs(b, c)) != MP_OK) return res;
    mp_int_zero(x); (void) mp_int_set_value(y, 1); return MP_OK;
  } 
  else if(cb == 0) {
    if((res = mp_int_abs(a, c)) != MP_OK) return res;
    (void) mp_int_set_value(x, 1); mp_int_zero(y); return MP_OK;
  }

  /* Initialize temporaries:
     A:0, B:1, C:2, D:3, u:4, v:5, ou:6, ov:7 */
  for(last = 0; last < 4; ++last) 
    mp_int_init(TEMP(last));
  TEMP(0)->digits[0] = 1;
  TEMP(3)->digits[0] = 1;

  SETUP(mp_int_init_copy(TEMP(4), a), last);
  SETUP(mp_int_init_copy(TEMP(5), b), last);

  /* We will work with absolute values here */
  MP_SIGN(TEMP(4)) = MP_ZPOS;
  MP_SIGN(TEMP(5)) = MP_ZPOS;

  { /* Divide out common factors of 2 from u and v */
    int  div2_u = s_dp2k(TEMP(4)), div2_v = s_dp2k(TEMP(5));
    
    k = MIN(div2_u, div2_v);
    s_qdiv(TEMP(4), k);
    s_qdiv(TEMP(5), k);
  }

  SETUP(mp_int_init_copy(TEMP(6), TEMP(4)), last);
  SETUP(mp_int_init_copy(TEMP(7), TEMP(5)), last);

  for(;;) {
    while(mp_int_is_even(TEMP(4))) {
      s_qdiv(TEMP(4), 1);
      
      if(mp_int_is_odd(TEMP(0)) || mp_int_is_odd(TEMP(1))) {
	if((res = mp_int_add(TEMP(0), TEMP(7), TEMP(0))) != MP_OK) 
	  goto CLEANUP;
	if((res = mp_int_sub(TEMP(1), TEMP(6), TEMP(1))) != MP_OK) 
	  goto CLEANUP;
      }

      s_qdiv(TEMP(0), 1);
      s_qdiv(TEMP(1), 1);
    }
    
    while(mp_int_is_even(TEMP(5))) {
      s_qdiv(TEMP(5), 1);

      if(mp_int_is_odd(TEMP(2)) || mp_int_is_odd(TEMP(3))) {
	if((res = mp_int_add(TEMP(2), TEMP(7), TEMP(2))) != MP_OK) 
	  goto CLEANUP;
	if((res = mp_int_sub(TEMP(3), TEMP(6), TEMP(3))) != MP_OK) 
	  goto CLEANUP;
      }

      s_qdiv(TEMP(2), 1);
      s_qdiv(TEMP(3), 1);
    }

    if(mp_int_compare(TEMP(4), TEMP(5)) >= 0) {
      if((res = mp_int_sub(TEMP(4), TEMP(5), TEMP(4))) != MP_OK) goto CLEANUP;
      if((res = mp_int_sub(TEMP(0), TEMP(2), TEMP(0))) != MP_OK) goto CLEANUP;
      if((res = mp_int_sub(TEMP(1), TEMP(3), TEMP(1))) != MP_OK) goto CLEANUP;
    } 
    else {
      if((res = mp_int_sub(TEMP(5), TEMP(4), TEMP(5))) != MP_OK) goto CLEANUP;
      if((res = mp_int_sub(TEMP(2), TEMP(0), TEMP(2))) != MP_OK) goto CLEANUP;
      if((res = mp_int_sub(TEMP(3), TEMP(1), TEMP(3))) != MP_OK) goto CLEANUP;
    }

    if(CMPZ(TEMP(4)) == 0) {
      if(x && (res = mp_int_copy(TEMP(2), x)) != MP_OK) goto CLEANUP;
      if(y && (res = mp_int_copy(TEMP(3), y)) != MP_OK) goto CLEANUP;
      if(c) {
	if(!s_qmul(TEMP(5), k)) {
	  res = MP_MEMORY;
	  goto CLEANUP;
	}
	 
	res = mp_int_copy(TEMP(5), c);
      }

      break;
    }
  }

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

  return res;
}

/* }}} */

/* {{{ mp_int_divisible_value(a, v) */

int       mp_int_divisible_value(mp_int a, int v)
{
  int       rem = 0;

  if(mp_int_div_value(a, v, NULL, &rem) != MP_OK)
    return 0;

  return rem == 0;
}

/* }}} */

/* {{{ mp_int_is_pow2(z) */

int       mp_int_is_pow2(mp_int z)
{
  CHECK(z != NULL);

  return s_isp2(z);
}

/* }}} */

/* {{{ mp_int_sqrt(a, c) */

mp_result mp_int_sqrt(mp_int a, mp_int c)
{
  mp_result  res = MP_OK;
  mpz_t      temp[2];
  int        last = 0;

  CHECK(a != NULL && c != NULL);

  /* The square root of a negative value does not exist in the integers. */
  if(MP_SIGN(a) == MP_NEG)
    return MP_UNDEF;

  SETUP(mp_int_init_copy(TEMP(last), a), last);
  SETUP(mp_int_init(TEMP(last)), last);

  for(;;) {
    if((res = mp_int_sqr(TEMP(0), TEMP(1))) != MP_OK)
      goto CLEANUP;

    if(mp_int_compare_unsigned(a, TEMP(1)) == 0) break;

    if((res = mp_int_copy(a, TEMP(1))) != MP_OK) 
      goto CLEANUP;
    if((res = mp_int_div(TEMP(1), TEMP(0), TEMP(1), NULL)) != MP_OK) 
      goto CLEANUP;
    if((res = mp_int_add(TEMP(0), TEMP(1), TEMP(1))) != MP_OK) 
      goto CLEANUP;
    if((res = mp_int_div_pow2(TEMP(1), 1, TEMP(1), NULL)) != MP_OK)
      goto CLEANUP;

    if(mp_int_compare_unsigned(TEMP(0), TEMP(1)) == 0) break;
    if((res = mp_int_sub_value(TEMP(0), 1, TEMP(0))) != MP_OK) goto CLEANUP;
    if(mp_int_compare_unsigned(TEMP(0), TEMP(1)) == 0) break;

    if((res = mp_int_copy(TEMP(1), TEMP(0))) != MP_OK) goto CLEANUP;
  }
  
  res = mp_int_copy(TEMP(0), c);

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

/* }}} */

/* {{{ mp_int_to_int(z, out) */

mp_result mp_int_to_int(mp_int z, int *out)
{
  unsigned int uv = 0;
  mp_size   uz;
  mp_digit *dz;
  mp_sign   sz;

  CHECK(z != NULL);

  /* Make sure the value is representable as an int */
  sz = MP_SIGN(z);
  if((sz == MP_ZPOS && mp_int_compare_value(z, INT_MAX) > 0) ||
     mp_int_compare_value(z, INT_MIN) < 0)
    return MP_RANGE;
     
  uz = MP_USED(z);
  dz = MP_DIGITS(z) + uz - 1;
  
  while(uz > 0) {
    uv <<= MP_DIGIT_BIT/2;
    uv = (uv << (MP_DIGIT_BIT/2)) | *dz--;
    --uz;
  }

  if(out)
    *out = (sz == MP_NEG) ? -(int)uv : (int)uv;

  return MP_OK;
}

/* }}} */

/* {{{ mp_int_to_string(z, radix, str, limit) */

mp_result mp_int_to_string(mp_int z, mp_size radix, 
			   char *str, int limit)
{
  mp_result res;
  int       cmp = 0;

  CHECK(z != NULL && str != NULL && limit >= 2);

  if(radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
    return MP_RANGE;

  if(CMPZ(z) == 0) {
    *str++ = s_val2ch(0, 1);
  } 
  else {
    mpz_t tmp;
    char  *h, *t;

    if((res = mp_int_init_copy(&tmp, z)) != MP_OK)
      return res;

    if(MP_SIGN(z) == MP_NEG) {
      *str++ = '-';
      --limit;
    }
    h = str;

    /* Generate digits in reverse order until finished or limit reached */
    for(/* */; limit > 0; --limit) {
      mp_digit d;

      if((cmp = CMPZ(&tmp)) == 0)
	break;

      d = s_ddiv(&tmp, (mp_digit)radix);
      *str++ = s_val2ch(d, 1);
    }
    t = str - 1;

    /* Put digits back in correct output order */
    while(h < t) {
      char tc = *h;
      *h++ = *t;
      *t-- = tc;
    }

    mp_int_clear(&tmp);
  }

  *str = '\0';
  if(cmp == 0)
    return MP_OK;
  else
    return MP_TRUNC;
}

/* }}} */

/* {{{ mp_int_string_len(z, radix) */

mp_result mp_int_string_len(mp_int z, mp_size radix)
{ 
  int  len;

  CHECK(z != NULL);

  if(radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
    return MP_RANGE;

  len = s_outlen(z, radix) + 1; /* for terminator */

  /* Allow for sign marker on negatives */
  if(MP_SIGN(z) == MP_NEG)
    len += 1;

  return len;
}

/* }}} */

/* {{{ mp_int_read_string(z, radix, *str) */

/* Read zero-terminated string into z */
mp_result mp_int_read_string(mp_int z, mp_size radix, const char *str)
{
  return mp_int_read_cstring(z, radix, str, NULL);

}

/* }}} */

/* {{{ mp_int_read_cstring(z, radix, *str, **end) */

mp_result mp_int_read_cstring(mp_int z, mp_size radix, const char *str, char **end)
{ 
  int       ch;

  CHECK(z != NULL && str != NULL);

  if(radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
    return MP_RANGE;

  /* Skip leading whitespace */
  while(isspace((unsigned char)*str))
    ++str;

  /* Handle leading sign tag (+/-, positive default) */
  switch(*str) {
  case '-':
    MP_SIGN(z) = MP_NEG;
    ++str;
    break;
  case '+':
    ++str; /* fallthrough */
  default:
    MP_SIGN(z) = MP_ZPOS;
    break;
  }

  /* Skip leading zeroes */
  while((ch = s_ch2val(*str, radix)) == 0) 
    ++str;

  /* Make sure there is enough space for the value */
  if(!s_pad(z, s_inlen(strlen(str), radix)))
    return MP_MEMORY;

  MP_USED(z) = 1; z->digits[0] = 0;

  while(*str != '\0' && ((ch = s_ch2val(*str, radix)) >= 0)) {
    s_dmul(z, (mp_digit)radix);
    s_dadd(z, (mp_digit)ch);
    ++str;
  }
  
  CLAMP(z);

  /* Override sign for zero, even if negative specified. */
  if(CMPZ(z) == 0)
    MP_SIGN(z) = MP_ZPOS;
  
  if(end != NULL)
    *end = (char *)str;

  /* Return a truncation error if the string has unprocessed
     characters remaining, so the caller can tell if the whole string
     was done */
  if(*str != '\0') 
    return MP_TRUNC;
  else
    return MP_OK;
}

/* }}} */

/* {{{ mp_int_count_bits(z) */

mp_result mp_int_count_bits(mp_int z)
{
  mp_size  nbits = 0, uz;
  mp_digit d;

  CHECK(z != NULL);

  uz = MP_USED(z);
  if(uz == 1 && z->digits[0] == 0)
    return 1;