数理解析研究所講究録 1456 巻 2005 年 27-34

27

1 変数代数方程式の一つの近接根クラスタに含まれる近接根の計算章 *佐々木建昭 \dag er

AKIRA TERUITATEAKI SASAKI照井筑波大学大学院数理物質科学研究科
GRADUATE SCHOOL OF PURE AND APPLIED SCIENCES, UNIVERSITY OF TSUKUBA

Abstract
複素係数 1 変数多項式 $P(x)$ は近接根クラスタ (複数でもよい) を持ち、 クラスタは他根から十分に分離できるとする。 このクラスタに含まれる全根を同時に、他の根を計算しないで、効率良く計算する算法を提案する。算法は DuraJld-KerneI 法を近接根用に改変したものであり、 非常に簡単である。 末尾に実験例を示す。

1

はじめに
古来より、代数方程式の解法として非常に多くの方法が提案されており、 近接根を含まない多項式に対しては、 たとえ悪条件であっても見事に解$\langle$

方法も発表されている (Fortune [2], 他)。

また、 与えられた多項式を、複素円盤の内部と外部にそれぞれ根を持つ二つの多項式の積に分解する高速算法により、計
しかし、 近接根を含む多項式に関してはまだ研究の余地がある。
算量も決定されつつある Pan[7], 他)
Yakoubsohn [11] は、 近接根に対する Newton 法の収束定理を証明したが、 その算法は効率的とは言えな$\langle$。

い。最近、 Sakurai らが近接根用の巧妙な算法を発表した [5]。その方法は近接根を桁落ちなく計算するが、算法の複雑さが難点である。
一方、 近接根を含まない多項式に関しては、 Durand-Kerner 法 ([1], [4]) などの簡潔かつ効率的な方法がある。 同様に、一つの近接根クラスタ内の全ての根を同時に、他の根を一切計算することなく、 しかも桁落ちを起こすことなく、効率良く計算することはできないものであろう力 17 本稿ではそのような方法を提案する。 この方法では、与多項式は変換され、近接根以外の根の影響は多項式の商として取り込まれる。その際、変換された多項式および反復の度に計算される商は、 (クラスタの近接度に依存する) 一定次数分だけ計算すればよく、結果として高次多項式に対しては著しい効率化を達成する。たとえば、 10 個の近接根を含むクラスタを持つ 1000 次の多項式の場合でも、変換後の多項式は 20\sim 25 華分を、商は 10\sim 15 次分を計算すればよい。我々の方法は近接根の近くの根の影響だけを取り込んで計算を行う方法である、 と言えよう。
以下、本稿では次の内容を述べる。2 では算法の基本的なアイデアについて述べる。3 では近接根用 DurandKerner 法について述べる。 4 では実際に算法を実行する上での注意点について述べ、計算例を掲げる。 5 では実験結果を示す。

$*\mathrm{t}$erui Qmath${}^{\overline{\mathrm{t}}}\mathrm{s}\mathrm{a}\mathrm{s}\mathrm{a}\mathrm{k}\mathrm{i}$. sukuba .
. suku ba .

$\mathrm{a}\mathrm{C}.$$\mathrm{t}$(Dmath$\mathrm{t}$JP$\mathrm{a}\mathrm{c}.\mathrm{j}\mathrm{P}$28
2

算法の基本的アイデア
我々の算法は次の三つのステップからなる。

1. 近接根クラスタの位置$\rho$

クラスタ中の根の個数 $m$ を決定する。 ここで、 近接度とは全体を含む円盤の半径] で定義する。

