half-GCD
code:halfgcd.sage.py
def pdivmod(u, v):
"""
polynomial version of divmod
"""
q = u // v
r = u - q*v
return (q, r)
def hgcd(u, v, min_degree=10):
"""
Calculate Half-GCD of (u, v)
f and g are univariate polynomial
"""
x = u.parent().gen()
if u.degree() < v.degree():
u, v = v, u
if 2*v.degree() < u.degree() or u.degree() < min_degree:
q = u // v
return matrix(1, -q], [0, 1)
m = u.degree() // 2
b0, c0 = pdivmod(u, x^m)
b1, c1 = pdivmod(v, x^m)
R = hgcd(b0, b1)
DE = R * matrix(u], [v)
q, f = pdivmod(d, e)
g0 = e // x^(m//2)
g1 = f // x^(m//2)
S = hgcd(g0, g1)
return S * matrix(0, 1], [1, -q) * R
def pgcd(u, v):
"""
fast implementation of polynomial GCD
using hgcd
"""
if u.degree() < v.degree():
u, v = v, u
if v == 0:
return u
if u % v == 0:
return v
if u.degree() < 10:
while v != 0:
u, v = v, u % v
return u
R = hgcd(u, v)
B = R * matrix(u], [v)
r = b0 % b1
if r == 0:
return b1
return pgcd(b1, r)
参考文献
を読もうー完ー
定理1
u,q,v,rを多項式とする。u=qv+r (deg(v)>deg(r)) のとき、qはu,vのうちxのdeg(u)冪から上位2(deg(u)-deg(v))+1桁で一意に定まる。
証明: 参考文献1を参照。
お気持ち:
$ u(x):=a_0+a_1x^1+\ldots+a_nx^n+\ldots+a_{n+m}x^{n+m}+\ldots+a_{n+2m}x^{n+2m}
$ v(x):=b_0+b_1x^1+\ldots+b_nx^n+\ldots+b_{n+m}x^{n+m}
とする。u-qvのn+m冪以上の項の係数を0にしたい。
筆算による除算によりqが求まる。
qは高々m次式、筆算の途中におけるu-qvの最上位の項はn+m冪以上である。
従って筆算の途中におけるu-qvの最上位の係数はvのn冪未満の係数の影響を受けない。
よって定理1が成り立つ（図を書こう！）。
この定理を用いて拡張ユークリッドの互除法を行うことで、GCD や $ v^{-1} \pmod u (但しGCD(u,v)=1) が高速に求まる。
N=deg(u)≧deg(v)とする。
u(x), v(x) に対して素朴に拡張ユークリッドの互除法を行うと一回の除算にO(N)掛かり、最悪ケースでは冪指数が1ずつしか減少しないため、全体で$ O(N^2)掛かる。一方、half-GCDでは$ O(N \log^2(N))となる。
拡張ユークリッドの互除法の各過程において
$ r_0=u,r_1=vと置いて
$ r_0=q_1r_1+r_2
$ r_1=q_2r_2+r_3
$ \vdots
$ r_{h-1}=q_hr_h
だとする。行列で書けば
$ \left( \begin{array}{c}\gcd(u,v) \\ 0 \end{array} \right) =\begin{pmatrix} 0 & 1 \\ 1 & -q_h \end{pmatrix}\dots\begin{pmatrix} 0 & 1 \\ 1 & -q_2 \end{pmatrix}\begin{pmatrix} 0 & 1 \\ 1 & -q_1 \end{pmatrix} \left( \begin{array}{c} u \\ v \end{array} \right) ・・・（＊）
となる。従って右辺の行列積が高速に求まればよい。
gcd(u,v)=1ならば
$ \begin{pmatrix} a & b \\ c & d \end{pmatrix} := \begin{pmatrix} 0 & 1 \\ 1 & -q_h \end{pmatrix}\dots\begin{pmatrix} 0 & 1 \\ 1 & -q_2 \end{pmatrix}\begin{pmatrix} 0 & 1 \\ 1 & -q_1 \end{pmatrix}
としたとき、$ au+bv=1だから$ b \equiv v^{-1}\bmod uが求まる。
$ \mathrm{HGCD}(r_0,r_1)＝一方の次数が$ \deg(r_0)/2-1以下になるまで（★）拡張互除法を適用した時の（＊）の右辺の行列積を返す関数
とする。
f div g ＝（ fをgで割ったときの商）とする。
以下のアルゴリズムを適用すると、次数の最大値が半減する。従って高々$ \log_2 ( \deg(u) )回適用すれば$ \gcd(u,v)および（＊）の右辺の行列積が求まる。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
$ R=\mathrm{HGCD}(u,v)
$ \left( \begin{array}{c} u \\ v \end{array} \right) \leftarrow R \left( \begin{array}{c} u \\ v \end{array} \right)
$ q \leftarrow u\ \mathrm{div}\ v
$ u \leftarrow u - qv
$ swap(u,v)
$ R \leftarrow \begin{pmatrix} 0 & 1 \\ 1 & -q \end{pmatrix} R
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
deg(u)-deg(v)<deg(v) を仮定してよい。さもなくば、既に★を満たしている。
m=ceil(deg(u)/2) とする。
以下のアルゴリズムを適用するとdeg(u)の次数が(3/4)倍になる。
一行目でHGCDが使えることは定理１とdeg(u)-deg(v)<deg(v)から従う。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
$ R_1=\mathrm{HGCD}(u\ \mathrm{div}\ x^m, v\ \mathrm{div}\ x^m)
$ \left( \begin{array}{c} u \\ v \end{array} \right) \leftarrow R_1 \left( \begin{array}{c} u \\ v \end{array} \right)
$ q=u\ \mathrm{div}\ v
$ u \leftarrow u - qv
$ swap(u,v)
$ R_1 \leftarrow \begin{pmatrix} 0 & 1 \\ 1 & -q \end{pmatrix} R_1
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
ここで一旦★を満たしていないか確認する。★を満たしていなければ次を実行する。
すると、★が満たされる。返り値は$ R_2R_1。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
$ R_2=\mathrm{HGCD}(u\ \mathrm{div}\ x^{m/2}, v\ \mathrm{div}\ x^{m/2})
$ \left( \begin{array}{c} u \\ v \end{array} \right) \leftarrow R_2 \left( \begin{array}{c} u \\ v \end{array} \right)
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
HGCDの計算量をT(N) （N=deg(u))とすると
T(N)=2T(N/2)+O(NlogN)
だからT(N)=N log^2(N)となる。
また、HGCDをlog_2(N)回適用する部分の計算量は
$ O(T(N)+T(N/2)+T(N/4)+\dots)=O(N \log^2(N))
である。