---

数論的関数
[多倍長整数]

モジュール
下位関数群
関数
void	ymp_euclid (mp_ref_t g, mp_cref_t a, mp_cref_t b)
	多倍長整数の最小公倍数を計算. より詳しく...
void	ymp_euclid_ex (mp_ref_t g, mp_ref_t s, mp_ref_t t, mp_cref_t a, mp_cref_t b)
	拡張Euclid互除法を行なう. より詳しく...
int	ymp_modinv (mp_ref_t result, mp_cref_t self, mp_cref_t modulus)
	剰余環Z/rZ上の逆元. より詳しく...
int	ymp_binary_sunzi (mp_ref_t result, const mp_cref_t remainders[2], const mp_cref_t moduli[2])
	孫子(中国)の剰余定理の解を求める. より詳しく...
void	ymp_multi_euclid (mp_ref_t result, size_t len, const struct multiprec *ary)
int	ymp_proot (mp_ref_t result, mp_cref_t self, size_t pow)
	Perfect powerの冪根を求める. より詳しく...
int	ymp_is_perfect_power (mp_cref_t self)
	Perfect powerかどうかを判定. より詳しく...

---

関数の解説

void ymp_euclid (  mp_ref_t    g, mp_cref_t    a, mp_cref_t    b )

多倍長整数の最小公倍数を計算. Euclid互除法による。 引数: g  結果を格納されるstruct multiprecへのポインタ a  第１項を表すstruct multiprecへのポインタ b  第２項を表すstruct multiprecへのポインタ 事後条件: 求められた最小公倍数は非負 euclid.c の 22 行で定義されています。 参照 digit_t, multiprec::digits, mp_cref_t, multiprec::used, ymp_assign_array, ymp_divmodabs, と YMP_TEMP_ALLOCATE. { struct { size_t used; digit_t *digits; } r0, r1; r0.used = a->used; r1.used = b->used; YMP_TEMP_ALLOCATE(r0.digits, digit_t, r0.used); YMP_TEMP_ALLOCATE(r1.digits, digit_t, r1.used); memcpy(r0.digits, a->digits, r0.used*sizeof(digit_t)); memcpy(r1.digits, b->digits, r1.used*sizeof(digit_t)); while (r1.used !=0 && (r1.used != 1 || r1.digits[0] != 0)) { register size_t tmp_used; register digit_t *tmp_digits; if (r0.used >= r1.used) { r0.used = ymp_divmodabs(NULL, r0.digits, r0.digits, r1.digits, r0.used, r1.used); } tmp_used = r0.used; r0.used = r1.used; r1.used = tmp_used; tmp_digits = r0.digits; r0.digits = r1.digits; r1.digits = tmp_digits; } ymp_assign_array(g, r0.used, r0.digits); }

void ymp_euclid_ex (  mp_ref_t    g, mp_ref_t    s, mp_ref_t    t, mp_cref_t    a, mp_cref_t    b )

拡張Euclid互除法を行なう. gcd(a,b) = g = as + btなるg, s, tを算出する. 引数: g  結果を格納されるstruct multiprecへのポインタ 不用なときはNULLを指定する s  結果を格納されるstruct multiprecへのポインタ 不用なときはNULLを指定する t  結果を格納されるstruct multiprecへのポインタ 不用なときはNULLを指定する a  第１項を表すstruct multiprecへのポインタ b  第２項を表すstruct multiprecへのポインタ 覚え書き: g, s, tの符号は不定 勿論、(s,t)=1 euclid.c の 68 行で定義されています。 参照 MAX, mp_cref_t, multiprec::used, ymp_alloca_initialize, ymp_assign, ymp_destroy, ymp_divmod, ymp_initialize, ymp_initialize_by_digit, ymp_is_zero, ymp_mul, ymp_sub, と ymp_swap. { struct multiprec r0, r1, s0, s1, t0, t1, q, tmp; ymp_alloca_initialize(&r0, a->used); ymp_alloca_initialize(&r1, b->used); ymp_alloca_initialize(&q, MAX(a->used, b->used)); ymp_assign(&r0, a); ymp_assign(&r1, b); ymp_initialize_by_digit(&s0, 1); ymp_initialize_by_digit(&s1, 0); ymp_initialize_by_digit(&t0, 0); ymp_initialize_by_digit(&t1, 1); ymp_initialize(&tmp); while (!ymp_is_zero(&r1)) { ymp_divmod(&q, &r0, &r0, &r1); ymp_mul(&tmp, &q, &s1); ymp_sub(&s0, &s0, &tmp); ymp_mul(&tmp, &q, &t1); ymp_sub(&t0, &t0, &tmp); ymp_swap(&r0, &r1); ymp_swap(&s0, &s1); ymp_swap(&t0, &t1); } if (g) ymp_assign(g, &r0); if (s) ymp_assign(s, &s0); if (t) ymp_assign(t, &t0); ymp_destroy(&s0); ymp_destroy(&s1); ymp_destroy(&t0); ymp_destroy(&t1); ymp_destroy(&tmp); }

