---

mul.c

解説を見る。

#include "multiprec.h"
#include <string.h>

inline size_t
ymp_mulabs(digit_t *result, const digit_t *lhs, const digit_t *rhs,
           size_t lhs_len, size_t rhs_len)
{
  double_digit_t tmp=0;
  size_t i=0, j=0;

  if (lhs_len == 0 || rhs_len == 0) return 0;
  if ((lhs_len==1 && lhs[0]==0) || (rhs_len==1 && rhs[0]==0)) return 0;

  memset(result, 0, (lhs_len+rhs_len)*sizeof(digit_t));
  for (i=0; i<lhs_len; ++i)
    {
      tmp = 0;
      for (j=0; j<rhs_len; ++j)
        {
          tmp += ((double_digit_t)lhs[i]) * rhs[j];
          tmp += result[i+j];

          result[i+j] = LOW_DIGIT(tmp);
          tmp = HIGH_DIGIT(tmp);
        }
      result[i+j] = LOW_DIGIT(tmp);
    }
  if (tmp) ++i;
  return i+j-1;
}





inline size_t
ymp_mulabs_digit(digit_t *result, 
                 const digit_t *lhs, digit_t rhs, size_t lhs_len)
{
  double_digit_t tmp = 0, tmp2, tmp3;
  size_t i;
  for (i=0; i<lhs_len; ++i)
    {
      tmp2 = ((double_digit_t)lhs[i]) * rhs;
      tmp3 = LOW_DIGIT(tmp);
      tmp3 += LOW_DIGIT(tmp2);

      result[i] = LOW_DIGIT(tmp3);

      tmp = HIGH_DIGIT(tmp);
      tmp += HIGH_DIGIT(tmp2);
      tmp += HIGH_DIGIT(tmp3);
    }
  if (LOW_DIGIT(tmp))
    {
      result[i] = LOW_DIGIT(tmp);
      ++i;
    }
  return i;
}


inline size_t
ymp_modmulabs(digit_t *result, 
              const digit_t *lhs, const digit_t *rhs, const digit_t *modulus,
              size_t lhs_len, size_t rhs_len, size_t modulus_len)
{
  digit_t *product_digits;
  size_t product_used;

  YMP_TEMP_ALLOCATE(product_digits, digit_t, lhs_len+rhs_len);

  product_used = ymp_mulabs(product_digits, lhs, rhs, lhs_len, rhs_len);

  if (product_used < modulus_len)
    {
      memcpy(result, product_digits, sizeof(digit_t)*product_used);
      return product_used;
    }
  else
    {
      return ymp_divmodabs(NULL, result, product_digits, modulus,
                           product_used, modulus_len);
    }
}



inline size_t
ymp_modmulabs_2exp(digit_t *result, const digit_t *lhs, const digit_t *rhs, 
                   size_t pow, size_t lhs_len, size_t rhs_len)
{
  size_t n_digits, n_bits;

  size_t max_len;
  digit_t *product;

  double_digit_t tmp;
  size_t i, j;

  if (lhs_len == 0 || rhs_len == 0) return 0;
  if ((lhs_len==1 && lhs[0]==0) || (rhs_len==1 && rhs[0]==0)) return 0;


  n_digits = pow / DIGIT_BIT; n_bits = pow % DIGIT_BIT;
  max_len = n_digits + (n_bits != 0);
  
  if (result == lhs || result == rhs)
    {
      YMP_TEMP_ALLOCATE(product, digit_t, max_len);
    }
  else
    {
      product = result;
    }

  memset(product, 0, max_len*sizeof(digit_t));

  if (lhs_len > max_len) lhs_len = max_len;
  for (i=0; i<lhs_len; ++i)
    {
      size_t bound = max_len - i;
      tmp = 0;
      for (j=0; j<rhs_len && j<bound ; ++j)
        {
          tmp += ((double_digit_t)lhs[i]) * rhs[j];
          tmp += product[i+j];

          product[i+j] = LOW_DIGIT(tmp);
          tmp = HIGH_DIGIT(tmp);
        }
      if (i+j < max_len) product[i+j] = LOW_DIGIT(tmp);
    }
  
  if (n_bits) product[n_digits] &= (1u << n_bits)-1;

  while (max_len-- > 0)
    {
      if (product[max_len]) break;
    }
  ++max_len;

  if (product != result) memcpy(result, product, max_len * sizeof(digit_t));

  return max_len;
}


void
ymp_mul(mp_ref_t result, mp_cref_t self, mp_cref_t other)
{
  digit_t *self_digits = self->digits;
  digit_t *other_digits = other->digits;
  if (result == self)
    {
      YMP_TEMP_ALLOCATE(self_digits, digit_t, self->used);
      memcpy(self_digits, self->digits, sizeof(digit_t)*self->used);
    }
  if (result == other)
    {
      YMP_TEMP_ALLOCATE(other_digits, digit_t, self->used);
      memcpy(other_digits, other->digits, sizeof(digit_t)*self->used);
    }

  ymp_reserve(result, self->used+other->used);
  result->used = ymp_mulabs(result->digits, self_digits, other_digits,
                            self->used, other->used);

  result->sign = self->sign != other->sign;
}



