Emacs + 暗号 + 数学 + プログラム
The Art and Mathematics of Genji-Kō – OranLooney.com を読んだ。
和書 (小林 2024) を一回さらっと読んだが、もう一度読んでみよう。
Skyline algorithm for packing 2D rectangles – Julien Vernay を読んだ。
(Gathen and Gerhard 2013) を参考にして、中国の剰余定理のアルゴリズムを実装する。
def alg10_20(cs, i_s, i_e, M, i_r0, j_r0): if i_e - i_s == 1: return cs[i_s] else: return M[i_r0][j_r0+1] * alg10_20(cs, i_s, i_s+(i_e-i_s)//2, M, i_r0-1, 2*j_r0) \ + M[i_r0][j_r0] * alg10_20(cs, i_s+(i_e-i_s)//2, i_e, M, i_r0-1, 2*j_r0+2)
ms = [2, 3, 5, 7, 11, 13, 17, 19] ms, r, k = pack_ms(ms) M = alg10_3(ms) cs = list(range(len(ms))) alg10_20(cs, 0, len(ms), M, -2, 0)
25524916
The stereographic projection of the Stern–Brocot tree を読んだ。
興味深いが後半は理解が覚束無い。記事に内に同一の筆者による記事：Some notes on the Stern–Brocot tree があるので読んでみたい。
Writes large correct programs を読んだ。
(Gathen and Gerhard 2013) を参考にして、中国の剰余定理のアルゴリズムを実装する。
このアルゴリムズは入力m1,…,mr-1で、法miでのm=m1…mi-1,mi+1,…,mr-1 の逆元を算出するものである。
拡張ユークリッド互除法法を使用するため拡張ユークリッド互除法を参考にした。
実装は以下のようになる。
def alg10_18(ms): m = functools.reduce(lambda a, b: a*b, ms, 1) ms2 = [mi**2 for mi in ms] # 1. call alg10_16 (alg10_16はalg10_3とalg10_14を併せたもの) M = alg10_3(ms2) m_rem_mi2 = alg10_14x(0, len(ms2), m, M, -2, 0) # 2. m/mi rem mi m_rem_mi = [m_rem2//mi for mi, m_rem2 in zip(ms, m_rem_mi2)] # 3. extended euclidean algo si = [] for mi, mi_rem in zip(ms, m_rem_mi): a, b = xgcd(mi, mi_rem) if 0 < b: si.append(b) else: si.append(b+mi) return si
期待する結果は以下である。
ms = [2, 3, 5, 7, 11, 13, 17, 19] m = functools.reduce(lambda a, b: a*b, ms, 1) si = [] for mi in ms: for i in range(1, mi): if 1 == (i * m//mi) % mi: si.append(i) break; si
| 1 | 1 | 2 | 6 | 7 | 5 | 16 | 18 |
実装した結果は以下になり、一致している。
ms = [2, 3, 5, 7, 11, 13, 17, 19] ms, r, k = pack_ms(ms) alg10_18(ms)
| 1 | 1 | 2 | 6 | 7 | 5 | 16 | 18 |
Traits are a Local Maxima · thunderseethe’s devlog を読んだ。
036.crt.html に引き続き、 (Gathen and Gerhard 2013) を参考に中国の剰余アルゴリズムを実装している。
“10.3 Fast Chinese remaindering” に6つの疑似コードがある。
本記事では “Algorithm.10.14” を調査、実装する。
このアルゴリズムは一つの整数(被乗数)に対して、複数の剰余を求めるものである。
入力を以下とする。( \(r=2^{k}\) とする.)
出力を以下とする。
このアルゴリズムは再帰呼び出しで定義されている。 ms=[2, 3, 5, 7, 11, 13, 17, 19], f=9699791 とすると以下のような呼び出しとなる。
実装は疑似コードを忠実に再現したものと、効率を考えたもので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
| 1 | 2 | 1 | 3 | 2 | 10 | 16 | 6 |
036.crt.html に引き続き、 (Gathen and Gerhard 2013) を参考に中国の剰余アルゴリズムを実装している。
“10.3 Fast Chinese remaindering” に6つの疑似コードがある。
本記事では “Algorithm.10.14” を調査、実装する。
このアルゴリズムは一つの整数(被乗数)に対して、複数の剰余を求めるものである。
入力を以下とする。( \(r=2^{k}\) とする.)
出力を以下とする。
このアルゴリズムは再帰呼び出しで定義されている。 ms=[2, 3, 5, 7, 11, 13, 17, 19], f=9699791 とすると以下のような呼び出しとなる。
実装は疑似コードを忠実に再現したものと、効率を考えたもので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
| 1 | 2 | 1 | 3 | 2 | 10 | 16 | 6 |
(Gathen and Gerhard 2013) を参考に中国の剰余アルゴリズムを実装している。
“10.3 Fast Chinese remaindering” に6つの疑似コードがある。
上記をPythonで実装しようと考えた。
このアルゴリズムは他のアルゴリズムから利用される補助的なものである。本文では多項式環を扱えるように一般化しているが、今回は \(\mathbb{Z}/p\mathbb{Z}\) だけを扱うことにする。
入力は互いに素な数である。入力のサイズはパラメータ \(k\) で決定する。 \(r=2^{k}\) 個の値を入力とすることになる。出力は \(r \times r\) の行列となる。
k=2として具体的な計算手順を実行する。入力は \(r=2^{k}=4\) の数 \(m_{1},m_{2},m_{3},m_{4}\) とする。出力となる行列 \(M\) を \(r\times r=4×4\) の大きさで初期化する。
\begin{pmatrix}
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
\end{pmatrix}
初期化した行列の1行目に入力の数を配置する。
\begin{pmatrix}
m_{1} & m_{2} & m_{3} & m_{4}\\
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
\end{pmatrix}
行列の2行目に1行目の隣同士の積を左詰めで設定する。
\begin{pmatrix}
m_{1} & m_{2} & m_{3} & m_{4}\\
m_{1}m_{2} & m_{3}m_{4} & 0 & 0\\
0 & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
\end{pmatrix}
行列の3行目に2行目の隣同士の積を左詰めで設定する。
\begin{pmatrix}
m_{1} & m_{2} & m_{3} & m_{4}\\
m_{1}m_{2} & m_{3}m_{4} & 0 & 0\\
m_{1}m_{2}m_{3}m_{4} & 0 & 0 & 0\\
0 & 0 & 0 & 0\\
\end{pmatrix}
上記が出力となる。
参考文献の疑似コードをPythonで実装した。
def alg10_3(ms): k=0 while 2**k < len(ms): k += 1 r = 2**k ms = ms + ([0] * (r-len(ms))) M = [[0]*r for _ in range(k+1)] for i in range(r): M[0][i] = ms[i] for i in range(1, k+1): for j in range(0, 2**(k-i)): M[i][j] = M[i-1][2*j] * M[i-1][2*j+1] return M
出力は以下のようなる。
alg10_3([2, 3, 5, 7, 11, 13, 17, 19])
| 2 | 3 | 5 | 7 | 11 | 13 | 17 | 19 |
| 6 | 35 | 143 | 323 | 0 | 0 | 0 | 0 |
| 210 | 46189 | 0 | 0 | 0 | 0 | 0 | 0 |
| 9699690 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |