東邦大学理学部情報科学科 2011 年度卒業論文 安定化手法のユークリッド互除法と スツルムのアルゴリズムへの適用 指導教員 : 白柳潔 提出日 :2012 年 2 月 24 日 提出者 :5508034 菊池裕也
2011 年度東邦大学理学部情報科学卒業研究 安定化手法のユークリッド互除法と スツルムのアルゴリズムへの適用 学籍番号 5508034 氏名菊池裕也 要旨 白柳らは 代数的アルゴリズムを近似計算したときに起こる不安定性に対し アルゴリズムの安定化手法を提案した 今日までに 安定化手法を数式処理システムで適用させる研究は行われてきているが 数式処理システム Maple14 を使って行った例はない 本研究では Maple14 で白柳らの安定化手法を適用させ その結果について考察を行う はじめに安定化手法とは 1. 入力の多項式に対し 係数の区間化を行う 2. アルゴリズムの構造は変えず 命令に従い区間演算を行いながら進行する 3.YES,NO の判定の際 区間に 0 が含まれる場合その区間を 0 とみなす ( ゼロ書き換え ) である 精度を上げながら再計算をすることにより 真の解を得るか またはそれに近づく解を得ることができる Maple14 では区間演算がサポートされていないので それらをできるようにした さらに その完成した区間演算をユークリッドの互除法とスツルムのアルゴリズムに適用させた 白柳研究室 i
目次 序論安定化手法について 第 1 章 第 2 章 第 3 章 第 4 章 1-1 ユークリッド互除法 (euclid,kukaneuclid) 1-2 区間演算の定義 1-3 係数の区間化 (kukan) 区間の定義直し (kukannaosi) ゼロ書き換え (zerohantei) 2-1 係数の区間演算 (tas,hiku,kakeru,waru) 2-2 多項式の区間演算 (tas2,hiku2,kakeru2,waru2,amari2) 2-3 ユークリッド互除法への適用結果 3-1 スツルムの定理 (suturumu) 3-2 符号の変わる値を求めるプログラム (hugou,kukanhugou) 3-3 区間多項式の常微分 (kukanbibun) 3-4 区間多項式のスツルムの定理 (kukansuturumu) 3-5 スツルムの定理への適用結果 まとめと考察 参考文献 ii
序論 コンピュータが数値計算を行うとき 分数は浮動小数点として計算されるこ とが多い ところが 浮動小数点で計算を行ったとき あるアルゴリズムでは 間違った出力を出してしまうことがある 今 割り切れるかという判定のアルゴリズムがあったとする 割り切れるか という条件に剰余が 0 になるかという定義を与える 入力 : [ 厳密計算 ] x 1 6 6 6x 1 6x 1 0 条件文では Yes という出力を返す しかし 浮動小数点計算では [ 近似値計算 ] 6 x 0.1667 6x 1 6 x 0 1..0002 0002 で No という間違った計算結果を出力してしまう 以上を安定化手法 ( 区間演算 ) を用いて計算すると [ 安定化計算 ] [6,6] [1,1] x [0.1666,0.1668] [6,6] x [1,1] [6,6] x [0.996,1.008] [ 0.004,0.008] という結果を出力する ここで白柳らの考えた手法は区間の中に 0 が含まれる 場合 0 とみなすもので この場合 区間 は 0 となる そして 剰余が 0 であるかという判定では Yes を返すことになる 以上の定義 が安定化手法の一例である 1
第 1 章 1-1 ユークリッド互除法 (euclid kukaneuclid) ユークリッド互除法とは任意の二つの式に対し 最大公約式を求める手法である 任意の二つの式 に対し f の g による剰余を r とすると f と g の最大公約式は g と r との最大公約式に等しいという性質が成り立つ この性質を利用して g を f r を g と置き換え計算を繰り返すと 剰余が 0 になったときの割る数が最大公約式となる 尚 プログラムでは f を a g を b と置いてある 下図の左図は基本となるユークリッド互除法のプログラム これを区間多項式のユークリッド互除法に移したときのプログラム ( 右図 ) を表している ユークリッド互除法のアルゴリズム 区間ユークリッドのアルゴリズム 以上のように基本となるユークリッド互除法のプログラムを作り 区間演算可能なプログラムへと拡張させる ここで必要なのは zerohantei と amari2 のプログラムである 剰余を求めるには区間多項式演算の和 差 積が必要となるのでこれらを求めるプログラムを構築し kukaneuclid に拡張させる 2
1-2 区間演算の定義 [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)] 区間演算を行う前にその他必要なプログラムを定義する 3
1-3 係数の区間化 (kukan) kukan プログラムは累乗根と分数に対して区間化を行うものである 例として という多項式に対して区間化を行う ここで重要な のは分数 または累乗根の近似値が区間の範囲として含まれていることである 実行結果 ( 整数の区間は手動で決めている ) kukan(f,n) で表されたプログラムは f は元になる数 n は精度桁の桁数の入力 を求められている このプログラムで重要なのは evalf で小数化を行った後 を足し引きし 上限と下限を決めているところにある 4
区間の定義直し (kukannaosi) Maple14 では区間演算のうち和に対してサポートを行ってくれることがわかっ た しかし 差に対しては次のような結果を出力してしまう [ 例 ] この出力は区間定義の a1<=a2 に反しているので正しい値が出力されたとは言 えず 又区間によっては差の計算が行われていないものもある よって区間の定義を直すプログラムが必要である 実行結果 このプログラムの利点は項ごとが + で分けられるところにもある このあとに 記述する多項式の区間演算でそれらは役立つ 5
ゼロ書き換え (zerohantei) 区間の中に 0 が含まれているか判別するプログラムである 主な判定方法は 下限か上限が 0 ならば 0 とみなす 又は下限が 0 以下かつ上限が 0 以上ならば 0 とみなすである 実行結果 zerohantei により 0 が含まれる項のみ消され 区間に 0 を含まない多項式が残った プログラム内では degree で最高次数を求め 高次数から while 文で下げていき必要な個所で判定をする方法を取ったが 後に記述する区間多項式の演算にも同様の方法で必要個所に応じ係数の区間演算を適用させる方法を採用する 6
第 2 章 2-1 係数の区間演算 (tas,hiku,kakeru,waru) 前述した区間演算の定義に習ったプログラムである まず係数間での定義を明 確にし 多項式へと拡張していく 和 差 積は以上の通りである しかし 商に至ってはより明確に定義する必要がある まず 分母に 0 が当てられる商は明らかに計算不可能なので 割る数の下限又は上限が 0 又は 0 を含む区間のとき計算不可能と定義する 次に 区間の商の定義に従うと区間の中に分数を含む結果が得られてしまう [ 例 ] この結果でも正しいのだが 後のプログラムで op の参照ができない不具合が生 7
じてしまうので 分数になった場合区間化を当てている また 精度桁は必ず 3 になるようにしている 実行結果 8
2-2 多項式の区間演算 (tas2,hiku2,kakeru2,waru2,amari2) 多項式の演算は全て最高次数から降順に見ていき 係数間で演算を行わせるものである 前述したように Maple14 では多項式の和についてはサポートされているのだが プログラムの基盤を考える上で一番取り組みやすいものが和と考えたので本研究では多項式の和についても考えることにした このプログラムでは まず二つの多項式のうち次数の高いものを a とする a の最高次数を t として t を while 文で下げていく その時 多項式 b の最高次 数 g と重なり合った箇所を始めとし 係数間で tas 演算を行い x に t 乗を掛け合 9
わせ answer に答えを代入していく 同時に g の値も下げていき 出力に answer を表示させる 又 係数に値が入っていない場合 Maple14 では coeff で 0 と表される為 if 文で係数 0 のときの場合分けも行った そのようにしなければ tas プロセスの際 op が参照できないエラーを生じてしまう 実行結果 [ 例 ] > 10
引かれる多項式を y=a と決めてしまう a の次数 t が b の次数 g よりも低い場 合 一回目の while 文に入り t=g となるまで処理を行う t>g の場合一回目の while 文には入らず 最下部の else まで飛び t=g となるまで処理を行う ( 処理 11
はプログラムを参照 ) 以下 t=g になった時は和と同様のプロセスで hiku 演算を 行わせている 実行結果 備考尚 例としての演算を行ったとする Maple14 での結果はが出力される しかしこの出力は正しくない なぜなら区間演算の差の定義では差 :[a1,a2]-[b1,b2]=[a1-b2,a2-b1] を表す よって正しい出力としては = が正しい出力にならなければいけない 今回 hiku 演算を作ったことで上記の定 義が満たされた 12
Maple14 の既存の関数 collect を使用したプログラム 2 つの多項式に対して collect 関数を用いると [ 例 ] 以上のように x の項でまとめてくれる これを使い y,z をでまとめ s についてプログラムを組む 注意点としては collect でまとめた際同じ係数では係数の累乗の形で表されるのと Maple14 では積の区間演算がサポー 13
トされていないので 例の第 2 項のような係数の積と和の形で表される点である 先の 係数の累乗で表されるものに対しては type(op(2,coeff(s,x,t)),`integer`) で二つ目の op が整数 すなわち区間でない係数の次数を表すものと定義をして kakeru 演算で y と z に同じ係数を定義して計算させる 第 2 項のような形で表されるものに対しては nops 関数でカッコ内の項数を数え 項ごとに kakeru 演算を行い answer に代入していき最後に再び collect 関数を使いまとめている 前述しているように Maple14 では和をサポートしているので係数どうしの積演算を行ってあげればまとめられる 実行結果 14
前述の kakeru2 hiku2 をプログラム内で用いている 注意は hiku2 での計算で例としてなど同じ係数の計算を行った際 [0,0] の区間として答えがでてしまう点である このプログラムでは剰余は必要としないので剰余を zerohantei で一掃している 実行結果 15
16
waru2 のプログラムをもとに answer 代入部を変えた剰余のプログラム waru2 では b3 に対し zerohantei で一掃を行ったが 今回は b3 の結果に同じ値を引いて消している [ 例 ] この演算では hiku2 を使っていない hiku2 ではが残ってしまう 又 この方法では多項式からある特定の項を消す方法が可能だが 一変数式の消去には 0 と出力され係数どうしの消去にはと残ってしまうので場合分けが必要である [ 例 ] Maple14 では degree(0) に対し を使い場合分け 実行結果 と定義している ( 他の整数は 0 と出力 ) これ 17
2-3 区間多項式ユークリッド互除法への適用結果 厳密計算 素朴な近似計算 安定化計算の比較を表した図である 実験 1 出力結果 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 実験 1は比較的簡素な多項式での実験結果 CPU 時間はどれも 0 でかからな かったが精度桁で変化が見られた 近似計算で 精度桁 2までは近似計算での正しい出力が見られたが 精度桁 3~6ではコンピュータによる桁落ちが生じ若干のずれが見られる 精度桁 7 以降では剰余が 0 という結果を出してしまい ユークリッド計算が一回目で終 了し厳密計算に近い結果が得られたが 不安定性を感じる 一方安定化計算では精度桁 2 以降で厳密計算に近い値を出力している 又 精度桁を上げることにより 真の値に近い出力結果を表している 18
実験 2 出力結果 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 実験 2は の形の多項式での実験結果 近似計算での精度桁 2~7では厳密計算とは明らかに違う値を示していて 正しい出力とは言えない 精度桁 8からは厳密計算に近い値を示しているが値 が収束せず不安定 安定化手法では精度桁 2から厳密計算に近い値を示している さらに 精度 桁を上げることにより 真の値に収束していることが実験からわかる CPU 時 間は多少かかるが気にならない程度である 19
実験 3 [, ] 厳密計算近似計算 出力結果 CPU 精 度 桁 係数膨張 0.062 +.. 0.015 3 +.. 0.031 10 +.. 0.031 100 安定化 計算 +.. 0.483 3 0.390 10 +.. 0.468 100 +.. 実験 3は 共通因子をもたずユークリッド互除法で係数膨張がおきてしまう例で行った 精度桁を増やすと CPU 時間も格段に増えてしまうと予想していたが結果は疎らで違った 近似計算では精度 10 桁と100 桁では変わらず この値が真の値であるかはわからない 安定化計算では低次で収束が見られる 20
第 3 章 3-1 スツルムの定理 (suturumu) スツルムの定理とは 重根を持たない実数係数の代数方程式の 任意に与えられた閉区間における根の個数を計算できるアルゴリズムを言う 1 多項式 y と閉区間 [w,z] を入力 2 多項式 y を f(0),y の常微分を f(1) とする 3 i:=0 として (rem は多項式 f(i),f(i+1) の x についての偏微分を表す ) 4 以下 剰余が 0 になるまでループ 5 f(i)(x) に閉区間の下限 w と上限 z を代入 6 f(i)(w) について符号の変わった回数を v1 f(i)(z) について符号の変わった回数を v2 7 が求める実数解の個数になる 以上の 1~4 で求めた f(i)(x) をスツルム列という 安定化手法を適用するため新たに区間の演算が必要な個所は2の微分化 6の符号の変わる個所の定義である まず基本となるスツルムの定理のプログラムと符号の変わる値を求めるプログラムを記述する Maple14 において diff(y,x) が y の x についての偏微分を表す 21
プログラム内の unapply は変数として扱っていた f(i) を x の関数 f(i)(x) に直す関数 である これで x に値を代入することができる 22
3-2 符号の変わる値を求めるプログラム (hugou kukanhugou) suturumu プログラムの後半で使うプログラム if 判定で符号の変わった場合に対 し +1 を加えている このプログラムを区間多項式上では以下のように書き換え る 判定方法は区間 y の上限と区間 z の下限で符号の変わったときの場合に対し +1 を加えている これは区間定義の a1<=a2 を満たすことから成り立っている 23
3-3 区間多項式の常微分 (kukanbibun) Maple14 では微分を行う際 diff 関数が備わっている しかし 区間多項式の微 分を行うには新たに関数の定義が必要である 実行結果 24
3-4 区間多項式のスツルムの定理 (kukansuturumu) 25
26
3-5 スツルムの定理への適用結果 実験 1 実験 2 27
28
第 4 章 まとめと考察 Maple14 への区間演算とそれらのアルゴリズムへの適用は成功したと言える また ユークリッド互除法への安定化手法について精度桁を変えることにより優位性を見られる結果を得られた しかし 係数膨張を生み出す不安定な多項式について実験を1つ行ったが 十分な結果は得られず 出力結果が真の値に近づいているのかがわからなかった スツルムのアルゴリズムへの適用には安定化手法の優位性を得られる明確な結果は得られなかった 実験研究期間からプログラムの準備に時間がかかりすぎて実験が満足に行えなかったのが悔やまれる また区間演算について 課題があるのは積の場合 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 年秋学期ゼミナールⅢ 29