void
ymp_modmul(mp_ref_t result, mp_cref_t self, mp_cref_t other, mp_cref_t modulus)
{
  digit_t *self_digits = self->digits;
  digit_t *other_digits = other->digits;

  if (result == self)
    {
      YMP_TEMP_ALLOCATE(self_digits, digit_t, self->used);
      memcpy(self_digits, self->digits, sizeof(digit_t)*self->used);
    }
  if (result == other)
    {
      YMP_TEMP_ALLOCATE(other_digits, digit_t, self->used);
      memcpy(other_digits, other->digits, sizeof(digit_t)*self->used);
    }

  ymp_reserve(result, modulus->used);

  result->used = ymp_modmulabs(result->digits, 
                               self_digits, other->digits, modulus->digits,
                               self->used, other->used, modulus->used);

  result->sign = self->sign != other->sign;
}


void
ymp_modmul_2exp(mp_ref_t result, mp_cref_t self, mp_cref_t other, size_t pow)
{
  ymp_reserve(result, pow/DIGIT_BIT + 1);

  result->used = ymp_modmulabs_2exp(result->digits, 
                                    self->digits, other->digits,
                                    pow, self->used, other->used);

  result->sign = self->sign != other->sign;
}


void
ymp_mul_digit(mp_ref_t result, mp_cref_t self, digit_t other)
{
  ymp_reserve(result, self->used+1);
  result->used = ymp_mulabs_digit(result->digits, 
                                  self->digits, other, self->used);
  result->sign = self->sign;
}


size_t
ymp_modpowabs_z(digit_t *result, const digit_t *source, const digit_t *modulus,
                size_t source_len, size_t modulus_len, size_t pow)
{
  /* 繰りかえし2乗法による */

  digit_t *powered_digits; /* self**(2の冪)を保持 */
  size_t powered_used;     /* powered_digitsの先に格納された有効要素数 */
  size_t result_used;

  /* result == modulus の場合、複製 */
  if (result == modulus)
    {
      digit_t *new_modulus;
      YMP_TEMP_ALLOCATE(new_modulus, digit_t, modulus_len);
      memcpy(new_modulus, modulus, sizeof(digit_t)*modulus_len);
      modulus = new_modulus;
    }

  /* <powered> = <self> ** (2**0) */
  YMP_TEMP_ALLOCATE(powered_digits, digit_t, modulus_len*2);
  if (source_len < modulus_len)
    {
      powered_used = source_len;
      memcpy(powered_digits, source, source_len*sizeof(digit_t));
    }
  else
    {
      powered_used = ymp_divmodabs(NULL, powered_digits,
                                   source, modulus, source_len, modulus_len);
                               
    }

  /* <result> = 1 */
  result_used = 1; result[0] = 1;

  for (;;)
    {
      if (pow & 0x01) 
        {
          result_used = ymp_modmulabs(result, 
                                      result, powered_digits, modulus,
                                      result_used, powered_used, modulus_len);
        }

      pow >>= 1;
      if (!pow) return result_used;

      powered_used = ymp_modmulabs(powered_digits, 
                                   powered_digits, powered_digits, modulus,
                                   powered_used, powered_used, modulus_len);
    }
}


size_t
ymp_modpowabs_2exp_z(digit_t *result, const digit_t *source,
                     size_t len, size_t emod, size_t pow)
{
  /* 繰りかえし2乗法による */

  digit_t *powered_digits; /* self**(2の冪)を保持 */
  size_t powered_used;     /* powered_digitsの先に格納された有効要素数 */
  size_t result_used;      /* resultに格納されている要素数 */
  size_t mod_len = emod / DIGIT_BIT + 1;

  /* <powered> = <self> ** (2**0) */
  YMP_TEMP_ALLOCATE(powered_digits, digit_t, mod_len);
  powered_used = ymp_modabs_2exp(powered_digits, source, len, emod);

  /* <result> = 1 */
  result_used = 1; result[0] = 1;

  for (;;)
    {
      if (pow & 0x01) 
        {
          result_used = ymp_modmulabs_2exp(result, result, powered_digits,
                                           emod, result_used, powered_used);
        }

      pow >>= 1;
      if (!pow) return result_used;

      powered_used = ymp_modmulabs_2exp(powered_digits, 
                                        powered_digits, powered_digits,
                                        emod, powered_used, powered_used);
    }
}



