---

root.c

解説を見る。

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

void
ymp_iroot(mp_ref_t result, mp_cref_t self, size_t n)
{
  /* Newton法による */

  size_t i;
  size_t ilog2, bit_len;
  struct multiprec x, x_prev; /* 近似途中の値を保持 */
  struct multiprec tmp;       /* 計算上の一時領域   */
  struct multiprec self_copy;

  if (!self->sign)
    { /* 非負ならば良い */
      result->sign = positive_sign;
    }
  else
    { /* そうでないとき */
      /* 偶数乗根はエラー */
      if ( (n & 0x01) == 0) ymp_fatal("負数の偶数乗根");

      /* 結果の符号が確定 */
      result->sign = self->sign;

      /* 以下、-<self> >= 0 の絶対値を算出すれば良い */
      memcpy(&self_copy, self, sizeof(struct multiprec));
      self_copy.sign = positive_sign;
      self = &self_copy;
    }

  /* bit長により、初期値を大まかに決定 */
  /* (初期値) >= (収束値) を保証して収束が遅くなるのを防止する */
  /* ilog2 := { 2**k > |<self>|なる最小のk } */
  ilog2 = ymp_ilog2(self); ++ilog2;

  bit_len = ilog2 / n;
  if (bit_len == 0)
    {
      if (ymp_is_zero(self)) ymp_assign_abs_digit(result, 0);
      else                   ymp_assign_abs_digit(result, 1);
      return;
    }
  if (ilog2 % n) ++bit_len;

  ymp_alloca_initialize(&x, self->used * n + 1);
  for (i=0; i<bit_len/DIGIT_BIT; ++i) x.digits[i] = 0;  
  x.digits[i] = 1u << (bit_len%DIGIT_BIT);
  x.used = ++i;

  /* 今、x は (2**k)**n > |<self>| なる最小の2**k */


  /* 作業領域確保 */
  ymp_alloca_initialize(&tmp, self->used * n + 1);
  ymp_alloca_initialize(&x_prev, self->used * n + 1);

  
  /* f(x) = x**n - |<self>| に対するNetwon法 */
  for (;;)
    {
    l_newton_iteration:
      ymp_assign(&x_prev, &x);
      
      /* x_{i+1} = ((n-1) * x_i**n + a) / (n * x_i**(n-1)) */
      ymp_pow(&tmp, &x_prev, n-1);  /* tmp = x_i ** (n-1)  */

      ymp_assign(&x, &tmp);         /* x = x_i ** n */
      ymp_mul(&x, &x, &x_prev);
      ymp_mul_digit(&x, &x, n-1);       /* x = (n-1) * x_i**n  */
      ymp_add(&x, &x, self);        /* (n-1) * x_i**n + a  */
      
      ymp_mul_digit(&tmp, &tmp, n); /* tmp = n*x_i**(n-1)  */
      
      ymp_div(&x, &x, &tmp);        /* x = x_{i+1} */


      if (x.used != x_prev.used) continue;  /* まだ誤差が大きい */

      for (i=1; i<x.used; ++i)
        {
          if (x.digits[i] != x_prev.digits[i]) /* まだ誤差が大きい */
            goto l_newton_iteration;           /* Newton法続行 */
        }

      /* かなり収束した。現在の誤差の絶対値が1以下なら、あとは微調整 */
      switch ((int)(x.digits[0] - x_prev.digits[0]))
        {
        case 0:
          ymp_assign_abs(result, &x);
          return;

        case 1:
          for (;;)
            {
              ymp_pow(&tmp, &x, n);
              if (ymp_compare_abs(&tmp, self) > 0) break;
              ymp_add_digit(&x, &x, 1);
            }
          ymp_sub_digit(&x, &x, 1);
          ymp_assign_abs(result, &x);
          return;
          
        case -1:
          for (;;)
            {
              ymp_pow(&tmp, &x, n);
              if (ymp_compare_abs(&tmp, self) <= 0) break;
              ymp_sub_digit(&x, &x, 1);
            }
          ymp_assign_abs(result, &x);
          return;
        }
    }
}

---