多 倍 長 計 算 手 法 平 成 6 年 度 第 4 四 半 期 目 次. はじめに. assless 計 算. vtx-. vtx-.3 box 3.sinc 求 積 法 によるHadaard 有 限 部 分 計 算 3. SE,DEの 計 算 方 法 3. ε- 算 法 3.3 計 算 結 果 精 度 比 較 4. 多 倍 長 減 算 桁 落 ちメモ 4. 3 倍 精 度 4. 5 倍 精 度 4.3 6 倍 精 度 4.4 7 倍 精 度 5. 仮 数 部 のビット 数 と 乗 算 のメモ 6. DD,D(long double,dqのインライン 化 メモ 6. FORTRAN 6. C
.はじめに ファインマンループ 積 分 のassless 計 算 に 関 してその 発 端 から 数 学 的 根 拠 などをまとめ ました これに 関 連 して, 端 点 と 内 部 に 同 時 に 特 異 点 を 持 つ 積 分 に 対 してHadaardの 有 限 部 分 を sinc 求 積 法 と 今 までのε- 算 法 で 実 行 した 時 の 結 果 の 精 度 の 比 較 をしています また 非 対 称 行 列 反 復 計 算 の 調 査 で 作 成 した メモをまとめて 記 述 しています ( 整 数 演 算 での 減 算 の 桁 落 ち 計 算 ( 仮 数 部 のビット 毎 の 乗 算 またDD,D(longdouble,DQソースのインライ ン 化, 並 列 化 のためのメモをまとめています
.assless 計 算. vtx- inf x をs 5 x ra vtx dydx D D sxy (x y 計 算 しますがの 値 が 小 さくなると, 反 復 回 数 が 多 くなり 計 算 に 時 間 がかかったり, 計 算 できない DD形 式 の4倍 精 度 など,.5で 計 算 する 場 合, D D iとして ( x y 事 が 発 生 するので, dydx D sxy (x y D D i, D log( の 値 が 充 分 小 さくない と 結 果 の 精 度 が 良 くな い 事 があり ます の 計 算 において ( 以 下 の 式 で 計 算 する 事 があります この 場 合 は など 倍 精 度, 5 5 5 5 5 実 測 値 と 解 析 式 の 相 対 誤 差 実 数 部 6.3 6.3 3 5 6.68.6 7 6.4 7 9 虚 数 部.38.59 5.58.6.49 4 6 7 8
3 x 3 3 x 4 J I A sq( sq( log ]dq sq( sq( log( sq( [log( dq sq( sq( ( dpdq sq( ( p p log( y (x sxy ( J 4 B 4 b,4ac,c b, sq( a c b b 4ac arctan b 4ac a B sq( sq( log A B A dpdq p ( p q( sp p s dydx y x ( y (x sxy I 誤 差 数 学 的 には 以 下 のよう になります
] 6 s ( log s log( [log( s ] 3 s ( [log s s log( log( s dq sq( sq( log dq sq( log( dq sq( sq( log A s A の 様 にして 導 く 事 がで きます の 場 合 の 解 析 近 似 解 は 以 下 の 計 算 から
S<では inf ra vtx s 5 n,.5, log( n 実 測 値 解 析 近 似 解.7948469985D+.79484573467959D+ 3.569367859386D+.56936759564598D+ 4.3598964879377D+.35989455857D+ 5.5448954737D+.544895364498D+ 6.35387733349354D+.3538773377444D+ 7.356859554779D+.356859553763563D+ 8.4784769938477D+.478476898779D+ 9.45883988435886D+.4588398758834D+.589866955856D+.5898669496D+.55987884549476D+.55987884485545D+.67776443659367D+.67776434598D+ 3.666758659498D+.66675869858D+ 4.757358547554D+.757358678337D+ 5.763477997958D+.763477993768D+ 6.84369958957787D+.843699578968684D+ 7.8656873687386D+.86568736435634D+ 8.9666395545356D+.96663955453854D+ 9.96764674345479D+.96764673647465D+.796839694D+3.79683953397D+3.688665797D+3.6886685938D+3.97598979453D+3.975989465434D+3 3.7657548383768D+3.76575487699D+3 4.555766974748D+3.555766677965D+3 5.745394855668394D+3.74539485833D+3 6.33357459965D+3.333573888948D+3 7.3745397538889D+3.374539494876D+3 8.4548643486386D+3.45486483D+3 9.47646835994494D+3.47646835976765D+3 3.5694557854759D+3.5694557836544D+3 と 非 常 に 良 い 結 果 が 得 られています
. vtx- さらに, 計 算 すると, I x x D D dydx ( ( ( ( x x 計 算 出 来 ています 結 果 は 積 分 変 数 の 変 換 区 間 (, での 最 小 値 に 依 存 し, か5かに 依 存 では dydx と その 精 度 は進 桁 と進 4桁 ( 反 復 法 を 使 用 せず の 差 があり, 前 者 では4回 の 反 復 により進 8桁 となって います このケースで はD ありませんでした D xy の 場 合 を となり, 二 重 指 数 関 数 型 積 分 で 容 易 に 5 465.( log( dx B(, ( ( ( ( 指 数 部 のビット 数 が.7474 に 相 等 で D iとする 操 作 は 必 要
.3 box 少 し 複 雑 なケース log( I D xz y( x y z では, x x y I ( 4 ( ( ( log ( x log( x ( x これは 以 下 の 事 からき ています I [ a から [ a I - ( ( q a ( ( q D q ( 次 ぎの 様 になります p dzdydx とすると,a exp[ ( q ( ( p ( ( ( k ( k 4,a ( dq] [ ( 4,a (k ( k ( p ( ( dp] ( dq 4 k k p k (k ( x k 5 4( より となるがより 詳 細 には k ] k dp] 5 3 4
( ( ( ( q k の 定 数 項 は ( 3 ( 4 4 のの 一 次 の 項 は, 負 の 数 の 和 の 絶 対 値 で 計 算 すると, 分 子 の 項. (k!(k 分 子 の 項 k k k q (k これから( ( q ( 全 体 では q ( ( オイラーの 定 数 内 のの 一 次 の 項 の 係 数 は k (k!.57756649 [ k をかけて. 取 った 場 合 の 誤 差 とな リます k r dq.4857384 3 k 3 3... n A li n n / n / / n log A.48754473 log(k (k [ dq (k![log(k ] (k!(k (k (k k log A n q ( [ ]dq ( q ( ( e 6 n 4 ( (...( k (k!(k ( t ( log( ] ( 8 k (k![log(k ] r log(k (k ( t t k log( ] 8 8.8473 (k (k dt.34895684.5533568 となりこれが 一 次 式 ま で
実 測 例 としては 以 下 のものがあります 解 析 近 似 式 を 使 用 して, ~ 45 で 計 算 しました ただしgaa 関 数 はサポートされていない 機 種 もあり, 展 開 式 より 計 算 しています a= 4.65949445499595 a= 4.7756665456379 a= -.449346676876644536454 変 数 変 換 区 間 SR6[ [ epsilon算 法 の 収 束 具 合 は x557 x67 sr6, 5 で 3 ] の4倍 精 度 では 値,.9433984457834D + 7.943398444577885D + 7.94339844458948Q + 7 ], x557とx67では 誤 差.46333337879D - 5.7953883447D - 7.56658999Q - 5-5 で x557.98959569564375d + 9 x67.989595598969d + 9 sr6.98959568658934q + 9.99533977D -.477849985379986D -.4785749989Q - -5 で x557.9879837876d +.779848538333458D + 4 x67.98796344433d +.838476838465D + 3 sr6.9878644738q +.9874965Q + ( 注 SR6では 他 に 変 化 の 小 さい 要 素 がある となっています
あとSR6でのあるの 範 囲 でのa,a, a は 以 下 の 様 になっていま す rada=^{-6}~~[-3} a= 4.34474677645844344675746 a= 3.9986476938755466874576558365 a= -.639453446 rada=^{-3}~~[-45} a= 4.99577869564344849358836 a= 3.99939593353839597347984385565 a= -.3699579899999999999997596 rada=^{-46}~~[-6} a= 4.46479864546475 a= 3.99966395945374357968579 a= -.3589688999999999999979646
rada=^{-6}~~[-75} a= 4.967665987446494575 a= 3.999794993938866879338657374 a= -.378474887484836 rada=^{-76}~~[-9} a= 4.445986999759656 a= 3.999834997653644966 a= -.38486934638 rabda=^{-9}~^{-5} a= 4.36947898573396664484 a= 3.9999649489457975597983644 a= -.45458733756 の 値 が 小 さくなるほど 5 a は 3 4に 近 付 いています
3sin c求 積 法 によるHadaard有 限 部 分 計 算 端 点 と 領 域 内 部 で 同 時 に 特 異 点 をもつ 場 合 の 計 算 を, 数 理 解 析 研 究 所 講 究 録 第 79巻 99年 6-9 Hadaard有 限 部 分 積 分 に 関 する DE公 式 緒 方 秀 教, 杉 原 正 顯, 森 正 武 ( 東 京 大 学 工 学 部 物 理 工 学 科 の 資 料 から 以 下 の 問 題 を 測 定 しました f.p. F(x (x dx, F(z 3 4 厳 密 値 ( ( N F(x / '(x k k 公 式 I h k N (x k ( z 5 ( {cot[ ( ]F'( h h 4 4 z sin 4 '( F( } [ ( / h]
DE : x tanh( sinh(t (x sinh x s log(, (x log(s x ( x '(x x ( log( x t x SE : x tanh( (x log( x F'(z 3 5 4 4 ( z ( z x, となるt DE SE h 3. SE,DEの 計 算 方 法 x ( log( x y y t log( ( y log( t ln( の 場 合, DE : t 3.855,SE : t 74.667 6.5で 測 定 s '(x x ここでsinhの 逆 関 数 はサポートされていませんので, logを 使 用 した 式 を 使 っています
3. εー 算 法 f.p で 計 算 しています 8 n h.5,., n ( でn 5まで 計 算. 今 回 は 領 域 を 分 割 しな い 場 合 と, 領 域 を [, ],[,] と 分 割 した 場 合 を 計 算 しました 変 数 変 換 区 間 [, ], 領 域 分 割 なしでは, N 976ですが, 領 域 分 割 する 場 合, F(x dx (x li.9でオーバーフローを 防 ぐ 処 理 が 必 要 になります (x x3(i * cnt cntで 発 生 する 丸 め 誤 差 の ため N 976の 場 合, x 計 算 するとします ま た 方 法 もあります F(x (x i dx の 場 合,.q.,.q x.の 場 合 のみ 4.5 6 としてN 96とする
3.3 計 算 結 果 精 度 比 較 4つの 方 法 で 精 度 を 検 証 した 結 果 は 以 下 の 様 になっています 変 数 変 換 区 間 は 全 て[, ] 6 としています 相 対 誤 算 一 覧 表 λ DE( 分 割 なし DE( 分 割 あり sinc (DE sinc (SE -.9.4D-4.38Q-6.8469D-.85D- -.8.6763D-6.5489Q-.645D-4.4563D- -.7.565D-6.3596Q-4.757D-4.8D- -.6.43D-6.5Q-5.855D-5.43D- -.5.6756D-7.3Q-5.966D-6.73D-3 -.4.359D-7.449Q-7.654D-6.57D-3 -.3.333D-7.588Q-6.4696D-6.374D-3 -..34D-7.4Q-6.5978D-6.85D-3 -..65D-7.7873Q-7.845D-6.53D-3..47D-8.4837Q-.95D-6.58D-3..4D-7.85Q-7.7D-5.68D-3.3.76D-6.65Q-7.36D-6.5D-3.4.74D-6.45Q-7.7D-6.934D-4.5.96D-7.8Q-6.D-6.8D-4.6.75D-6.356Q-5.99D-6.7D-4.7.89D-6.59Q-4.585D-5.699D-4.8.44D-6.3936Q-.86D-4.563D-4.9.5647D-6.73Q-6.583D-.557D-4 今 回 の 結 果 ではsinc(DE 分 点 数 6で 最 も 適 し ているという 結 果 になりました
4. 多 倍 長 減 算 桁 落 ちメモ 整 数 演 算 方 式 では 加 減 算 において 内 部 処 理 が 減 算 の 場 合 の 桁 落 ちと 丸 め 処 理 が 特 に 問 題 となりますので, 作 成 時 使 用 しました メモを 紹 介 します 多 倍 長 整 数 減 算 IZ IX IY (IX IY の 結 果 が64ビット 整 数 変 数 の6ビットのみ 使 用 するとします 4. 3倍 精 度 配 列 数 IZ( IZ( IZ( IZ( IZ( IZ(ikk IZ( IZ( IZ(ikk 丸 め 位 置 IZ((: 桁 落 ち 8 (ikk 6 ik 丸 めなし ( ikk, ik 59 ik 桁 落 ちなし. 丸 め 位 置 IZ(( : 丸 め 後 IZ( 桁 落 ちなし. IZ( 桁 落 ち IZ( IZ( IZ( 桁 落 ち 丸 めなし 4. 5倍 精 度 5 4 配 列 数 3 IZ(3 IZ(3 4 6 6 IZ(3 丸 め 後 5 IZ(3 IZ(3 IZ(3 IZ(3.IZ(3 IZ(3 IZ(ikk IZ(3 IZ(3 6 6 IZ(3 IZ(ikk 丸 め 位 置 IZ((: 桁 落 ちなし 桁 落 ち 桁 落 ちなし. 桁 落 ち 46 (ikk 6 ik 丸 めなし ( ikk 3, ik 59 6 5 4 ik 丸 め 位 置 IZ(( : 桁 落 ち 丸 めなし
4.3 6倍 精 度 57 56 配 列 数 3 IZ(3 IZ(3 56 58 58 IZ(3 丸 め 後 57 IZ(3 IZ(3 IZ(3 IZ(3 58 IZ(3 桁 落 ち 56 IZ(3 IZ(3 IZ(ikk IZ(3 IZ(3 IZ(ikk 桁 落 ちなし. 桁 落 ち 78 (ikk 6 ik 丸 めなし ( ikk 3, ik 59 58 58 桁 落 ちなし. 丸 め 位 置 IZ((: 57 ik 丸 め 位 置 IZ(( : 桁 落 ち 丸 めなし 4.4 7倍 精 度 9 8 配 列 数 4 IZ(4 IZ(4 IZ(4 IZ(4 8 3 3 9 IZ(4 IZ(4 丸 め 後 IZ(4.IZ(4 IZ(4 IZ(ikk IZ(4 IZ(4 3 3 IZ(4 IZ(ikk 丸 め 位 置 IZ((: 桁 落 ちなし 桁 落 ち 桁 落 ちなし. 桁 落 ち (ikk 6 ik 丸 めなし ( ikk 4, ik 59 3 9 8 ik 丸 め 位 置 IZ(( : 桁 落 ち 丸 めなし
5. 仮 数 部 のビット 数 と 乗 算 のメモ 非 対 称 行 列 の 反 復 解 法 で 使 用 した 仮 数 部 の ビット 数 を 調 整 した 乗 算 のメモを 紹 介 します 仮 数 部 のビット 数 における 乗 算 の 場 合 の 桁 上 がりと 丸 め 多 倍 長 整 数 乗 算 IZ IX IY (IX が3ビット 整 数 変 数 の3ビットのみ 使 用 するとします 符 号 部 ビット, 指 数 部 5ビットを 想 定, IY の 結 果 (6ビット 配 列 数 5 IZ(5 IZ(5 (64ビット 配 列 数 5 IZ(5 IZ(5 (368ビット 配 列 数 5 IZ(5 IZ(5 (47ビット 配 列 数 5 IZ(5 IZ(5 (576ビット 配 列 数 6 IZ(6 IZ(6 桁 上 がりあり 桁 上 がりなし 9 9 7 7 5 5 3 3 桁 上 がりあり 桁 上 がりなし 桁 上 がりあり 桁 上 がりなし 桁 上 がりあり 桁 上 がりなし 桁 上 がりあり 桁 上 がりなし 丸 め 位 置 IZ(3( : 丸 め 位 置 IZ((9 : 9 丸 め 位 置 IZ(3(4 : 4 丸 め 位 置 IZ(3(3 : 3 丸 め 位 置 IZ(3(8 : 8 丸 め 位 置 IZ(3(7 : 7 丸 め 位 置 IZ(3( : 丸 め 位 置 IZ(3(: 丸 め 位 置 IZ(3(6 :6 丸 め 位 置 IZ(3(5 :5
6 DD,D(long double,dqのインライン 化 メモ E5-67,E5-66,Phi5P 等 でつの 変 数 の 和 で 表 す 演 算 を 使 用 した 時, 並 列 化 効 果 を 出 す ために 作 成 したソースのメモを 紹 介 します 6. FORTRAN ( 加 算 c=a+b ( 減 算 c=a-b p=a( p=a( q=b( q=b( t=p+q t=t-p t3=(p-(t-t+(q-t t4=p+q t5=t4-p t6=(p-(t4-t5+(q-t5 t7=t3+t4+t6 t8=t+t7 t9=t8-t t=(t-(t8-t9+(t7-t9 c(=t8 c(=t p=a( p=a( q=b( q=b( t=p-q t=t-p t3=(p-(t-t-(q+t t4=p-q t5=t4-p t6=(p-(t4-t5-(q+t5 t7=t3+t4+t6 t8=t+t7 t9=t8-t t=(t-(t8-t9+(t7-t9 c(=t8 c(=t
(3 乗 算 c=a*b p=a( p=a( q=b( q=b( t=p*q t=p*r a=t-(t-p a=p-a t3=q*r b=t3-(t3- b=q-b t4=((a*b-t+a*b+a*b+a*b t5=t4+p*q t6=t5+p*q t7=t6+p*q t8=t+t7 t9=t8-t t=(t-(t8-t9+(t7-t9 c(=t8 c(=t r r 34779 57 7 (4倍 精 度 ( 倍 精 度
(4 除 算 c=a/b p3=a( p4=a( q3=b( q4=b( ts=p3/q3 q=ts*q3 q=ts*q4 p=p3 p=p4 t=p-q t=t-p t3=(p-(t-t-(q+t t4=p-q t5=t4-p t6=(p-(t4-t5-(q+t5 t7=t3+t4+t6 t8=t+t7 t=ts t7=t8/q3 t8=t+t7 t9=t8-t t=(t-(t8-t9+(t7-t9 c(=t8 c(=t
6. C ( 加 算 c=work+work ( 減 算 work=one-xx p=work[] ; p=work[] ; q=work[] ; q=work[] ; t=p+q ; t=t-p ; t3=(p-(t-t+(q-t ; t4=p+q ; t5=t4-p ; t6=(p-(t4-t5+(q-t5 ; t7=t3+t4+t6 ; t8=t+t7 ; t9=t8-t ; t=(t-(t8-t9+(t7-t9 ; c[]=t8 ; c[]=t ; p=one[] ; p=one[] ; q=xx[] ; q=xx[] ; t=p-q ; t=t-p ; t3=(p-(t-t-(q+t ; t4=p-q ; t5=t4-p ; t6=(p-(t4-t5-(q+t5 ; t7=t3+t4+t6 ; t8=t+t7 ; t9=t8-t ; t=(t-(t8-t9+(t7-t9 ; work[]=t8 ; work[]=t ;
(3 乗 算 work3=work*work p=work[] ; p=work[] ; q=work[] ; q=work[] ; t=p*q ; t=p*r ; a=t-(t-p ; a=p-a ; t3=q*r ; b=t3-(t3- ; b=q-b ; t4=((a*b-t+a*b+a*b+a*b ; t5=t4+p*q ; t6=t5+p*q ; t7=t6+p*q ; t8=t+t7 ; t9=t8-t ; t=(t-(t8-t9+(t7-t9 ; work3[]=t8 ; work3[]=t ; r r 34779 8589934593 7 33 ( 倍 精 度 ( 拡 張 倍 精 度
(4 除 算 w3=work8/work p3=work8[] ; p4=work8[] ; q3=work[] ; q4=work[] ; ts=p3/q3 ; q=ts*q3 ; q=ts*q4 ; p=p3 ; p=p4 ; t=p-q ; t=t-p ; t3=(p-(t-t-(q+t ; t4=p-q ; t5=t4-p ; t6=(p-(t4-t5-(q+t5 ; t7=t3+t4+t6 ; t8=t+t7 ; t=ts ; t7=t8/q3 ; t8=t+t7 ; t9=t8-t ; t=(t-(t8-t9+(t7-t9 ; w3[]=t8 ; w3[]=t ;