４０００桁の円周率を求める

前田稔(Maeda Minoru)の超初心者のプログラム入門

円周率の話

pai=
    1: 3.1415 9265 3589 7932 3846 2643 3832 7950 2884 1971 
   41:   6939 9375 1058 2097 4944 5923 0781 6406 2862 0899 
   81:   8628 0348 2534 2117 0679 8214 8086 5132 8230 6647 
  121:   0938 4460 9550 5822 3172 5359 4081 2848 1117 4502 
  161:   8410 2701 9385 2110 5559 6446 2294 8954 9303 8196 
  201:   4428 8109 7566 5933 4461 2847 5648 2337 8678 3165 
  241:   2712 0190 9145 6485 6692 3460 3486 1045 4326 6482 
  281:   1339 3607 2602 4914 1273 7245 8700 6606 3155 8817 
  321:   4881 5209 2096 2829 2540 9171 5364 3678 9259 0360 
  361:   0113 3053 0548 8204 6652 1384 1469 5194 1511 6094 

プロジェクトの作成

1. ソースプログラムです。
   

   ファイル名	説明
   Dos_PAI.cpp	円周率

2. 新規プロジェクトから、空の Console Application を作成して下さい。
   
3. 計算には次の Ｍａｃｈｉｎ の公式を使います。
   
   
4. ここで問題となるのが Arctan() 関数ですが、幸いにして収束が速く高精度で計算できる公式があります。
   この公式を利用して、自前で高精度の Arctan() 関数を作成するのが最大の難問です。
   
   
5. 円周率を計算する領域は、求める桁数に合わせて定義します。
   領域は short タイプで一個に４桁ずつ格納するので「1000＊4＝4000 桁」分確保しています。
   小数点は「０個目 pai[0] と１個目 pai[1] の間」にあるものとして計算します。
   

   #define S 1000 //配列のサイズ（Ｓ＊４桁） #define C1 3000 //第一項のループ回数 #define C2 900 //第二項のループ回数 short pai[S],t[S]; short a[S],w[S];

6. Arctan() 関数の展開式を見て下さい。
   4000 桁同士の加算や減算は簡単でしょう。
   y[S] と z[S] を足して x[S] に格納します。
   

   //★ x[]= y[] + z[] void add(short x[], short y[], short z[])

   y[S] から z[S] を引いて x[S] に格納します。
   

   //★ x[]= y[] - z[] void sub(short x[], short y[], short z[])

7. 4000 桁を一桁の整数で掛けるのは比較的簡単です。
   y[S] を一桁の n で掛けて x[S] に格納します。
   

   //★ x[]= y[] * n void mlt(short x[], short y[], short n)

8. 最もやっかいなのが次の計算です。
   

   X**3, X**5, X**7, X**9, ...

   但し、X の値は 1/5 と 1/239 なので、4000 桁を三桁以内の整数で割る関数を作成します。
   y[S] を三桁以内の整数 n で割って x[S] に格納します。
   

   //★ x[]= y[] / n void dev(short x[], short y[], short n)

9. X が 1/5 のときは、4000 桁のデータを 25 で割って次の項を求めます。
   X が 1/239 のときは、4000 桁のデータを 239 で二回割って次の項を求めます。
   