int ymp_modinv (  mp_ref_t    result, mp_cref_t    self, mp_cref_t    modulus )

剰余環Z/rZ上の逆元. selfのmodulusを法とした乗法逆元を計算する 引数: result  結果を格納されるstruct multiprecへのポインタ self  元の数を表すstruct multiprecへのポインタ modulus  法を表すstruct multiprecへのポインタ 戻り値: 0  逆元が正しく計算された 1  gcd(self, modulus)が1でなかったために、正しく計算されなかった euclid.c の 114 行で定義されています。 参照 multiprec::digits, MAX, mp_cref_t, multiprec::sign, multiprec::used, ymp_alloca_initialize, ymp_assign, ymp_assign_digit, ymp_divmod, ymp_is_zero, ymp_mod, ymp_mul, ymp_sub, と ymp_swap. { struct multiprec r0, r1, s0, s1, q, tmp; register size_t len = MAX(self->used, modulus->used); ymp_alloca_initialize(&r0, self->used); /* |r0| <= |self| */ ymp_alloca_initialize(&r1, modulus->used); /* |r1| <= |modulus| */ ymp_alloca_initialize(&s0, modulus->used+len+1); /* |s0| <= |r0|*|q|+|s1| */ ymp_alloca_initialize(&s1, modulus->used+len+1); /* |s1| <= |r1|*|q|+|s0| */ ymp_alloca_initialize(&q, len); /* |q| <= max(|r0|,|r1|)*/ ymp_alloca_initialize(&tmp, modulus->used+len); /* |tmp| <= |r_i|*|q| */ ymp_assign(&r0, self); ymp_assign(&r1, modulus); ymp_assign_digit(&s0, 1); ymp_assign_digit(&s1, 0); while (!ymp_is_zero(&r1)) { ymp_divmod(&q, &r0, &r0, &r1); ymp_mul(&tmp, &s1, &q); ymp_sub(&s0, &s0, &tmp); ymp_mod(&s0, &s0, modulus); ymp_swap(&r0, &r1); ymp_swap(&s0, &s1); } if (r0.used == 1 && r0.digits[0] == 1) { ymp_assign(result, &s0); if (r0.sign) result->sign = !result->sign; return 0; } else { return 1; } }

int ymp_binary_sunzi (  mp_ref_t    result, const mp_cref_t    remainders[2], const mp_cref_t    moduli[2] )

孫子(中国)の剰余定理の解を求める. u = u_0 (mod m_0), u = u_1 (mod m_1) なる u を、 m_0 * m_1を法として求める 引数: result  結果uを格納するstruct multiprecへのポインタ remainders  u_0, u_1をそれぞれ表すstruct multiprecへのポインタの配列 moduli  m_0, m_1をそれぞれ表すstruct multiprecへのポインタの配列 事前条件: <m_0>と<m_1>は互いに素であること。そうでない場合、結果は正しくない 戻り値: 0  正しく計算された 1  <m_0>と<m_1>は互いに素でないので正しく計算されなかった euclid.c の 166 行で定義されています。 参照 MAX, mp_cref_t, ymp_add, ymp_alloca_initialize, ymp_mod, ymp_modinv, ymp_modmul, ymp_mul, と ymp_sub. { /* * H. L. Garnerの高速解法 * c_{0,1} := m_0^{-1} (mod m_1) * * v_0 := u_0 mod m_0 * v_1 := (u_1 - u_0) * c_{0,1} mod m_1 * * u = v_1 * m_0 + v_0 */ struct multiprec m_0_inv; /* m_0^{-1} (mod m_1) */ struct multiprec v[2]; register size_t len; ymp_alloca_initialize(&m_0_inv, moduli[1]->used); if (ymp_modinv(&m_0_inv, moduli[0], moduli[1])) return 1; ymp_alloca_initialize(&v[0], moduli[0]->used); ymp_mod(&v[0], remainders[0], moduli[0]); len = MAX(remainders[1]->used, moduli[0]->used + moduli[1]->used); ymp_alloca_initialize(&v[1], len); ymp_sub(&v[1], remainders[1], &v[0]); ymp_modmul(&v[1], &v[1], &m_0_inv, moduli[1]); ymp_mul(&v[1], &v[1], moduli[0]); ymp_add(result, &v[1], &v[0]); return 0; }

