π(円周率)を小数以下10,000桁まで計算しなさい。

(1)

プログラミング言語 I 第 6 回 πの計算 1

埼玉大学工学部電気電子システム工学科伊藤和人

課題

π(円周率)を小数以下10,000桁まで計算しなさい。

マーチンの公式を利用

⎟ ⎠

⎜ ⎞

⎝ + ⎛

⎟ ⎠

⎜ ⎞

⎝

− ⎛

⎟ ⎠

⎜ ⎞

⎝ + ⎛

⎟ ⎠

⎜ ⎞

⎝

= ⎛

110443 arctan 1 239 48

arctan 1 57 20

arctan 1 49 128

arctan 1 π 48

⎟ ⎠

⎜ ⎞

⎝

− ⎛

⎟ ⎠

⎜ ⎞

⎝

= ⎛

239 arctan 1 5 4

arctan 1 π 16

他にも、以下のような公式あり

⎟ ⎠

⎜ ⎞

⎝ + ⎛

⎟ ⎠

⎜ ⎞

⎝

− ⎛

⎟ ⎠

⎜ ⎞

⎝ + ⎛

⎟ ⎠

⎜ ⎞

⎝

= ⎛

12943 arctan 1 682 96

arctan 1 239 48

arctan 1 57 28

arctan 1 π 176

πの計算

arctan の計算 ⇒ テイラー展開を利用

+ L

− +

−

− =

= ^∞ − ⁻

=

∑ ⁽ ₂ ¹ ⁾ ₁ + ₃ ₅ ₇

arctan

7 5 3 1

2 1

1 x x x

x i x

x ⁱ

i i

= S

⋅ +

⋅ −

⋅ +

−

⎟ =

⎠

⎜ ⎞

⎝

⎛ 3 5 7 5 7 L

16 5 5

16 5 3

16 5 16 5 arctan 1 16

⎟ の場合

⎠

⎜ ⎞

⎝

⎛ 5 arctan 1 16

テイラー展開による arctan 計算

5 16

1 =

T S = T 1

25 1 3

T = T

3 T 3

S S = − 25

3 5

T = T

5 T 5

S S = + 25

5 7

T = T

7 T 7

S

S = − T _n が十分小さくなるまで繰り返す除算、加算、減算を実行

= S

⋅ +

⋅ −

⋅ +

−

⎟ =

⎠

⎜ ⎞

⎝

⎛ 3 5 7 5 7 L

16 5 5

16 5 3

16 5 16 5 arctan 1 16

高精度演算の考え方

小数以下10,000桁 = 10 ^-10000 の精度

細かな数値を表すための C 言語データ型

±1.175×10 ^-38 ～±1.175×10 ⁺³⁸ (32ビット、有効桁数7)

float double long double

±2.225×10 ^-308 ～±2.225×10 ⁺³⁰⁸ (64ビット、有効桁数15)

±1.084×10 ^-4932 ～±1.084×10 ⁺⁴⁹³² (80ビット、有効桁数19)

long doubleを使用しても精度が不足しかも、19桁分しか値を記憶できない

別の方法が必要

桁数の長い数=多倍長数

多数のデータ(各桁の値)をまとめて処理

整数型配列の各要素が桁の値を記憶

要素数を増やせば長い桁に対応可能配列を利用

3.141592…

3 1 4 1 5 9 2 ・・・・ ..

int pi[] pi[0]

pi[1] pi[2]

pi[3] pi[4]

pi[5]

.

pi[6] pi[10000]

(2)

長い桁数の演算

C 言語では int 型データが基本演算単位

1桁の演算と2桁の演算はどちらもint型

⇒ 計算時間は同じ ( 同じ 1 回の演算処理 )

1 回の int 型演算でなるべく多くの桁を計算させたい

1つの要素がM桁を表す場合、

乗算の過程で桁数が最大で 2M 桁になる

int型で2M桁の整数が扱えるような桁数Mを

選択

長い桁数の演算

int 型ではなく long 型を指定その理由は …