10. １９８９年１２月に PC-9801UV(10MHZ) で実行したときは１７分間以上かかったのですが今では？ (^_^;
    ウインドウに表示された実行の結果です。
    計算精度の関係で最後の桁は誤差が出るので、何桁まで合っているか確認して下さい。
    

    pai= 1: 3.1415 9265 3589 7932 3846 2643 3832 7950 2884 1971 41: 6939 9375 1058 2097 4944 5923 0781 6406 2862 0899 81: 8628 0348 2534 2117 0679 8214 8086 5132 8230 6647 121: 0938 4460 9550 5822 3172 5359 4081 2848 1117 4502 161: 8410 2701 9385 2110 5559 6446 2294 8954 9303 8196 201: 4428 8109 7566 5933 4461 2847 5648 2337 8678 3165 241: 2712 0190 9145 6485 6692 3460 3486 1045 4326 6482 281: 1339 3607 2602 4914 1273 7245 8700 6606 3155 8817 321: 4881 5209 2096 2829 2540 9171 5364 3678 9259 0360 361: 0113 3053 0548 8204 6652 1384 1469 5194 1511 6094 401: 3305 7270 3657 5959 1953 0921 8611 7381 9326 1179 441: 3105 1185 4807 4462 3799 6274 9567 3518 8575 2724 481: 8912 2793 8183 0119 4912 9833 6733 6244 0656 6430 521: 8602 1394 9463 9522 4737 1907 0217 9860 9437 0277 561: 0539 2171 7629 3176 7523 8467 4818 4676 6940 5132 601: 0005 6812 7145 2635 6082 7785 7713 4275 7789 6091 641: 7363 7178 7214 6844 0901 2249 5343 0146 5495 8537 681: 1050 7922 7968 9258 9235 4201 9956 1121 2902 1960 721: 8640 3441 8159 8136 2977 4771 3099 6051 8707 2113 761: 4999 9998 3729 7804 9951 0597 3173 2816 0963 1859 801: 5024 4594 5534 6908 3026 4252 2308 2533 4468 5035 841: 2619 3118 8171 0100 0313 7838 7528 8658 7533 2083 881: 8142 0617 1776 6914 7303 5982 5349 0428 7554 6873 921: 1159 5628 6388 2353 7875 9375 1957 7818 5778 0532 961: 1712 2680 6613 0019 2787 6611 1959 0921 6420 1989 1001: 3809 5257 2010 6548 5863 2788 6593 6153 3818 2796 1041: 8230 3019 5203 5301 8529 6899 5773 6225 9941 3891 1081: 2497 2177 5283 4791 3151 5574 8572 4245 4150 6959 1121: 5082 9533 1168 6172 7855 8890 7509 8381 7546 3746 1161: 4939 3192 5506 0400 9277 0167 1139 0098 4882 4012 1201: 8583 6160 3563 7076 6010 4710 1819 4295 5596 1989 1241: 4676 7837 4494 4825 5379 7747 2684 7104 0475 3464 1281: 6208 0466 8425 9069 4912 9331 3677 0289 8915 2104 1321: 7521 6205 6966 0240 5803 8150 1935 1125 3382 4300 1361: 3558 7640 2474 9647 3263 9141 9927 2604 2699 2279 1401: 6782 3547 8163 6009 3417 2164 1219 9245 8631 5030 1441: 2861 8297 4555 7067 4983 8505 4945 8858 6926 9956 1481: 9092 7210 7975 0930 2955 3211 6534 4987 2027 5596 1521: 0236 4806 6549 9119 8818 3479 7753 5663 6980 7426 1561: 5425 2786 2551 8184 1757 4672 8909 7777 2793 8000 1601: 8164 7060 0161 4524 9192 1732 1721 4772 3501 4144 1641: 1973 5685 4816 1361 1573 5255 2133 4757 4184 9468 1681: 4385 2332 3907 3941 4333 4547 7624 1686 2518 9835 1721: 6948 5562 0992 1922 2184 2725 5025 4256 8876 7179 1761: 0494 6016 5346 6804 9886 2723 2791 7860 8578 4383 1801: 8279 6797 6681 4541 0095 3883 7863 6095 0680 0642 1841: 2512 5205 1173 9298 4896 0841 2848 8626 9456 0424 1881: 1965 2850 2221 0661 1863 0674 4278 6220 3919 4945 1921: 0471 2371 3786 9609 5636 4371 9172 8746 7764 6575 1961: 7396 2413 8908 6583 2645 9958 1339 0478 0275 9009 2001: 9465 7640 7895 1269 4683 9835 2595 7098 2582 2620 2041: 5224 8940 7726 7194 7826 8482 6014 7699 0902 6401 2081: 3639 4437 4553 0506 8203 4962 5245 1749 3996 5143 2121: 1429 8091 9065 9250 9372 2169 6461 5157 0985 8387 2161: 4105 9788 5959 7729 7549 8930 1617 5392 8468 1382 2201: 6868 3868 9427 7415 5991 8559 2524 5953 9594 3104 2241: 9972 5246 8084 5987 2736 4469 5848 6538 3673 6222 2281: 6260 9912 4608 0512 4388 4390 4512 4413 6549 7627 2321: 8079 7715 6914 3599 7700 1296 1608 9441 6948 6855 2361: 5848 4063 5342 2072 2258 2848 8648 1584 5602 8506 2401: 0168 4273 9452 2674 6767 8895 2521 3852 2549 9546 2441: 6672 7823 9864 5659 6116 3548 8623 0577 4564 9803 2481: 5593 6345 6817 4324 1125 1507 6069 4794 5109 6596 2521: 0940 2522 8879 7108 9314 5669 1368 6722 8748 9405 2561: 6010 1503 3086 1792 8680 9208 7476 0917 8249 3858 2601: 9009 7149 0967 5985 2613 6554 9781 8931 2978 4821 2641: 6829 9894 8722 6588 0485 7564 0142 7047 7555 1323 2681: 7964 1451 5237 4623 4364 5428 5844 4795 2658 6782 2721: 1051 1413 5473 5739 5231 1342 7166 1021 3596 9536 2761: 2314 4295 2484 9371 8711 0145 7654 0359 0279 9344 2801: 0374 2007 3105 7853 9062 1983 8744 7808 4784 8968 2841: 3321 4457 1386 8751 9435 0643 0218 4531 9104 8481 2881: 0053 7061 4680 6749 1927 8191 1979 3995 2061 4196 2921: 6342 8754 4406 4374 5123 7181 9217 9998 3910 1591 2961: 9561 8146 7514 2691 2397 4894 0907 1864 9423 1961 3001: 5679 4520 8095 1465 5022 5231 6038 8193 0142 0937 3041: 6213 7855 9566 3893 7787 0830 3906 9792 0773 4672 3081: 2182 5625 9966 1501 4215 0306 8038 4477 3454 9202 3121: 6054 1466 5925 2014 9744 2850 7325 1866 6002 1324 3161: 3408 8190 7104 8633 1734 6496 5145 3905 7962 6856 3201: 1005 5081 0665 8796 9981 6357 4736 3840 5257 1459 3241: 1028 9706 4140 1109 7120 6280 4390 3975 9515 6771 3281: 5770 0420 3378 6993 6007 2305 5876 3176 3594 2187 3321: 3125 1471 2053 2928 1918 2618 6125 8673 2157 9198 3361: 4148 4882 9164 4706 0957 5270 6957 2209 1756 7116 3401: 7229 1098 1690 9152 8017 3506 7127 4858 3222 8718 3441: 3520 9353 9657 2512 1083 5791 5136 9882 0914 4421 3481: 0067 5103 3467 1103 1412 6711 1369 9086 5851 6398 3521: 3150 1970 1651 5116 8517 1437 6576 1835 1556 5088 3561: 4909 9898 5998 2387 3455 2833 1635 5076 4791 8535 3601: 8932 2618 5489 6321 3293 3089 8570 6420 4675 2590 3641: 7091 5481 4165 4985 9461 6371 8027 0981 9943 0992 3681: 4488 9575 7128 2890 5923 2332 6097 2997 1208 4433 3721: 5732 6548 9382 3911 9325 9746 3667 3058 3604 1428 3761: 1388 3032 0382 4903 7589 8524 3744 1702 9132 7656 3801: 1809 3773 4440 3070 7469 2112 0191 3020 3303 8019 3841: 7621 1011 0044 9293 2151 6084 2444 8596 3766 9838 3881: 9522 8684 7831 2355 2658 2131 4495 7685 7262 4334 3921: 4189 3039 6864 2624 3410 7732 2697 8028 0731 8915 3961: 4411 0104 4682 3252 7162 0105 2652 2721 0892

超初心者のプログラム入門(C/C++)