数理解析研究所講究録
第 1843 巻 2013 年 101-113

101

SNAP package and Improved QRGCD algorithm
増井貴明
TAKAAKI MASUI*

神戸大学人間発達環境学研究科
GRADUATE SCHOOL OF HUMAN DEVELOPMENT AND ENVIRONMENT, KOBE UNIVERSITY

長坂耕作
KOSAKU NAGASAKA \dag er

神戸大学人間発達環境学研究科
GRADUATE SCHOOL OF HUMAN DEVELOPMENT AND ENVIRONMENT, KOBE UNIVERSITY

1

はじめに

数式処理において，多項式の最大公約因子 ( $GCD$ ) 計算は重要なトピックのひとつであり，多くのアルゴ
リズムが知られている (詳しくは [4] を参照されたい).
ど

しかしながら，実際の問題 (制御理論や画像処理なにおいて，多くの場合，入力多項式の係数は浮動小数点数による数値計算の結果や実験値などから表されており，一般的に誤差を含んでいる．このような多項式に対しては通常のアルゴリズムでは $GCD$ を計算することは難しい．このような問題に対して近似 $GCD$ 計算という考え方が提唱され，これまでに多くの$)$研究が行われている．
QRGCD [2] (以下，[2] を元論文と呼ぶ) は一変数多項式の近似 $GCD$ 計算を行うアルゴリズムであり，数式処理ンフトウェア (CAS) Maple の近似計算ライブラリである SNAP パッヶにも実装されている．
$-\backslash ^{\backslash }\nearrow^{\backslash }$従って多くの新しく提唱されるアルゴリズムに対してペンチマークとしての役割を担っていると言える．しかしながら，QRGCD の理論的な背景はあまり明らかにされておらず，元論文の主定理の証明に誤りがあること，元論文のアルゴリズムと SNAP で実装されているアルゴリズムとの間にはいくつかの相違点があることが分かった．例えば，Bini, Boito [1] らは微小主係数多項式に対して QRGCD は正しい次数の近似$GCD$を計算することは難しいと述べているが，このことは実際には SNAP の実装における前処理が関係しているものの元論文などでは明記されていない．本講演では，これらのことについて速報する．
本稿では，2 章で元論文の QRGCD の概要について簡単に説明を行い，元論文で挙げられている主定理についての修正を行う．3 章では，元論文のアルゴリズムと SNAP で実装されているアルゴリズムとの相違点や問題点について述べ，4 章では，修正された主定理に基づく改善版 QRGCD アルゴリズム (ExQRGCD)を提案する．5 章でこれらのアルゴリズムに対していくつかの数値実験を行う．
なお，本講演で扱う一変数多項式の近似 $GCD$ は以下の定義にょるものとする．

inasui@main.h.kob6u.ac.jp\dagger nagasaka@main.h.kob . ac.jp$\triangleright u$102
定義 1(一変数多項式の近似 $GCD$ )$f(x),g(x)\in \mathbb{R}[x]$に対して，$d(x)\in \mathbb{R}[x]$が摂動の許容度$\epsilon>0$

における $f(x),g(x)$ の近似 $GCD$ である$d(x)$ が次を満たすときである．
とは，.

$f(x)+\Delta_{f}(x)=f_{1}(x)d(x), g(x)+\Delta_{g}(x)=g_{1}(x)d(x)$

ここで，$\Delta_{f}(x),$ $\Delta_{g}(x),$$fi(x),g_{1}(x)\in R[x]$ は$\deg(\Delta_{f})\leq\deg(f),\deg(\Delta_{g})\leq\deg(g), ||\Delta_{f}(x)||_{2}<\epsilon||f(x)||_{2}, ||\Delta_{9}(x)||_{2}<\epsilon||g(x)\Vert_{2}$を満たすものである ( $fi(x),$ $g_{1}(x)$ はそれぞれ $f(x),$ $g(x)$ の余因子と呼ぷ).

2

理論的な概要2.1

QRGCD の復習$\triangleleft$

本章では，QRGCD の元論文 [2] に基づく QRGCD の概要について述べる．次数が $m,n(m\geq n)$ の二つの一変数多項式 $f(x),g(x)\in R[x]$ を$f(x)=f_{m}x^{m}+f_{m-1}x^{m-1}+\cdots+f_{1}x+f_{0},$(1)$g(x)=g_{n}x^{n}+g_{n-1}x^{n-1}+\cdots+g_{1}x+g_{0}$.

$||f(x)||_{2}=||g(x)||_{2}=1)$$f(x)$ と $g(x)$ はそれぞれ単位 2- ノルムを持つとする すなわち，
このとき，
とし，
$($

$f(\cdot x)$と $g(x)$ の Sylvester 行列を次で定義する．

窟義 2(Sylvester 行列)$Syl(f,g)=[f_{m}g_{n} f_{m-i}g_{n-1}f_{m}g_{n} f_{m-1}g_{n-1} f_{m}g_{n}g_{1}f_{1}.\cdot\cdot f_{m,.-1}g_{n-1}g_{0}f_{0}.

\ldots g_{1}f_{1}.

g_{0}f_{0}]$.

(2)$\triangleleft$

以下特に断らない限り，ベクトルおよび行列に対して$||\cdot||_{2}$を}ノルム，$||\cdot||_{F}$をフロベニウスノルムとする．

Sylvester 行列の $QR$ 分解と $f(x)$ と $g(x)$ の $GCD$ との関係についてよく知られている定理を以下に示す．
定理 3$R\in \mathbb{R}^{(m+n)x(m+n)}$(2) 式の $QR$ 分解を $Syl(f,g)=QR$ とする．ここで $Q\in \mathbb{R}^{(m+n)x(m+n)}$ は直交行列，$R$ の 6 でない最後の行は $f(x)$ と $g(x)$ の (exact な) $GCD$ の係数ベクトは上三角行列である．このとき，ノレである$\triangleleft$

