多倍長整数演算,円周率


多倍長整数演算ルーチン

公開鍵暗号に使うために作成したもので,それに使う関数しか作っていません。筆算の方法で計算していますので,乗算と除算の実行 速度はO(n2) ですが,2000bit程度ならば十分な速さだと思います。自由に改造してご利用下さい。基数は216で,例えば unsigned long A[1000]; に多倍長整数を入れる場合,A[1]以降の下位16bitにデータを入れ,A[0]にデータ長を入れる構造になっています。
除算のアルゴリズム  multiple.pdf
多倍長整数演算用関数  multi.c

ニュートン法と類似のアルゴリズムによる円周率の計算法

円 周率の計算を,ニュートン法と類似のアルゴリズムで求める方法で,どこまで高速化できるのか考えてみました。
初期値を θ≒π/2としてθ:=θ+cosθを繰り返すと,θは3次収束でπ/2に近づく。BASICのプログラムは次 のようになります。
なお,初期値を θ≒πとしてθ:=θ+sinθを繰り返しても,θは3次収束でπに近づきます。

LET pi2=3/2
FOR i=1 TO 3
   LET pi2=pi2+COS(pi2)
   PRINT 2*pi2
NEXT i
END

証明pi.pdf

次のものは,π/2=x1/B+x2/B3+x3/B9+…+xn/B3^(n-1)と 分解し加法定理を用いて高速に計算する方法で,DRM法という計算法の一部のアイデアを利用したものです。高速フーリエ変換を使うような大きな桁でも,効 果があります。計算量は,π/2の値がわかっているときにcos(π/2)の値を1回だけ求める場合と同じになります。以下に,考え方を示すためのプログ ラムを書きます。

LET pi2=3/2
LET d=pi2
LET s=0
LET c=1
LET n=10
FOR i =1 TO 3
   LET s2=s*COS(d)+c*SIN(d)
   LET c=c*COS(d)-s*SIN(d)
   LET s=s2
   LET n=n^3
   LET d=INT(c*n)/n
   LET pi2=pi2+d
   PRINT 2*pi2
NEXT i
END

最速の円周率の計算法

60桁の円周率の種から1000桁の円周率を求める,1000桁精度の十進BASICのプログラムを書くと,以下のようになります。
p0を60桁精度のπ/2とします。
1 cos(p0/2^50)を計算し,2倍角の公式を用いてcos(p0)を求めます。
2 p0+arcsin(cos(p0))を求めます。これが1000桁精度の円周率と成ります。
証明
p0+arcsin(cos(p0))において,p0=π/2-εとおくと
    cos(p0+arcsin(cos(p0)))=cos(π/2-ε+arcsin(cos(π/2-ε)))=cos(π/2-ε+arcsin (sinε))=cos(π/2-ε+ε)=cos(π/2)=0
p0+arcsin(cos(p0))はπ/2に十分近いから,
    p0+arcsin(cos(p0))=π/2

また,cosθを求めるときにθが小さいと精度落ちが起きるため,途中の計算では,cosθの代わりにcosθ-1を求めています。

三角関数を用いた円周率の計算法

! 三角関数法
OPTION ARITHMETIC DECIMAL_HIGH
DECLARE EXTERNAL SUB SC
DECLARE EXTERNAL FUNCTION arcsin

LET t0=TIME
LET m=10^60
LET p0=INT(PI/2*m)/m  !円周率の種
LET kurikaesi=66
LET p=p0/(2^kurikaesi)
LET n=10^(kurikaesi/3*2)
LET d=INT(p*n)/n
CALL SC(d,s,c)
LET p=p-d
DO
   LET n=n^2
   LET d=INT(p*n)/n
   CALL SC(d,s1,c1)
   LET p=p-d
   IF p=0 THEN
      LET c=c*(c1+1)+c1-s*s1   ! 加法定理
      EXIT DO
   END IF
   LET s2=s*(c1+1)+(c+1)*s1  ! 加法定理
   LET c=c*(c1+1)+c1-s*s1
   LET s=s2    
LOOP
FOR k=1 TO kurikaesi
   LET c=2*c*(c+2)  ! 2倍角の公式
NEXT k
LET pai=(p0+arcsin(c+1))*2
PRINT TIME-t0
PRINT pai-PI
END