現在の 32 ビットプロセッサでは int 型と long 型と同じ

⇒long型演算を最短時間(1クロック)で実行

もし 16 ビットプロセッサを使用したとしても、

long型=32ビットとして正しく実行できる

要素の桁数の決定

long 型で扱える値の範囲

long 型の最大桁数と M の関係

M として 4 を選択

⇒ 1つの要素が0から9999(4桁)を表す符号付 : -2,147,483,648 ～ 2,147,483,647 符号なし : 0 ～ 4,294,967,295(=2 ³² -1)

(2 ³¹ -1) の十進数桁数 = 9 ≧ 2M (-2 ³¹ ～ 2 ³¹ -1)

多倍長数と配列

小数点以下 10000 桁を表すには 10000/4=2500 要素の配列が必要

桁と配列要素の対応少なくとも 2 通りあり

3. 1415 9265 3589 7932 3846 … 1989 3

2 4 249

1 0 配列添え字

246 247 245 0

248 249 配列添え字

π = 方法1 方法2

方法2を採用除算(後述)の過程で、整数部分に拡張する場合あり方法1では拡張不可能

以下は250要素の場合

多倍長加算

下位桁から順に処理する void LVadd( long a[], long b[] ){

long i, c=0;

for( i=0 ; i<N ; i++ ){

a[i] += (b[i]+c);

if( a[i] >= D ){

a[i] -= D;

c = 1;

}else{

c = 0;

} } }

#define N 2500

#define D 10000 として、

加算結果がD以上ならば桁上げ

a ← a+b

多倍長減算

下位桁から順に処理する void LVsub( long a[], long b[] ){

long i, c=0;

for(i=0 ; i<N ; i++ ){

a[i] -= (b[i]+c);

if( a[i] < 0 ){

a[i] += D;

c = 1;

}else{

c = 0;

} } }

減算結果が負ならば桁借り

a ← a-b

#define N 2500

(3)

多倍長精度除算

被除数(割られる数): a = a _s D ^s + a _s

₋₁

D ^s

⁻¹

+ L + a

₀

除数(割る数): b = b _t D ^t + b _t

₋₁

D ^t

⁻¹

+ L + b

₀

a<bの場合は除算は不要(q=0)なので a≧bとするすなわち

商: q = q _u D ^u + q _u

₋₁

D ^u

⁻¹

+ L + q

₀

余り: r = r _t D ^t + r _t

₋₁

D ^t

⁻¹

+ L + r

₀

r qb

a = + ( 0 ≤ r < b )

0 ,

0 ≠

≠ _t

s b

a

多倍長精度除算

0 1

1

D b

b D b

kb = _t ′ ^t + _t ′

₋

^t

⁻

+ L + ′

bをk倍した数を考える

D D b

t ′ <

2 ≤

ただし k はを満たす整数

計算を効率化するための工夫

a も k 倍すれば ( ) ( ) ( ) ka = q kb + kr であるから、

z商は不変

z 余りは除算後に 1/k 倍すればよい以降では、予め被除数、除数とも k 倍し、

D D b

t <

2 ≤ が成り立つと仮定

多倍長精度除算 (3)

商の最上位桁位置uおよび桁値q _u を求める

の場合

のとき

t

s b

a ≥

t

bD s

a ≥

⁻

D a D D b

s

t < <

≤ ,

2 であるから、 < 2 ⇒ < 2

t t s

s b

b a a さらに ≥ 1

t s

b

a であり、商の最上位桁は 1 に定まるよって u = s − t , q _u = 1

多倍長精度除算 (4)

の場合(続き)

^s ^t

のとき

b a ≥

t

bD s

a <

⁻

1

1 −

−

< _t

s b

a とすれば a > bD ^s

⁻

^t − D ^s

よって u = s − t − 1 , q _u = D − 2 または D − 1

j t j s t s t

s b a b a b

a = ,

₋₁

=

₋₁

, K ,

₋

<

₋

( j ≥ 1 )

t t

t D D

D b

b ≥ ≥ 2 より 2 b ≥ D ^t