103$QR$ 分解は後退誤差の面では安定していることが知られているが，前進誤差の面では不数値計算において，安定になることがある．つまり，$S=Syt(f,g)$ に対して数値計算による $QR$ 分解の結果がだったとき，
後退誤差の $\Delta S=\hat{Q}\hat{R}-S$ が十分小さくても前進誤差である $\Delta R=\hat{R}-R$ が小さいことは保証されない．
$\hat{Q}\hat{R}$よって元論文では近似 $GCD$ に対応する根の検出に関して，以下の重要な二つの定理が示されている．

定理 4(Theorem 2, [2])$\omega$

き，
$R$

の

$\vec{0}$でない (近似的に$\vec{0}$.

$|\omega|<1)$上の) 単位円内にある共通根とする すなゎち，
このとでない) 最後の行は を含む $f(x)$ と $g(x)$ の近似 $GCD$ の近似因子の係数を $f(x)+\Delta_{f}(x),g(x)+\Delta_{g}(x)$ の ($\mathbb{C}$$($$\omega$ベクトルである$\triangleleft$

定理 5 (Theorem 3, [2])$f(x)$と $g(x)$ の Sylvester 行列の $QR$ 分解では，$f(x)$ と $g(x)$ の ($\mathbb{C}$上の) 単位円外にある共通根を検出することは難しい．

また，$\triangleleft$$p(x)\in \mathbb{R}[x]$に対して，reversal polynomial rev$(p)\in \mathbb{R}[x]$をrev $(p)=x^{\deg(p)}\cdot p(1/x)$で定義する．これらの定理および定義により，アルゴリズムは次のようになる．
アルゴリズム 1 (QRGCD [2])Input
OutputStep 1$u(x),v(x),$ $d(x)\in \mathbb{R}[x]^{1)}(d(x)$ が $f(x)$ と $g(x)$ の近似 $GCD)$

Step 2$R$$f(x),g(x)\in \mathbb{R}[x]$, 摂動の許容度$\epsilon>0$$Syl(f,g)=QR$ を計算する．
の右下から $(k+1)\cross(k+1)$ 次部分行列を$R_{22}^{(k)}$とし，$\Vert R_{22}^{(k)}\Vert_{2}>\epsilon$に対して以下のギャップを満たすを満たすとする．この(Case 1, 近似的に互いに素):$R_{22}^{(k)}$$k$

かっ $||R_{22}^{(k-1)}||_{2}<\epsilon$

を探す．

$\Vert R_{22}^{(0)}||_{2}>\epsilon$(i) $d_{1}(x);=1,$ $u(x),v(x)$ をの一番下の行から計算する．
(Case 2, 大きなギャップ検出):(i) $d_{1}(x)$ $:=R$ の下から 番目の行から計算する．
(Case 3, ギャップ検出):(i) $d_{1}(x)$ $:=R$ の下から番目の行から計算する．
(Case 4, Difficult case): それ以外
(i) $f(x),g(x)$ に対してアルゴリズム “Split” を行い，
$d_{1}(x)$ として $f(x),g(x)$ の単位円内にある共通根を計算する．
$f(x)$ と $g(x)$ の $d_{1}(x)$ に関する余因子の reversal polynomial に対して，Step $1-2$ を行$d_{2}(x)$ を計算する．
い，$Q^{T}$$\Vert R_{22}^{(k-1)}||_{2}<10\epsilon\Vert R_{22}^{(k)}\Vert_{2}$$k$$\exists k_{1},$$||R_{22}^{(k_{1}-1)}\Vert_{2}<10\epsilon\Vert R_{22}^{(k_{1})}\Vert_{2}$$k_{1}$

Step 3
Step 4$f(x)$ と $g(x)$ の $d(x)=d_{1}(x)d_{2}(x)$ に関する余因子に対して，StepStep5$u(x),$ $v(x),$ $d(x)$1-2 を行い，B\’ezout

係数 $u(x),v(x)$ を計算する．

を返す．

1 元論文では $u(x),$ $v(x),$ $d(x)$ の条件として $||f(x)-fi(x)d(X)||<\epsilon,$$)$$\Vert g(x)-g_{1}(x)d(x\rangle||<\epsilon, ||u(x)f(x)+v(x)g(x)-d(x)|\uparrow<\epsilon$が記されているが，証明も与えられておらず保証されないことから，ここでは省略する．

104

..
.

簡単にまとめると，アルゴリズムの流れは次のようになる．
$Syl(f,g)$ の $QR$ 分解により求める (normal side と呼ぶ).
単位円内にある近似 $GCD$ の因子 $d_{1}(x)$ は，単位円外にある近似 $GCD$ の因子晦 (x) は，reversal polynomial を用いて $Syl(rev(f/d_{1}),rev(g/d_{1}))$の $QR$ 分解により求める (rever 謝 side と呼ぶ).
$f(x)$$d(x)=d_{1}(x)\cdot d_{Q}(x)$ を求める．
と $g(x)$ の近似 $GCD$ として，注意 1 (Split” の概要)
‘sp 働” は入力多項式に対して，単位円内にある根と単位円外にある根に分離するアルゴリズムであり，

Graeffe’s root-sq uaring, ntour integration, Newton’s formula, Nevvton iteration, lifting steps に基づいている．Split” の目的については後ほど述べる$s$$\infty$$\triangleleft$

2.2

踏定理の修正前述の定理 4 および定理 5 に関して，元論文で与えられている証明にはいくつかの問題点があることが$S=Syl(f,g),$ $S+\Delta_{S}=Syl(f+\Delta_{f},g+\Delta_{g})$ とし (すな
分かった．以下にその概要を示す．簡単のため，わち，$\Delta_{S}=Syl(\Delta_{f}, \Delta_{9}))$,$QR$ 分解 $S=QR,$ $S+\Delta_{S}=\hat{Q}\hat{R}$ および $N,\hat{N}$をそれぞれ $S,$ $S+\Delta_{S}$ の零空間を張るベクトルを列に持つ行列とする．また，$d(x)$ を単位塾ノルムを持つ次数 の $f(x)$ と $g(x)$ の近似$f(x)+\Delta_{f}(x),g+\Delta_{g}(x)$ の exact な $GCD$ ) とする．このとき，$GCD$ (すなわち，$k$$SN=(S+\Delta_{S})\hat{N}=\vec{0},$(3)$(S+\Delta_{S})\hat{N}=(QR+\Delta_{S})\hat{N}\Rightarrow-R\hat{N}=Q^{T}\Delta_{S}\hat{N}$が成り立つ．

定瑠 4 の証明 ([2]) $r(x)$ を $R$ の行

るとき，ある$N_{k}\subset$$\vec{r}$を係数にもつ多項式とし， を$\omega$

抑があって，$|r(\omega)|\leq\Vert\Delta S||_{2}||N_{k}||_{2}$$|\omega|<1$

であるような $d(x)$ の根とすが成り立つことから $d(x)$ と (2- ノルムの意味で) 近 $V)r(x)$ が検出できるという流れになっている．しかしながら，これは $r(x)$ の根 における絶対評価であり，$d(x)$ と $r(x)$ が係数ベースで相対的に十分近いことは言っていない (例えば，$\omega$$10^{-16}(x-0.5)(x-0.8)$ と $10^{-16}(x+0.5)(x+0.8))$