と近接度[近接根クラスタの広がり] [$/$$\ovalbox{\tt\small REJECT} \mathrm{E}$$2\gamma_{\text{、}}$

2. ターゲットとするクラスタを原点に移動し、 $x\mapsto cx$ なるスケール変換で、 クラスタを半径 1 程度の大きさに拡大する (このステップは多倍長計算を必要とする)。

3. 一つのクラスタ中の全根だけを扱うように Durand-Kerner 法を改変し、 反復的に近似解を計算する(このステップは固定長計算で十分である)。

この算法は、 3) の反復算法が不明な点を除けば、 実に平凡かつ簡単である。 3) の反復算法も、 我々が以前に発表した論文 (実根用 Durand-Kerner 法 [10] : 実多項式の実根だけを同時反復法で計算する方法) と全く同じである。 こんな簡単な算法がうまく働くのか? と思われようが、後で見るように、 実にうまく働くのである。

まず、 ステップ 1 について。我々の算法では、 ターゲットとする近接根クラスタの位置

近接度$2\gamma_{\text{、}}$

およびクラスタに含まれる根の個数 m、が必要である。 これらの情報は、近似無平方分解 ([6], [S] を参照)

で$\rho_{\text{、}}$

得られるので、本稿ではこれらの情報は既知として議論を進める。
近接根がクラスタをなすと言っても、 そのクラスタが他の根から区別できなければ、実際的にクラスタをなすとは言えず、 我々の算法は破綻する。近接根クラスタが他の根から明確に区別できるかどうかについて、 著者らは数年前に簡潔な定理を得ており

次に、 ステップ 2 について。与式を$P(\bm{x})$$[9]_{\text{、}}$

その定理で区別できることを確認しておく。

とする:$P(x)=c_{n}x^{n}+\cdots$ 十 $c_{m}x^{m}+\cdots+c_{0}$.

(1)

説明の便宜上、 多項式 $P(x+\rho)$ の根は原点を中心とする半径 2 の複素円盤内に含まれるとする (この仮定となる) Terui-Sasaki のは実際の計算では必要ないが、 この仮定の下では近接根クラスタの広がりは約$2\gamma$

近接根分離定理によれば、 $P(x+\rho)$ の $m$ 次未満の項の係数は、 次数が一つ下がれば約し、 $\gamma<1/9$ であるは$\tilde{P}(x)$このことを念頭に、 スケール変換のノルム規格化のための数値である)$[9]_{\text{。}}$:$\tilde{P}(x)=\overline{c}_{n}x^{n}+\cdots+$

近接根に対応する$P(x)\mapsto\tilde{P}\{x$$\overline{P}(x)$)

。

$\gamma$倍減少する。 ただ$=\eta P(\gamma(x+\rho))$を考えよう

輪 $x^{m}+\cdots+\tilde{c}_{0}$ .