EXTERNAL SUB SC(x,s,c)
OPTION ARITHMETIC DECIMAL_HIGH
LET c=0  ! cはcos(x)-1
LET s=x
LET kou=-x^2/2
LET k=3
DO
   LET c=c+kou
   LET kou=kou*x/k
   LET s=s+kou
   LET kou=-kou*x/(k+1)
   LET k=k+2
LOOP WHILE ABS(kou)>=EPS(0)
END SUB

EXTERNAL FUNCTION arcsin(x)
OPTION ARITHMETIC DECIMAL_HIGH
LET y=x
LET xx=x*x
LET kou=x*xx/2
LET k=3
DO
   LET y=y+kou/k
   LET kou=kou*xx*k/(k+1)
   LET k=k+2
LOOP WHILE ABS(kou)>EPS(0)
LET arcsin=y
END function

十進BASICの1000桁モードでは,上記のプログラムは最速と言われているChudnovskyの公式によるプログラムとより高速です。なお,上記の計算法は,十進とは相性が悪く,プログラムを二進用に直すとさらに速く成ります。

他の計算法
比較のためにRamanujan系公式と算術幾何平均法と逆正接公式によるプログラムを掲載します。

次のプログラムは,Chudnovskyの公式によるものです。

!Chudnovskyの公式
OPTION ARITHMETIC DECIMAL_HIGH
LET t0=TIME
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 x=10005 ! SQR(10005)をコンパイラに事前に計算させないため
LET s=426880*SQR(x)/s
PRINT TIME-t0
PRINT s-PI
END

次のプログラムは,Ramanujanの公式によるものです。

!Ramanujanの公式
OPTION ARITHMETIC DECIMAL_HIGH
LET t0=TIME
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 x=8 ! SQR(8)をコンパイラに事前に計算させないため
LET s=9801/(s*SQR(x))
PRINT TIME-t0
PRINT s-PI
END

次のプログラムは,ガウス・ルジャンドルの算術幾何平均法によるもの です。

! ガウス・ルジャンドルの算術幾何平均法
OPTION ARITHMETIC DECIMAL_HIGH
LET t0=TIME
LET a=1
LET b=SQR(1/2)
LET t=1/4
LET k=1
FOR i=0 TO 8
   LET aa=(a+b)/2
   LET b=SQR(a*b)
   LET t=t-k*(a-aa)^2
   LET a=aa
   LET k=k*2
NEXT i
LET pai=(a+b)^2/(t*4)
PRINT TIME-t0
PRINT pai-pi
END

次のプログラムは,逆正接関数の6項最良の公式によるものです。

!6項最良の公式
OPTION ARITHMETIC DECIMAL_HIGH
DECLARE EXTERNAL FUNCTION arccot
LET t0=TIME
LET pai=83*arccot(107)+17*arccot(1710)-44*arccot(225443)-68*arccot(2513489)+22*arccot(42483057)+34*arccot(7939642926390344818)
PRINT TIME-t0
PRINT pai*4-pi
END

EXTERNAL FUNCTION arccot(x)
OPTION ARITHMETIC DECIMAL_HIGH
LET xx=-x*x
LET kou=1/(x*xx)
LET s=1/x
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

加法定理による高速版  addition_theorem.pdf

三角関数を用いた計算法は,加法定理と,下記のBSA法を用い,さらに高速フーリエ変換を用いると,円周率をn桁まで計算するときの計算量はO(n (logn)3) 以下であり,円周率の世界記録に挑戦するときの計算量です。
ただし,大 きな桁を計算するときには,2倍角の公式を使う回数を増やすことはできますが,それでも加法定理を使 う回数が増えるため,Chudnovskyの公式に対して低速であると考えられます。
なお,sinθ,cosθをテーラー展開したときの係数1/n!の速度への寄与は,2倍角の公式を使うことによる速度への寄与に比べて僅かです。

BSA法 (Binary Splitting Algorithm)BSA.pdf
BSA法(Binary Splitting Algorithm)は,級数の部分和や連分数などを高 速に計算する方法です。合成関 数を用いて独自な証明をしました。自分の理解のために解説を書きました。

モンゴメ リ乗算  Montgomery.pdf
モンゴメリ乗算は,時間のかかる剰余算を,乗算,加減算,シフト演算,論理積,条件分岐のみで求めるアルゴリズムです。 RSA暗号やElGamal暗号などの高速化のために使われます。自分の理解のために解説を書きました。

最初のページに戻る