多倍長整数演算，世界最速の円周率の計算法
円周率を4種類の方法で10億桁求めまし
た
円周率を求めるアルキメデスの方法
の加速
アルキメデスは，円に内接する正n角形と外接する正n角形の周の長さを用いて，
円周率の近似値を求める公式を幾何学的に導きましたが，ここでは，三角関数を用いました。さらに，円弧を放物線で近似して面積を求めることにより，半径1
の円に内接
する正n角形の周の長さの半分をpnとするときに(4pn-pn/2)/3が円周率の近似の加速することなどを見つけ， これらを用いて15°までの
三角比で円周率がほぼ3.14であることを示しました。
さらに発展させて，円周率を1000桁計算するのに僅か254角形まででできる加速法を見つけましたが，その加速法は，建部賢弘が用いた加速法と同等で，リチャードソン補外のようです。次の十進BASICのプログラムを実行してみて下さい。ただし，アルキメデスの方法とその加速版は，さらに大きな桁数を計算するには不向きです。
!アルキメデスの方法
OPTION ARITHMETIC DECIMAL_HIGH
LET p=2*SQR(2)
LET q=4
FOR i=1 TO 1660
LET q=2*p*q/(p+q)
LET p=SQR(p*q)
PRINT i,PI-p
NEXT i
END
!アルキメデスの方法の加速版
OPTION ARITHMETIC DECIMAL_HIGH
OPTION BASE 0
LET N=52
DIM p(N)
LET p(0)=2*SQR(2)
LET q=4
FOR i=1 TO N
LET temp1=p(0)
LET q=2*p(0)*q/(p(0)+q)
LET p(0)=SQR(p(0)*q)
LET d=1
FOR j=1 TO i
LET d=d*4
LET temp2=p(j)
LET p(j)=p(j-1)+(p(j-1)-temp1)/(d-1)
LET temp1=temp2
NEXT J
PRINT i,PI-p(i)
NEXT i
END
多倍長整数演算ルーチン
公開鍵暗号に使うために作成したもので，それに使う関数しか作っていません。筆算の方法で計算していますので，乗算と除算の実行
速度はO(n2)
ですが，2000bit程度ならば十分な速さだと思います。自由に改造してご利用下さい。基数は216で，例えば
unsigned long A[1000];
に多倍長整数を入れる場合，A[1]以降の下位16bitにデータを入れ，A[0]にデータ長を入れる構造になっています。
除算のアルゴリズム
multiple.pdf
多倍長整数演算用関数 multi.c
ニュートン法と類似のアルゴリズムによる円周率の計算法
世界記録を出した高速な円周率の計算法は，大きく分けて「ラマヌジャン型公式
(Chudnovskyの公式)」，「逆正接(arctan)公式」，「ガウス-ルジャンドルの算術幾何平均」の3種類があります。最近の世界記録の更新はどれもChudnovskyの公式によるもので，2022年6月には、Googleが100兆桁を計算し世界記録を更新したと発表しました。
ここでは，新たな円
周率の計算法として，ニュートン法と類似のアルゴリズム（以下，三角関数法）
を提案します。この方法がどこまで高速化できるかを工夫してみました。
円周率を
三角関数法でn桁まで求めるときの計算量は，O(n(log n)3)
程度までなります。
私が作ったシングルタスクのプログラムでは、メモリ8GBで10億桁の計算ができ，速い順に
Chudnovskyの公式(ラマヌジャン型公式)， ガウス-ルジャンドルの算術幾何平均， 逆正接(arctan)公式， 三角関数法
となり，三角関数法の計算時間は，Chudnovskyの公式の4.5倍，算術幾何平均の2倍，逆正接公式の1.2倍程度でした。
Chudnovskyの公式と逆正接公式は，分母分子の共通因数を約分することにより高速化できるようですが，三角関数法はそれができないため，速度差はさら
に拡がります。
また，ガウス-ルジャンドルの算術幾何平均の計算量はO(n(log n)2)
程度で，他のアルゴリズムはO(n(log n)3)
程度なので，n→∞の極限では，ガウス-ルジャンドルの算術幾何平均が最速になるようです。
他より遅い三角関数法ですが、利点は使用メモリがO(n)であることと，分散処理（並列処理）がやり易いことです。しかし，並列処理をするには、より多くのメモリが必要です。
三角関数法は，
他より遅いけれどメモリ8GBで10億桁が計算できるので，円周率の高速な計算法の1つに加えていただければ嬉しいです。
三角関数法は，次の補題においてcosxとarcsinxのテーラー展開を利用することにより円周率を求める方法です。
補題 0≦θ≦π とするとき，θ＋arcsin(cosθ)＝π/2
証明
0≦θ≦πより，-π/2≦π/2－θ≦π/2 であるから
arcsin(cosθ)＝arcsin(sin(π/2－θ))＝π/2－θ
よって θ+arcsin(cosθ)＝π/2
証明終
同様に，
π/2≦θ≦3π/2 とするとき，arcsin(sinθ)＝arcsin(sin(π-θ))＝π-θ より θ+arcsin(sinθ)＝π
0<θ<π とするとき，arctan(cotθ)＝arctan(tan(π/2-θ))＝π/2-θ より θ+arctan(cotθ)＝π/2
が成り立つ。
例 円周率の近似値を上記の補題を用いて計算してみます。π/2≒3/2ですから
PRINT (3/2+ASIN(COS(3/2)))*2
で求まりますが，ここではテーラー展開を用いて，四則演算だけで求めてみます。
cos 3/2 は，テーラー展開の項数を少なくするために，まず cos 3/64 を求め，2倍角の公式を5回用いて cos 3/2 を求めます。
LET pi2=3/2
LET x=pi2/32
LET c=1-x^2/2+x^4/24-x^6/720+x^8/40320 ! cos x のテーラー展開
FOR i=1 TO 5
LET c=2*c^2-1 ! 2倍角の公式
next I
PRINT 2*(pi2+c+c^3/6+c^5*3/40+c^7*5/112+c^9*35/1152) ! 2*(pi2+arcsin c のテーラー展開)
END
3.14159265358984
定理 初期値を θ≒π/2 として
θ：＝θ＋cosθ
を繰り返すと，θは3次収束でπ/2に近づく。
証明
-1＜x＜1のとき，
arcsin(x)をテーラー展開すると，
arcsin(x)=x＋1/6・x^3＋3/40・x^5＋5/112・x^7＋……
よって，0<θ<πとするとき，-1＜cosθ＜1であるから，補題により次の公式が得られる。
π/2＝θ＋cosθ＋1/6・(cosθ)^3＋3/40・(cosθ)^5＋5/112・(cosθ)^7＋……
この公式において，ε＝π/2－θ，ε'＝π/2－(θ＋cosθ) とおき，整理すると
ε'＝1/6・(sinε)^3＋3/40・(sinε)^5＋5/112・(sinε)^7＋……
ここで，|sinε|≦|ε|であり，右辺の係数は正で単調に減少するから
|ε'|≦1/6・|ε|^3＋3/40・|ε|^5＋5/112・|ε|^7＋……≦1/6・(|ε|^3＋|ε|^5＋|ε|^7＋……)
よって，|ε|<1のとき
|ε'|≦1/6・|ε|^3/(1-|ε|^2)
したがって，|ε|≦1/2であれば
|ε'|≦2/9・|ε|^3
初期値を θ≒π/2 とすると，ε≒0であるから
θ：＝θ＋cosθ
を繰り返すと，εは3次収束で0に近づき，θは3次収束でπ/2に近づく。
証明終
同様に
θ：＝θ＋cosθ＋1/6・(cos θ)^3
を繰り返すと，θは5次収束で
θ：＝θ＋cosθ＋1/6・(cosθ)^3＋3/40・(cosθ)^5
を繰り返すと，θは7次収束でπ/2に近づく。
また，
初期値を θ≒πとして
θ：＝θ＋sinθ
を繰り返すと，θは3次収束でπに近づく。(これをBeelerの公式と言うようです)
なお，ニュートン法によって
方程式 cos x=0 (0<x<π) の解 x=π/2 を求める漸化式は，xn+1=xn+cos xn/sin xn
となってしまいますが
方程式 cos x/(1+sin x)=0 (-π/2<x<3π/2) の解 x=π/2 を求める漸化式は，xn+1=xn+cos xn
となります。
円周率を求めるBASICのプログラムは，次のようになります。
LET pi2=3/2
FOR i=1 TO 3
LET pi2=pi2+COS(pi2)
PRINT 2*pi2
NEXT i
END
三角関数法pi.pdf
漸化式の加法定理による変形
漸化式 xn+1=xn+cos xn ではcos xnを計算することが必須ですが，10億桁を求めるプログラムにおいてcos xn を高速に計算するために，この漸化式を加法定理を用いて変形します。
sin xn+1=sin (xn+xn+1-xn)=sin xn cos(xn+1-xn)+cos xn sin(xn+1-xn)
cos xn+1=cos (xn+xn+1-xn)=cos xn cos(xn+1-xn)-sin xn sin(xn+1-xn)
ここで，xn+1-xn=cos xn であるから
sin xn+1=sin xn cos(cos xn)+cos xn sin(cos xn)
cos xn+1=cos xn cos(cos xn)-sin xn sin(cos xn)
初期値を x1 として sin x1，cos x1 は直接計算する。このとき，
xn+1=x1+Σcos xk
である。
10億桁を計算するとき，後に解説するBSA法では，cos(1/√自然数) は高速に計算することができます。下記のプログラムは
π/2=1+1/√3-1/√23281+……
と分解されて π/2 の近似値が求まります。
LET pi2=1
LET s=SIN(pi2)
LET c=COS(pi2)
DO
LET QQ=INT(1/c^2+1/6)
IF c>0 THEN LET e=1/SQR(QQ) ELSE LET e=-1/SQR(QQ)
LET pi2=pi2+e
PRINT pi2*2
LET s1=SIN(e)
LET c1=COS(e)
LET c2=c*c1-s*s1
IF ABS(c2)<EPS(1) THEN EXIT DO
LET s=s*c1+c*s1
LET c=c2
LOOP
END
三角関数法の実装
円周率の世界記録に挑戦するときには，高速フーリエ変換を用いた乗算の計算が必須ですが，これは，多倍長任意精度の算術ライブラリに組み込まれているので，知る必要はありません。
Chudnovskyの公式や三角関数法では，さらにBSA法（Binary Splitting Algorithm）が必須です。三角関数法では，ニュートン法も必須です。
これらを用いたとき，三角関数法で円周率をn桁まで求めるときの計算量はO(n (log
n)3) 程度です。世界記録達成に用いられるChudnovsky
の公式によるの計算量もO(n (log
n)3)で，オーダーとしては同じですが，Chudnovskyの公式は n (log
n)3 に掛かる係数が小さいため，この計
算法は計算量では敵いません。
BSA法は，級数の部分和や連分数などを高
速に計算する方法です。色々なサイトを参考にして、BSAについて解説（BSA.pdf）を書きました。合成関数を用いて独自な証明をしましたので，読んで頂ければ嬉しいです。また，
BSA法を用いるときに効率が良いように，次のような工夫をしてみました。
① π/2を，符号付きの分子が1で，分母が自然数の平方根である分数（エジプト式
分数の変形）に分
解します。(符号付き強欲算法)
π/2=1+1/√3-1/√23281-1/√372082905075224+1/√2579681866751687411045599324489688909153448939+……
そして，qを自然数としてcos(1/√q)を計算していきます。分母が自然数の場合は2次収束となりますが，この場合3次収束し，項の数が少なくて済み
ます。
cos(1/√q)の計算は下で解説するBinary Splitting
Algorithmと非常に相性が良く，高速かつ省メモリに計算することができます。
このとき，1つのcos(1/√q)の計算量はO(n (log
n)2) 程度です。qの値が log ｎ個程度あるため，全体の計算量はO(n (log
n)3) 程度となります。
② sin(1/√q)は公式sinθ=√(1-cos2θ)を用いてcos(1/√q)から計算します。なお，余弦の加法定理を用いる部分は，次の式変形から分かるように，多倍長の積は２回，平方根は１回で済みます。
cos(α+β)=cosαcosβ-sinαsinβ=cosαcosβ±√(1-cos2α)(1-cos2β)=cosαcosβ±√(1-cosα)(1-cosβ)(1+cosα)(1+cosβ)=cosαcosβ±√{1-(cosα+cosβ)+cosαcosβ}{1+(cosα+cosβ)+cosαcosβ}
ここで，s=cosαcosβ，t=1-(cosα+cosβ)+cosαcosβ とおくと，次のようになる。
cos(α+β)=s±√t(2s-t+2)
! 三角関数法
OPTION ARITHMETIC DECIMAL_HIGH
LET pi2=1
LET n=INT(kosuu(1))
CALL BSA(P,Q,1,0,n-1)
LET y=1-P/Q
DO
LET QQ=INT(1/y^2+1/6)
LET n=INT(kosuu(QQ))
CALL BSA(P,Q,QQ,0,n-1)
LET z=P/Q
LET e=y*z
LET s=y-e
LET t=z-e
IF y>0 THEN
LET y=s-SQR(t*(2*s-t+2))
LET pi2=pi2+1/SQR(QQ)
PRINT QQ
ELSE
LET y=s+SQR(t*(2*s-t+2))
LET pi2=pi2-1/SQR(QQ)
PRINT -QQ
END IF
LOOP WHILE ABS(y)>=1E-333
LET pi2=pi2+y
PRINT pi2*2-PI
END
EXTERNAL FUNCTION kosuu(QQ)
OPTION ARITHMETIC DECIMAL_HIGH
LET log_Q=LOG(QQ)/2
LET x=1000
FOR i=1 TO 5
LET x=(x+LOG(10)*1000)/(LOG(x)+log_Q)
NEXT I
LET kosuu=x/2
END FUNCTION
EXTERNAL SUB BSA(P,Q,QQ,left,right)
OPTION ARITHMETIC DECIMAL_HIGH
IF right-left>0 THEN
LET middle=INT((left+right)/2)
CALL BSA(P0,Q0,QQ,left,MIDDLE)
CALL BSA(P1,Q1,QQ,MIDDLE+1,right)
LET P=P1+P0*Q1
LET Q=Q0*Q1
ELSE
LET P=-1
LET Q=-(2*left+1)*(2*left+2)*QQ
END IF
END SUB
高速フーリエ変換を用いて円周率を10億桁求めるには，GMP（The GNU Multiple Precision Arithmetic Library）などの高速フーリエ変換が予め組み込んである，多倍長任意精度の算術ライブラリが必須です。GMPは，符号付き整数，有理数，浮動小数点数を扱うことができます。
三角関数法で10億桁求める実装では，cos(1/√q)を計算する回数を減らすために，次のような工夫を追加してみました。
これらによって計算量は少しだけ減少します。
① π/2の，エジプト式分数への分解は，最初の6項を
π/2=2^25(1/√456310044896370+1/√1191109328952345107005661159230615032546365914207 +1/√……)
とし，2倍角の公式を25回用いました。
② エジプト式分数への分解を途中で止めて，補題 θ＋arcsin(cosθ)＝π/2 を用いてπ/2の値を計算しました。
BSA法とGMPを用いたCで実装したところ，8GBのメモリで10億桁が計算でき，Chudnovskyの公式の単純な実装と比べて4.5倍の
計
算時間がかかりました。
また，ガウス・ルジャンドルの算術幾何平均法と比べて約2倍の計
算時間がかかりました。
Chudnovskyの公式には，計算途中で約分するという高速化の方法がありますが，三角関数法ではそれが使えません。したがって，三角関数法は
Chudnovskyの公式と比べてかなり遅いと考えられます。
なお，GMPには 1/√a を求める関数がなく，√a を計算した後に逆数を取るのはあまりにも非効率ですから，GMPのサイトにあった
gmp-chudnovsky.c を参考にして関数 myf_1_div_sqrt_z() を作りました。
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <math.h>
#include <time.h>
void mpf_set_prec_raw_clear(mpf_t x,unsigned long prec){
mpf_set_prec_raw(x,prec);
mpf_clear(x);
}
void f_1_div_sqrt_z(mpf_t x,mpz_t a,mpf_t x0){ // x=1/SQRT(a)
unsigned long prec,prec0,bits=0;
mpf_t temp,a_t;
long prec_a;
mpz_get_d_2exp(&prec_a,a);
prec0=mpf_get_prec(x);
mpf_init2(temp,prec0);
mpf_init2(a_t,prec_a);
mpf_set_z(a_t,a);
mpf_set_prec_raw(x,prec_a);
mpf_set(x,x0);
mpf_abs(x,x);
for (prec=prec0; prec>prec_a;){
bits=(bits<<1)|(prec&1);
prec=(prec+1)/2;//prec=prec/2では精度不足になる
}
while (prec<prec0){
prec=prec*2-(bits&1);
bits=bits>>1;
// x=x+x*(1-a*x*x)/2;
mpf_set_prec_raw(temp,prec);
mpf_mul(temp,x,x);
mpf_mul(temp,temp,a_t);
mpf_ui_sub(temp,1,temp);
mpf_set_prec_raw(temp,(prec+1)/2);
mpf_div_2exp(temp,temp,1);
mpf_mul(temp,temp,x);
mpf_set_prec_raw(x,prec);
mpf_add(x,x,temp);
}
mpf_set_prec_raw_clear(temp,prec0);
mpf_clear(a_t);
}
void BSA(mpz_t P,mpz_t Q,mpz_t QQ,long left,long right) {
if (right>left){
long middle=(left+right)/2;
mpz_t P1,Q1;
mpz_init(P1);
mpz_init(Q1);
BSA(P,Q,QQ,left,middle);
BSA(P1,Q1,QQ,middle+1,right);
mpz_mul(P,P,Q1);
mpz_add(P,P,P1);
mpz_clear(P1);
mpz_mul(Q,Q,Q1);
mpz_clear(Q1);
}else{
mpz_set_si(P,-1);
mpz_mul_si(Q,QQ,-(2*left+1));
mpz_mul_ui(Q,Q,2*left+2);
}
}
long kousuu(mpz_t QQ,unsigned long digits){ // BSAに必要な項数を、ニュートン法により求める
double x=digits,a,log_Q;
long exp;
int i;
a=mpz_get_d_2exp(&exp,QQ);
log_Q=(log2(a)+exp)*log(2)/2;
for (i=0;i<5;i++){
x=(x+log(10)*digits)/(log(x)+log_Q);
}
return x/2;
}
void Round(mpz_t N,mpf_t F,double d){// 正の実数を丸めて自然数に
mpf_t temp;
mpf_init2(temp,10);
mpf_set_d(temp,d);
mpf_add(F,F,temp);
mpz_set_f(N,F);
mpf_clear(temp);
}
void COS_ADD_Q(mpf_t cos,mpz_t QQ,unsigned long digits,int sgn){
// cos=cos(1/√QQ)
mpz_t P,Q;
mpf_t p,q;
long n=kousuu(QQ,digits);
printf("%ld ",n);
mpz_init(P);
mpz_init(Q);
BSA(P,Q,QQ,0,n);
mpf_init(p);
mpf_set_z(p,P);
mpz_clear(P);
mpf_init(q);
mpf_set_z(q,Q);
mpz_clear(Q);
mpf_div(p,p,q);
if (mpf_sgn (cos)==0){
mpf_ui_sub(cos,1,p);
}else{
//余弦の加法定理
mpf_mul(q,cos,p);
mpf_sub(cos,cos,q);
mpf_sub(p,p,q);
mpf_mul_2exp(q,cos,1);
mpf_sub(q,q,p);
mpf_add_ui(q,q,2);
mpf_mul(q,q,p);
mpf_sqrt(p,q);
if (sgn>0){
mpf_sub(cos,cos,p);
}else{
mpf_add(cos,cos,p);
}
}
mpf_clear(q);
mpf_clear(p);
}
#define TIME printf("%.3f s\n",(double)(clock()-start)/CLOCKS_PER_SEC);
#define PI2n 25
void Pi2(mpf_t pi2pre,mpf_t pi2,unsigned long digits,clock_t start) {
long exp;
int i,sgn;
mpf_div_2exp(pi2pre,pi2pre,PI2n);
unsigned long prec_t,prec=(digits+10)*log2(10);
mpf_set_default_prec(prec);
mpf_t cos,temp,temp2;
mpz_t QQ;
mpf_init(cos);
mpf_init(pi2);
mpz_init(QQ);
unsigned long precpre=mpf_get_prec (pi2pre);
mpf_init2(temp,precpre);
do{
mpf_mul(temp,pi2pre,pi2pre);
mpf_ui_div(temp,1,temp);
Round(QQ,temp,0.5);// tempを四捨五入してQQに代入
mpf_init(temp2);
f_1_div_sqrt_z(temp2,QQ,pi2pre);
sgn=mpf_sgn(pi2pre);
if (sgn>0){
mpf_add(pi2,pi2,temp2);
mpf_sub(pi2pre,pi2pre,temp2);
}else{
mpf_sub(pi2,pi2,temp2);
mpf_add(pi2pre,pi2pre,temp2);
}
mpf_clear(temp2);
COS_ADD_Q(cos,QQ,digits,sgn);
TIME
mpf_get_d_2exp(&exp,pi2pre);
}while (-3*exp<precpre);
mpf_clear(pi2pre);
mpf_clear(temp);
for (i=0;i<PI2n;i++){
mpf_mul(cos,cos,cos);
mpf_mul_2exp(cos,cos,1);
mpf_sub_ui(cos,cos,1);
}
mpf_mul_2exp(pi2,pi2,PI2n);
TIME
do {
mpf_get_d_2exp(&exp,cos);
prec_t=-exp;
mpf_set_default_prec(prec-prec_t);
mpf_set_prec_raw(cos,prec-prec_t);
mpf_init2(temp,prec_t*2+10);
mpf_mul(temp,cos,cos);
mpf_ui_div(temp,1,temp);
Round(QQ,temp,1.0/6.0);
mpf_clear(temp);
mpf_init(temp2);
f_1_div_sqrt_z(temp2,QQ,cos);
sgn=mpf_sgn(cos);
if (sgn>0){
mpf_add(pi2,pi2,temp2);
}else{
mpf_sub(pi2,pi2,temp2);
}
mpf_clear(temp2);
COS_ADD_Q(cos,QQ,digits,sgn);
TIME
}while(prec_t<prec/200);
mpz_clear(QQ);
mpf_get_d_2exp(&exp,cos);
unsigned long prec_t0=-exp;
mpf_set_prec_raw(cos,prec-prec_t0);
prec_t=prec_t0*3;
unsigned long prec1=prec-prec_t;
mpf_t xx,kou;
mpf_init2(xx,prec1);
mpf_init2(kou,prec1);
mpf_mul(xx,cos,cos);
mpf_set(kou,cos);
i=1;
while(1){
mpf_mul(kou,kou,xx);
if (i!=1) mpf_mul_ui(kou,kou,(2*i-1)*(2*i-1));
mpf_div_ui(kou,kou,2*i*(2*i+1));
mpf_add(cos,cos,kou);
prec_t=prec_t+prec_t0*2;
if (prec_t>=prec) break;
mpf_set_prec_raw(xx,prec-prec_t);
mpf_set_prec_raw(kou,prec-prec_t);
i++;
}
mpf_set_prec_raw_clear(xx,prec1);
mpf_set_prec_raw_clear(kou,prec1);
mpf_add(pi2,pi2,cos);
mpf_set_prec_raw_clear(cos,prec);
TIME
}
int main(){
mpf_t pi0,pi1,pi2;
unsigned long digits=1000000000;
mpf_init2(pi0,53);//doubleの精度から53とした
mpf_set_d(pi0,1.5707963267948966);//π/2の近似値をセット
clock_t start=clock();
Pi2(pi0,pi1,6000,start);//π/2を6000桁求める
Pi2(pi1,pi2,digits,start);//π/2をdigits桁求める
mpf_mul_2exp(pi2,pi2,1);
FILE* fp;
fp=fopen("picos.txt","w");
mpf_out_str(fp,10,digits+1,pi2);
fclose(fp);
mpf_clear(pi2);
TIME
return 0;
}
三角関数法は、Chudnovskyの公式と比較して、さらに分散処理(並列処理)がしやすいという利点があります。
(Chudnovskyの公式と逆正接の公式は
Binary Splitting
Algorithmを用いているので並列化ができ，ガウス・ルジャンドルの算術幾何平均法は並列化がしづらい)
Chudnovskyの公式\[
\frac{1}{\pi}=\frac{1}{426880\sqrt{10005}}{ \sum_{n=0}^{\infty}}\frac{{
(-1)}^n (6n)!~(545140134n+13591409)}{(n!)^3(3n)!~640320^{3n}}
\]
十進BASICの単純なプログラムは次のようになります。これが世界
最速の円周率の計算法への基本となるものです。
Chudnovskyの公式によりπを計算するプログラムは次のようになりま
す。なお，作成に当たっては
http://www.kk62526.server-shared.com/pi/BinSplit.html を参考にさせて頂きました。
!Chudnovskyの公式
OPTION ARITHMETIC DECIMAL_HIGH
LET kou=1
LET s=kou*13591409
LET i=1
DO
LET
kou=-kou*(2*i-1)*(6*i-5)*(6*i-1)/(10939058860032000*i^3)
LET s=s+kou*(545140134*i+13591409)
LET i=i+1
LOOP WHILE ABS(kou)>EPS(0)
LET s=426880*SQR(10005)/s
PRINT s-PI
END
上記のプログラムをBAS法を用いて書くと次のようになります。
! Chudnovskyの公式 Binary Splitting
OPTION ARITHMETIC DECIMAL_HIGH
LET n=51 ! n個の和
CALL BSA(P,Q,R,0,n-1)
LET pai=SQR(10005)*R/Q
PRINT pai-PI
END
EXTERNAL SUB BSA(P,Q,R,left,right)
OPTION ARITHMETIC DECIMAL_HIGH
IF right-left>0 THEN
LET middle=INT((left+right)/2)
CALL BSA(P0,Q0,R0,left,MIDDLE)
CALL BSA(P1,Q1,R1,middle+1,right)
LET P=P0*P1
LET Q=P0*Q1+Q0*R1
LET R=R0*R1
ELSE
LET P=-(2*left+1)*(6*left+1)*(6*left+5)
LET Q=13591409+545140134*left
IF left=0 THEN
LET R=426880
ELSE
LET R=10939058860032000*left^3
END IF
END IF
END SUB
上記の十進BASICのプログラムをGMPを使ってCに移植したところ，8GBのメモ
リで高速に10億桁が計算できました。BSA
法による計算量はO(n
(log n)3)
程度ですが，n
(log n)3に掛かる係数が小さいため，世界最速の公式と
言われています。
10億桁の出力ファイルと
三角関数法の出力ファイルは，SHA-256によってハッシュ値が一致することを確かめました。ただし，10億桁目は
10億1桁目が四捨五入さ
れ，かつ，末尾の0は出力されないようです。
なお，作成に当たっては
https://itchyny.hatenablog.com/entry/20120304/1330870932
を参考にさせて頂きました。そのサイトのプログラムとの速度差はほとんどありません。
しかし，「y-cruncher」という円周率の計算ソフトは，以下
のプログラムと比べて，10億桁を計算したときに4.5倍速く，「y-
cruncher」の凄さを実感しました。
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <math.h>
#include <time.h>
void BSA(mpz_t T,mpz_t P,mpz_t Q,long left,long right,long end) {
if (right>left){
long middle=(left+right)/2;
mpz_t T1,P1,Q1;
mpz_init(T1);
mpz_init(P1);
mpz_init(Q1);
BSA(T,P,Q,left,middle,end);
BSA(T1,P1,Q1,middle+1,right,end);
mpz_mul(P,P,Q1);
mpz_addmul(P,T,P1);
mpz_clear(P1);
if (right==end) {
mpz_clear(T);
mpz_init(T);
}else{
mpz_mul(T,T,T1);
}
mpz_clear(T1);
mpz_mul(Q,Q,Q1);
mpz_clear(Q1);
}else{
mpz_set_si(T,-(2*left+1));
mpz_mul_ui(T,T,6*left+1);
mpz_mul_ui(T,T,6*left+5);
mpz_set_ui(P,545140134);
mpz_mul_ui(P,P,left);
mpz_add_ui(P,P,13591409);
if (left==0){
mpz_set_ui(Q,426880);
}else{
mpz_set_str(Q,"10939058860032000",10);
mpz_mul_ui(Q,Q,left);
mpz_mul_ui(Q,Q,left);
mpz_mul_ui(Q,Q,left);
}
}
}
int main () {
FILE* fp;
unsigned long digits=1000000000;
unsigned long prec=(digits+5)*log2(10);
long n=ceil(digits/14.18);
mpf_t temp,pi;
mpz_t T,P,Q;
clock_t start=clock(),end;
mpf_set_default_prec(prec);
mpz_init(T);
mpz_init(P);
mpz_init(Q);
BSA(T,P,Q,0,n,n);
mpz_clear(T);
mpf_init(pi);
mpf_sqrt_ui(pi,10005);
mpf_init(temp);
mpf_set_z(temp,Q);
mpz_clear(Q);
mpf_mul(pi,pi,temp);
mpf_set_z(temp,P);
mpz_clear(P);
mpf_div(pi,pi,temp);
mpf_clear(temp);
end=clock();
printf("%.3f s\n",(double)(end-start)/CLOCKS_PER_SEC);
fp=fopen("pi.txt","w");
mpf_out_str(fp,10,digits+1,pi);
fclose(fp);
mpf_clear(pi);
end=clock();
printf("%.3f s\n",(double)(end-start)/CLOCKS_PER_SEC);
return 0;
}
高速なπの計算法
参考のために，Ramanujanの公式と算術幾何平均法と逆正接公式と連分数による十進BASICの単純なプログラムを掲載します。
さらに，Binary Splitting Algorithmを用いてCとGMPへの移植したプログラムを示します。
次のプログラムは，Ramanujanの公式\[\frac{1}{\pi}
= \frac{\sqrt{8}}{9801} \sum_{n=0}^{\infty}
\frac{(4n)!~(26390n+1103)}{(n!)^4 ~ 396^{4n}}\]によるもので，BSA法
による計算量はO(n
(log n)3)
程度ですが，数兆桁程度を計算するときの計算量はn
(log n)2程度で，Chudnovskyの公式の次に速い公式と言われています。
!Ramanujanの公式
OPTION ARITHMETIC DECIMAL_HIGH
LET kou=1
LET s=kou*1103
LET i=1
DO
LET kou=kou*(4*i-1)*(4*i-2)*(4*i-3)/(6147814464*i^3)
LET s=s+kou*(26390*i+1103)
LET i=i+1
LOOP WHILE ABS(kou)>EPS(0)
LET s=9801/(s*SQR(8))
PRINT s-PI
END
次のプログラムは，ガウス・ルジャンドルの算術幾何平均法\[a_0=1,~b_0=\frac{1}{\sqrt{2}},~t_0=\frac{1}{4},~k_0=1
\\ a_{n+1}=\frac{a_n+b_n}{2} \\ b_{n+1}=\sqrt{a_n b_n} \\
t_{n+1}=t_n-k_n (a_n-a_{n+1})^2 \\ k_{n+1}=2k_n \\ \pi \approx
\frac{(a+b)^2}{4t} \]によるもの
です。計算量はO(n
(log n)2)
程度ですが，並列処理化は難しそうです。
! ガウス・ルジャンドルの算術幾何平均法
OPTION ARITHMETIC DECIMAL_HIGH
LET a=1
LET b=SQR(1/2)
LET t=1/4
LET k=1
DO
LET aa=(a+b)/2
LET b=SQR(a*b)
LET t=t-k*(a-aa)^2
LET a=aa
LET k=k*2
LOOP WHILE a-b>=10^(-500)
LET pai=(a+b)^2/(t*4)
PRINT pai-PI
END
次のプログラムは，私が発見した逆正接関数の6項最良の公式\[\frac
{\pi}
{4}=83\arctan \frac{1}{107} +17\arctan \frac{1}{1710} -44\arctan
\frac{1}{225443} -68\arctan \frac{1}{2513489} +22\arctan
\frac{1}{42483057}+34\arctan \frac{1}{7939642926390344818} \]によるも
ので，BSA
法による計算量はO(n
(log n)3)
程度です。
!6項最良の公式
OPTION ARITHMETIC DECIMAL_HIGH
LET
pai=83*arccot(107)+17*arccot(1710)-44*arccot(225443)-68*arccot(2513489)+22*arccot(42483057)+34*arccot(7939642926390344818)
PRINT pai*4-PI
END
EXTERNAL FUNCTION arccot(x)
OPTION ARITHMETIC DECIMAL_HIGH
LET s=1/x
LET xx=-s*s
LET kou=s*xx
LET k=1
DO
LET s=s+kou/(2*k+1)
LET kou=kou*xx
LET k=k+1
LOOP WHILE ABS(kou)>EPS(0)
LET arccot=s
END FUNCTION
次のプログラムは，連分数\[\pi=4 \arctan
1=\frac{4}{1+\cfrac{1^2}{3+\cfrac{2^2}{5+\cfrac{3^2}{7+\cfrac{4^2}{9+\ddots}}}}}
\]によるもので，BSA 法による計算量はO(n (log n)3) 程度です。
! 連分数
OPTION ARITHMETIC DECIMAL_HIGH
LET x=0
FOR i=1310 TO 1 STEP -1
LET x=i^2/(2*i+1+x)
NEXT i
PRINT 4/(1+x)-PI
END
また，次の式が任意の自然数nについて成り立つようで，BSA 法で実装したときにこちらの方が高速です。\[{\pi}=4\arctan 1
=
4\left(\frac{1}{1}-\frac{1}{3}+\frac{1}{5}-{\cdots}+
\frac{(-1)^{n-1}}{2n-1}\right)+\frac{2 \cdot
(-1)^{n}}{2n+\cfrac{1^2}{2n+\cfrac{2^2}{2n+\cfrac{3^2}{2n+\cfrac{4^2}{2n+\ddots}}}}}
\]
! 連分数
OPTION ARITHMETIC DECIMAL_HIGH
LET n=654
LET s=0
LET t=1
FOR i=1 TO n
LET s=s+4*t/(2*i-1)
LET t=-t
next i
PRINT s-pi
LET x=0
FOR i=n TO 1 STEP -1
LET x=i^2/(2*n+x)
NEXT i
LET x=2*t/(2*n+x)
PRINT s+x-pi
END
Ramanujanの公式による計算法のCとGMPへの移植
Ramanujanの公式による
プログラムをBinary Splitting
Algorithmを用いてGMPを使ってCで書いたところ，8GBのメモリではメモリがギリギリの状態で10億桁を計算できました。
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <math.h>
#include <time.h>
void BSA(mpz_t T,mpz_t P,mpz_t Q,long left,long right,long end) {
if (right>left){
long middle=(left+right)/2;
mpz_t T1,P1,Q1;
mpz_init(T1);
mpz_init(P1);
mpz_init(Q1);
BSA(T,P,Q,left,middle,end);
BSA(T1,P1,Q1,middle+1,right,end);
mpz_mul(P,P,Q1);
mpz_addmul(P,T,P1);
mpz_clear(P1);
if (right==end) {
mpz_clear(T);
mpz_init(T);
}else{
mpz_mul(T,T,T1);
}
mpz_clear(T1);
mpz_mul(Q,Q,Q1);
mpz_clear(Q1);
}else{
mpz_set_ui(T,2*left+1);
mpz_mul_ui(T,T,4*left+1);
mpz_mul_ui(T,T,4*left+3);
mpz_set_ui(P,26390);
mpz_mul_ui(P,P,left);
mpz_add_ui(P,P,1103);
if (left==0){
mpz_set_ui(Q,9801);
}else{
mpz_set_ui(Q,3073907232);
mpz_mul_ui(Q,Q,left);
mpz_mul_ui(Q,Q,left);
mpz_mul_ui(Q,Q,left);
}
}
}
int main () {
FILE* fp;
unsigned long digits=1000000000;
unsigned long prec=(digits+5)*log2(10);
long n=ceil(digits/7.98);
mpf_t temp,pi;
mpz_t T,P,Q;
clock_t start=clock(),end;
mpf_set_default_prec(prec);
mpz_init(T);
mpz_init(P);
mpz_init(Q);
BSA(T,P,Q,0,n,n);
mpz_clear(T);
mpf_init(pi);
mpf_sqrt_ui(pi,8);
mpf_init(temp);
mpf_set_z(temp,P);
mpz_clear(P);
mpf_mul(pi,pi,temp);
mpf_set_z(temp,Q);
mpz_clear(Q);
mpf_div(pi,temp,pi);
mpf_clear(temp);
end=clock();
printf("%.3f s\n",(double)(end-start)/CLOCKS_PER_SEC);
fp=fopen("pi.txt","w");
mpf_out_str(fp,10,digits+1,pi);
fclose(fp);
mpf_clear(pi);
end=clock();
printf("%.3f s\n",(double)(end-start)/CLOCKS_PER_SEC);
return 0;
}
ガウス・ルジャンドルの算術幾何平均法のCとGMPへの移植
十進BASICの算術幾何平均法のプログラムをThe GNU Multiple Precision Arithmetic Library
(GMP)を使ってCに移植したところ，8GBのメモリで10億桁が計算できました。
なお，古典的な円周率計算ソフトである「スーパーπ」も，ガウス・ルジャンドルの算術幾何平均法を用いて計算しています。
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <math.h>
#include <time.h>
int main(){
mpf_t a,b,t,aa,temp;
clock_t start=clock(),end;
long exp;
unsigned long digits=1000000000;
unsigned long prec=(digits+10)*log2(10);
mpf_set_default_prec (prec);
mpf_init(a);
mpf_init(b);
mpf_init(t);
mpf_init(aa);
mpf_init(temp);
mpf_set_ui(a,1);
mpf_sqrt_ui(b,2);
mpf_div_2exp(b,b,1);
mpf_div_2exp(t,a,2);
int k=0;
printf("%d ",k);
end=clock();
printf("%.3f s\n",(double)(end-start)/CLOCKS_PER_SEC);
do{
mpf_add(aa,a,b);
mpf_div_2exp(aa,aa,1);
mpf_mul(temp,a,b);
mpf_sqrt(b,temp);
mpf_sub(temp,a,aa);
mpf_get_d_2exp(&exp,temp);
mpf_mul(temp,temp,temp);
mpf_mul_2exp(temp,temp,k);
mpf_sub(t,t,temp);
mpf_set(a,aa);
k++;
printf("%d ",k);
end=clock();
printf("%.3f s\n",(double)(end-start)/CLOCKS_PER_SEC);
}while (-exp*4<=prec);
mpf_add(temp,a,b);
mpf_mul(temp,temp,temp);
mpf_div(temp,temp,t);
mpf_div_2exp(temp,temp,2);
mpf_clear(a);
mpf_clear(b);
mpf_clear(t);
mpf_clear(aa);
end=clock();
printf("%.3f s\n",(double)(end-start)/CLOCKS_PER_SEC);
FILE* fp;
fp=fopen("pi.txt","w");
mpf_out_str(fp,10,digits+1,temp);
fclose(fp);
mpf_clear(temp);
end=clock();
printf("%.3f s\n",(double)(end-start)/CLOCKS_PER_SEC);
return 0;
}
逆正接(arctan)公式による計算法のCとGMPへの移植
逆正接公式のプログラムをBinary Splitting
Algorithmを用いてGMPを使ってCで書いたところ，8GBのメモリで10億桁が計算できました。
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <math.h>
#include <time.h>
void BSA(mpz_t T,mpz_t P,mpz_t Q,long p,mpz_t QZ,mpz_t QQ,long left,long right,long end) {
if (right>left){
long middle=(left+right)/2;
mpz_t T1,P1,Q1;
mpz_init(T1);
mpz_init(P1);
mpz_init(Q1);
BSA(T,P,Q,p,QZ,QQ,left,middle,end);
BSA(T1,P1,Q1,p,QZ,QQ,middle+1,right,end);
mpz_mul(P,P,Q1);
mpz_addmul(P,T,P1);
mpz_clear(P1);
if (right==end) {
mpz_clear(T);
mpz_init(T);
}else{
mpz_mul(T,T,T1);
}
mpz_clear(T1);
mpz_mul(Q,Q,Q1);
mpz_clear(Q1);
}else{
if (left==0){
mpz_set_si(T,-4*p);
mpz_set_si(P,4*p);
mpz_set(Q,QZ);
}else{
mpz_set_si(T,-2*left-1);
mpz_set_ui(P,1);
mpz_mul_ui(Q,QQ,2*left+1);
}
}
}
void ARCTAN0(mpf_t arctan,long p,char *q,unsigned long digits){
mpz_t P,Q,T,QZ,QQ;
mpf_t temp;
mpz_init(QZ);
mpz_init(QQ);
mpz_set_str(QZ,q,10);
mpz_mul(QQ,QZ,QZ);
long exp;
double a=mpz_get_d_2exp(&exp,QQ);
long n=digits/((exp+log2(a))*log10(2));
mpz_init(P);
mpz_init(Q);
mpz_init(T);
BSA(T,P,Q,p,QZ,QQ,0,n,n);
mpz_clear(QZ);
mpz_clear(QQ);
mpz_clear(T);
mpf_init(arctan);
mpf_set_z(arctan,P);
mpz_clear(P);
mpf_init(temp);
mpf_set_z(temp,Q);
mpz_clear(Q);
mpf_div(arctan,arctan,temp);
mpf_clear(temp);
}
void ARCTAN(mpf_t pi,long p,char *q,unsigned long digits){
mpf_t arctan;
ARCTAN0(arctan,p,q,digits);
mpf_add(pi,pi,arctan);
mpf_clear(arctan);
}
void TIME(clock_t start){
clock_t end=clock();
printf("%.3f s\n",(double)(end-start)/CLOCKS_PER_SEC);
}
int main () {
FILE* fp;
unsigned long digits=1000000000;
unsigned long prec=(digits+5)*log2(10);
mpf_t pi;
clock_t start=clock();
mpf_set_default_prec(prec);
ARCTAN0(pi,83,"107",digits);
TIME(start);
ARCTAN(pi,17,"1710",digits);
TIME(start);
ARCTAN(pi,-44,"225443",digits);
TIME(start);
ARCTAN(pi,-68,"2513489",digits);
TIME(start);
ARCTAN(pi,22,"42483057",digits);
TIME(start);
ARCTAN(pi,34,"7939642926390344818",digits);
TIME(start);
fp=fopen("pi.txt","w");
mpf_out_str(fp,10,digits+1,pi);
fclose(fp);
mpf_clear(pi);
TIME(start);
return 0;
}
連分数による計算法のCとGMPへの移植
連分数のプログラムをBinary Splitting
Algorithmを用いてGMPを使ってCで書いたところ，8GBのメモリで1億桁が計算できましたが，10億桁はメモリ不足のため計算できませんでし
た。
以下のプログラムは2つの連分数のうちの後者の方のものです。この計算法は，メモリを食い，かつ，遅くて，世界記録の挑戦には使われていません。
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <math.h>
#include <time.h>
void BSA0(mpz_t T,mpz_t P,mpz_t Q,long left,long right,long end) {
if (right>left){
long middle=(left+right)/2;
mpz_t T1,P1,Q1;
mpz_init(T1);
mpz_init(P1);
mpz_init(Q1);
BSA0(T,P,Q,left,middle,end);
BSA0(T1,P1,Q1,middle+1,right,end);
mpz_mul(P,P,Q1);
mpz_addmul(P,T,P1);
mpz_clear(P1);
if (right==end) {
mpz_clear(T);
mpz_init(T);
}else{
mpz_mul(T,T,T1);
}
mpz_clear(T1);
mpz_mul(Q,Q,Q1);
mpz_clear(Q1);
}else{
if (left==0){
mpz_set_si(T,-4);
mpz_set_si(P,4);
mpz_set_ui(Q,1);
}else{
mpz_set_si(T,-2*left-1);
mpz_set_ui(P,1);
mpz_set_ui(Q,2*left+1);
}
}
}
void BSA(mpz_t T,mpz_t U,mpz_t P,mpz_t Q,long left,long right,long end) {
if (right>left){
long middle=(left+right)/2;
mpz_t T1,U1,P1,Q1,TEMP;
mpz_init(T1);
mpz_init(U1);
mpz_init(P1);
mpz_init(Q1);
BSA(T,U,P,Q,left,middle,end);
BSA(T1,U1,P1,Q1,middle+1,right,end);
if (right==end){
mpz_clear(T1);
mpz_clear(U1);
mpz_mul(P,P,Q1);
mpz_addmul(P,T,P1);
mpz_clear(T);
mpz_init(T);
mpz_mul(Q,Q,Q1);
mpz_clear(Q1);
mpz_addmul(Q,U,P1);
mpz_clear(P1);
mpz_clear(U);
mpz_init(U);
}else{
mpz_init(TEMP);
mpz_mul(TEMP,T,T1);
mpz_addmul(TEMP,P,U1);
mpz_mul(U1,Q,U1);
mpz_addmul(U1,U,T1);
mpz_clear(T1);
mpz_mul(P,P,Q1);
mpz_addmul(P,T,P1);
mpz_set(T,TEMP);
mpz_clear(TEMP);
mpz_mul(Q,Q,Q1);
mpz_clear(Q1);
mpz_addmul(Q,U,P1);
mpz_clear(P1);
mpz_set(U,U1);
mpz_clear(U1);
}
}else{
mpz_set_ui(T,0);
mpz_set_ui(U,1);
mpz_set_ui(Q,2*end);
if (left==0){
mpz_set_ui(P,2);
}else{
mpz_set_ui(P,left);
mpz_mul_ui(P,P,left);
}
}
}
int main () {
FILE* fp;
unsigned long digits=100000000;
unsigned long prec=(digits+5)*log2(10);
mpf_set_default_prec(prec);
long n=ceil(digits*0.654);
mpf_t temp,pi0,pi;
mpz_t T,U,P,Q;
clock_t start=clock(),end;
mpz_init(T);
mpz_init(P);
mpz_init(Q);
BSA0(T,P,Q,0,n-1,n-1);
mpz_clear(T);
mpf_init(pi0);
mpf_set_z(pi0,P);
mpz_clear(P);
mpf_init(temp);
mpf_set_z(temp,Q);
mpz_clear(Q);
mpf_div(pi0,pi0,temp);
mpf_clear(temp);
end=clock();
printf("%.3f s\n",(double)(end-start)/CLOCKS_PER_SEC);
mpz_init(T);
mpz_init(U);
mpz_init(P);
mpz_init(Q);
BSA(T,U,P,Q,0,n,n);
mpz_clear(T);
mpz_clear(U);
mpf_init(pi);
mpf_set_z(pi,P);
mpz_clear(P);
mpf_init(temp);
mpf_set_z(temp,Q);
mpz_clear(Q);
mpf_div(pi,pi,temp);
if ((n & 1)==1) mpf_neg(pi,pi);
mpf_clear(temp);
mpf_add(pi,pi,pi0);
mpf_clear(pi0);
end=clock();
printf("%.3f s\n",(double)(end-start)/CLOCKS_PER_SEC);
fp=fopen("pi.txt","w");
mpf_out_str(fp,10,digits+1,pi);
fclose(fp);
mpf_clear(pi);
end=clock();
printf("%.3f s\n",(double)(end-start)/CLOCKS_PER_SEC);
return 0;
}
モンゴメ
リ乗算
Montgomery.pdf
モンゴメリ乗算は，時間のかかる剰余算を，乗算，加減算，シフト演算，論理積，条件分岐のみで求めるアルゴリズムです。
RSA暗号やElGamal暗号などの高速化のために使われます。自分の理解のために解説を書きました。