以前、大学でプログラミング実習を教えたときの課題だった、複数のアルゴリズムで円周率を求めるプログラムが出てきたのでここに書いておく。実行時間も測定している。これはCで書いたけど、使用言語はFortranでも可だったはず。
出力結果:
Machin's formula: 3.14159265358979 Time: 0.437000s (1000000 times) Gauss-Legendre : 3.14159265358979 Time: 0.219000s (1000000 times) Count points : 3.14552000000000 Time: 0.265000s (10 times) Monte Carlo : 3.14394000000000 Time: 0.594000s (10 times) Circumference : 3.14159265358979 Time: 3.047000s (1000000 times) Integral : 3.14159265393433 Time: 0.375000s (10 times)
Machin's formula, Gauss-Legendre, Circumferenceは100万回計算、Count points, Monte Carlo, Integralは10回計算して時間を測定している。それにしても、算術幾何平均を使ったガウス-ルジャンドルのアルゴリズムは速いなぁ。
以下、ソースコード。
/*======================================================= ratio of the circumference of a circle to its diameter =======================================================*/ #include <stdio.h> #include <stdlib.h> #include <math.h> #include <time.h> /*------------------------------------------------------- Machin's formula: PI/4 = 4*arctan(1/5) - arctan(1/239) --------------------------------------------------------*/ double pi_machin(int dummy) { int k; double p, temp, last; p = 0.0; /* 4*arctan(1/5) */ k = 1; temp = 16.0 / 5.0; /* arctan(x) = x - x^3/3 + x^5/5 - x^7/7 + ... */ do { last = p; p += temp / k; temp /= -5.0 * 5.0; k += 2; } while (p != last); /* arctan(1/239) */ k = 1; temp = 4.0 / 239.0; do { last = p; p -= temp / k; temp /= -239.0 * 239.0; k += 2; } while (p != last); return p; } /*------------------------------------------------------- Gauss-Legendre algorithm a[1] = 1, b[1]= 1/sqrt(2) arithmetic means: a[n+1] = (a[n] + b[n]) / 2 geometric means: b[n+1] = sqrt(a[n] * b[n]) PI[n] = 2*a[n+1]^2 / (1-SUM^n_k=0{2^k*(a[k]^2-b[k]^2)}) = 2*a[n+1]^2 / (1-SUM^n_k=0{2^k*(a[k]-a[k-1])^2}) PI[n->inf] -> PI --------------------------------------------------------*/ double pi_gauss_legendre(int dummy) { int i; double a, b, s, temp, last; a = 1.0; b = 1.0 / sqrt(2.0); s = 1.0; temp = 4; do { last = a; a = (a + b) / 2.0; /* arithmetic means */ b = sqrt(last * b); /* geometric means */ s -= temp * (a - last) * (a - last); temp *= 2.0; } while (a != last); return (a + b) * (a + b) / s; } /*------------------------------------------------------- Using count points in circle num: length of side -------------------------------------------------------*/ double pi_points(int num) { int i, j, hit; double x, y, p; hit = 0; for (i = 0; i < num; i++) { for (j = 0; j < num; j++) { x = (double)i / num; y = (double)j / num; if (x * x + y * y < 1.0) hit++; } } p = (double)hit / (num * num); return p * 4.0; } /*------------------------------------------------------- Using Monte Carlo method num: number of random points -------------------------------------------------------*/ double pi_monte(int num) { int i, hit; double x, y, p; hit = 0; for (i = 0; i < num; i++) { x = (double)rand() / RAND_MAX; y = (double)rand() / RAND_MAX; if (x * x + y * y < 1.0) hit++; } p = (double)hit / num; return p * 4.0; } /*------------------------------------------------------- 2*PI <- length of side * N div: number of division on a side of hexagon -------------------------------------------------------*/ double pi_circumference(int div) { /* 2*PI = sin(theta) * 360.0 / theta */ int i; double a, b, multi; a = 1.0; multi = 1.0; for (i = 0; i < div; i++) { b = 1.0 - sqrt(1.0 - a * a / 4.0); a = sqrt(b * b + a * a / 4.0); multi *= 2.0; } return a * 6.0 * multi / 2.0; } /*------------------------------------------------------- x^2 + y^2 = 1, y = sqrt(1 - x^2) PI/4 = INT{y}dx div: number of division on integrals -------------------------------------------------------*/ double pi_integral(int div) { int i; double p, delta, y, a; y = 0.0; for (i = 0; i < div; i++) { a = (i + 0.5) / div; y += sqrt(1.0 - a * a); } return y / div * 4.0; } /*------------------------------------------------------- Measurement of time func: function of calculation of pi arg: argument of function num: number of repeats msg: message of output -------------------------------------------------------*/ void run(double (*func)(int), int arg, int num, char *msg) { clock_t start, finish; int i; start = clock(); for (i = 0; i < num - 1; i++) (*func)(arg); printf("%-16s: %16.14f\n", msg, (*func)(arg)); finish = clock(); printf("Time: %14.6fs (%d times)\n",(double)(finish - start) / CLOCKS_PER_SEC, num); printf("\n"); } int main() { srand((unsigned int)time(NULL)); run(pi_machin, 0, 1000000, "Machin's formula"); run(pi_gauss_legendre, 0, 1000000, "Gauss-Legendre"); run(pi_points, 1000, 10, "Count points"); run(pi_monte, 1000000, 10, "Monte Carlo"); run(pi_circumference, 50, 1000000, "Circumference"); run(pi_integral, 1000000, 10, "Integral"); return 0; }
出力結果:
Machin's formula: 3.14159265358979 Time: 0.437000s (1000000 times) Gauss-Legendre : 3.14159265358979 Time: 0.219000s (1000000 times) Count points : 3.14552000000000 Time: 0.265000s (10 times) Monte Carlo : 3.14394000000000 Time: 0.594000s (10 times) Circumference : 3.14159265358979 Time: 3.047000s (1000000 times) Integral : 3.14159265393433 Time: 0.375000s (10 times)
Machin's formula, Gauss-Legendre, Circumferenceは100万回計算、Count points, Monte Carlo, Integralは10回計算して時間を測定している。それにしても、算術幾何平均を使ったガウス-ルジャンドルのアルゴリズムは速いなぁ。
以下、ソースコード。
/*======================================================= ratio of the circumference of a circle to its diameter =======================================================*/ #include <stdio.h> #include <stdlib.h> #include <math.h> #include <time.h> /*------------------------------------------------------- Machin's formula: PI/4 = 4*arctan(1/5) - arctan(1/239) --------------------------------------------------------*/ double pi_machin(int dummy) { int k; double p, temp, last; p = 0.0; /* 4*arctan(1/5) */ k = 1; temp = 16.0 / 5.0; /* arctan(x) = x - x^3/3 + x^5/5 - x^7/7 + ... */ do { last = p; p += temp / k; temp /= -5.0 * 5.0; k += 2; } while (p != last); /* arctan(1/239) */ k = 1; temp = 4.0 / 239.0; do { last = p; p -= temp / k; temp /= -239.0 * 239.0; k += 2; } while (p != last); return p; } /*------------------------------------------------------- Gauss-Legendre algorithm a[1] = 1, b[1]= 1/sqrt(2) arithmetic means: a[n+1] = (a[n] + b[n]) / 2 geometric means: b[n+1] = sqrt(a[n] * b[n]) PI[n] = 2*a[n+1]^2 / (1-SUM^n_k=0{2^k*(a[k]^2-b[k]^2)}) = 2*a[n+1]^2 / (1-SUM^n_k=0{2^k*(a[k]-a[k-1])^2}) PI[n->inf] -> PI --------------------------------------------------------*/ double pi_gauss_legendre(int dummy) { int i; double a, b, s, temp, last; a = 1.0; b = 1.0 / sqrt(2.0); s = 1.0; temp = 4; do { last = a; a = (a + b) / 2.0; /* arithmetic means */ b = sqrt(last * b); /* geometric means */ s -= temp * (a - last) * (a - last); temp *= 2.0; } while (a != last); return (a + b) * (a + b) / s; } /*------------------------------------------------------- Using count points in circle num: length of side -------------------------------------------------------*/ double pi_points(int num) { int i, j, hit; double x, y, p; hit = 0; for (i = 0; i < num; i++) { for (j = 0; j < num; j++) { x = (double)i / num; y = (double)j / num; if (x * x + y * y < 1.0) hit++; } } p = (double)hit / (num * num); return p * 4.0; } /*------------------------------------------------------- Using Monte Carlo method num: number of random points -------------------------------------------------------*/ double pi_monte(int num) { int i, hit; double x, y, p; hit = 0; for (i = 0; i < num; i++) { x = (double)rand() / RAND_MAX; y = (double)rand() / RAND_MAX; if (x * x + y * y < 1.0) hit++; } p = (double)hit / num; return p * 4.0; } /*------------------------------------------------------- 2*PI <- length of side * N div: number of division on a side of hexagon -------------------------------------------------------*/ double pi_circumference(int div) { /* 2*PI = sin(theta) * 360.0 / theta */ int i; double a, b, multi; a = 1.0; multi = 1.0; for (i = 0; i < div; i++) { b = 1.0 - sqrt(1.0 - a * a / 4.0); a = sqrt(b * b + a * a / 4.0); multi *= 2.0; } return a * 6.0 * multi / 2.0; } /*------------------------------------------------------- x^2 + y^2 = 1, y = sqrt(1 - x^2) PI/4 = INT{y}dx div: number of division on integrals -------------------------------------------------------*/ double pi_integral(int div) { int i; double p, delta, y, a; y = 0.0; for (i = 0; i < div; i++) { a = (i + 0.5) / div; y += sqrt(1.0 - a * a); } return y / div * 4.0; } /*------------------------------------------------------- Measurement of time func: function of calculation of pi arg: argument of function num: number of repeats msg: message of output -------------------------------------------------------*/ void run(double (*func)(int), int arg, int num, char *msg) { clock_t start, finish; int i; start = clock(); for (i = 0; i < num - 1; i++) (*func)(arg); printf("%-16s: %16.14f\n", msg, (*func)(arg)); finish = clock(); printf("Time: %14.6fs (%d times)\n",(double)(finish - start) / CLOCKS_PER_SEC, num); printf("\n"); } int main() { srand((unsigned int)time(NULL)); run(pi_machin, 0, 1000000, "Machin's formula"); run(pi_gauss_legendre, 0, 1000000, "Gauss-Legendre"); run(pi_points, 1000, 10, "Count points"); run(pi_monte, 1000000, 10, "Monte Carlo"); run(pi_circumference, 50, 1000000, "Circumference"); run(pi_integral, 1000000, 10, "Integral"); return 0; }
コメント