文献あり

安定化理論によるEuclidの互除法の安定化

はじめに

この記事は日曜数学 Advent Calendar 2025 の13日目の記事です．
この記事では，Euclidの互除法による多項式GCD計算を例に，代数計算に近似を入れる恐ろしさ，安定化理論の概要について紹介します．
主に1を参考に執筆しましたが，他にも日本語の文献がRIMS講究録に多数あがっています（安定化理論提案者の一人が日本人だからです，強い恩恵ですね）．

Euclidの互除法と係数膨張

整数のGCDを計算するアルゴリズムとして，Euclidの互除法をどこかで習った・習うと思います．多項式GCDでもまったく同じように互除法で計算ができます．
疑似コード(アルゴリズムを架空のプログラミング言語で記述したもの)は以下の通りです．読み方を解説しておくと、while 条件 do ～ end whileで条件が真の間～内の処理を繰り返す、:=で左辺を右辺で更新する、returnはアルゴリズムの出力を意味しています。

入力：$f,g\in \mathbb{R}[x]\backslash\{0\}$ with $\deg f\geq \deg g$
出力：$\gcd(f,g)\in \mathbb{R}[x]$

$a:=f,b:=g$
while $b\neq 0$ do
$r:=$(aをbで割った余り)
$a:=b,b:=r$
end while
return $a$

整数の場合，Euclidの互除法はかなり効率の良いアルゴリズムですが，多項式の場合はそうではありません．以下のような現象が起こるからです．

$f=x^9-11 x^7+26 x^6-30 x^5+19 x^4+5 x^3-24 x^2+19 x-5=(x - 1)^5 (x^4 + 5 x^3 + 4 x^2 + 6 x + 5),$
$g=5 x^9-28 x^8+67 x^7-96 x^6+106 x^5-105 x^4+83 x^3-42 x^2+11 x-1=(x - 1)^5 (5 x^4 - 3 x^3 + 2 x^2 - 6 x + 1)$
と定め，Euclidの互除法で$\gcd(f,g)$を計算しよう．手計算は地獄なので，余りの計算はMathematicaのPolynomialRemainderを使用した（他の数式処理システムにも同様の関数がある）．
$Out[i]$が互除法の過程で現れる割り算の余りを意味する．

$ In[20]:= PolynomialRemainder[-5 + 19 x - 24 x^2 + 5 x^3 + 19 x^4 - 30 x^5 + 26 x^6 - 11 x^7 + x^9, -1 + 11 x - 42 x^2 + 83 x^3 - 105 x^4 + 106 x^5 - 96 x^6 + 67 x^7 - 28 x^8 + 5 x^9, x]$
$ Out[20]= \frac{28 x^8}{5}-\frac{122 x^7}{5}+\frac{226 x^6}{5}-\frac{256 x^5}{5}+40 x^4-\frac{58 x^3}{5}-\frac{78 x^2}{5}+\frac{84 x}{5}-\frac{24}{5} $
$In[21]:= PolynomialRemainder[-1 + 11 x - 42 x^2 + 83 x^3 - 105 x^4 + 106 x^5 - 96 x^6 + 67 x^7 - 28 x^8 + 5 x^9, -(24/5) + (84 x)/5 - ( 78 x^2)/5 - (58 x^3)/5 + 40 x^4 - (256 x^5)/5 + (226 x^6)/5 - ( 122 x^7)/5 + (28 x^8)/5, x]$
$ Out[21]= -\frac{85 x^7}{196}-\frac{25 x^6}{196}+\frac{660 x^5}{49}-\frac{4925 x^4}{98}+\frac{16475 x^3}{196}-\frac{14565 x^2}{196}+\frac{475 x}{14}-\frac{310}{49}$
$ In[22]:= PolynomialRemainder[-(24/5) + (84 x)/5 - (78 x^2)/5 - ( 58 x^3)/5 + 40 x^4 - (256 x^5)/5 + (226 x^6)/5 - (122 x^7)/5 + ( 28 x^8)/5, -(310/49) + (475 x)/14 - (14565 x^2)/196 + (16475 x^3)/ 196 - (4925 x^4)/98 + (660 x^5)/49 - (25 x^6)/196 - (85 x^7)/196, x]$
$ Out[22]= \frac{327712 x^6}{1445}-\frac{2180696 x^5}{1445}+\frac{1197560 x^4}{289}-\frac{1739696 x^3}{289}+\frac{1411984 x^2}{289}-\frac{3038392 x}{1445}+\frac{542136}{1445}$
$ In[23]:= PolynomialRemainder[-(310/49) + (475 x)/14 - (14565 x^2)/ 196 + (16475 x^3)/196 - (4925 x^4)/98 + (660 x^5)/49 - (25 x^6)/ 196 - (85 x^7)/196, 542136/1445 - (3038392 x)/1445 + (1411984 x^2)/289 - (1739696 x^3)/ 289 + (1197560 x^4)/289 - (2180696 x^5)/1445 + (327712 x^6)/1445, x]$
$ Out[23]= \frac{183767875 x^5}{136983616}-\frac{918839375 x^4}{136983616}+\frac{918839375 x^3}{68491808}-\frac{918839375 x^2}{68491808}+\frac{918839375 x}{136983616}-\frac{183767875}{136983616}$
$ In[24]:= PolynomialRemainder[ 542136/1445 - (3038392 x)/1445 + (1411984 x^2)/289 - (1739696 x^3)/ 289 + (1197560 x^4)/289 - (2180696 x^5)/1445 + (327712 x^6)/1445 , -(183767875/136983616) + (918839375 x)/136983616 - (918839375 x^2)/ 68491808 + (918839375 x^3)/68491808 - (918839375 x^4)/136983616 + ( 183767875 x^5)/136983616, x] $
$Out[24]= 0$