$(\eta$(2)の $m$ 個の根は原点を中心とする半径 2 程度の複素円盤内にあるので、$\tilde{c}_{m}=1$,$\max\{|\tilde{c}_{m-1}|, \cdots, |\tilde{c}_{0}|\}\approx 1$,(3)となるように決めることができる。 このとき、高次項の係数は次のようになる。
$\overline{c}_{m+j}=o(\gamma^{j-1})$$(j=1,2, \ldots,n-m)$ .

この性質は実際の算法において非常に重要な役割をはたす。

(4)

29
近接根用 Durand-Kerner 法

3

式 (2) の$\overline{P}(x)$に対して、 論文 [10] のアイデアに従い、 Durand-Kerner 法を改変する。

$\tilde{F}(x)$に対し、方程式$\overline{P}(x)=0$の全根を Durand-Kerner の 2 次法で求める場合、根の初期値を与え、次の反復式によって、 根$\zeta_{1},$$\zeta_{n}$$\ldots,$に対応する$z_{j}^{(\nu+1)}=z_{f}^{\langle\nu)}- \frac{\tilde{P}(z_{j}^{(l^{J}\rangle})}{\tilde{c}_{n}\tilde{Q}’(z_{j}^{(\nu)})}$

ここで、 求める近接根を$\zeta_{1},$

,

$\nu$

回反復後の近似値$\tilde{Q}(x)=\prod_{j=1}^{n}(x-z_{j}^{(\nu)})$

その近似値をそれぞれ$\zeta_{m\text{、}}$$\ldots,$,$z_{n}^{(\nu)}$$z_{1}^{(\nu)},$$\ldots,$$j=1,$$\ldots,$$z_{m}^{(\iota\prime)}(m<n)$$z_{1}^{(\nu)},$$\ldots,$$m$を用いないため、 式 (5) の$\overline{Q}(x)$$z_{j}^{(\nu)}$$\text{接}\ovalbox{\tt\small REJECT} \mathrm{f}\mathrm{l}$の近 -‘\lambda ffi\llcorner からなる
$f1$$\mathrm{E}\mp\prod_{j=1}^{m}(x-z_{f}^{(\nu)})$で除しを計算する。

(5)とする。本稿の算法は、を構成することができない。
に近づく (係数同士が互いにほぼ等し$\tilde{Q}(x)$$\tilde{P}(x)$$\tilde{Q}(x)$を、 $P(x)$ を近$\ldots,$.

,...,
これ以外の根の近似値はところで、すべての近似値
が根らに十分近づくと、
くなる) 。そこで、本稿の算法では,のうち, 求める近接根以外の近似値からなる因子$z_{n}^{(\nu)}$$z_{m+1}^{(\nu)}$$z_{n}^{(0)}$$z_{1}^{(0)},$$\prod_{j=m+1}^{n}(x-z_{j}^{(\nu)})$$_{arrow}^{\sim} \Re-\mathrm{q}\mathrm{u}\mathrm{o}(\tilde{P}(x), \prod_{j=1}^{m}(x-z_{j}^{(_{V})}))$で近似する。すなわち、本稿における Durand-Kerner 法の反復式として、 次式を採用する。
$z_{j}^{(\iota/+1)}=z_{j}^{(\nu)}- \frac{\tilde{P}(z_{j}^{\langle\nu)})}{Q_{\nu}’(z_{j}^{(\nu)})}$,$Q_{\nu}(x)= \prod_{j=1}^{m}(x-z_{j}^{(\nu)})\cdot \mathrm{q}\mathrm{u}\mathrm{o}(\tilde{P}(x),\prod_{j=1}^{m}(x-z_{j}^{(\nu)}))$

反復式 (6) は次の収束定理をみたす,$j=1_{7}\ldots,$ $n$.

(6)$[10]_{\text{。}}$

定理 1
をそれぞれ
を方程式 $\tilde{P}(x)=0$ の根とし、 互いに異なるものとする。
はに+分近いものとする。 このとき、 反復式 (6) は 2 次#疎する。
近似値とし、$z_{m}^{(\nu)}$$z_{1}^{(\nu)},$$\zeta_{m}$$\zeta_{1},$$\ldots,$$\ldots,$$z_{j}^{(\nu)}$$—\ovalbox{\tt\small REJECT}-\Xi B$