⁺¹

⇒ 2 bD ^s

⁻

^t

⁻¹

≥ D ^s すなわち a > bD ^s

⁻

^t − D ^s ≥ bD ^s

⁻

^t − 2 bD ^s

⁻

^t

⁻¹

= bD ^s

⁻

^t

⁻¹

( D − 2 )

( ^j ⁼ ¹ )

多倍長精度除算(5-1)

の場合

まず、簡単化のための場合を考える

t

s b

a <

r Qb a D

a _s + _s

₋₁

= _t + とすれば ⎥

⎦

⎢ ⎥

⎣

⎢ +

≤

⁻

t s s

b a D

Q a

¹

b Q a b R

a D Q a

t s s

0 0

1

0

⎥ , = −

⎦

⎢ ⎥

⎣

⎢ +

=

⁻

とおく

+ 1

= t s

Q ₀ は商の候補であるが、除数を小さく見積もっているため、 Q ₀ が実際の商 q より大きい、

すなわち ( b _t D R ^t ₀ ≤ が負になっている場合も存在する b )

商の候補を 1 増やしながら、余りが非負に

多倍長精度除算(5-2)

a _s < b _t の場合(続き)

0

, r R

Q

q ← ← とし、

を繰り返す b r r qb a = + , <

b Q a b R

a D Q a

t s

s

1 0 0

0

⎥ , = −

⎦

⎢ ⎥

⎣

⎢ +

=

⁻

とおく

≥ 0

r となるまで b

r r q

q ← − 1 , ← +

だから

− 1

>

⇒ +

< b

q a b qb a このとき

+ 1

= t

s の場合

繰り返し回数

⎟ ⎠

⎜ ⎞

⎝ ⎛ −

−

⎟ ≤

⎠

⎜ ⎞

⎝ ⎛ − + −

<

−

⁻¹

1

₋₁

1

0

b

a D b

a b

a D q a

Q _s

t t

s s

※1

(4)

多倍長精度除算 (5-3)

a _s < b _t の場合 ( 続き ) s = t + 1 の場合

すなわち※1の操作は、高々2回繰り返せば正しい商が得られる

1

₁

1

0 1

⎟⎟ +

⎠

⎜⎜ ⎞

⎝

= ⎛ −

⎟ ⎠

⎜ ⎞

⎝ ⎛ −

−

<

−

₋

_s

₋⁻

t s t s

t b D

D b b b a b

a D b q a

Q (A)

+ 1

= t

s および b − b _t D ^t = b _t

₋₁

D ^t

⁻¹

+ L + b

₀

< D ^t より 3 1 1

1

⎟⎟ ⎠ + = + ≤

⎜⎜ ⎞

⎝

< ⎛

⎟⎟ +

⎠

⎜⎜ ⎞

⎝

⎛ −

−

t t

t t s

t s t

b D D

b D D D

b D b b b

a (B)

式(A),(B)より、繰り返し回数 Q

₀

− q < 3

よって u = s − t − 1 = 0 , q

₀

= Q

₀

または Q

₀

− 1 または Q

₀

− 2

多倍長精度除算 (5-4)

a _s < b _t の場合 ( 一般の場合 )

0 1

0

D , r R

Q

q ′ ← ^s

⁻

^t

⁻

′ ← とし、

を繰り返す r

b QD r b q

a = ′ + ′ = ^s

⁻

^t

⁻¹

+ ′ b D Q a b R

a D

Q a ^s ^t

t s

s

1

0 0

1

0 ⁻

⎥ , = −

⁻⁻

⎦

⎢ ⎥

⎣

⎢ +

= とおく

≥ 0

r ′ となるまで b

D r r D q

q ′ ← ′ − ^s

⁻

^t

⁻¹

, ′ ← ′ + ^s

⁻

^t

⁻¹

だから

1

+ ⇒ > −

<

⁻⁻ ⁻⁻ ₋₋

b D Q a b D b QD

a ^s ^t ^s ^t _s _t

このとき

繰り返し回数

⎟ ⎠