以上の計算結果から，$\gcd(f,g)=\frac{183767875 x^5}{136983616}-\frac{918839375 x^4}{136983616}+\frac{918839375 x^3}{68491808}-\frac{918839375 x^2}{68491808}+\frac{918839375 x}{136983616}-\frac{183767875}{136983616}$となる．係数をうまく取っ払うと，$\gcd(f,g)=(x-1)^5$となる．

$f,g$の係数は高々三桁なのに，互除法の計算過程で現れる多項式の係数の分子，分母は巨大な桁数となっています．このような現象を係数膨張といいます．整数/有理数演算は桁数が大きくなるほど処理が重いため，Euclidの互除法を素朴に実装してもあまり効率的ではありません．

係数膨張をどうやって回避する？近似の恐ろしさ

有限体上のGCDから元のGCDを復元する，部分終結式剰余列を利用するなど，係数膨張を回避する方法は色々知られています．今回は，浮動小数を用いて係数膨張を回避できるか検討してみましょう．そもそも係数膨張が起こる原因は整数/有理数をそのまま扱っているからなので，浮動小数で置き換えればうまくいくんじゃないか，というわけです．
浮動小数で置き換えようと思うと，精度の問題は避けられません．例えば$\frac13$の場合，$0.333\ldots$と循環小数のため，どれだけ精度を上げたところで浮動小数で正確に表現することはできません．試しに，このような数を浮動小数で近似した時，互除法が機能するか具体例で試してみましょう．

$f=3x-1,g=x-\frac13$とする．明らかに$\gcd(f,g)=x-\frac13$だが，$\frac13$を浮動小数で近似して互除法を実行しよう．
$\frac13\approx0.333$のとき：
$In[49]:= PolynomialRemainder[3 x - 1, x - 0.333, x]$
$Out[49]= -0.001$
$PolynomialRemainder[x - 0.333, -0.001, x]$
$Out[50]= 0$
ゆえに$\gcd(f,g)=1$と出力されるが，これは誤った結果である．
$\frac13\approx0.33333333333$のとき：
$In[75]:= PolynomialRemainder[3 x - 1, x - 0.33333333333, x]$
$Out[75]= -1.\times 10^{-11}$
$In[76]:= PolynomialRemainder[ x - 0.33333333333, -1.* 10^{-11}, x]$
$Out[76]= 0$
先ほど同様に$\gcd(f,g)=1$となり，誤った結果が得られた．

すぐ分かることですが，どれだけ精度を上げても，$\gcd(f,g)=1$という誤った結果しか得られません．簡単な例ではありますが，浮動小数で置き換えるだけではうまくいかない，ということが分かってしまいました．
失敗の原因について考えてみましょう．疑似コードを再掲します．

入力：$f,g\in \mathbb{R}[x]\backslash\{0\}$ with $\deg f\geq \deg g$
出力：$\gcd(f,g)\in \mathbb{R}[x]$

$a:=f,b:=g$
while $b\neq 0$ do
$r:=$(aをbで割った余り)
$a:=b,b:=r$
end while
return $a$

4行目で，$b$は直前で計算した余りに更新されています．whileループは4行目が末尾なので，whileループを抜けるかどうか，判定が発生します．whileループを抜けるための条件は$b=0$，いわゆるゼロ判定というものになっています．
失敗の原因は，最初の余りの計算結果がゼロにはならず(浮動小数の近似精度をどんなに上げても，ぴったりゼロにはなりえない)，ゼロ判定をすり抜けてしまい，次の余りの計算に進んでしまったためです．互除法に限らず，代数計算（多項式の整除性判定，グレブナー基底，etc...）ではゼロ判定が頻出するため，代数計算に浮動小数による近似を導入すると，ゼロ判定が正常に機能せず正しい結果が得られません．このようなアルゴリズムを不安定と呼ぶことにします．

安定化理論

ゼロ判定を誤ることが不安定性をもたらしていることが分かったため，ゼロ判定がうまく機能するようにアルゴリズムを改良できればよいわけです．