定理 5 の証明 ([2])$\vec{\omega}_{2}=(\omega^{k}, \cdots,\omega^{1},\omega^{0})^{T}\in N_{k}$.

.

で (3) 式に対して $(S+\Delta_{S})\tilde{\omega}_{l}=\hat{Q}_{1}肪_{}*=\vec{0}$ が成り立つ

ことから， の $QR$ 分解で近似 $GCD$ を検出するためには， s のノルムが十分小さくなる必要がであれば近似 $GCD$ の検出は難しいという流れにある．よってもし単位円外にある共通根$\Delta$$S$$\omega$$(|\omega|>1)$

なっている．この議諭は近似 $GCD$ を根ベースで考える場合には正しいが，本稿および元論文の係数ペースで考える場合には根がずれてしまう可能性があるため正しいとは言えない．

以上の理由から，本章では定理 4 および定理 5 に類似した形で新しい定理を修正版として示すこととする．
そのために，いくつかの補題を以下に示す．
福題 6$S=QR,$ $S+\Delta_{S}=\hat{Q}\hat{R}$ とし， を $S+\Delta_{S}$ の零空間を張るベクトルを列に持つ行列とする．このとき，
が成り立つ
に対して
$\hat{N}$$\forall\vec{\omega}_{*}\in\hat{N}$$\triangleleft$$|r(\omega)|\leq||\Delta_{S}||_{2}||\vec{\omega}_{*}\Vert_{2}$

証明 (3) 式より，$-R\vec{\omega}_{*}=Q^{T}\Delta_{S}\vec{\omega}_{*}.$

105
したがって，$\Rightarrow\Vert R\vec{\omega}_{*}\Vert_{2}=\Vert Q^{T}\Delta_{S}\vec{\omega}_{*}\Vert_{2}$$\Rightarrow\Vert R\vec{\omega}_{*}\Vert_{2}\leq\Vert\Delta_{S}\vec{\omega}_{*}\Vert_{2}$$\Rightarrow\Vert r(\omega)|\leq\Vert R\vec{w}_{*}\Vert_{2}\leq\Vert\Delta_{S}\vec{\omega}_{*}\Vert_{2}\leq\Vert\Delta_{S}||_{2}||\vec{\omega}_{*}\Vert_{2}.$$\iota$また，次の補題に関する準備として，$p(x)= \sum_{i=0}^{k}p_{i}x^{i},$ $p_{i}\in \mathbb{C}$に対して，多項式集合$\mathcal{P}_{\epsilon}(p)$を$\mathcal{P}_{e}(p)=\{\tilde{p}(x)\in \mathbb{C}[x]|||\tilde{p}-\vec{\tilde{p}}||_{2}\leq\epsilon||p\neg|_{2}, \deg(\tilde{p})\leq\deg(p)\}$で定義する．ここで，$\vec{p}=$$(Pk, \ldots,Po)$,$\vec{\tilde{p}}=(\tilde{p}_{k}, \ldots,\tilde{p}0)$である．

補題 7 (Corollary 7, [5])に対して 2),$\omega\in \mathbb{C}$

ヨ$\tilde{p}(x)\in \mathcal{P}_{c}(p),\tilde{p}(\omega)=0\Leftrightarrow|p(\omega)|\leq\epsilon||p(x)||_{2}||\vec{\omega}_{*}||_{2}$$\triangleleft$

これらの補題により前述の定理 4 を次のように修正する．
定理 8(修正版定瑠 4)$r(x)$ を
これまでの記法に従うものとし，

とき，$\ovalbox{\tt\small REJECT}_{2}^{\Delta_{\mathcal{S}2}}rx\ll 1$

証明

ならば，

$r(x)$

補題 6 より，$R$の下から $(k+1)$ 番目の行

は $f(x)$ と $g(x)$ の近似 $GCD$$\vec{r}$である．
$\triangleleft$$|r(\omega)|\leq\Vert\Delta_{S}||_{2}\Vert\vec{\omega}_{*}||_{2}=_{rx}^{As}\ovalbox{\tt\small REJECT}_{2}\cdot||r(x)||_{2}\Vert||\vec{\omega}_{*}\Vert_{2}$

ならば，な $GCD$ とすると，$|r(\omega)|\leq\epsilon||r(x)||_{2}||\vec{\omega}_{*}\Vert$

