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

(1)

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

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

(2)

課題

π(円周率)を小数以下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

(3)

πの計算

arctan

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

+L

− +

−

− =

= ^∞ − ⁻

=

∑

⁽₂¹⁾ ₁+ ₃ ₅ ₇

arctan

7 5

1 3 2 1

1 x x x

x i x

x ⁱ

i

= S

⋅ +

⋅ −

⋅ +

−

⎟ =

⎠

⎜ ⎞

⎝

⎛ ₃ ₅ ₇ L

5 7

16 5

5 16 5

3 16 5

16 5

arctan 1 16

⎟ の場合

⎠

⎜ ⎞

⎝

⎛ 5 arctan 1 16

(4)

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

5 16

1 =

T S = T₁

25

3 1

T = T

3 T3

S S = −

25

5 3

T = T

5 T5

S S = +

25

7 5

T = T

7 T7

S

S = − T_n

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

= S

⋅ +

⋅ −

⋅ +

−

⎟ =

⎠

⎜ ⎞

⎝

⎛ ₃ ₅ ₇ L

5 7

16 5

5 16 5

3 16 5

16 5

arctan 1 16

(5)

高精度演算の考え方

小数以下

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桁分しか値を記憶できない

別の方法が必要

(6)

桁数の長い数 = 多倍長数

多数のデータ

(

各桁の値

)

をまとめて処理

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

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

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]

多倍長数という

(7)

長い桁数の演算

C

言語では

int

型データが基本演算単位

1

桁の演算と

2

桁の演算はどちらも

int

型

⇒ 計算時間は同じ

(

同じ

1

回の演算処理

)

1

回の

int

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

1

つの要素が

M

桁を表す場合、

乗算の過程で桁数が最大で

2M

桁になる

int

型で

2M

桁の整数が扱えるような桁数

M

を

選択

(8)

長い桁数の演算

int

型ではなく

long

型を指定その理由は

…

現在の

32

ビットプロセッサでは

int

型と

long

型と同じ

⇒

long

型演算を最短時間

(1

クロック

)

で実行

もし

16

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

long

型

=32

ビットとして正しく実行できる

(9)

要素の桁数の決定

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)

(10)

多倍長数と配列

小数点以下

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要素の場合

(11)

多倍長加算

下位桁から順に処理する

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

(12)

多倍長減算

下位桁から順に処理する

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

(13)

多倍長精度除算

被除数(割られる数): a = a_sD^s + a_s₋₁D^s⁻¹ +L+ a₀ 除数(割る数): b = b_tD^t + b_t₋₁D^t⁻¹ +L+ b₀

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

商: q = q_uD^u + q_u₋₁D^u⁻¹ +L+ q₀ 余り: r = r_tD^t + r_t₋₁D^t⁻¹ +L+ r₀

r qb

a = +

(

⁰ ^≤ ^r ^< ^b

)

0 ,

0 ≠

≠ _t

s b

a

(14)

多倍長精度除算

0 1

1D 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 ≤ ^{が成り立つと仮定}

(15)

多倍長精度除算 (3)

商の最上位桁位置

u

および桁値

q_u

を求める

の場合

のとき

t

s b

a ≥

t

bDs

a ≥ ⁻

D a

D D b

s

t < <

≤ ,

2 ^{であるから、} < 2 ⇒ < 2

t s t

s b

b a a

さらに ≥ 1

t s

b

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

(16)

多倍長精度除算 (4)

の場合(続き)

のとき

t

s b

a ≥

t

bDs

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 b a b a b

a = , ₋₁ = ₋₁,K, ₋ < ₋

(

^j ^≥¹

)

t t

t D D

D b

b ≥ ≥ 2 ^より 2b ≥ D^t⁺¹ ⇒ 2bD^s⁻^t⁻¹ ≥ D^s すなわち^a ^> ^bD^s⁻^t ⁻ ^D^s ^≥ ^bD^s⁻^t ⁻ ²^bD^s⁻^t⁻¹ ⁼ ^bD^s⁻^t⁻¹

(

^D ⁻ ²

)

(

^j ⁼¹

)

(17)