⎜ ⎞

⎝

⎛ −

+ −

<

−

⁻¹ ₋₋₁

1

0

D b

a b

a D Q a

Q _s _t

t s s

※1

多倍長精度除算 (5-5)

a _s < b _t の場合 ( 一般の場合 ) ( 続き ) 1

1

₁ ₁

1

⎟ ≤ − +

⎠

⎜ ⎞

⎝

⎛ −

+ −

−

b D

a D

b a b

D a b

a D a

t s s t t

s t

s s

1

₁

1

⎟⎟ ⎠ +

⎜⎜ ⎞

⎝

= ⎛ − +

−

=

₋₋ ₋₋ ₋₋

_t

t t t t

s t

s t s t

t b D

D b b b D

a b

D a D

D b

a

3 1 1 = + ≤

⎟⎟ +

⎠

⎜⎜ ⎞

⎝

< ⎛

t t

b D D

よって、繰り返し回数 Q

₀

− Q < 3

2 1

,

1 =

₀ ₀

−

₀

−

= s t q Q Q Q

u _u またはまたは

多倍長精度除算 (6)

商の最上位桁位置uおよび桁値q _u を求める

の場合

のとき

の場合

t

s b

a ≥

t

bD s

a ≥

⁻

1 , =

−

= s t q _u u

t

bD s

a <

⁻

1 2

,

1 = − −

−

= s t q D D

u _u または

t

s b

a <

2 1

,

1 =

₀ ₀

−

₀

−

= s t q Q Q Q

u _u またはまたは

⎥ ⎦

⎢ ⎥

⎣

⎢ +

=

⁻

t s s

b a D

Q

₀

a

¹

ただし

多倍長精度除算(7)

商の最上位桁位置 u および桁値 q _u が得られたら

b D q a r = − _u ^u a ← r として再度除算を実行

中間余りを計算

a < b になったとき、除算完了

多倍長除算の例

D=10 として 41567÷32 を行う

0. D b D になるように被除数a、除数bをk倍する

t <

2 ≤

a=83134=8×10 ⁴ +3×10 ³ +1×10 ² +3×10+4 k=2 として

b= 64= 6×10+4 1. 83134÷64 ⇒

t

s b

a ≥ , a ≥ bD ^s

⁻

^t より ,

= 83134 – 64000=19134

= 3

−

= s t

u q _u = q

₃

= 1

3 3

× 10

×

−

← a b q a

6 , 1 , 8 ,

4 = = =

= a _s t b _t

s

(5)

多倍長除算の例 ( 続き 1)

2. 19134÷64 ⇒ より,

= 19134 – 12800=6334

t

s b

a < 3

6

1

19

0

=

⎥⎦ ⎥

⎢⎣ ⎢

⎥ =

⎦

⎢ ⎥

⎣

⎢ +

=

⁻

t s s

b a D Q a

64×3×10 ² = 19200 > 19134

64×9×10 = 5760 ≦ 6334

・・・NG

・・・OK 3. 6334÷64 ⇒

t

s b

a ≥ , a < bD ^s

⁻

^t より

= 6334 – 5760=574 64×2×10 ² = 12800 ≦ 19134 ・・・OK

6 , 1 , 1 ,

4 = = =

= a _s t b _t

s 2 1 =

−

= s t

u ,

6 , 1 , 6 ,

3 = = =

= a _s t b _t

s

2 2

× 10

×

−

← a b q a

1 1 =

−

= s t

u , q _u = q

₁

= 9 :

1

= 9 q

:

2

= 3 q

:

2

= 2 q

1 1

× 10

×

−

← a b q a

多倍長除算の例 ( 続き 2)

4. 574÷64 ⇒ より

= 574 – 512=62

t

s b

a < 9

6

1

57

0

=

⎥⎦ ⎥

⎢⎣ ⎢

⎥ =

⎦

⎢ ⎥

⎣

⎢ +

=

⁻

t s s

b a D Q a

64×9 = 576 > 574 ・・・NG 64×8 = 512 ≦ 574 ・・・OK