を係数に持つ多項式とする．このが成り立つ．もし$\ovalbox{\tt\small REJECT}_{2}^{2}rx=\epsilon\Delta\ll 1$

である．よって補題 7 から，$d(x)$ を $f(x)+\Delta_{f}(x),g(x)+\Delta_{g}(x)$ の exactが成り立つ．よって $r(x)$ は $f(x)$ と $g(x)$ の近似 $GCD$ である．

$\Vert\tilde{r}-d]|_{2}\leq\epsilon||r\urcorner|_{2}$$I$

次に，単位円外にある共通根の検出性に関して，次の補題を述べる．
補題 9

Sylvester 行列の $QR$ 分解により得られた$R$$r(x)$ を
に対して，

$R$

のある行

$\vec{r}$を係数に持っ多項式とする．こ$r(x)$ に含まれる根は相対的に近いものから検出されるわけではない．
のとき，

証明

証明については，本稿では省略するが，あらすじとしては Sylvester’s single sum [3] に基づいて$r(x)$ に含まれる根の大きさと $PRS$の関係について調べることができる．発表では割愛したものとして，以下に具体例を挙げる．

6

例 1 (artificial common inside roots)
次の $f(x)$ と $g(x)$ に対して Sylvester 行列の $QR$ 分解を考える．

$f_{-}(x)=(x+0.01)(x-0.01)(x+100)(x-100)$$g(x)=(x+0.010000001)(x-0.09)(x+100.00001)(x-200)$

ここで，摂動の許容度$\epsilon=10^{-12}$に対して $f(x)$ と $g(x)$ の近似 $GCDd(x)$ は，$d(x)=9.99998\cross 10^{-4}x^{2}+0.999999x+9.9999S\cross 10^{-4}$$\approx 9.99998\cross 10^{-4}(x+0$2 [5] では双対ノルムについて示されているが，ここでは 3/$)$.00100000 $)(x+1000.00)$

ルムの相対評価に書き換えている．

106となる．しかしながら， の下から 2 番目，3 番目および 4 番目の行ベクトルを係数に持つ多項式は次のよ
$R$

うになる．

4 番目: 0. $00143003x^{3}+1.39855x^{2}-0.00559397x-6.99252\cross 10^{-6}$$\approx o.\omega 1430oe(x-0.00499981)(x+$0.001 $)(x+977.993)$ ,3 番目: 0. $20978x^{2}-0.0W839074x-1.04885x10^{-6}$$\approx 0.20978(x-0.00499978)(x+$ 0.001 ,$)$2 番目: $o.\alpha\}565685x+5.65685x10^{-6}$$\approx 0$.00565685 $(x+$ 0.001 .
$)$3 番目の多項式は $d(x)$ とは無関係な $f(x)$ と $g(x)$ の共通根でないものを含んでいる．この影響により本来の単位円外にある共通根が正しく検出されていない．以後このような単位円外にある共通根のかわりに検出されてしまう単位円内の根を artfflcial $\infty mmon$ inside roots と呼ぶ$\triangleleft$

以上の理由から前述の定理 5 を次のように修正する．

窟瑠 10(燭日阪建瑠 5)$Syl$

一行列の $QR$ 分解では，単位円外の共通根を含む近似 $GCD$ を計算することは難しい．正しく計算を行うには，多くの場合に normal side と reveoed side の両方から計算を行い，それらを合わせる必要がある．

$\triangleleft$

3SNAP の実装
QRGCD は Maple の SNAP パッケージとして実装されている．本章では，まず SNAP の実装を示し (アルゴリズム 2), 元論文のアルゴリズムとの間にあるいくつかの相違点について述べる．なお，行列 に対して
$A$$||A||s$は $A$ の各行ベクトルの》ノルムの和を表すとする $( すなわち，||A||_{S}=\sum||\vec{a}_{i}\Vert_{2}, A=(\vec{a}_{1}, \ldots,\vec{a}_{n})^{T})$.

注慮 2 (Split” の目的)元論文でも用いられているアルゴリズム ‘iSp 廊 “ の単位円内の根と単位円外の根に分離する目的に関しては

元論文では明記されていない．分離する理由として考えられることは次の二つである，一つ目は低次の多項式に対する $QR$ 分解の方が高次のものに比べて前進農差がより小さくなることで精度が保たれることにある．二つ目は前述の a 鵬 iHcial $\infty mmon$ inside $m$詔の影響により，単位円外にある共通根が正しく検出できない可能性があることである．

3.1$\triangleleft$

元論文との相違点

Maple の SNAP における実装が元論文のアルゴリズムと異なる点として，以下の 4 点について述べる．
用いるノルム SNAP の実装では元論文で用いられている 3 ノルムではなく，$||\cdot||s$が用いられているが，$||A||_{2}\leq||A||_{F}\leq\Vert A||s$ であるこ
そのことは元論文では明記されていない．任意の行列 に対して，$||A||s$とが分かるが，は緩和しすぎではないかと考えられる．実験により計算される近似 $GCD$ に大きな違いが見受けられないことは確認したが，後述の改善版 QRGCD では計算コストとの兼ね合い$A$からもフロベニウスノルムを採用することとする．

107
アルゴリズム 2 (QRGCD in the SNAP package)

Input$f(x),g(x)\in \mathbb{R}[x]$Output$u(x),v(x),d(x)\in \mathbb{R}[x]$, 摂動の許容度$\epsilon>0$

sucb that $\Vert f(x)-f_{1}(x)d(x)\Vert_{2}+||g(x)-g_{1}(x)d(x)||_{2}<10\epsilon$, or

failure”
Step 1
Step 2
Step 3

閾値を$\epsilon/100$

として’“find non-zero term” を $f(x),g(x)$ に対して行う．