近似値$z_{j}^{(\nu)}$$\zeta_{m}$$\zeta_{1},$$\ldots,$の$(_{j}$の誤差を $\epsilon_{j}^{(\nu)}=z_{j}^{(\nu)}-\zeta_{j}(j=1, \ldots, m)$ とおくと、 次式を得る。

(7)$\epsilon_{j}=z_{j}-(\nu+1)(\nu+1)\zeta_{j}=\epsilon_{j}-(\nu)\frac{\tilde{P}(z_{j}^{(\nu)})}{Q_{\nu}’(z_{j}^{(\nu)})}=\in_{j}-(\nu)\epsilon_{j}^{(\nu\rangle}\frac{\prod_{k=1,\neq j}^{n}(z_{j}^{(\nu)}-\zeta_{k}\rangle}{Q_{\nu}’(z_{j}^{(\nu)})}$

一方, $R(x)= \prod_{k=m+1}^{n}$ ($x$-\mbox{\boldmath $\zeta$}k) \in (\mbox{\boldmath $\nu$})=max{夏 l\mbox{\boldmath $\nu$})l, . . . ,$)$$|\epsilon_{m}^{(\nu)}|$} とおくと,$Q_{\nu}’(z_{j}^{(\nu)})$は次式で表される.

.

$Q_{\nu}’(z_{j}^{(\nu)})= \prod_{k=1,\neq j}^{m}(z_{j}^{(\nu)}-z_{k}^{(\nu)})\prod_{k=m+1}^{n}(z_{j}^{(\nu)}-\zeta_{k})+(\epsilon_{1}^{(/)}+\cdots+\epsilon_{m}^{(\nu)})|R’(z_{j}^{(\nu)})+O((\epsilon^{(\nu)})^{2})$(8)

式 (8) を用いることにより, 式 (7) の最岱明の有理式は次の通り変形できる,$\frac{\prod_{k=1,\neq j}^{n}(z_{j}^{(\nu)}-\zeta_{k})}{Q_{\nu}(z_{j}^{(\nu)})},=\prod_{k=1,\neq j}^{m}(\frac{z_{j}^{(\nu)}-\zeta_{k}}{z_{j}^{(\nu)}-z_{k}^{(\nu)}})+O(\epsilon^{(\nu)})=\prod_{k=1,\neq j}^{m}(1+\frac{\epsilon_{k}^{(\nu)}}{\zeta_{j}+\epsilon_{j}^{(\nu)}-z_{k}^{(\nu)}})+O(\epsilon^{(\nu)})$.

$= \prod_{k=1,\neq j}^{m}(1+\frac{\epsilon_{k}^{(\nu)}}{\zeta_{j}-z_{k}^{(\nu)}}+O(\in_{k}(\nu)\epsilon_{j}^{(1J)}))+O(\epsilon^{(\nu)})=1+\sum_{k=1,\neq j}^{m}\frac{\epsilon_{k}^{(\nu)}}{\zeta_{j}-z_{k}^{(\nu)}}+O(\epsilon^{(\nu)})$$(9\rangle$

よって, 式 (9) を式 (7) に代入することにより)$\epsilon_{j}^{(\nu+1)}=O((\epsilon^{(Ir)})^{2})$

を得る.

30
4

算法の実際と計算例
具体的に算法を構成するには、 次の 2 点に注意する必要がある。

a) 近接根に対応する

b)$\tilde{P}(z_{j}^{(\nu\rangle})$と$\tilde{P}(x)$$Q_{\nu}’(z_{j}^{\langle\nu)})$a) に関しては、$\overline{P}(x)$の初期値の選択。

の計算の効率化。

の $m$ 次以下の項からなる式(10)$\tilde{c}_{m}x^{m}+\cdots+\tilde{c}_{0}$の根を (通常の) Durand-Kerner 法で計算し、 それらを近接雑用 Durand-Kerner 法の初期値に設定する。
(仮定より、求める近接根以外の根は複素平面上の単位円から十分離れており、式 (10) の根は求める近接根に近い位置に存在する。)

