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

(1)

プログラミング言語 I 第 7 回 πの計算 2

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

課題

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

マーチンの公式を利用

計算の高速化

無駄な計算を省く

プログラムを見直し

⎟ ⎠

⎜ ⎞

⎝

− ⎛

⎟ ⎠

⎜ ⎞

⎝

= ⎛

239 arctan 1 5 4

arctan 1 π 16

πの計算

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

πの計算プログラム (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

を求める

πの計算プログラム(2)

d2[0] = 25;

d[0] = 3;

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

LVdiv( t, d2, t, r ); // t = t/25 if( GetMSD(t) < 1 ) 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

多倍長除算

long LVdiv( long a0[], long b0[], long q[] ) { …

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;

s = GetMSD( a );

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

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

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

q[u]--;

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

} }

…

これらすべて多倍長数の最上位桁位置を必要とする

(2)

最上位桁を求める処理

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

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

if( a[k] ) break;

} return k;

}

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

値が小さい ( 最上位桁が添え字 0 に近い ) と時間がかかる

同じ値に対して何度も繰り返し実行している最上位桁を求める処理を節約する

最上位桁を得る

最上位桁を毎回求めるのではなく、変数に記憶しておく

1 つの多倍長数に対して

多倍長数を記憶する配列v[N2]

最上位桁位置を記憶する変数nMSD

意味の異なるデータをまとめて扱いたい構造体を利用

構造体の宣言

上に示したのは名前付き構造体

他に匿名(anonymous)の構造体もあり

⇒ ここでは省略

struct

名前

{

メンバ宣言

};

書式

1. 1 つの構造体としてまとめるデータ項目 ( メンバ ) の名前とデータ型を宣言 2. 構造体の型の名前を宣言

構造体型宣言の例

多倍長数の場合の例 struct LongValue_data {

long v[N2];

int nMSD;

};

例

名前LongValue_data の構造体の型を宣言

2つのメンバの

　データ型と名前を宣言

通常の変数宣言と同じ書式

構造体名の規則は変数や関数と同じ

(アルファベット、数字、‘_’の組合せ)

変数/関数/配列と名前が重なっても問題ない

構造体型の変数宣言

通常のデータ型と同様に変数や配列を宣言可能

struct LongValue_data x, y;

例

pは構造体personal_dataを指すポインタ struct LongValue_data A[100];

struct LongValue_data

*p;

2つの構造体変数xとyを宣言

要素数100の構造体配列Aを宣言例

例

(後述)

構造体の型宣言と変数宣言

同時に宣言可能

構造体の型は1度だけ宣言すればよい struct st {

int n;

double f;

} x, y;

例名前stの構造体の型を宣言同時に構造体stの

変数xとyを宣言

struct st { int n;

double f;

} x, y;

struct st z;

例

宣言済みの構造体stの変数zを宣言

名前stの構造体の型と構造体stの変数xとyを宣言

(3)

構造体の代入

同じ型の構造体は代入が可能

たとえメンバが全く同じでも、名前が異なる構造体は代入できない

struct LongValue_data x, y;

x = y;

例

struct st1 { int a,b; } x;

struct st2 { int a,b; } y;

x = y;

異なる名前の構造体の代入

⇒コンパイルエラー例

メンバの参照

メンバのデータを読み書きするには struct st1 { int a,b; } x, y[10];

x.a = 10;

y[1].b += x.a;

例

メンバは通常の変数とまったく同じように使用できる構造体変数名とメンバ名をピリオドでつなぐ

構造体の初期化

配列の初期化と同じ文法で行う struct st1 {

int a,b;

char s[10];

};

struct st1 x = {1, 10, “Hello” };

例

メンバaの初期化データ

メンバbの初期化データ

メンバsの初期化データ初期化

メンバより初期化データが少ない場合は、

初期化データとして0やNULLを補う

型別名の宣言

typedef 宣言子を利用 typedef

型別名

;

書式型に別名を与える