inline size_t
ymp_modpowabs(digit_t *result, 
              const digit_t *source, const digit_t *modulus, 
              const digit_t *pow,
              size_t source_len, size_t modulus_len, size_t pow_len)
{
  /* 繰りかえし2乗法による */

  size_t i;                /* pow->digits[]内のindex */
  digit_t bit_mask;        /* pow->digits[i]内のbitテスト用 */

  digit_t *powered_digits; /* self**(2の冪)を保持 */
  size_t powered_used;     /* powered_digitsの先に格納された有効要素数 */

  size_t result_used;      /* resultの先に格納された有効要素数 */

  /* result == modulus の場合、複製 */
  if (result == modulus)
    {
      digit_t *new_modulus;
      YMP_TEMP_ALLOCATE(new_modulus, digit_t, modulus_len);
      memcpy(new_modulus, modulus, sizeof(digit_t)*modulus_len);
      modulus = new_modulus;
    }

  /* result == pow の場合、複製 */
  if (result == pow)
    {
      digit_t *new_pow;
      YMP_TEMP_ALLOCATE(new_pow, digit_t, pow_len);
      memcpy(new_pow, pow, sizeof(digit_t)*pow_len);
      pow = new_pow;
    }

  /* <powered> = <self> ** (2**0) */
  YMP_TEMP_ALLOCATE(powered_digits, digit_t, modulus_len*2);
  if (source_len < modulus_len)
    {
      powered_used = source_len;
      memcpy(powered_digits, source, sizeof(digit_t)*source_len);
    }
  else
    {
      powered_used = ymp_divmodabs(NULL, powered_digits,
                                   source, modulus, source_len, modulus_len);
                               
    }

  /* <result> = 1 */
  result_used = 1; result[0] = 1;

  i=0; bit_mask=1;
  for (;;)
    {
      if (pow[i] & bit_mask) 
        {
          result_used = ymp_modmulabs(result, 
                                      result, powered_digits, modulus,
                                      result_used, powered_used, modulus_len);
        }

      if (bit_mask == DIGIT_HIGHEST_BIT)
        { 
          bit_mask = 1;
          ++i;
          if (i==pow_len) return result_used;
        }
      else
        {
          bit_mask <<= 1;
          if (i==pow_len-1 && pow[i] < bit_mask) return result_used;
        }

      powered_used = ymp_modmulabs(powered_digits, 
                                   powered_digits, powered_digits, modulus,
                                   powered_used, powered_used, modulus_len);
    }
}



inline size_t
ymp_powabs(digit_t *result, const digit_t *source, 
           size_t source_len, size_t pow) 
{
  /* 繰りかえし2乗法による */
  digit_t *tmp_digits;     /* 一時的にdigit_t列を待避する為の領域をポイント */
  digit_t *powered_digits; /* self**(2の冪)を保持 */
  size_t powered_used;     /* powered_digitsの先に格納された有効要素数 */
  size_t result_used;

  if (pow == 0)
    {
      result[0] = 1;
      return 1;
    }
  /* 以下、<result> >= <source>を仮定できる */


  /* 一時領域を確保 */
  YMP_TEMP_ALLOCATE(tmp_digits, digit_t, source_len*pow);

  /* <powered> = <self> ** (2**0) */
  YMP_TEMP_ALLOCATE(powered_digits, digit_t, source_len*pow);
  powered_used = source_len;
  memcpy(powered_digits, source, sizeof(digit_t)*source_len);

  /* <result> = 1 */
  result_used = 1; result[0] = 1;



  /* 繰りかえし2乗法 */
  for (;;)
    {
      if (pow & 0x01) 
        {
          memcpy(tmp_digits, result, sizeof(digit_t)*result_used);
          result_used = ymp_mulabs(result, tmp_digits, powered_digits,
                                    result_used, powered_used);
        }

      pow >>= 1;
      if (!pow) return result_used;

      memcpy(tmp_digits, powered_digits, sizeof(digit_t)*powered_used);
      powered_used = ymp_mulabs(powered_digits, tmp_digits, tmp_digits,
                                powered_used, powered_used);
    }
}
void
ymp_pow(mp_ref_t result, mp_cref_t self, size_t pow)
{
  /* 結果の格納領域を確保 */
  ymp_reserve(result, self->used*MAX(pow, 1));

  result->sign = self->sign && (pow & 0x01);
  result->used = ymp_powabs(result->digits, self->digits, self->used, pow);
}

void
ymp_modpow_z(mp_ref_t result, mp_cref_t self, mp_cref_t modulus, size_t pow)
{
  /* 十分な格納領域を確保 */
  ymp_reserve(result, modulus->used);

  result->sign = self->sign && (pow & 0x01);
  result->used = ymp_modpowabs_z(result->digits, self->digits, modulus->digits,
                                 self->used, modulus->used, pow);
}


void
ymp_modpow_2exp_z(mp_ref_t result, mp_cref_t self, size_t emod, size_t pow)
{
  /* 十分な格納領域を確保 */
  ymp_reserve(result, emod / DIGIT_BIT + 1);

  result->sign = self->sign && (pow & 0x01);
  result->used = ymp_modpowabs_2exp_z(result->digits, self->digits,
                                      self->used, emod, pow);
}


void
ymp_modpow(mp_ref_t result, mp_cref_t self, mp_cref_t modulus, mp_cref_t pow)
{
  /* 十分な格納領域を確保 */
  ymp_reserve(result, modulus->used);

  result->sign = self->sign && (pow->used != 0 || (pow->digits[0] & 0x01));
  result->used = ymp_modpowabs(result->digits, 
                               self->digits, modulus->digits, pow->digits,
                               self->used, modulus->used, pow->used);
}

---