数に関する
(4) に
b) に関しては
半径 2 程度の複素円盤内にあることから、$\prime \mathrm{r}4^{\hslash}\mathrm{g}$$\mathrm{f}\mathrm{f}+_{\backslash }$$\ovalbox{\tt\small REJECT}^{\backslash }\Xi$$\text{、}$$P_{\nu}^{\mathrm{A}}(x)$する。 $\hat{P}_{\nu}(x)=\prod_{j=1}^{m}(x-z_{j}^{(\nu)})$ とおくと、$z_{1}^{(\nu)},$$\cdots$, z#ゝがの各項の係数は 1 よりそれほど大きくはない。 したがって、(11)$qk\mathit{2}_{m}x^{n-m}+\cdots+q_{1}^{(\nu)}x+q_{0}^{(\nu)}$$Q_{\nu}(x)=\mathrm{q}\mathrm{u}\mathrm{o}(\overline{P}(x),\hat{P}_{\nu}(x))=$とおくと、 $|q_{j1}^{(\nu)}||=o(\gamma^{j-1})(j=1,2, \ldots, n-m)$ と見倣せる.
反$\text{復}/$$Q_{\nu}’(x)$式 (6) では、$\overline{P}(z_{j}^{(\nu)})$と$Q_{\nu}’(z_{j}^{(\nu)})$さえ F.f\doteqdot \not\in よく$\equiv-+\mathrm{p}\text{算}$できればよいが、の高次項は計算結果にはほとんど影響しない。 したがって、$\tilde{P}(x)$$|z_{j}^{(\nu)}|\sim<2$なので、$\tilde{P}(x)$と $Q_{\nu}(x)$ はそれぞれ $\mu+m$ 次、次以下の項を計算しておけば十分である。 ここで、 は数値根に必要な精度ととえば、 倍精度計算で $\gamma=0.1$ の場合、 $\mu=20$ 程度で十分である。
次に計算例を示す。 (この計算は数式処理システム GAL を用いて行った。)$\mu$$\gamma$と$\mu$で決まる整数である。 た