多倍長精度除算 (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より大きい、

すなわちR₀が負になっている場合も存在する

(

^bt ^D^t ≤ ^b

)

商の候補を1増やしながら、余りが非負になる真の商を見つける

(18)

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

a_s < b_t

の場合

(

続き

)

0

0, r R

Q

q ← ← ^とし、

を繰り返す b

r r qb

a = + , <

b Q a

b R a D

Q a

t s s

0 0

1

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 b

a D

q a

Q _s

t t

s s

※1

(19)

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

a_s < b_t

の場合

(

続き

) s = t +1 ^の場合

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

1

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_tD^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 ₋₁ ₌ ₀_, _q₀ ₌ _Q₀ または _Q₀ ₋₁または _Q₀ ₋ ₂

(20)

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

a_s < b_t

の場合

(

一般の場合

)

0 1

0D , 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 1 0

0 ⁻ ⎥, = − ⁻ ⁻

⎦

⎢ ⎥

⎣

⎢ +

= とおく

≥ 0

r′ ^{となるまで} b

D r

r D

q

q′ ← ′ − ^s⁻^t⁻¹, ′ ← ′ + ^s⁻^t⁻¹

だから

1 1

1

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

(21)

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

a_s < b_t

の場合

(

一般の場合

) (

続き

)

1

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 ₁

1

1 ⎟⎟ +

⎠

⎜⎜ ⎞

⎝

= ⎛ − +

−

= ₋ ₋ ₋ ₋ ₋ ₋ _t

t

t t 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 またはまたは

(22)

多倍長精度除算 (6)

商の最上位桁位置

u

および桁値

q_u

を求める

の場合

のとき

の場合

t

s b

a ≥

t

bDs

a ≥ ⁻

1

, =

−

= s t q_u u

t

bDs

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 ¹

ただし

(23)

多倍長精度除算 (7)

商の最上位桁位置

u

および桁値

q_u

が得られたら

b D q a

r = − _u ^u

a ← r として再度除算を実行

中間余りを計算

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

(24)

多倍長除算の例

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

(25)

多倍長除算の例 ( 続き 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

(26)

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

4. 574÷64 ⇒ より

= 574 – 512=62

t

s b

a < ₉

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

(27)

多倍長除算の準備 (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

(28)

整数部を考慮

π

は整数部

(=3)

あり

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

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

3. 1415 9265 3589 7932 3846 … 1989

246

247 245 0

248 249

拡張した部分

π =

N=250要素の場合

250 251

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

#define N2 (N+2)

(29)

多倍長除算の準備 (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ならば桁上げ下位からの桁上げを加算

(30)

多倍長除算の準備 (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以上にするために必要な最小係数

(31)

多倍長除算 (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に初期化

(32)

多倍長除算 (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 を検査

bDs

a ≥ ⁻

(33)

多倍長除算 (3)

while( 1 ){

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

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

…

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

q[u]--;

}LVsubOffset( a, m, u );

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

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

u, q[u]の候補を求める(省略)

(34)

多倍長除算の準備 (4)

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

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

(35)

多倍長除算の準備 (5)

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

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]が一致

(36)

多倍長除算の準備 (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を引く

次のスライドへ続く

(37)

多倍長除算の準備 (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

を比較

前のスライドから続き

(38)

被除数と除数の桁あわせ (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

小数点位置

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

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

課題

πの計算

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

∑

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

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

高精度演算の考え方

小数以下

桁

の精度

細かな数値を表すための

言語データ型

桁数の長い数 = 多倍長数

多数のデータ

各桁の値

をまとめて処理

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

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

多倍長数という

長い桁数の演算

言語では

型データが基本演算単位

桁の演算と

桁の演算はどちらも

型

⇒ 計算時間は同じ

同じ

回の演算処理

回の

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

つの要素が

桁を表す場合、

乗算の過程で桁数が最大で

桁になる

型で

桁の整数が扱えるような桁数

を

選択

長い桁数の演算

型ではなく

型を指定 その理由は

現在の

ビットプロセッサでは

型と

型 と同じ

⇒

型演算を最短時間

クロック

で実行

もし

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

型

ビットとして正しく実行できる

要素の桁数の決定

型で扱える値の範囲

型の最大桁数と

の関係

として

を選択

⇒

つの要素が

から

桁

を表す 符号付

～

符号なし

～

の十進数桁数

≧

～

多倍長数と配列

小数点以下

桁を表すには

要素の配列が必要

桁と配列要素の対応 少なくとも

通りあり

多倍長加算

下位桁から順に処理する

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

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

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

型を指定その理由は

型と同じ

を表す符号付

桁と配列要素の対応少なくとも

が得られたら

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

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