$Syl(f,g)=QR$ を計算する．
$R$ の右下から $(k+1)\cross(k+1)$
次部分行列を$R_{22}^{(k)}$とし，$\Vert R_{22}^{(k)}$Isに対して以下のギャップを満たすを満たすとする．この(Case 1, 近似的に互いに素):$\hslash_{2}^{k)}$$>\epsilon$

$k$

かつ

$\Vert R_{22}^{(k-1)}\Vert s<\epsilon$

を探す．

$||R_{22}^{(0)}||s>\epsilon$(i) $d_{1}(x):=1,$ $u(x),v(x)$ をの一番下の行から計算する．
(Case 2, 大きなギャッ 7 検出):$:=R$ の下から
(i)
番目の行から計算する．
(Case 3, ギャップ検出):$;=R$ の下から
(i)
番目の行から計算する．
(Case 4, Difflcult case): それ以外)
(i) $R$ の行ベクトルで a ノルムが以上であるものを係数に持つ多項式を$Q^{T}$$\Vert R_{22}^{(k-1)}\Vert s<10\epsilon||R_{22}^{(k)}\Vert s$$d_{1}(x)$$k$$\exists k_{1},$$d_{1}(x)$$||R_{22}^{(k_{1}-1)}||s<10\epsilon||R_{22}^{(k_{1})}\Vert s$$k_{1}$$5\alpha$$p_{1}(x),p_{2}(x),$ $(\deg(p_{1})>\deg(p_{2}))$$0\epsilon$とする．

(ii) アルゴリズム “Split” を用いて $p_{1}(x)$ の単位円内にある根からなる多項式を計算する．
(iii)

Step 4$p_{1_{in}}(x)$$p_{1_{n}}:(x)$$d_{1}(x)$ を計算する．
と働 (x) に対して Step 2-3 を行い，

$d_{1}(x)=1$ ならば rev $(f),rev(g)$に対して再帰呼び出しを行い，結果の reversal$polyn\infty$

mial をとる．
Step 5$f(x),g(x)$ の $d_{1}(x)$ に関する余因子の reversal polynomial に対して Step 2-3 を行い，
$d_{2}(x)$を計算する．

Step 6$f(x),g(x)$ の d(x) $=$ dl(x) 吻 $(x)$ に関する余因子に対して Step2-3 を行い，
$u(x),v(x)$Step7$u(x),v(x),d(x)$ を返す．

$\phi$

を計算する．

Split” を適用する多項式元論文では入力多項式に対して “Split” を適用するが，SNAP の実装では，上三角行列 $R$ に現れる多項式 (すなわち，$PRS$ ) のーつに対してのみ “Split” を適用する．
$PRS$ の Split

は入力多項式の Split に比べて次数が低いので，計算コストが小さくて済むと考えられるが，$5000\epsilon$

以上の 2- ノルムを持つ多項式を選ぶ根拠については示されておらず，実験的な値であると考えられる．

fail-safe ループ元論文のアルゴリズムでは，常に単位円内にある共通根から検出を試み，その後で単位円外にある共通根の検出を行う．しかし，もし単位円外にある共通根の個数が単位円内のものより非常に多い場合，このような単位円内にある共通根は検出が難しい可能性があるため，単位円外にある共通根から検出を試みる方がよい．SNAP では normal side から計算を行って，近似 $GCD$ が計算されない場合に reversal side から計算し直す実装になっている．

ここでは，特に重要な相違点として，SNAP の前処理のーつである “find non-zero terms” について以下で詳しく述べることとする．
‘喰 nd non-zero terms”

SNAP の実装では前処理の一つとして “find non-zero terms” が実行される．このアルゴリズムは，入力多項式に対して主係数および定数を$\epsilon/100$

以上となるように書き換えるもの

108である．たとえば，$f(x)= \sum_{1=0}^{m}f_{i}x^{i}$に対して，$|f_{m}|, |f_{m-1}|, \ldots, |f_{m’+1}|<\epsilon/100, |f_{m’}|\geq\epsilon/100,$$|f_{m}|\geq\epsilon/100,$$|$

漏 n-ll?

とすると，“find non-zero terms” によって$\cdots$, $|fo|<\epsilon/100$$f(x)= \sum_{=0}^{\dot{m}’m^{n}}f_{m^{n}+i^{X^{i}}}$となる．この前処理によって微小主係数多項式に対して QRGCD を計算すると，いくつかの間題が生じることが分かった．Bini, Boito[1] らは微小主係数多項式が QRGCD の -conditioned であると述べているが，このことは次の実験から，“find non-zero terms” を除いた SNAP の実装と比較することで微小主係数多項式が QRGCDのnditiond ではないことが分かる．
$i\mathbb{I}$$ill-\infty$

