(19) 中国の剰余定理アルゴリズムの高速版(要加筆)

互いに素な数で割った余りがそれぞれ、であるようなnはを法として一意に存在します。(中国の剰余定理)

解nは以下の式で求めることができます。

ここで、、、はを法とするの逆元です。すなわち、を満たします。

これをMaximaで素直に書き下してみましょう。

p : [3, 5, 7]$
u : [2, 3, 2]$
k : length(p)$
pi : product(p[i], i, 1, k)$
c : makelist(pi/p[i], i, 1, k)$
mod(sum(c[i]*inv_mod(c[i], p[i])*u[i], i, 1, k), pi);

使用した機能

• inv_mod
• product

このnを求めるのに、kが2のべき乗の時には高速なアルゴリズムがあるそうです。
これをMaximaで天下り的に載せておきます。

/* 4.4.11m */
matchdeclare(i, integerp)$
defrule(exponent_of_2, 2^i, i)$

p : [3, 7, 11, 101, 103, 257, 701, 7001]$
u : [1, 6, 5, 67, 87, 230, 321, 5999]$

k : length(p)$
t : exponent_of_2(factor(k))$
if t = false then error("not power of 2")$

pp : product(p[m], m, 1, k)$
d : makelist(inv_mod(pp / p[n], p[n]), n, 1, k)$
/* d[n]*(product of ps except for p[n]) = 1 modulo p[n] */
q : zeromatrix(k, t + 1)$
s : zeromatrix(k, t + 1)$
for j : 1 thru t + 1 do
  for i : 1 next i + 2^(j - 1) thru k do (
    q[i][j] : product(p[m], m, i, i + 2^(j - 1) - 1))$
for b : 1 thru k do s[b][1] : d[b] * u[b]$
for j : 2 thru t + 1 do
  for i : 1 next i + 2^(j - 1) thru k do (
    s[i][j] : s[i][j - 1] * q[i + 2^(j - 2)][j - 1] + s[i + 2^(j - 2)][j - 1] * q[i][j - 1])$
answer : mod(s[1][t + 1], q[1][t + 1]);

使用した機能

• matchdeclare
• integerp
• defrule

参考文献

1. Crandall, Mathematica―理工系ツールとしての (アジソン ウェスレイ・トッパン情報科学シリーズ) p.75-p.78
2. ウィキペディア, 中国の剰余定理