例 1$\delta=1.0\mathrm{x}10^{-4}$に対し、$P(x)=(x-1.2+\delta)(x-1.2)(x-1.2-\delta)(x-1.2-2\delta)$(12)$\cross(x+1)\{(x+2)^{2}+2\}\{(x-5)^{2}+1\}\{(x+5\rangle^{2}+3\}\{(x-7)^{2}+3\}\{x-10)$とおき、 $x=1.2$ 付近の近接根を計算する。多倍長浮動 」 数演算の精度は 10 進 34 桁 (IEEE 倍精度の約$’$$\backslash$2 倍) とする。 多倍長演算で近似無平方分解を行うと、 x–L20004997253175025157497235308794450 が$P(x)$ の近似 4 重因子として計算される。そこで、 多倍長演算を用いて、 $x=1.2000499725317502515\cdots$ の
位置に $P(x)$ の原点を移動し、 4 個の近接根が原点を中心に半径 2 の円盤内に含まれるよう、 スケール変換と規格化を行う。 このことにより、 式 (2) に対応する多項式として次式を得る。
$\overline{P}(x)$$\tilde{P}(x)=-1.7578537183681$ $\mathrm{x}10^{-45}x^{14}+\cdots$ +4.0748102580441 $\mathrm{x}10^{-22}x^{9}$$+1.0683056800151\mathrm{x}10^{-17}x^{8}-$ $2.318794575^{\cdot}3916\mathrm{x}10^{-13}x^{7}-5.0470051845535\mathrm{x}10^{-9}x^{6}$(13)+5599830886856 $\mathrm{x}10^{-5}x^{5}+x^{4}-0.000750\mathrm{S}9624770413x^{3}-1.0x^{2}$$+0.0003524888550432x+0.090000002786884$ .
式 (13) のれる。

$\tilde{P}(x)$$\tilde{P}(x)$の 5 次以上の係数は、次数力$[searrow]\backslash ^{\backslash }1$増加するたびに約 10倍ずつ小さくなっていくのが見てとの 9 次以上の項は、係数の大きさは $10^{-20}$ より小さく、 反復式 (6) におけるは影響しない。

$\tilde{P}(z_{j}^{(\nu)})$

の計算に

31
$\overline{P}(x)$の 4 次以下の項からなる式の根を通常の Durand-Kerner 法で求めると、 反復計算は 8 回で収束し、根として$x=0.31640116584382-2.9387358770557\mathrm{x}10^{-39}\mathrm{i}$,$x=-0.94848132719237$ ,(14)$x=-0.31605441686571$ ,$x=0.94888547446197$を得る ( は虚数単位を表す)$\mathrm{i}$。

(これらの根は、 与多項式において、 それぞれ

$x=1.2000999999445-4.6465490204057\mathrm{x}10^{-43}\mathrm{i}$ ,$x=1.1999000044786$ ,(15)$x=1.1999999999448$ ,$x=1.2002000044862$
を表す。)

としては、
次に、式 (14) を初期値として近接根元 Durand-Kemer 法を実行する。反復計算に用いる式 (13) の 8 次以下の項のみを用いる。 同様に、 反復式 (6) の初期値に対し、 式 (11) において $\nu=0$ とおい$\overline{P}(x)$た多項式$Q\mathrm{o}(x)$として、 次に掲げる 4 次以下の項のみからなる式を用いる。
$Q\mathrm{o}(x)=1.0683056800151$$+$$\mathrm{x}10^{-17}x^{4}$$(-2.318794495173\cross 10^{-13}- 3 1394682295228 \mathrm{x}10^{-56}\mathrm{i})x^{3}$$+(-0.0000000050470053479879+6.8142252413433\mathrm{x}10^{-52}i)x^{2}$
$+$

(16)

$(0.000055998304846903+1 4832031290705 \mathrm{x}10^{-47}\mathrm{i})x$$+(1.0000000370019-1.645595346359\mathrm{x}10^{-43}\mathrm{i})$ .

この結果、 近接根用 Durand-Kemer 法は反復回数 3 回で収束し、 近接根として$x=1.2001000000000000-2.9323846931128\mathrm{x}10^{-72}\mathrm{i}$ ,$x=1.1999000000000000-7.3309617327821$$\cross 10^{-73}\mathrm{i}$

,

(17)

$x=1.2000000000000000$ ,$x=1.2002000000000000+2.9323846931128\mathrm{x}10^{-72}\mathrm{i}$

を得る。

5

実験

以上で述べた算法の効率を確かめるための実験を行った。 実験には SPARCStation 5 (microSPARC
SunOS 4 垣 4) 上で数式処理システム GAL を用 た。
, RAM
本実験では、 1 つの近接根クラスタをもつ多項式を係数を乱数で与え、 それらの近接根を計算した。多項式の次数は 10 次から 50 次まで 10 次おきに定め、それぞれの次数において多項式を 10 個ずつ生成した。近接根の多重度は 4、近接度は 10 程度とし、 近接根クラスタの中心は複素平面上の単位円盤内 (かつその虚部は実軸から近接度程度の範囲内) に乱数で与えた。
近接根用 Durand-Kemer 法の実行にあたっては、 4 で述べた初期値の選択 (a) および算法の効率化 (b)$\mathrm{I}\mathrm{I}$$\mathrm{t}_{\mathit{4}^{\mathrm{a}}}$$85\mathrm{M}\mathrm{H}\mathrm{z}$$32\mathrm{M}\mathrm{B},$を行った場合と行わなかった場合 (以下ではそれぞれ「効率化あり」「効率化なし」 と呼ぶ) の計算時間を

32

表 1: 効率化なし/ありのそれぞれの場合の, 各次数における 1 多項式あたりの平均計算時問 (単位は$\mathrm{m}\mathrm{s}$, 以下同様).

図 1: 表 1 のグラフ (実線/破線は, それぞれ効率化なし/ありの場合を表す. 以下同様) .

表 2: 効率化なしの場合における近接根用 D-K 法の平均反復回数と, 効率化ありの場合における初期値計算の低次項の平均用/近接油鼠 D-K 法の平均反復回数 (それぞれ$n2$ ), D-K 反復に用いた$n_{1},$

次数 (それぞれ$d_{1},$$d_{2}$$\overline{P}(x),$$Q_{\nu}’(x)$).

比較した。計算時間は、 近接根クラスタの中心を計算する近似無平方分解の時間を除き、 原点移動と近接
根回 Durand-Kerner 法の実行時間を測定した。 ガーベジコレクション等の時問を除いた実計算時間を測定し、各次数毎に多項式毎の計算時間の平均値を算出、 比較した。

5.1

実験 1

効率化なし/ありのそれぞれの場合の各次数における 1 多項式あたりの平均計算時間 (単位 :mse.c.) を
表 1 および図 1 に示す。
次に、効率化なし/ありのそれぞれの場合における Durand-Kerner 法の平均反復回数を表 2 に示す。表中、効率化なしの場合における反復回数は、 近接根用 Durand-Kerner 法の平均反復回数を表す。一方、 効率化ありの場合においては、 以下の記号を用いている。
$n_{1}$$n_{2}$$d_{1}$

初期値計算用 Durand-Kemer 法の平均反復回数
近接根用 Durand-Kerner 法の平均反復回数

初期値計算用 /近接根用 Durand-Kerner 法において、 反復公式の分子の計算に用いた P-(x). の低次項の平均次数$d_{2}$

近接根用 Durand-Kem er 法において、 反復公式の分母の計算に用いた$Q_{\nu}’(x)$の低次項の平均次数

効率化なしの場合は、入力多項式の次数が変化しても近接根用 Durand-Kerner 法の反復回数がほぼ一定である。一方、効率化ありの場合は、 入力多項式の次数の増大にかかわらず、 初期値計算用 / 近接根用

33

表 3: 効率化あり /なしのそれぞれの場合の, 各次
図 2: 表 3 のグラフ.

数における D-K 法の平均計算時間.

$\grave{\prime}J:\Re$$\# X_{}^{\ovalbox{\tt\small REJECT}}4\mathrm{b}t_{I}$$\llcorner$$ffl^{r}\ovalbox{\tt\small REJECT}\backslash \not\in \mathrm{b}h$

2.8
13.2
33.2
71.8
112.8

10
20
30
40
50

$\mp’\backslash \mathrm{J}\S ffi \mathrm{t}_{\sim}$$\mathfrak{h}$

3.8
14.0
34.8
69.8
116.8

3.3
13.6
34.0
70.8
114.8

表 4: 各次数におけ 6 原点移動時間. 効率化あり /なしのそれぞれの場合および両者の平均値を表す

図 3: 図 1 に表 4 の平均値を重ねたグラ 7. 灰色の部分は表 4 の平均値を示す.

Durand-Kerner 法の反復公式に用いられる多項式の次数が分母、 分子ともほぼ一定で、 かつ反復回数もほぼ一定であるので、 Durand-Kerner 法の計算量は、 入力多項式の次数にかかわらずほとんど変化しないは
図 1 の通り、全体の計算時間は入力多項式の次数の増加に伴って増ずである。 それにもかかわらず、 表$1_{\text{、}}$えている。
そこで、 Durand-Kerner 法の実行時間のみを取り出したところ、表 3 および図 2 の通りになった。 Durand-

Kerner 法の計算時問は表 2 の反復回数をほぼ反映したものである。
さらに計算時間の内訳を調べたところ、近接根クラスタを原点付近に移動する原点移動の時間が、計算時間全体において大きな割合を占めていることがわかった。結果を表 4 および図 3 に示す (図 3 の灰色で塗りつぶされた部分が、効率化なし/ありのそれぞれの場合の原点移動時間の平均値を示す) 。多項式の原点移動は GAL の多倍長浮動小数演算 (精度 10 進 34 桁) で行っており、多項式の次数が増大するほど原点移動の時間が計算時間全体に影響することがわかる。

52

実験 2:

$\mathrm{Q}\mathrm{D}$Library を用いた原点移動の効率化

そこで、原点移動における多倍長浮動小数演算に外部のライブラリを用いることにより、原点移動の効率化を図った。

Library は IEEE 倍精度
Library [3] と呼ばれるライブラリを用いた。
倍精度と呼ぶ)の浮動小数演算を C++で実装したの約 2 倍と約 4 倍の精度 (以下、 それぞれ 4 倍精度、 8

多倍長浮動小数演算には$\mathrm{Q}\mathrm{D}$$\mathrm{Q}\mathrm{D}$の他に Fortran-90 用のインターフェースが用意されている。 ソースコードは公開されており、 (著者らによる調査の限りライセンス条項は不明であるが) 無償で使用できる。
ライブラリで、$\mathrm{C}++$

Library から
本実験では、要求される浮動小数演算の精度が IEEE 倍精度の約 2 倍の精度であるため、
Library
4 倍精度のライブラリを用いた。 CAL は外部ライブラリの関数を呼び出すことができるので、$\mathrm{Q}\mathrm{D}$$\mathrm{Q}\mathrm{D}$34

表 5: 効率化なし/ありのそれぞれの場合の, 各次数における 1 多項式あたりの平均計算時間と原点
Library を使用) .
移動の平均時間 ($\mathrm{Q}\mathrm{D}$

図 4: 表 5 のグラフ. 灰色の部分は原点移動の平均時間を表す.

を用いて原点移動のライブラリを作成し、 GAL から同ライブラリを呼び出すことにより、 原点移動を行った。 そして、 実験 1 と同様の実験を行い、 計算時事を測定、 比較した。
実験結果を表 5 および図 4 に示す。本実験結果の通り、原点移動を効率化することにより、 4 で述べた効率化の本来の効果が現われたことがわかる。

参考文献[1] E. Durand. Solutions Numiriques des \’Equations Alg\’ebriques, Tome I. Masson, Paris, 1960.

[2] S. Fortune. An iterated eigenvalue algorithm for approximating roots of univariate polynomials.
Symbolic Comput., Vol. 33, No. 5, pp. 627-646, 2002.

.

$J$[3] Y. Hida, X. S. Li, and D. H. Bailey. Algorithms for quad-double precision floating point arithmetic. In
Proc. the 15th IEEE Symposium on Computer Arithrnetic (ARITH’01), PP. 155-162, 2001. Software.berkeley.
can be obtained from Hida’s home page: http$://\mathrm{w}\mathrm{w}\mathrm{w}.\mathrm{c}\mathrm{s}$$\mathrm{e}\mathrm{d}\mathrm{u}/-\mathrm{y}\mathrm{o}\mathrm{z}\mathrm{o}/$[4] I. O. Kerner. Ein Gesamtschrittverfahren zur Berechnung der Nullstellen von Polynomen. Nu$mer$. Math., Vol. 8, PP. 290-294, 1966.
[5] X.-M. Niu and T. Sakurai. A method for finding the zeros of polynomials using a companion matrix.
Japan J. Indust. AppL Math., Vol. 20, pp. 239-256, 2003.
[6] M.-T. Noda and T. Sasaki. Approximate GCD and its application to ill-conditioned algebraic equations. J. Comput. Appl. Math., Vol. 38, No. 1-3, PP. 335-351, 1991.
[7] V. Y. Pan, Univariate polynomials: nearly optimal algorithms for numerical factorization and rootfinding. J. Symbolic Comput., Vol. 33, No. 5, pP. 701-733, 2002.
[8] T. Sasaki and M.-T. Noda. Approximate square-free decomposition and root-finding of ill-conditioned
algebraic equations. J. Inform. Process., Vol. 12, No. 2, pp. 159-168, 1989.

[9] A. Terui and T. Sasaki. “Approximate zero-points” of real univariate poiynomial with large error
terms. 情報処理学会論文誌 Vol. 41, No. 4, pp. 974-989, April 2000.
[10] A. Terui and T. Sasaki. Durand-Kerner method for the real roots. Japan J. Indus. Appl. Math.,. 19 No. 1, . 19-38, January 2002.
$\mathrm{V}\mathrm{o}\mathrm{l}$$\mathrm{P}\mathrm{P}$[11] J.-C. Yakoubsohn. Finding a cluster of zeros of univariate polynomials. J. Compleity, Vol. 16,
No. 3, pp. 603-638, 2000.