白柳,Sweedlerは，区間演算とゼロ書き換えという操作を利用することで，不安定なアルゴリズムを安定化する，安定化理論を提案しました．
安定化理論の中身に入る前に，区間演算について解説しましょう．

区間演算

まず，区間演算からです．真の値$x$の代わりに，$x$を含む区間で表現し，区間どうしの加減乗除を定義します．区間演算の結果には，真の値どうしの演算結果が含まれるように演算を定義します．

$X=[a,b],Y=[c,d],a,b,c,d\in \mathbb{R}$として，以下のように演算を定める．
$X+Y=[a+c,b+d],$
$X-Y=[a-d,b-c],$
$X\cdot Y=[\min\{ab,ad,bc ,bd\},\max\{ab,ad,bc ,bd\}]$
$X/Y=X\cdot [1/c,1/d]$ (ただし，$0\notin Y$).

$\pi$を$X=[3.14159265,3.14159266]$，$e$を$Y=[2.71828182,2.71828183]$と区間で表現すると，
$X+Y=[5.85987447,5.85987449],X-Y=[0.42331082,0.42331084].$
$\pi+e=5.859874482048838...\in X+Y,$
$\pi-e=0.423310825130748...\in X-Y$と，真の値どうしの演算結果が含まれていることが分かる．

区間演算を定めたことで，区間を係数とする多項式も扱えるようになりました．

ゼロ書き換え

ひとまず，互除法に区間演算を導入してみましょう．
互除法のアルゴリズムはそのままで，入力多項式の係数を区間で表現してアルゴリズムを実行します．整数は幅なしの区間で表せることに注意してください．

$f=3x-1,g=x-\frac13$の係数を区間で表現し，$F=[3,3]x-[1,1],G=[1,1]x-[0.333333,0.333334] $と定める．
$F$の$G$による余りは$F-3G=[-2\times 10^{-6},1.0\times 10^{-6}]．$

次に，余りのゼロ判定が実行されますが，結局余りはゼロ（区間だと$[0,0]$）ではないので，whileループを抜け出せません．ここで，ゼロ書き換えという操作を行います．
ゼロ書き換えとは，ゼロ判定の直前にゼロを含む区間係数を$[0,0]$に置き換えること！
今回の場合，$[-2.0\times 10^{-6},1.0\times 10^{-6}]$はゼロを含むため，$[0,0]$に置き換えられ，whileループを抜け出すことができ，$[1,1]x-[0.333333,0.333334] $が出力されます（区間係数の範囲から，この出力は$\gcd(f,g)=x-\frac13$を含んでいることに注意します）．
改めて，安定化されたEuclidの互除法の疑似コードを記載しておきましょう．ZeroWritingはゼロ書き換えを施す関数です．

入力：区間多項式$f,g$ with $\deg f\geq \deg g$
出力：区間多項式$\gcd(f,g)$

$a:=f,b:=g$
while $\text{ZeroWriting}(b)\neq [0,0]$ do
$r:=$(aをbで割った余り)
$a:=b,b:=r$
end while
return $a$

出力が区間多項式となっていて，有理数係数のGCDが求められておらず，あくまで近似でしかないのでは？と思うかもしれない．実際その通りだが，実は区間に更に情報を付加し，計算履歴をとっておくことで，厳密係数のGCDを計算することができる．詳しくは，2を参照せよ．

ゼロを含むだけで$[0,0]$に書き換えてしまう，というかなりヤバそうな操作をしているわけですが，実際，入力多項式の係数を区間に変換する際に精度よく（=区間幅を狭く）変換しないと，ゼロ判定にかからないはずの区間もゼロ書き換えしてしまうことがあります．これが発生するとアルゴリズムは破綻します．
それじゃダメじゃないか！と思われるかもしれませんが，なんと，精度を十分上げることで，ゼロ書き換えは正しく行われること，アルゴリズムが安定化される，すなわち，出力が元の問題の解に収束することが示されています．
そのため，ある精度でアルゴリズムを実行してうまくいかなければ，さらに精度を上げて再度実行する，という操作を繰り返すことで，元の問題の解を近似することができるのです．

安定化理論の課題

十分精度を上げれば，アルゴリズムを安定化することができる，とのことでした．では，十分ってどこから？というのが気になってきます．例えばとあるアルゴリズムを安定化するために，係数の区間への変換の際，精度を$1$桁ずつ上げて実験するとします．$50$桁あたりから計算結果が安定してきて，$100$桁まで確認しても同様の計算結果だったとしましょう．よっしゃ，これが答えや！で終わって良いのでしょうか？もしかすると$101$桁目から急に結果が変わるかもしれません，$300$桁目からまた変わるかもしれません．このように，"十分な精度"がどこからか分からないと，計算結果がいまいち信頼しかねます（因数分解のように，後から入力と照らし合わせることができる問題ならまた話は別ですが）．
今でも必要精度の評価は難しいらしく（アルゴリズムだけでなく，入力にも依存するらしい），計算機実験に頼るしかないようです．