中国の剰余定理のアルゴリズム(Algorithm.10.14)

036.crt.html に引き続き、 (Gathen and Gerhard 2013) を参考に中国の剰余アルゴリズムを実装している。

“10.3 Fast Chinese remaindering” に6つの疑似コードがある。

• 済: Algorithm.10.3
• Algorithm.10.14
• Algorithm.10.16
• Algorithm.10.18
• Algorithm.10.20
• Algorithm.10.22

本記事では “Algorithm.10.14” を調査、実装する。

このアルゴリズムは一つの整数(被乗数)に対して、複数の剰余を求めるものである。

入力を以下とする。( \(r=2^{k}\) とする.)

• \(ms=m_{1},\cdots,m_{r-1}\) : 互いに素な整数(除数)のリスト
• \(M\) : alg10_3 に \(ms\) を入力して得られる行列
• \(f\) : 整数(被乗数)

出力を以下とする。

• \(f \mod m_{1}, \cdots ,f \mod m_{r-1}\) : 複数の剰余

このアルゴリズムは再帰呼び出しで定義されている。 ms=[2, 3, 5, 7, 11, 13, 17, 19], f=9699791 とすると以下のような呼び出しとなる。

• alg10_14(ms=[2, 3, 5, 7, 11, 13, 17, 19], M=algo10_3(ms), f=9699791);
  • alg10_14(ms=[2, 3, 5, 7], M=algo10_3(ms), f=101)
    • alg10_14(ms=[2, 3], M=algo10_3(ms), f=5)
      • alg10_14(ms=[2], M=algo10_3(ms), f=1)
      • alg10_14(ms=[3], M=algo10_3(ms), f=2)
    • alg10_14(ms=[5, 7], M=algo10_3(ms), f=31)
      • alg10_14(ms=[5], M=algo10_3(ms), f=1)
      • alg10_14(ms=[7], M=algo10_3(ms), f=3)
  • alg10_14(ms=[11, 13, 17, 19], M=algo10_3(ms), f=101)
    • alg10_14(ms=[11, 13], M=algo10_3(ms), f=101)
      • alg10_14(ms=[11], M=algo10_3(ms), f=2)
      • alg10_14(ms=[13], M=algo10_3(ms), f=10)
    • alg10_14(ms=[17, 19], M=algo10_3(ms), f=101)
      • alg10_14(ms=[17], M=algo10_3(ms), f=16)
      • alg10_14(ms=[19], M=algo10_3(ms), f=6)

実装は疑似コードを忠実に再現したものと、効率を考えたもので2つ作成した。

def alg10_14(ms, f):
    if len(ms) == 1:
        return [f]
    else:
        M = alg10_3(ms)
        f0 = f % M[-2][0]
        f1 = f % M[-2][1]
        return alg10_14(ms[:len(ms)//2], f0) + alg10_14(ms[len(ms)//2:], f1)

def alg10_14x(i_ms_s, i_ms_e, f, M, i_f0, j_f0):
    if i_ms_e - i_ms_s == 1:
        return [f]
    else:
        f0 = f % M[i_f0][j_f0]
        f1 = f % M[i_f0][j_f0+1]
        return alg10_14x(i_ms_s, i_ms_s+(i_ms_e-i_ms_s)//2, f0, M, i_f0-1, 2*j_f0) \
             + alg10_14x(i_ms_s+(i_ms_e-i_ms_s)//2, i_ms_e, f1, M, i_f0-1, 2*j_f0+2)

ms = [2, 3, 5, 7, 11, 13, 17, 19]
f = 9699791

res0 = [f % m for m in ms]

res1 = alg10_14(ms=ms, f=f)

M = alg10_3(ms)
res2 = alg10_14x(0, len(ms), f, M, -2, 0)

assert all([x==y for x, y in zip(res0, res1)])
assert all([x==y for x, y in zip(res1, res2)])
res1