---

div.c

解説を見る。

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

digit_t
ymp_divmodabs_digit(digit_t *quotient, const digit_t *lhs, const digit_t rhs,
                    size_t lhs_len)
{
  double_digit_t tmp = 0;
  register size_t i = lhs_len;
  
  while (i-- > 0)
    {
      tmp = LOW_DIGIT_TO_HIGH_DIGIT(tmp);
      tmp += lhs[i];
      if (quotient) quotient[i] = tmp/rhs;
      tmp %= rhs;
    }
  return LOW_DIGIT(tmp);
}




inline static void
do_divmodabs(digit_t *quotient, digit_t *lhs, const digit_t *rhs, 
             size_t lhs_len, size_t rhs_len)

{
  double_digit_t lhs_top;
  size_t i;

  i = lhs_len; lhs_top = lhs[--i];
  while (i-- >= rhs_len)
    {
      size_t j;
      digit_t q;
      double_digit_t tmp;

      /* 最上位桁による仮商を定める */
      lhs_top = LOW_DIGIT_TO_HIGH_DIGIT(lhs_top);
      lhs_top += lhs[i];
      q = lhs_top / rhs[rhs_len-1];

      /* (被除数)*(仮商) を被除数から引いてみる*/
      tmp = 0;
      for (j=0; j<rhs_len; ++j)
        {
          double_digit_t tmp2 = MULTIPREC_RADIX;
          tmp += ((double_digit_t)rhs[j]) * q;

          tmp2 += lhs[i-rhs_len+j+1];
          tmp2 -= LOW_DIGIT(tmp);
          lhs[i-rhs_len+j+1] = LOW_DIGIT(tmp2);

          tmp = HIGH_DIGIT(tmp);
          tmp += !HIGH_DIGIT(tmp2);
        }
      tmp = MULTIPREC_RADIX - tmp;
      tmp += lhs[i+1];
      lhs[i+1] = LOW_DIGIT(tmp);

      /* 引き過ぎていたら、最上位桁の最上位bitは1の筈 */
      while (lhs[i+1] & DIGIT_HIGHEST_BIT) 
        {
          /* 商を減らし */
          --q;

          /* その分だけ足し戻す */
          tmp = 0;
          for (j=0; j<rhs_len; ++j)
            {
              tmp += lhs[i-rhs_len+j+1];
              tmp += rhs[j];
              lhs[i-rhs_len+j+1] = LOW_DIGIT(tmp);
              tmp = HIGH_DIGIT(tmp);
            }
          lhs[i+1] += LOW_DIGIT(tmp);
        }

      /* 得られた商を書き込む */
      if (quotient) quotient[i+1-rhs_len] = q;
      
      lhs_top = lhs[i];
    }
}



size_t
ymp_divmodabs(digit_t *quotient, digit_t *remainder,
              const digit_t *lhs, const digit_t *rhs,
              size_t lhs_len, size_t rhs_len)
{
  int shift_width; /* シフト幅 */
  digit_t hd;      /* rhsの最上位桁 */
  size_t i;
  digit_t *new_lhs, *new_rhs; /* 被除数、除数を各々シフトしたもの */
  double_digit_t tmp;

  /* 0除算を弾く */
  if (rhs_len==0 || (rhs_len==1 && rhs[0]==0))
    {
      ymp_fatal("0による除算");
      return 0;
    }

  /* シフトした被除数、除数を格納するための領域を確保 */
  YMP_TEMP_ALLOCATE(new_lhs, digit_t, lhs_len+1);
  YMP_TEMP_ALLOCATE(new_rhs, digit_t, rhs_len);

  /* 適切なシフト幅を算出 */
  shift_width = 0;
  hd = rhs[rhs_len-1];
  while ( !(hd & DIGIT_HIGHEST_BIT) )
    {
      ++shift_width;
      hd <<= 1;
    }

  /* 除数をshift_widthだけシフト */
  i = rhs_len; tmp = rhs[--i];
  while (i-- > 0)
    {
      tmp = LOW_DIGIT_TO_HIGH_DIGIT(tmp);
      tmp += rhs[i];
      tmp <<= shift_width;
      new_rhs[i+1] = HIGH_DIGIT(tmp);
      tmp = rhs[i];
    }
  new_rhs[0] = rhs[0] << shift_width;

  /* 被除数をshift_widthだけシフト */
  tmp = 0; i=lhs_len;
  while(i-- > 0)
    {
      tmp += lhs[i];
      tmp <<= shift_width;
      new_lhs[i+1] = HIGH_DIGIT(tmp);

      tmp = lhs[i];
      tmp = LOW_DIGIT_TO_HIGH_DIGIT(tmp);
    }
  new_lhs[0] = lhs[0] << shift_width;

  /* シフトした状態で、除算を実行 */
  do_divmodabs(quotient, new_lhs, new_rhs, lhs_len+1, rhs_len);

  /* 今、シフトされた剰余がnew_lhsに入っている */

  if (remainder)
    {
      /* 剰余を逆シフトして復元しながら、非0の最高桁数を探す */
      size_t rem_used = 0;
      tmp = 0; i=rhs_len;
      while (i-- > 0)
        {
          tmp += new_lhs[i];
          remainder[i] = LOW_DIGIT(tmp >> shift_width);

          if (remainder[i] && !rem_used) rem_used = i+1;

          tmp = new_lhs[i];
          tmp = LOW_DIGIT_TO_HIGH_DIGIT(tmp);
        }
      return rem_used;
    }
  else
    {
      return 0;
    }
}
void
ymp_divmod(mp_ref_t quotient, mp_ref_t remainder, 
           mp_cref_t self, mp_cref_t other)
{
  enum mp_sign_t rem_sign;
  size_t rem_used;

  /* 0除算を弾く */
  if (ymp_is_zero(other))
    {
      ymp_fatal("0による除算");
      return;
    }

  /* 明らかに (被除数) < (除数) と分かる場合:
     商は0、剰余は被除数に等しい */
  if (self->used < other->used)
    {
      if (quotient) quotient->used = 0;
      if (self != remainder && remainder) ymp_assign(remainder, self);
      return;
    }

  /* 商・剰余の格納に十分な領域を確保 */
  if (quotient) ymp_reserve(quotient, self->used-other->used+1);
  if (remainder) ymp_reserve(remainder, other->used);

  /* 絶対値同士の除算を行う */
  rem_sign = self->sign;
  rem_used = ymp_divmodabs(quotient?quotient->digits:NULL, 
                           remainder?remainder->digits:NULL, 
                           self->digits, other->digits,
                           self->used, other->used);
  
  /*
   *  remainder==selfまたはremainder==otherのときにも
   *  商のの符号と桁数の調整に影響を与えないように、
   *  剰余の符号と桁数は直ぐには書き込まない
   */

  if (quotient)
    {
      /* 商の符号と桁数の調整 */
      quotient->sign = self->sign != other->sign;
      quotient->used = self->used - other->used + 1;
      if (!quotient->digits[quotient->used-1]) --quotient->used;
    }

  if (remainder)
    {
      /* 剰余の符号と桁数を書き込む */
      remainder->sign = rem_sign;
      remainder->used = rem_used;
    }
}