例 2(微小主係数多項式)
この実験は Binf, Boito [11 により行われており，計算される近似 $GCD$ の次数が問題となる．ここでは 10000 ペアの以下の多項式を生成した．
$f(x)=(\alpha x^{3}+2x^{2}-x+5)(x^{4}+7x^{2}-x+1)$ ,$g(x)=(\alpha x^{3}+2x^{2}-x+5)(x^{3}-x^{2}+4x-2)$

上で，$\alpha\in[10^{-10},10^{-5}]$$\epsilon=10^{-6}$.

はランダムに選ぶ．これら $f(x),g(x)$ に対して Maple 16 上で摂動の許容度

で実験を行った (D\’igits$:\Rightarrow 16$). 表 1 で，$E_{K}QRGCD$ ” , “SNAP” および “SNAP $(w/0$

表 1: 微小主係数多項式: 実験結果

F.n.t)” はそれぞれ我々の改善版 QRGCD アルゴリズム，SNAP の実装，‘ffind non-zero terms” を除いた SNAP の実装を表す．結果として，‘i 伽 non-zero tems” の影響により近似 $GCD$ として正し$d$

い次数のものが検出されていないことが分かる．一方で，SNAP $(w/0 f.n.t)$ と ExQRGCD は全ての多項式に対して次数 3 のものを検出できている

4$\triangleleft$

改魯版 QRGCD アルゴリズム

本章では，2.2 章で述べた修正版定理に基づいた，改善版 QRGCD アルゴリズム (ExQRGCD) について述べる．ExQRGCD (アルゴリズム 3) は SNAP や元論文より計算速度は遅いが，数学的根拠に基づいており: 与えた摂動の許容度 を満たすことを保証している．ExQRGCD では，行列のノルム計算に関して計算$\epsilon$

コストを抑えるために，2- ノルムのかわりにフロベニウスノルムを採用している．また，artificial common
inside roots の影響を抑えるために ExQRGCD では，理想的な近似 $GCD$ に十分近いものから少しずつ計算を行う．従って，normal side と reversal side の計算を何度も繰り返す必要があり，計算コストは元論文やSNAP に比べると多くかかる．以下では元諭文を書き換えた Case l, 2 について詳しく述べることとする．

109
アルゴリズム 3 (ExQRGCD)

Input
OutputStep 1

Step 2
Step 3, 摂動の許容度$u(x),v(x),d(x)\in \mathbb{R}[x]$ such $that||f(x)-fi(x)d(x)\Vert_{2}<\epsilon,$ $\Vert g(x)-g_{1}(x)d(x)\Vert_{2}<\epsilon$
frst $ide$ (norm 1 か revers 1) を決定する．
$i:=1,$
Hrst side が $noma1\Rightarrow f_{1}(x)$ $:=f(x),g_{1}(x)$ $:=g(x)$ ,first side が $revelsal\Rightarrow f_{1}(x)$ $:=rev(f),g_{1}(x)$ $:=rev(g)$$f(x),g(x)\in \mathbb{R}[x]$$\epsilon>0$$Syl(f,g)=QR$ を計算する．
$R$の右下から $(k+1)\cross(k+1)$ 次部分行列を$R_{22}^{(k)}$

とし， を

$R$

$\vec{r}_{k}$の下から $(k+1)$ 番目の行ベクトルとする．
(Case 1, 近似的に互いに棄):( ) $d_{i}(x):=1.$(Case 2, 試し創り):(i)を満たす全ての $k\geq 1$ に対して，を計算する．
$||R_{22}^{(0)}||_{F}>\epsilon\sqrt{m+n}$$i$$||R_{22}^{(0)}||_{F}\leq\epsilon\sqrt{m+n}$$||R_{22}^{(k-1)}||_{F}\leq\epsilon\sqrt{m+n}$(ii) 三つの最小の$er_{k}\ovalbox{\tt\small REJECT}$$er_{k}:=\Vert R_{22}^{(k-1)}\Vert_{F}/\Vert\vec{r}_{k}\Vert_{2}$

こ対して，$d_{i}(x);=r_{k}(x)$$d(x);= \prod d_{T}(x)$ に対してとし，$\exists f_{i}(x),g_{1}(x)$such that

かつ $\Vert g(x)-g_{1}(x)d(x)\Vert_{2}$ ならば，Step 4 へ進む (以下の (v) における再帰呼び出しにおいて，$fi(x),g_{1}(x)$ が見つからなけれ$||f(x)-fi(x)d(x)\Vert_{2}<\epsilon$

ば $d_{i}(x):=1$ として Step 4 へ進む).

(iii)

$R$

の行ベクトルで 2- ノルムが 1.0 以上となる 2 つのものをとり，それぞれ$p_{1}(x),p_{2}(x)(\deg(p_{1})>\deg($ 勉 $))$とする．

(iv) $p_{1}(x),p_{2}(x)$ に対してアルゴリズム “Split” を行い，単位円内にある根のみからな$(x)$ を求める．
る多項式$p_{1_{n}},(x),p_{2_{i}}$。

(v)に対して Step 2-3 を行い，を計算する．
$d(x):= \prod d_{i}(x)$ に関する余因子$i:=i+1$ とし，normal $rightarrow reversd$ を入れ替え，$p_{1_{\ln}}(x),p_{2_{In}}(x)$Step 4$fi(x),g_{1}(x)$(reversal side ならば reversal polynomial をとったもの) に対して Step 2-

3 を行い，あるStep 54.1$d_{i}(x)$$i$$fi(x),g_{1}(x),d(x)$に対して $d_{i}(x)=d_{i+1}(x)=1$ となるまで，$d_{i+1}(x)$ を計算する．

を返す．

Case 1 の判定 : 近似的に互いに素QRGCD と ExQRGCD は，与えられた二つの多項式 $f(x),g(x)$ の Sylvester 行列 $S=Syl(f,g)$ を $QR$
分解することで $R$ の行ベクトルを係数に持っ多項式として近似因子を検出するアルゴリズムであることは既に述べたとおりである．従って，$S+\Delta_{S}=\hat{Q}\hat{R}$としたとき瓦から得られる $f(x)+\Delta_{f}(x),g(x)+\Delta_{9}(x)$に対して，QRGCD と ExQRGCD は $S=QR$ の から得られる $r(x)$ がど
れだけ $d(x)$ に近いかという観点に立つものである．よって $d(x)$ そのものをから直接的に計算することの exact な $GCD$ である $d(x)$

$R$

$R$

はできない．

,
ここで，$S=QR$ により から得られた次数 の $r(x)$ を exact な $GCD$ として持つ多項式をを考える．このとき，Sylvester 行列の $QR$とし，その Sylvester 行列の $QR$ 分解$R$$9(x)+\delta_{g}(x)$$k$$f(x)+\delta_{f}(x)$$S+\delta_{S}=\tilde{Q}\tilde{R}$分解の性質から$S+.U, ||U||_{F}=||R_{22}^{(k-1)}||_{F}\Rightarrow rank(S+U)=m+n-k$110であることが分かる．ここで$R_{22}^{(k-1)}$

は unstmcturedである．一方，本来$S$に加えるべき $U$ は structuredである必要があるが，多くの場合に stmctured な行列のノルムは unstructured な行列のノルムより大きい
ため，$\Vert R_{22}^{(k-1)}||_{F}\leq||\delta_{S}||_{F}$が成り立つ．従って，$||R_{22}^{(k-1)}\Vert_{F}>\epsilon\sqrt{m+n}$であれば，多くの場合に$||\delta_{S}||_{F}>\epsilon\sqrt{m+n}$が成り立つので，近似的に素である，と判定できる．

4.2

Case 2 の判定 : 試し割り

これまでの議諭により，一度の $QR$ 分解では全ての単位円内にある共通根を検出できない可能性がある．
ExQRGCD では，ギャップの判定条件としてを計算するが，$r(x)$ としてできる限り高次数のものを検出するために，alfficial common inside roots の影響を最小限にする必要がある．従っでに関して $r(x)$ を計算しようと考えた．実験により，多くの場合に最小のて，できる限り小さな$erk:=\Vert R_{22}^{(k-1)}\Vert_{F}/||\vec{r}_{k}\Vert_{2}$$er_{k}$摂動の許容度$\epsilon$$er_{k}$を十分に満たす近似 $GCD$ の因子が計算できることを確認したが，念のために最小のもの三つを計算することとした．また，“Split” を適用する多項式に対しては，次のように変更した．

“Split” 奄適用する多項式 $p_{1}(x),p_{2}(x)$ を $f(x),g(x)$ の $PRS,$ $r(x)$ を近似 $GCD$ とする．このとき，$f(x)+\delta_{f}(x)=f_{1}(x)r(x),g(x)+\delta_{g}(x)=g_{1}(x)r(x)$,$u_{1}(x)f(x)+v_{1}(x)g(x)=p_{1}(x), u_{2}(x)f(x)+v_{2}g(x)=p_{2}(x)$$u_{1}(x),v_{1}(x),$は単位 2- ノルムを持ち，$Syl(f,g)$ の $QR$ 分解で得らが成り立つ．ここで，れた直行行列のある列ベクトルを係数にもつ多項式である．ここで，$u_{2}(x),$ $v_{2}(x)$$Q$$p_{1}(x)=(u_{i}(x)fi(x)+v_{1}(x)g_{1}(x))r(x)-(u_{1}(x)\delta_{f}(x)+v_{1}(x)\delta_{g}(x))$,$p_{2}(x)=(u_{2}(x)fi(x)+v_{2}(x)g_{1}(x))r(x)-(u_{2}(x)\delta_{f}(x)+v_{2}(x)\delta_{g}(x))$$r(x)$ は摂動の許容度
より，$(\Vert\delta_{f}(x)||_{2}+||\delta_{9}(x)||_{2})/\dot{m}n\{||p_{1}(x)||_{2}, ||p_{2}(x)||_{2}\}$

近似 $GCD$ である．従って，”Split” を適用する多項式として $p_{1}(x)$ , 勉 (x) の
い

$($

例えば，

$||p_{1}(x)||_{2},$ $||p_{2}(x)||_{2}\geq 1.0)$における $p_{1}(x),p_{2}(x)$ の$2-\nearrow$

ルムが十分小さくなものを取ることとする (このような $p_{1}(x),p_{2}(x)$ は $f(x),g(x)$の共通根以外のものを含んでいる可能性がある).

5

数値実験

以下に，3 つの実験に関してその結果を示す．なお，表記として，ExQRGCD は改善版 QRGCD アルゴ
リズム，SNAP は SNAP パッケージに実装されている QRGCD アルゴリズム，QRGCD は元論文のアルの対数グラフになっていること
ゴリズムを示している．図 2 において摂動量の平均に関して，縦軸がに注意されたい．なお，発表においては実験 1,3 のみを扱い，本稿では新たに実験 2 を追加した．
$log_{10}$実験 1, (exact な根を持つランダム多項式)$i=1,$$\ldots,$$10$に対して，以下を満たす 100 個の多項式ペア $f(x),g(x)$ を生成する．
$f(x)=f_{1}(x)d(x),g(x)=g_{1}(x)d(x)$ ,$d(x)= \sum_{j=0}^{5i}d_{j}\dot{d}, f_{1}(x)=\sum_{j=0}^{\epsilon:}f_{1,j}\dot{\theta},g_{1}(x)=j\wedge\sum_{-}^{\epsilon:}g_{1,j}i.$

111
上で，

$f_{i,j,g_{1,j},d_{j}}\in[-99,99]\subset \mathbb{Z}$$||g(x)||_{2}=1)$を満たすものをランダムに選び，$f(x),g(x)$ は単位 2- ノルム $(||f(x)||_{2}=$を持つように正規化を行い，Digits: $=10$ で丸めた．これらに対して摂動の許容度$\epsilon=10^{-5}$で計算した (Digits $:=16$ ).

図 1 より，検出次数に関して ExQRGCD と SNAP の間には大きな違いはないが，QRGCD は他に比べて良くない．SNAP と QRGCD は表 2 から “failure” や与えた摂動の許容度 を満たさないものがいくっか検出されている (ExQRGCD は全ての多項式に対してこのような問題は発生していない). 計算速度に関しては，ExQRGCD と QRGCD は SNAP と比べてそれぞれ，1.59 倍，2.08 倍遅い結果となった．“伽 lure”の数を除けば，摂動量に関しては SNAP が最も良い結果となった．
$\epsilon$

実験 2. (摂動多項式を加えたランダム多項式)$i=1,$$\ldots,$$10$に対して，以下を満たす 100 個の多項式ペアノ (x), $g(x)$ を生成する．
$f(x)=f_{1}(x)d(x)/||fi(x)d(x)||_{2}+10^{-8}\Delta_{f}(x)/||\Delta_{f}(x)\Vert_{2},$$g(x)=g_{1}(x)d\langle x)/||g_{1}(x)d(x)\Vert_{2}+10^{-8}\Delta_{g}(x)/\Vert\Delta_{g}(x)\Vert_{2},$$\Delta_{f}(x)=\sum_{j=0}^{10i}\Delta_{f_{j}}\dot{d}, \Delta_{g}(x)=\sum_{j=0}^{10i}\Delta_{g_{j}}\dot{d}.$上で，を満たすものをランダムに選び，$fi(x),g_{1}(x),d(x)$ は実験 1. の多項式とで計算した (Digits $:=16$ ).
し，Digits: $=10$ で丸めた．これらに対して摂動の許容度$\Delta_{f_{j}},$$\Delta_{g_{j}}\in-.[-99,99]\subset \mathbb{Z}$$\epsilon=10^{-5}$