void ymp_multi_euclid (  mp_ref_t    result, size_t    len, const struct multiprec *    ary )

int ymp_proot (  mp_ref_t    result, mp_cref_t    self, size_t    pow )

Perfect powerの冪根を求める. <self>が完全pow乗数であればそのpow乗根を<result>に格納 アルゴリズムについては剰余による冪根の高速計算(dvi)を参照 引数: result  冪根を格納するstruct multiprecへのポインタ self  対象の整数を表すstruct multiprecへのポインタ pow  冪指数 戻り値: <self>が完全pow乗数であれば非0。そうでなければ0 事後条件: <self>がpow乗数でなければ<result>は変化しない。 perfect_power.c の 254 行で定義されています。 { struct multiprec tmp; size_t n_square; /* powの素因数分解の中の2の個数 */ size_t e2_digits, e2_bits; /* 下位の0の並んだ部分の、digit数, bit数 */ register digit_t ld; /* 正数以外は除く */ if (self->sign || ymp_is_zero(self)) return 0; /* e2 := {2**e が <self> を割る最大の e} */ for (e2_digits=0; e2_digits<self->used; ++e2_digits) { if (self->digits[e2_digits] != 0) break; } e2_bits = 0; for (ld = self->digits[e2_digits]; (ld & 0x1) == 0; ld>>=1) { ++e2_bits; } /* 素因数分解の中の2の冪指数をチェック */ if ((e2_digits * DIGIT_BIT + e2_bits) % pow != 0) return 0; /* 以下、2の冪は考えなくて良い */ n_square = 0; while ( (pow & 0x01) == 0 ) { pow >>= 1; ++n_square; } ymp_alloca_initialize(&tmp, self->used - e2_digits); ymp_div_2exp(&tmp, self, e2_digits * DIGIT_BIT + e2_bits); if (!ymp_odd_root(&tmp, &tmp, pow)) { return 0; } /* 以下、冪指数のうち2冪の部分のみ調べる */ while (n_square-- > 0) { if (!ymp_odd_sqrt(&tmp, &tmp)) return 0; } ymp_assign(result, &tmp); return 1; }

int ymp_is_perfect_power (  mp_cref_t    self )

Perfect powerかどうかを判定. <self> = n**e なるn, e in N が存在するかどうかを判定。 x アルゴリズムについては剰余による冪根の高速計算(dvi)を参照 引数: self  対象の整数を表すstruct multiprecへのポインタ 戻り値: <self>がperfect powerであれば非0。そうでなければ0 perfect_power.c の 321 行で定義されています。 { struct multiprec tmp; size_t e2_digits, e2_bits, e2; /* 下位の0の並んだ部分の、digit数, bit数 */ size_t pow1, pow2, max_pow; register digit_t ld; /* 正数以外は除く */ if (self->sign || ymp_is_zero(self)) return 0; /* e2 := {2**e が <self> を割る最大の e} */ for (e2_digits=0; e2_digits<self->used; ++e2_digits) { if (self->digits[e2_digits] != 0) break; } e2_bits = 0; for (ld = self->digits[e2_digits]; (ld & 0x1) == 0; ld>>=1) { ++e2_bits; } e2 = e2_digits * DIGIT_BIT + e2_bits; ymp_alloca_initialize(&tmp, self->used - e2_digits); ymp_div_2exp(&tmp, self, e2); /* 完全平方か? */ if (e2 % 2 == 0 && ymp_odd_sqrt(NULL, &tmp)) return 1; /* 完全立方か? */ if (e2 % 3 == 0 && ymp_odd_root(NULL, &tmp, 3)) return 1; /* 完全奇素数乗か ? */ max_pow = ymp_abs_clog2(tmp.digits, tmp.used); for (pow1=1, pow2=5; 1; pow1 += 6, pow2 += 6) { if (pow1 > max_pow) break; if (e2 % pow1 == 0 && ymp_odd_root(NULL, &tmp, pow1)) return 1; if (pow2 > max_pow) break; if (e2 % pow2 == 0 && ymp_odd_root(NULL, &tmp, pow2)) return 1; } return 0; }

---