digit_t
ymp_divmod_digit(mp_ref_t quotient, mp_cref_t self, digit_t other)
{
  digit_t rem; /* 剰余の絶対値 */

  /* 0除算を弾く */
  if (other == 0)
    {
      ymp_fatal("0による除算");
      return 0;
    }

  /* 被除数が0の場合、商と剰余も0 */
  if (self->used == 0)
    {
      if (quotient) quotient->used = 0;
      return 0;
    }

  /* 商を格納するのに十分な領域を確保 */
  if (quotient) ymp_reserve(quotient, self->used);

  /* 絶対値の除算を行う */
  rem = ymp_divmodabs_digit(quotient?quotient->digits:NULL,
                            self->digits,
                            other,
                            self->used);

  /* 桁数と符号の設定 */
  if (quotient)
    {
      if (quotient->digits[quotient->used-1]) --quotient->used;
      quotient->sign = self->sign;
    }

  if (self->sign && rem)
    {
      /* 剰余が正になるように調整 */
      if (quotient)
        {
          quotient->used = ymp_addabs(quotient->digits, quotient->digits, 
                                      &other, quotient->len, 1);
        }
      return other - rem;
    }
  else
    {
      return rem;
    }
}

int
ymp_div_if_divisible(mp_ref_t quotient, mp_cref_t self, mp_cref_t other)
{
  size_t rem_used; /* 剰余の桁数 */
  digit_t *quot_buf=NULL, *rem_buf=NULL; /* 商・剰余の一時格納用 */

  /* 0除算を弾く */
  if (ymp_is_zero(other))
    {
      ymp_fatal("0による除算");
      return 2;
    }

  /* 明らかな場合 */
  if (self->used < other->used)
    {
      if (ymp_is_zero(self))
        {
          if (quotient) quotient->used = 0;
          return 1;
        }
      else
        {
          return 0;
        }
    }

  /* 格納領域を確保 */
  if (quotient) YMP_TEMP_ALLOCATE(quot_buf, digit_t, self->used-other->used+1);
  YMP_TEMP_ALLOCATE(rem_buf, digit_t, other->used);

  /* 絶対値の除算を実行 */
  rem_used = ymp_divmodabs(quot_buf, rem_buf, self->digits, other->digits,
                           self->used, other->used);

  /* 剰余が0なら (そのときrem_usedは必ず0) */
  if (rem_used==0)
    {
      if (quotient)
        {
          /* *quotientに商のデータをコピー */
          size_t quot_used;
          quotient->sign = self->sign != other->sign;
          
          quot_used = self->used - other->used + 1;
          if (!quot_buf[quot_used-1]) --quot_used;
          quotient->used = quot_used;
          ymp_reserve(quotient, quot_used);
          memcpy(quotient->digits, quot_buf, quot_used);
        }
      return 1;
    }
  else
    {
      return 0;
    }
}



int
ymp_div_if_divisible_digit(mp_ref_t quotient, mp_cref_t self, digit_t other)
{
  digit_t rem;
  digit_t *quot_buf=NULL;

  /* 0除算を弾く */
  if (other == 0)
    {
      ymp_fatal("0による除算");
      return 2;
    }

  /* 明らかな場合 */
  if (ymp_is_zero(self))
    {
      if (quotient) quotient->used = 0;
      return 1;
    }

  /* 格納領域を確保 */
  if (quotient) YMP_TEMP_ALLOCATE(quot_buf, digit_t, self->used);

  /* 絶対値の除算を実行 */
  rem = ymp_divmodabs_digit(quot_buf, self->digits, other, self->used);

  /* 剰余が0なら */
  if (rem == 0)
    {
      if (quotient)
        {
          /* *quotientに商のデータをコピー */
          size_t quot_used;
          quotient->sign = self->sign;
          
          quot_used = self->used;
          if (!quot_buf[quot_used-1]) --quot_used;
          quotient->used = quot_used;
          ymp_reserve(quotient, quot_used);
          memcpy(quotient->digits, quot_buf, quot_used);
        }
      return 1;
    }
  else
    {
      return 0;
    }
}

---