図 1 と表 2 より，ExQRGCD が SNAP と QRGCD に比べて非常に良いことが分かる．計算速度に関しては，EXQRGCD と QRGCD は SNAP と比べてそれぞれ，1.99 倍，2.26 倍遅い結果となった．“failure”の数を除けば，摂動量に関しては SNAP が最も良い結果となった．

$2356910 35678910 3456789 10$
Case

Case

Cw

図 1: 検出次数 : 実験結果 (failure” は$0$として換算)

Random Polynornlals without Perturbation

Random PolynomiaJs with Perturbation

Artificlrl Common Hoots

CR

Case

$C\varpi$

図 2: 摂動量

$($

すなわち，

$(\Vert\Delta_{f}||_{2}+\Vert\Delta_{9}||_{2})/2)$: 実験結果 (failure” は除く)

112

表 2: SNAP, QRGCD における “failure” と適切でない摂動量の数 : 実験結果

契験 3. (artificial common roots を含む多項式)$i=1,$$\ldots,$$10$に対して，以下を満たす 100 個の多項式ペア $f(x),g(x)$ を生成する．
$f(x)=d(x) \prod_{j=1}^{2i}(x-\omega_{f,j})\prod_{j=1}^{2i}(x- 。f,j)$,,$g(x)=d(x) \prod_{J=1}^{2i}(x-\omega_{g,j})\prod_{j=1}^{2i}(x-\hat{\omega}_{g,j})$$d(x)= \prod_{j=1}^{31}(x-\omega_{d,j})\prod_{j=1}^{3i}(x-\hat{\omega}_{d,j})$.

