---

root.c の解説

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

ソースコードを見る。

関数
void	ymp_iroot (mp_ref_t result, mp_cref_t self, size_t n)
	多倍長整数の冪乗根の整数部分. より詳しく...

---

関数の解説

void ymp_iroot (  mp_ref_t    result, mp_cref_t    self, size_t    n )

多倍長整数の冪乗根の整数部分. result := floor( <self> ** (1/n) ) 引数: result  結果を格納するstruct multiprecへのポインタ self  元の数を表すstruct multiprecへのポインタ n  冪指数 事前条件: nが偶数ならば、selfは非負であること。 そうでないときはymp_fatalが呼ばれる。 nは0でないこと. そうでないとき、結果は不定 root.c の 16 行で定義されています。 参照 DIGIT_BIT, multiprec::digits, mp_cref_t, positive_sign, multiprec::sign, multiprec::used, ymp_add, ymp_add_digit, ymp_alloca_initialize, ymp_assign, ymp_assign_abs, ymp_assign_abs_digit, ymp_compare_abs, ymp_div, ymp_fatal, ymp_ilog2, ymp_is_zero, ymp_mul, ymp_mul_digit, ymp_pow, と ymp_sub_digit. { /* 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; } } }

---