typedef struct LongValue_data { long v[N2];

int nMSD;

} LongValue;

LongValue a, x[5];

例

構造体の型を宣言するとともに、別名

LongValueを与える

別名を使って

変数や配列を宣言変数a, 配列xはLongValue型として宣言される

(= struct LongValue_data型)

多倍長加算

加数bの最上位桁まで加算

void LVadd( LongValue pa, LongValue pb ) { long i, c=0;

for( i=0 ; i<=pb->nMSD ; i++ ){

pa->v[i] += (pb->v[i]+c);

if( pa->v[i] >= D ){

pa->v[i] -= D;

c = 1;

}else{

c = 0;

} }

…

#define D 10000

として、

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

a ← a+b

変数の前の*や->は何?

多倍長数の構造体を引数に

関数の呼び出しについて

定義した関数を呼び出して、処理を行う

#include <stdio.h>

#include <math.h>

double

solve2p( double a, double b, double c ) {

double

bunshi, bunbo;

bunshi = -b+sqrt(bb-4a*c);

bunbo = 2*a;

return

bunshi / bunbo;

}

int

main()

{

double

xp = solve2p( 1, 4, 3 );

printf( “xp=%f¥n”, xp );

}

関数の定義

関数の呼び出し

(4)

関数の呼び出し ( その 2)

実行例

プログラム中で関数を呼び出す部分

引数に実際の値を指定して関数を呼び出し

⇒実引数という

引数aに1、bに4、cに3が対応付けられる

実引数に対して、関数定義側の引数を仮引数という

xp=-1.000000

xp = solve2p( 1, 4, 3 );

2次方程式 x

²

+ 4 x + 3 = 0

の解-1, -3

doublesolve2p( doublea, doubleb, doublec )

xp = solve2p( 1, 4, 3 );

^並び順に_よる対応

関数呼び出しと引数

仮引数は変数

実引数の値が初期値として代入される int n=1, m=5, y;

y = func( n, m );

... int func( int a, int b ) { int x;

a = a*2;

... return x;

}

実引数(1, 5) 仮引数(a, b) 初期値として

aに1, bに5が与えら

れて関数が開始する

仮引数の値を変更しても実引数は変化しない以降a, bは変数として使用

C 言語の引数は ‘ 値渡し ’

値渡し : 実引数の値が仮引数にコピーされる

問題点が 2 つ

関数側で仮引数を書き換えても、実引数は変化しない

引数のデータ量が多い場合に、コピーのための時間とメモリの両方が無駄遣い

int n=1, m=5, y;

y = func( n, m );

... void func( int a, int b ) { a = a*2;

} ...

仮引数aを変更しても実引数nは変化しない

ポインタ

計算機のメモリ上に記憶された変数の位置 ( アドレス ) をポインタという

変数を名前(変数名)で指し示すだけでなく、

ポインタで指し示すこともできる

ポインタ変数

ここにデータを記憶

変数xの

‘ポインタ’

変数x

変数名

‘x’

変数a 変数w

ポインタに関する操作

変数のポインタを取得 ⇒ 単項演算子&

ポインタを記憶する変数: ポインタ変数

ポインタが指す変数を読み書き

⇒単項演算子 * int i=10, j, *p;

p = &i;

j = *p;

*p = 100;

例

int i;

例

&i

変数iのポインタ(アドレス)

int *p;

例

int型変数のポインタを記憶する

ポインタ変数pを宣言

pに変数iのポインタを代入 pが指す値(i=10)を読み出し pが指す変数(i)に代入

ポインタを引数に用いる

ポインタも関数引数に使用可能

実引数として与えたポインタが関数の仮引数に代入される

int m=5, y;

/ここではm=5/

func( &m );

/ここではm=10/

...

void func( int

*p )

{

*p = 10;

}

実引き数として変数mのポインタ

仮引き数pに変数mのポインタが代入される

*pは関数の外の変数mを参照

(5)

ポインタを引数に用いると

ポインタを通して、関数外の変数を変更可能

大域変数だけでなく局所変数も変更可能

ポインタは位置情報 ( アドレス ) であり、実引数から仮引数へのコピーは簡単

元のデータがどんなに巨大でもポインタは高々

32ビットないし64ビット

(pが指す(->)構造体のメンバa)

構造体とポインタ

構造体もポインタを利用可能

ポインタを通してメンバを読み書き struct st1 { int a,b; } x, y, *p;

p = &x;

y = *p;

例

struct st1 { int a,b; } x, *p;

p = &x;

(*p).a = 10;

例

ポインタの定義どおりの方法

p->a = 10;

等価な

記述演算子‘->’を利用構造体変数xのポインタをポインタ変数pに代入

多倍長加算

加数 b の最上位桁まで加算

void LVadd( LongValue pa, LongValue pb ) { long i, c=0;

for( i=0 ; i<=pb->nMSD ; i++ ){

pa->v[i] += (pb->v[i]+c);

if( pa->v[i] >= D ){

pa->v[i] -= D;

c = 1;

}else{

c = 0;

} }

…

ポインタと‘->’を用いて構造体のメンバを読み書き

a ← a+b

多倍長数構造体のポインタを引数

typedef struct LongValue_data { long v[N2]; /*

多倍長数データ

*/

int nMSD; /*

最上位桁位置

*/

} LongValue;

多倍長加算 ( 続き )

上位桁への桁上げを考慮

void LVadd( LongValue pa, LongValue pb ) {

… if( c ){

while( 1 ){

pa->v[i]++;

if( pa->v[i] >= D ) pa->v[i++] = 0;

else break;

} }

if( pa->nMSD < i ) pa->nMSD = i;

}

bの最上位桁まで加算(前ページ)

残りの上位桁への桁上げ処理

a ← a+b

上位桁への桁上げがなくなれば処理終了

加算後にaの最上位桁位置修正

多倍長減算

加数bの最上位桁まで減算

void LVsub( LongValue pa, LongValue pb ) { long i, c=0;

for( i=0 ; i<=pb->nMSD ; i++ ){

pa->v[i] -= (pb->v[i]+c);

if( pa->v[i] < 0 ){

pa->v[i] += D;

c = 1;

}else{

c = 0;

} }

…

引算結果が負ならば上位から桁借りして

Dを加えて正にする

a ← a-b

多倍長減算(続き1)

上位桁への桁借りを考慮

void LVsub( LongValue pa, LongValue pb ) {

… if( c ){

while( 1 ){

pa->v[i]--;

if( pa->v[i] < 0 ) pa->v[i++] = D-1;

else break;

} } } …

bの最上位桁まで減算(前ページ)

残りの上位桁への桁借り処理

a ← a-b

上位桁への桁借りがなくなれば処理終了

aの最上位桁位置修正(次ページ)

(6)

多倍長減算 ( 続き 2)

減算の結果、最上位桁が 0 になる場合あり void LVsub( LongValue pa, LongValue pb ) { …

if( pa->nMSD < i ) pa->nMSD = i;

while( pa->v[pa->nMSD] == 0 ){

pa->nMSD--;

} }

桁借り処理完了(前ページ)

a ← a-b

最上位桁が0ならば最上位桁位置を下げる最も上位にある非ゼロの桁が最上位桁

多倍長除算 (1)

void LVdiv(LongValue pa0, LongValue pb0, LongValue pq, LongValue pr)

{ long s, t, k, u, nk;

LongValue a, b;

nk = Normalize( pb0, &b );

LVmulS( pa0, nk, &a );

t = b.nMSD;

for( k=0 ; k<N2 ; k++ ) pq->v[k] = 0;

pq->nMSD = 0;

while(1){

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

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

…

q ← a0/b0

bを正規化(nk倍) aもnk倍

商を0に初期化

多倍長除算 (2)

s = a.nMSD;

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

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

u = s-t; pq->v[u] = 1;

}else{

u = s-t-1; pq->v[u] = D-1;

}else{ }

u = s-t-1; pq->v[u] = (a.v[s]*D+a.v[s-1])/b.v[t];

} LVmulS( &b, pq->v[u], &m );

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

pq->v[u]--;

LVmulS( &b, pq->v[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, pq->v[u], &m );

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

pq->v[u]--;

LVmulS( &b, pq->v[u], &m );

} LVsubOffset( &a, &m, u );

}

pq->nMSD = GetMSD( pq );

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

^u

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

商qの最上位桁位置を求める

πの計算 (1)

void CalcPi( void )

{ LongValue s, t, s5, s239, q, r, d, d2;

int pm, k;

/* 定数初期化 */

for( k=0 ; k<N2 ; k++ ) t.v[k] = 0;

t.nMSD = 0;

LVcopy( &t, &d );

LVcopy( &t, &d2 );

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

t.v[N] = 16; t.nMSD = N;

d2.v[0] = 5; d2.nMSD = 0;

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

LVcopy( &s5, &t ); // t = 16/5;

t=d=d2=0に初期化

πの計算 (2)

d2.v[0] = 25; d2.nMSD = 0;

d.v[0] = 3; d.nMSD = 0;

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

LVdiv( &t, &d2, &t, &r ); // t = t/25 if( t.nMSD < 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

(7)

πの計算 (3)

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

t.v[N] = 4; t.nMSD = N;

d2.v[0] = 239; d2.nMSD = 0;

LVdiv( &t, &d2, &s239, &r );

LVcopy( &s239, &t );

d2.v[1] = 5; d2.v[0] = 7121;

d2.nMSD = 1; // 239*239

d.v[2] = 0; d.v[1] = 0; d.v[0] = 3;

d.nMSD = 0;

…

⎟ ⎠

⎜ ⎞

⎝

⎛ 239 arctan 1

4

_を求める

πの計算 (4)

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

LVdiv( &t, &d2, &t, &r ); // t = t/(239*239) if( t.nMSD < 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 5256375633

正しくは78

N=2500 ⇒ 小数点以下10000桁まで計算

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

3.

1415926535 8979323846 2643383279 5028841971 6939937510

…

4610126483 6999892256 9596881592 0560010165 5256375678 5628

プログラムの実行時間

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

0.01 0.10 1.00 10.00 100.00 1000.00

100 1000 10000

N

秒

PentiumIV 2.4GHz

オリジナル

最上位桁位置記録

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

プログラミング言語 I 第 7 回 πの計算 2

課題

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

マーチンの公式を利用

計算の高速化

⎟ ⎠

⎜ ⎞

⎝

− ⎛

⎟ ⎠

⎜ ⎞

⎝

= ⎛

239 arctan 1 5 4

arctan 1 π 16

πの計算

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

− +

−

− =

=

−

∑ ( 2 1 ) 1

3 5 7

arctan

x x x

x i x

x

= S

⋅ +

⋅ −

⋅ +

−

⎟ =

⎠

⎜ ⎞

⎝

⎛

7 5

L

16 5 5

16 5 3

16 5 16 5 arctan 1 16

⎟

⎠

⎜ ⎞

⎝

⎛ 5 arctan 1 16

πの計算プログラム (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;

/* 16*arctan(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

πの計算プログラム(2)

d2[0] = 25;

d[0] = 3;

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

LVdiv( t, d2, t, r ); // t = t/25 if( GetMSD(t) < 1 ) 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

多倍長除算

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

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

∑ ⁽ ₂ ¹ ⁾ ₁

₃ ₅ ₇

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

値が小さい ( 最上位桁が添え字 0 に近い ) と時間がかかる

同じ値に対して何度も繰り返し実行している最上位桁を求める処理を節約する

最上位桁を毎回求めるのではなく、変数に記憶しておく

意味の異なるデータをまとめて扱いたい構造体を利用

通常のデータ型と同様に変数や配列を宣言可能

たとえメンバが全く同じでも、名前が異なる構造体は代入できない