$\omega.,j=O(10^{-2}),\hat{\omega}.,j=0(10^{2})$$f(x),g(x)$ は単位》ノルム上で，を満たすものをランダムに選び，$||g(x)||_{2}=1)$$(\Vert f(x)\Vert_{2}=$を持つように正規化を行い，Digits: $=10$ で丸めた．これらに対して摂動の許容度$\epsilon=10^{-6}$で計算した (Digits: $=16$ ).

図 1 と表 2 より， QRGCD が SNAP と QRGCD に比べて非常に良いことが分かる．計算速度に関しては，ExQRGCD と QRGCD は SNAP と比べてそれぞれ，39.8 倍，52.6 倍遅い結果となった (SNAP は実験で用いた多項式のうち，約 62%に対して ‘繊皿 uxe” を返しており，SNAP の計算速度は ‘繊皿 ure” の出ない
簡単な問題に対してのみの平均となっている). 摂動量に関しては，ExQRGCD が最も良い結果となった．
$E$

6

まとめ

本講演では，QRGCD アルゴリズムの元論文に関して，主定理の修正を行い，また SNAP の実装と元論文のアルゴリズムの間にあるいくつかの相違点について速報した．特に SNAP の前処理の一つである“find non-zero tems” の影響により，本来 QRGCD の皿- nditioned ではない微小主係数多項式に対して，
SNAP が正しく計算できないことを確認した．我々は，修正した主定理に基づき改善版 QRGCD アルゴリ$\omega$113
ズム (ExQRGCD)を提案した．実験結果からは，摂動を含む多項式および arti 丘 cial common roots を含む多項式に対しては，ExQRGCD が SNAP や QRGCD と比べて非常に良い結果となった．特に artificial
common roots を含む多項式に対して，SNAP が多くの場合 (failure” を返す所を，ExQRGCD では与えた摂動の許容度を満たすもので高次数の近似 $GCD$ を計算できていることを確認した．今後は，Modern なア
ルゴリズム (Faetgcd, UVGCD, GPGCD など) とも比較を行いたい．

Acknowledgement
This work was supported in part by Japanese Ministry of Education, Culture, Sports, Science andechnology under Grant-in-Aid for Young Scientists, NEXT KAKENHI $(227\alpha)011)$ .

$I^{1}$

参考文献[1] Dario A. Bini and Paola Boito. $A$ fast algorithm for approximate polynomial $GCD$ based on struc-tured matrix computations. In Numerical methods for structured matrices and apphcations, Vol. 199
of Oper. Theory $Adv$ . Appl., pp. 155-173. Birkh\"auser Verlag, Basel, 2010.

[2] Robert M. Corless, Stephen M. Watt, and Lihong Zhi. $QR$ factoring to compute the $GCD$ of
univariate approximate polynomials. IEEE Trans. Signal Process., Vol. 52, No. 12, pp. 3394-3402,

2004.

[3] Carlos D’Andrea, Hoon Hong, Teresa Krick, and Agnes Szanto. An elementary proof of sylvester’s
double sums for subresultants. J. Symb. Comput., Vol. 42, No. 3, pp. 290-297, March 2007.
[4] Joachim Von Zur Gathen and Jurgen Gerhard. Modem Computer Algebra. Cambridge University
Press, New York, $NY$ , USA, 2 edition, 2003.

[5] Hans J. Stetter. The nearest polynomial with a given zero, and similar problems. SIGSAM Bull.,Vol. 33, No. 4, pp. 2-4, December 1999.