東邦大学理学部情報科学科 2011 年度卒業論文安定化手法のユークリッド互除法とスツルムのアルゴリズムへの適用指導教員 : 白柳潔提出日 :2012 年 2 月 24 日提出者 : 菊池裕也

(1)

東邦大学理学部情報科学科

2011 年度卒業論文

安定化手法のユークリッド互除法と

スツルムのアルゴリズムへの適用

指導教員：白柳潔

提出日：２０１２年

2 月 24 日

提出者：５５０８０３４菊池裕也

(2)

i ２０１１年度東邦大学理学部情報科学卒業研究

安定化手法のユークリッド互除法と

スツルムのアルゴリズムへの適用

学籍番号５５０８０３４氏名菊池裕也要旨白柳らは、代数的アルゴリズムを近似計算したときに起こる不安定性に対し、アルゴリズムの安定化手法を提案した。今日までに、安定化手法を数式処理システムで適用させる研究は行われてきているが、数式処理システムMaple14 を使って行った例はない。本研究では、Maple14 で白柳らの安定化手法を適用させ、その結果について考察を行う。はじめに安定化手法とは１．入力の多項式に対し、係数の区間化を行う。２．アルゴリズムの構造は変えず、命令に従い区間演算を行いながら進行する。３．YES,NO の判定の際、区間に 0 が含まれる場合その区間を 0 とみなす。（ゼロ書き換え）である。精度を上げながら再計算をすることにより、真の解を得るか、またはそれに近づく解を得ることができる。 Maple14 では区間演算がサポートされていないので、それらをできるようにした。さらに、その完成した区間演算をユークリッドの互除法とスツルムのアルゴリズムに適用させた。白柳研究室

(3)

ii