q ₃ = 1, q ₂ = 2, q ₁ = 9, q ₀ = 8 より

41567÷32 = 1298 あまり 62÷2=31 6 , 1 , 5 ,

2 = = =

= a _s t b _t

s

0 1 =

−

= s t

u ,

0 0

× 10

×

−

← a b q a

:

0

= 9 q

:

0

= 8 q

多倍長除算の準備 (1)

良く使う関数を準備

void LVcopy( long a[], long b[] ) { int k;

for( k=N2-1 ; k>=0 ; k-- ) b[k] = a[k];

} long GetMSD( long a[] ) { int k;

for( k=N2-1 ; k>=0 ; k-- ){

if( a[k] ) break;

} return k;

}

多倍長数aの最上位桁位置を求める多倍長数aを多倍長数bへコピー

上位桁から順に調べ非ゼロならば最上位

多倍長数がゼロならば戻り値は-1

整数部を考慮

π は整数部 (=3) あり

演算の過程で、整数部分に桁上げする場合あり

整数部も統一的に扱うため、整数部まで桁を拡張する

3. 1415 9265 3589 7932 3846 … 1989 246

247 245 0

248 249 拡張した部分

π =

N=250要素の場合 250

251 余裕をもって、配列2要素分(8桁)拡張

#define N2 (N+2)

多倍長除算の準備(2)

多倍長数に単精度数を乗じる