コンピュータが数値計算を行うとき、分数は浮動小数点として計算されることが多い。ところが、浮動小数点で計算を行ったとき、あるアルゴリズムでは間違った出力を出してしまうことがある。今、割り切れるかという判定のアルゴリズムがあったとする。割り切れるかという条件に剰余が0 になるかという定義を与える。入力： [厳密計算] 0 1 6 6 1 6 6 1    x x x 条件文ではYes という出力を返す。しかし、浮動小数点計算では [近似値計算] 0002 . 0.0002 1 6 6 1 6 1667 . 0   x x x でNo という間違った計算結果を出力してしまう。以上を安定化手法（区間演算）を用いて計算すると [安定化計算] ] 008 . 0 , 004 . 0 [ ] 008 . 1 , 996 . 0 [ ] 6 , 6 [ ] 6 , 6 [ ] 1 , 1 [ ] 6 , 6 [ ] 1668 . 0 , 1666 . 0 [ ] 1 , 1 [    x x x という結果を出力する。ここで白柳らの考えた手法は区間の中に 0 が含まれる場合0 とみなすもので、この場合、区間は 0 となる。そして、剰余が 0 であるかという判定では Yes を返すことになる。以上の定義が安定化手法の一例である。

(5)

2

第１章

１－１ユークリッド互除法（euclid、kukaneuclid）

ユークリッド互除法とは任意の二つの式に対し、最大公約式を求める手法である。任意の二つの式、に対し f の g による剰余を r とすると、f と g の最大公約式は g と r との最大公約式に等しいという性質が成り立つ。この性質を利用して、g を f、r を g と置き換え計算を繰り返すと、剰余が 0 になったときの割る数が最大公約式となる。尚、プログラムではf を a、g を b と置いてある。下図の左図は基本となるユークリッド互除法のプログラム。これを区間多項式のユークリッド互除法に移したときのプログラム（右図）を表している。ユークリッド互除法のアルゴリズム区間ユークリッドのアルゴリズム以上のように基本となるユークリッド互除法のプログラムを作り、区間演算可能なプログラムへと拡張させる。ここで必要なのは、zerohantei と amari2 のプログラムである。剰余を求めるには区間多項式演算の和、差、積が必要となるのでこれらを求めるプログラムを構築し、kukaneuclid に拡張させる。

(6)

3

１－２区間演算の定義

[a1,a2]で表記されたものを区間といい、a1 を下限、a2 を上限と表し、区間は常にa1<=a2 を満たすものとする。区間演算とは下記で表せた計算のことを言う。以下は区間多項式の係数間の計算を表したものである。和：[a1,a2]+[b1,b2]=[a1+b1,a2+b2] 差：[a1,a2]-[b1,b2]=[a1-b2,a2-b1] 積：[a1,a2]×[b1,b2] =[max(a1*b1,a1*b2,a2*b1,a2*b2) ,min(a1*b1,a1*b2,a2*b1,a2*b2)] 商：[a1,a2]÷[b1,b2] =[max(a1*1/b1,a1*1/b2,a2*1/b1,a2*1/b2) ,min(a1*1/b1,a1*1/b2,a2*1/b1,a2*1/b2)] 区間演算を行う前にその他必要なプログラムを定義する。

(7)

4

１－３係数の区間化(kukan)

kukanプログラムは累乗根と分数に対して区間化を行うものである。例としてという多項式に対して区間化を行う。ここで重要なのは分数、または累乗根の近似値が区間の範囲として含まれていることである。実行結果（整数の区間は手動で決めている。） kukan(f,n)で表されたプログラムは、f は元になる数、n は精度桁の桁数の入力を求められている。このプログラムで重要なのはevalf で小数化を行った後、 _{を足し引きし、上限と下限を決めているところにある。}

(8)

5

区間の定義直し（kukannaosi）

Maple14 では区間演算のうち和に対してサポートを行ってくれることがわかった。しかし、差に対しては次のような結果を出力してしまう。 [例] – この出力は区間定義の a1<=a2 に反しているので正しい値が出力されたとは言えず、又区間によっては差の計算が行われていないものもある。よって区間の定義を直すプログラムが必要である。実行結果 – このプログラムの利点は項ごとが＋で分けられるところにもある。このあとに記述する多項式の区間演算でそれらは役立つ。

(9)

6

ゼロ書き換え（zerohantei）

区間の中に 0 が含まれているか判別するプログラムである。主な判定方法は、下限か上限が0 ならば 0 とみなす。又は下限が 0 以下かつ上限が 0 以上ならば 0 とみなすである。実行結果 zerohantei により 0 が含まれる項のみ消され、区間に 0 を含まない多項式が残った。プログラム内ではdegree で最高次数を求め、高次数から while 文で下げていき必要な個所で判定をする方法を取ったが、後に記述する区間多項式の演算にも同様の方法で必要個所に応じ係数の区間演算を適用させる方法を採用する。

(10)

7

第２章

２－１係数の区間演算（tas,hiku,kakeru,waru）

前述した区間演算の定義に習ったプログラムである。まず係数間での定義を明確にし、多項式へと拡張していく。和、差、積は以上の通りである。しかし、商に至ってはより明確に定義する必要がある。まず、分母に 0 が当てられる商は明らかに計算不可能なので、割る数の下限又は上限が0、又は 0 を含む区間のとき計算不可能と定義する。次に、区間の商の定義に従うと区間の中に分数を含む結果が得られてしまう。 [例] この結果でも正しいのだが、後のプログラムでop の参照ができない不具合が生

(11)

8 じてしまうので、分数になった場合区間化を当てている。また、精度桁は必ず３になるようにしている。実行結果

(12)

9

２－２多項式の区間演算

（tas2,hiku2,kakeru2,waru2,amari2）

多項式の演算は全て最高次数から降順に見ていき、係数間で演算を行わせるものである。前述したようにMaple14 では多項式の和についてはサポートされているのだが、プログラムの基盤を考える上で一番取り組みやすいものが和と考えたので本研究では多項式の和についても考えることにした。このプログラムでは、まず二つの多項式のうち次数の高いものを a とする。a の最高次数をt として、t を while 文で下げていく。その時、多項式 b の最高次数g と重なり合った箇所を始めとし、係数間で tas 演算を行い x に t 乗を掛け合

(13)

10 わせ、answer に答えを代入していく。同時に g の値も下げていき、出力に answer を表示させる。又、係数に値が入っていない場合、Maple14 では coeff で 0 と表される為、if 文で係数 0 のときの場合分けも行った。そのようにしなければ tas プロセスの際op が参照できないエラーを生じてしまう。実行結果 [例] >

(14)

11

引かれる多項式をy=a と決めてしまう。a の次数 t が b の次数 g よりも低い場

合、一回目の while 文に入り t=g となるまで処理を行う。t>g の場合一回目の

(15)

12 はプログラムを参照）以下t=g になった時は和と同様のプロセスで hiku 演算を行わせている。実行結果備考尚、例としての演算を行ったとする。 Maple14 での結果はが出力される。しかしこの出力は正しくない。なぜなら区間演算の差の定義では差：[a1,a2]-[b1,b2]=[a1-b2,a2-b1] を表す。よって正しい出力としては = が正しい出力にならなければいけない。今回hiku 演算を作ったことで上記の定義が満たされた。

(16)

13 Maple14 の既存の関数 collect を使用したプログラム。２つの多項式に対して collect 関数を用いると [例] 以上のように x の項でまとめてくれる。これを使い y,z をでまとめs についてプログラムを組む。注意点としては collect でまとめた際同じ係数では係数の累乗の形で表されるのと、Maple14 では積の区間演算がサポー

(17)

14 トされていないので、例の第２項のような係数の積と和の形で表される点である。先の、係数の累乗で表されるものに対しては、type(op(2,coeff(s,x,t)),`integer`) で二つ目のop が整数、すなわち区間でない係数の次数を表すものと定義をして kakeru 演算で y と z に同じ係数を定義して計算させる。第２項のような形で表されるものに対してはnops 関数でカッコ内の項数を数

え、項ごとにkakeru 演算を行い answer に代入していき最後に再び collect 関

数を使いまとめている。前述しているようにMaple14 では和をサポートしているので係数どうしの積演算を行ってあげればまとめられる。実行結果

(18)

15 前述のkakeru2、hiku2 をプログラム内で用いている。注意は hiku2 での計算で例としてなど同じ係数の計算を行った際[0,0]の区間として答えがでてしまう点である。このプログラムでは剰余は必要としないので剰余をzerohantei で一掃している。実行結果

(19)

(20)

17 waru2 のプログラムをもとに answer 代入部を変えた剰余のプログラム。waru2

ではb3 に対し zerohantei で一掃を行ったが、今回は b3 の結果に同じ値を引いて消している。 [例] この演算では hiku2 を使っていない。hiku2 ではが残ってしまう。又、この方法では多項式からある特定の項を消す方法が可能だが、一変数式の消去には 0 と出力され係数どうしの消去にはと残ってしまうので場合分けが必要である。 [例] Maple14 では degree(0)に対しと定義している（他の整数は 0 と出力）。これを使い場合分け。実行結果

(21)

18

２－３区間多項式ユークリッド互除法への適用結果

厳密計算、素朴な近似計算、安定化計算の比較を表した図である。実験１出力結果 CPU 精度桁厳密計算 0 近似計算 0.02000000000 0 2 0.001999999998 0 3 0.000199999998 0 4 0.000019999998 0 5 0.000001999998 0 6 x-0.1666667 0 7 x-0.16666667 0 8 x-0.166666667 0 9 x-0.1666666667 0 10 安定化計算 0 2 0 3 0 5 0 10 実験１は比較的簡素な多項式での実験結果。CPU 時間はどれも 0 でかからなかったが精度桁で変化が見られた。近似計算で、精度桁２までは近似計算での正しい出力が見られたが、精度桁３～６ではコンピュータによる桁落ちが生じ若干のずれが見られる。精度桁７以降では剰余が 0 という結果を出してしまい、ユークリッド計算が一回目で終了し厳密計算に近い結果が得られたが、不安定性を感じる。一方安定化計算では精度桁２以降で厳密計算に近い値を出力している。又、精度桁を上げることにより、真の値に近い出力結果を表している。

(22)

19 実験２出力結果 CPU 精度桁厳密計算 0 近似計算 0 2 0 3 0 4 0 5 0.000002718853583 0 6 0 7 0 8 0 9 0 10 安定化計算 0.016 2 0.016 3 0.016 5 0.078 10 実験２は、の形の多項式での実験結果。近似計算での精度桁２～７では厳密計算とは明らかに違う値を示していて、正しい出力とは言えない。精度桁８からは厳密計算に近い値を示しているが値が収束せず不安定。安定化手法では精度桁２から厳密計算に近い値を示している。さらに、精度桁を上げることにより、真の値に収束していることが実験からわかる。CPU 時間は多少かかるが気にならない程度である。

(23)

20 実験３ [ _, ] 出力結果 CPU 精度桁厳密計算係数膨張 0.062 近似計算 +….. 0.015 3 +….. 0.031 10 +….. 0.031 100 安定化計算 +….. 0.483 3 +….. 0.390 10 +….. 0.468 100 実験３は、共通因子をもたずユークリッド互除法で係数膨張がおきてしまう例で行った。精度桁を増やすとCPU 時間も格段に増えてしまうと予想していたが結果は疎らで違った。近似計算では精度１０桁と１００桁では変わらず、この値が真の値であるかはわからない。安定化計算では低次で収束が見られる。

(24)

21

第３章

３－１スツルムの定理（suturumu）

スツルムの定理とは、重根を持たない実数係数の代数方程式の、任意に与えられた閉区間における根の個数を計算できるアルゴリズムを言う。１、多項式y と閉区間[w,z]を入力。２、多項式y を f(0),y の常微分を f(1)とする。

３、 i:=0 として（rem は多項式 f(i),f(i+1)の x についての偏微分を表す。）４、以下、剰余が0 になるまでループ。５、 f(i)(x)に閉区間の下限 w と上限 z を代入６、 f(i)(w)について符号の変わった回数を v1 f(i)(z)について符号の変わった回数を v2 ７、が求める実数解の個数になる。以上の1～4 で求めた f(i)(x)をスツルム列という。安定化手法を適用するため新たに区間の演算が必要な個所は２の微分化、６の符号の変わる個所の定義である。まず基本となるスツルムの定理のプログラムと符号の変わる値を求めるプログラムを記述する。Maple14 において diff（y,x）が y の x についての偏微分を表す。

(25)

22 プログラム内のunapplyは変数として扱っていたf(i)をxの関数f(i)(x)に直す関数である。これでxに値を代入することができる。

(26)

23

３－２符号の変わる値を求めるプログラム

(hugou、kukanhugou）

suturumuプログラムの後半で使うプログラム。if判定で符号の変わった場合に対し+1を加えている。このプログラムを区間多項式上では以下のように書き換える。判定方法は区間y の上限と区間 z の下限で符号の変わったときの場合に対し+1 を加えている。これは区間定義のa1<=a2 を満たすことから成り立っている。

(27)

24

３－３区間多項式の常微分（kukanbibun）

Maple14 では微分を行う際 diff 関数が備わっている。しかし、区間多項式の微分を行うには新たに関数の定義が必要である。実行結果

(28)

25

(29)

(30)

27

３－５スツルムの定理への適用結果

実験１実験２

(31)

28

(32)

29

第４章

まとめと考察

Maple14 への区間演算とそれらのアルゴリズムへの適用は成功したと言える。また、ユークリッド互除法への安定化手法について精度桁を変えることにより優位性を見られる結果を得られた。しかし、係数膨張を生み出す不安定な多項式について実験を１つ行ったが、十分な結果は得られず、出力結果が真の値に近づいているのかがわからなかった。スツルムのアルゴリズムへの適用には安定化手法の優位性を得られる明確な結果は得られなかった。実験研究期間からプログラムの準備に時間がかかりすぎて実験が満足に行えなかったのが悔やまれる。また区間演算について、課題があるのは積の場合collect 関数を使ってまとめたが全ての場合について網羅できたわけではなく、特別な条件下でエラーが起こってしまう場合も否定できない。他のwaru2、amari2 についても同様である。今回、高次から while 文で下して係数を見ていく手法を採用したが、他にも手法は無数にあると考えられる。例として、商を求める場合quo 関数を使うなど。最後に、今回の研究にあたり白柳先生また同研究室メンバーに力添えをしてもらったことを感謝いたします。

参考文献

・ K.Shirayanagi and M.Sweedler:A Theory of Stabilizing Algebraic Algorithms, Technical Report 95-28,Mathematical Sciences Institute,Cornell University,pp.1-92,1995. ・白柳潔(2003)『コンピュータのカオスをおさえる、新しい「安定」計算術』 NTT コミュニケーション科学基礎研究所監修・吉沢香織(2008)「東海大学理学部情報数理学科 2007 年度ゼミナールⅢレポート」・鮎田哲郎(2009)「東海大学理学部情報数理学科 2009 年秋学期ゼミナールⅢ」

東邦大学理学部情報科学科 2011 年度卒業論文安定化手法のユークリッド互除法とスツルムのアルゴリズムへの適用指導教員 : 白柳潔提出日 :2012 年 2 月 24 日提出者 : 菊池裕也