void LVmulS( long a[], long nb, long m[] ) { int k;

long w, nCarry = 0;

long nMSDa = GetMSD(a);

for( k=0 ; k<=nMSDa ; k++ ){

w = a[k]*nb + nCarry;

m[k] = w % D;

nCarry = w/D;

} m[k] = nCarry;

for( k++ ; k<N2 ; k++ ) m[k] = 0;

aの最上位桁位置を求める m ← a×nb

w≧Dならば桁上げ下位からの桁上げを加算

多倍長除算の準備(3)

除数の最上位桁をD/2以上、D未満にする long Normalize( long b[], long c[] )

{ long ns, v;

long nMSD = GetMSD(b);

v = b[nMSD];

if( v == 0 ) exit(1);

ns = D/(v+1);

if( ns > 1 ) LVmulS( b, ns, c );

else LVcopy( b, c );

return ns;

}

bの最上位桁位置を求める bを正規化しcに代入

最上位桁が0なら終了

nsはb[nMSD]をD/2以上にするために

(6)

多倍長除算 (1)

long LVdiv( long a0[], long b0[], long q[] ) { long s, t, k, u, nk;

long a[N2], b[N2];

nk = Normalize( b0, b );

LVmulS( a0, nk, a );

t = GetMSD( b );

for( k=0 ; k<N2 ; k++ ) q[k] = 0;

while(1){

if( Compare( a, b ) < 0 ) break;

// a < b ならば除算終了

… } …

q ← a0/b0 bを正規化(nk倍) aもnk倍

商を0に初期化

多倍長除算 (2)

s = GetMSD( a );

if( a[s] >= b[t] ){

if( CompareOffset( a, b, s-t ) >= 0 ){

u = s-t; q[u] = 1;

}else{

u = s-t-1; q[u] = D-1;

} }else{

u = s-t-1; q[u] = (a[s]*D+a[s-1])/b[t];

}

LVmulS( b, q[u], m );

while( CompareOffset( a, m, u ) < 0 ){

q[u]--;

LVmulS( b, q[u], m );

}

商の

最上位桁位置u および値q[u]の候補を求める

a≧m=b×q[u]を満たす最大のq[u]を求める

t を検査

bD s

a ≥

⁻

多倍長除算 (3)

while( 1 ){

if( Compare( a, b ) < 0 ) break;

// a < b ならば除算終了

…

LVmulS( b, q[u], m );

while( CompareOffset( a, m, u ) < 0 ){

q[u]--;

LVmulS( b, q[u], m );

} LVsubOffset( a, m, u );

} a ← a – b×q[u]*D ^u

a≧m=b×q[u]を満たす最大のq[u]を求める u, q[u]の候補を求める(省略)

多倍長除算の準備 (4)

int Compare( long a[], long b[] ) { int k;

long nMSDa = GetMSD(a);

long nMSDb = GetMSD(b);

if( nMSDa > nMSDb ) return 1;

if( nMSDa < nMSDb ) return -1;

for( k=nMSDa ; k>= 0 ; k-- ){

if( a[k] > b[k] ) return 1;

if( a[k] < b[k] ) return -1;

} return 0;

}

戻り値 a=b ⇒ 0 a>b ⇒ 1 a<b ⇒ -1

最上位桁が同じ場合

多倍長数 a と b を比較

すべて等しい ⇒ a=b

aの最上位桁が bの最上位桁より上なら a>b

多倍長除算の準備(5)

int CompareOffset( long a[], long b[], long No ) { int k;

long nMSDa = GetMSD(a);

long nMSDb = GetMSD(b)+No;

if( nMSDa > nMSDb ) return 1;

if( nMSDa < nMSDb ) return -1;

for( k=nMSDa ; k>=No ; k-- ){

if( a[k] > b[k-No] ) return 1;

if( a[k] < b[k-No] ) return -1;

} for( ; k>=0 ; k-- ) if( a[k] ) return 1;

return 0;

}

多倍長数aとbを桁をずらして比較

a[No-1]から a[0]の中に非ゼロのものがあれば aの方が大きい

ここに来たとき

a[nMSDa],...,a[No]

と

b[nMSDa-No],...,b[0]

が一致

多倍長除算の準備(6-1)

void LVsubOffset( long a[], long b[], long No ) { int k;

long v = 0;

long nMSDb = GetMSD(b);

for( k=0 ; k<=nMSDb ; k++ ){

if( a[k+No] >= b[k]+v ){

a[k+No] –= (b[k]+v); v = 0;

}else{

a[k+No] += D;

a[k+No] –= (b[k]+v); v = 1;

} }

a ← a – b×D ^No

桁借りがある場合

多倍長数aからbを桁をずらして引き算

最上位桁までbを引く

次のスライドへ続く

(7)

多倍長除算の準備 (6-2)

void LVsubOffset( long a[], long b[], long No ) { …

for( ; v ; k++ ){

if( a[k+No] ){

a[k+No]--; v = 0;

}else{

a[k+No] = D-1; v = 1;

} } }

上位桁への桁借りを処理

多倍長数 a と b を比較

前のスライドから続き

被除数と除数の桁あわせ (1)

多倍長数が表す数

0 N-2

N-1

A ’ =

N N+1

(D=10000)

∑ ⁺

=

= ¹

0 D ] [ a ' ^N

k

k k

A 配列添え字

やはり整数多倍長除算結果(商) ∑ ⁺

=

= ¹

0 D ] [ ' q

' ' ^N

k

k k

B Q A

多倍長除算は、多倍長数aが整数値を表すと想定多倍長数a

1 2345 6789 … 9012 要素の値

整数値

6789 2345 1

0 9012

小数点位置

被除数と除数の桁あわせ (2)

多倍長数による小数の表現 _(D=10000)

N N k

N

k A D

k

A ⁻

+

=

− =

= ∑ ¹ ^a ^[ ^] ^D ^'

0 (多倍長数aが表す整数値 A’ と実際値 A の関係) 小数点位置を移動

0 N-2

N-1

A =

N N+1 配列添え字

1.2345 6789 … 9012 要素の値

実際値

6789 2345 1

0 9012

被除数と除数の桁あわせ (3)

多倍長数による小数の除算被除数 A を除数 B で除算 ⇒

' ' D '

D '

B A B

A B A

N N =

=

= ₋ ⁻

商

必要な商は ^N _N

N

k

N k

BD D A

B k A

Q = ⁺ = ⁻ =

=

∑ ¹ − 0

D ] [ q つまり

小数どうしの除算は実際値 A ( B )を整数値 A ’( B ’)とみなして多倍長除算を行うただし、

商は整数値形式で得られる ∑ ⁺

=

= ¹

0 D ] [ ' q ' ^N

k

k k

B A B A

除数はD ^N 倍して除算を行えばよい

除算の実行

long q[N2], r[N2];

long a[N2] = {0}, b[N2] = {0};

a[N] = 16;

b[0] = 5;

LVdiv( a, b, q, r ); // q = 16/5

16/5を実行

16×D ^N-N =16 5×D ^0-N = 5×D ^-N

0 N-2

N-1 N N+1 添え字

a 0 16 0 0 0

b 0 0 0 0 5

q 0 3 2 0 0

A =16 Q =3.2 B ’=5

πの計算 (1)

void CalcPi( void )

{ long s[N+2], s5[N+2], s239[N+2];

long t[N+2] = {0}, q[N2], r[N2];

long d2[N+2] = {0}, d[N+2] = {0};

int pm;

/* 16arctan(1/5) を計算 /

t[N] = 16;

d2[0] = 5;

LVdiv( t, d2, s5, r ); // s5 = 16/5 LVcopy( s5, t ); // t = 16/5

…

⎟ ⎠

⎜ ⎞

⎝

⎛ 5 arctan 1

16 を求める

(8)

πの計算 (2)

d2[0] = 25;

d[0] = 3;

pm = 0; // 減算 while(1){

LVdiv( t, d2, t, r ); // t = t/25 if( GetMSD(t) < 0 ) break;

LVdiv( t, d, s, r ); // s = t/d if( pm ){ // 加算

LVadd( s5, s ); pm = 0;

}else{ // 減算

LVsub( s5, s ); pm = 1;

} LVaddS( d, 2 );

}

tが0になれば終了ここでは t=16/5

s5 ← s5 + s s5 ← s5 - s d ← d+2

πの計算 (3)

… /* 4arctan(1/239) を計算 /

t[N] = 4;

d2[0] = 239;

LVdiv( t, d2, s239, r ); // s239 = 4/239 LVcopy( s239, t ); // t = 4/239

…

⎟ ⎠

⎜ ⎞

⎝

⎛ 239 arctan 1

4 _を求める

πの計算 (4)

d2[1] = 5; d2[0] = 7121; // 239*239 d[2] = 0; d[1] = 0; d[0] = 3;

pm = 0; // 減算 while(1){

LVdiv( t, d2, t, r ); // t = t/(239*239) if( GetMSD(t) < 0 ) break;

LVdiv( t, d, s, r ); // s = t/d if( pm ){ // 加算

LVadd( s239, s ); pm = 0;

}else{ // 減算

LVsub( s239, s ); pm = 1;

} LVaddS( d, 2 );

}

tが0になれば終了

ここでは t=4/239

s239 += s s239 -= s d ← d+2

πの計算 (5)

void CalcPi( void )

{ /* 16arctan(1/5) を計算 /

… /* 4arctan(1/239) を計算 /

…

LVsub( s5, s239 );

} s5-s239 = π

実行結果

3. 1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164 0628620899 8628034825 3421170679

…

2162484391 4035998953 5394590944 0704691209 1409387001 2645600162 3742880210 9276457931 0657922955 2498872758 4610126483 6999892256 9596881592 0560010165 5256375662

正しくは78 N=2500 ⇒ 小数点以下10000桁まで計算

正しく求めるには、N=2501 とする

3. 1415926535 8979323846 2643383279 5028841971 6939937510

…

4610126483 6999892256 9596881592 0560010165 5256375678 5646

プログラムの実行時間

桁数Nとプログラム実行時間の関係

0.01 0.10 1.00 10.00 100.00 1000.00

100 1000 10000

N

秒

O(n ³ )

PentiumIV

2.4GHz

Celeron M

1.5GHz

(9)

まとめ

πの計算

整数型の配列を利用して多倍長数を演算

改善の余地あり

π(円周率)を小数以下10,000桁まで計算 しなさい。

プログラミング言語 I 第 6 回 πの計算 1

課題