複素数モードをサポートする十進BASICでは、あまり必要性はないと思いますが、
アルゴリズムの勉強として参考にしてください。
サンプルは1000桁モードですが、収束判定の精度を落とすことで、
10進や2進モードでも実行できます。その場合は、複素数モードと同じことになりますが、、、
サンプル リーマンのゼータ関数の零点
!複素数の計算
LET t0=TIME
!リーマンのゼータ関数
! ζ(s)=1/(1-2^(1-s))*∑[m=0,∞]{2^(-(m+1))*∑[j=0,m]{(-1)^j*comb(m,j)*(j+1)^(-s)}}
LET M=100 !精度20桁程度 ※精度1000桁程度は、3000
DIM C(0 TO M) !∑2^(-(M+1))*∑comb(M,J)
LET C(0)=1
FOR K=1 TO M !パスカルの三角形のM段より、二項係数comb(m,j)を求める
FOR J=K TO 1 STEP -1
LET C(J)=(C(J)+C(J-1))/2 !※2^(-(m+1))も加味する
NEXT J
NEXT K
PRINT C(M) !debug
SUB fnZeta(ReS,ImS, ReZ,ImZ) !リーマンのゼータ関数
local J
local ReU,ImU,ReT1,ImT1,ReT2,ImT2
LET ReU=0 !LET U=0
LET ImU=0
FOR J=M TO 0 STEP -1
CALL CompPow(J+1,0, ReS,ImS, ReT1,ImT1) !LET U=U+(-1)^J*C(J)/(J+1)^S
CALL CompDiv((-1)^J*C(J),0, ReT1,ImT1, ReT2,ImT2)
CALL CompAdd(ReU,ImU, ReT2,ImT2, ReU,ImU)
NEXT J
CALL CompPow(2,0, 1-ReS,-ImS, ReT1,ImT1) !LET fnZeta=U/(1-2^(1-S))
CALL CompDiv(ReU,ImU, 1-ReT1,ImT1, ReZ,ImZ)
END SUB
!ゼータ関数の零点の計算（Riemann zeta zeros）
LET ReS=1/2 !非自明な零点 1/2+14*i
LET ImS=14 !1/2+t*i t=14, 21, 25, 30, 33, 38, 41, 43, 48, 50, … 付近
!ニュートン法で零点を求める
! sn+1 = sn-ζ(sn)/ζ'(sn)、ζ'(sn)≒(ζ(sn+h)-ζ(sn))/h hは微小な値より
! sn+1 = sn-ζ(sn)*h/(ζ(sn+h)-ζ(sn)) = sn+h/(1-ζ(sn+h)/ζ(sn)) となる。
! h=ζ(sn)/sn として、収束の判定にはζ(s)の絶対値を使用する。
DO
CALL fnZeta(ReS,ImS, ReZ,ImZ) !LET Z=fnZeta(S)
CALL CompDiv(ReZ,ImZ,ReS,ImS, ReH,ImH) !LET H=Z/S
CALL fnZeta(ReS+ReH,ImS+ImH, ReW,ImW) !LET W=fnZeta(S+H)
CALL CompDiv(ReW,ImW, ReZ,ImZ, ReT1,ImT1) !LET S=S+H/(1-W/Z)
CALL CompDiv(ReH,ImH, 1-ReT1,-ImT1, ReT2,ImT2)
CALL CompAdd(ReS,ImS, ReT2,ImT2, ReS,ImS)
PRINT USING "####.#################### ####.####################": ReS,ImS !PRINT S
LOOP UNTIL CompABS(ReZ,ImZ)<1/10^18 !LOOP UNTIL ABS(Z)<1/10^18
!収束するまで ※精度は調整が必要である。10進モード、2進モードでは1/10^8程度
PRINT "計算時間=";TIME-t0
END
!複素数の計算
EXTERNAL FUNCTION CompRe(x,y) !実部 (x+y*i)、x,yは実数、iは虚数単位
LET CompRe=x
END FUNCTION
EXTERNAL FUNCTION CompIm(x,y) !虚部 Im(x+y*i)
LET CompIm=y
END FUNCTION
EXTERNAL SUB CompConj(x1,y1, x,y) !共役複素数
LET xx=x1
LET yy=-y1
LET x=xx
LET y=yy
END SUB
EXTERNAL FUNCTION CompABS(x,y) !絶対値 |x+y*i|
!a'=MAX(ABS(a),ABS(b))、b'=MIN(ABS(a),ABS(b))として、SQR(a*a+b*b)=a'*SQR(1+(b'/a')^2)を求める。
IF x=0 THEN
LET r=ABS(y)
ELSEIF y=0 THEN
LET r=ABS(x)
ELSEIF ABS(y)>ABS(x) THEN
LET t=x/y
LET r=ABS(y)*SQR(1+t*t)
ELSE
LET t=y/x
LET r=ABS(x)*SQR(1+t*t)
END IF
LET CompABS=r
END FUNCTION
EXTERNAL FUNCTION CompARG(x,y) !偏角θ (-π,π]
IF x>0 THEN !第１象限、第４象限
LET r=ATN2(y/x)
ELSE
IF x<0 THEN
IF y>=0 THEN !第２象限
LET r=ATN2(y/x)+PI
ELSE !第３象限
LET r=ATN2(y/x)-PI
END IF
ELSE !x=0
IF y>0 THEN !∞
LET r=PI/2
ELSE
IF y<0 THEN !-∞
LET r=-PI/2
ELSE !y=0
PRINT "ARG関数は不定です。"
STOP
END IF
END IF
END IF
END IF
LET CompARG=r
END FUNCTION
!演算関連 べき乗
EXTERNAL SUB CompPow(x1,y1, x2,y2, x,y) !べき乗 (x1+y1*i) ^ (x2+y2*i) → x+y*i
CALL CompLOG(x1,y1, s,t) !z1^z2=exp(log(z1)*z2)より
CALL CompMul(s,t, x2,y2, a,b)
CALL CompEXP(a,b, x,y)
END SUB
!演算関連 四則演算、関数など
EXTERNAL SUB CompAdd(x1,y1, x2,y2, x,y) !加算 (x1+y1*i) + (x2+y2*i) → x+y*i
LET xx=x1+x2 !実部
LET yy=y1+y2 !虚部
LET x=xx
LET y=yy
END SUB
EXTERNAL SUB CompSub(x1,y1, x2,y2, x,y) !減算 (x1+y1*i) - (x2+y2*i) → x+y*i
LET xx=x1-x2
LET yy=y1-y2
LET x=xx
LET y=yy
END SUB
EXTERNAL SUB CompMul(x1,y1, x2,y2, x,y) !乗算 (x1+y1*i) * (x2+y2*i) → x+y*i
LET xx=x1*x2-y1*y2 !丸め誤差対策
LET yy=x1*y2+y1*x2
LET x=xx
LET y=yy
END SUB
EXTERNAL SUB CompDiv(x1,y1, x2,y2, x,y) !除算 (x1+y1*i) / (x2+y2*i) → x+y*i
IF x2=0 AND y2=0 THEN
PRINT "０では割れません。"
STOP
END IF
IF ABS(x2)>=ABS(y2) THEN !上位桁あふれ対策
LET w=y2/x2
LET tt=x2+y2*w
LET xx=(x1+y1*w)/tt
LET yy=(y1-x1*w)/tt
ELSE
LET w=x2/y2
LET tt=x2*w+y2
LET xx=(x1*w+y1)/tt
LET yy=(y1*w-x1)/tt
END IF
LET x=xx
LET y=yy
END SUB
EXTERNAL SUB CompSQR(x1,y1, x,y) !平方根
LET SQRT05=SQR(1/2) !0.707106781186547524
LET r=CompABS(x1,y1) !r=SQR(x1*x1+y1*y1)
LET w=SQR(r+ABS(x1))
IF x1>=0 THEN
LET x=SQRT05*w
LET y=SQRT05*y1/w
ELSE
LET x=SQRT05*ABS(y1)/w
IF y1>=0 THEN LET y=SQRT05*w ELSE LET y=-SQRT05*w
END IF
END SUB
EXTERNAL SUB CompEXP(x1,y1, x,y) !指数関数
DECLARE EXTERNAL FUNCTION EXP,SIN,COS !関数のオーバーロード
LET t=EXP(x1) !EXP(x+i*y)=EXP(x)*EXP(i*y)=EXP(x)*(COS(y)+i*SIN(y))
LET xx=t*COS(y1)
LET yy=t*SIN(y1)
LET x=xx
LET y=yy
END SUB
EXTERNAL SUB CompLOG(x1,y1, x,y) !対数関数
DECLARE EXTERNAL FUNCTION LOG !関数のオーバーロード
LET xx=0.5*LOG(x1*x1+y1*y1) !LOG(r)+i*θ、-π<θ<=π
LET yy=CompARG(x1,y1)
LET x=xx
LET y=yy
END SUB
!補助ルーチン ※多桁の関数
EXTERNAL FUNCTION ATN2(x) !アークタンジェント (-π/2,π/2)
IF x>1 THEN
LET cSGN=1
LET x=1/x
ELSEIF x<-1 THEN
LET cSGN=-1
LET x=1/x
ELSE
LET cSGN=0
END IF
LET a=0
FOR i=1500 TO 1 STEP -1 !※繰り返し回数は調整が必要である
LET a0=a
LET a=(i*i*x*x)/(2*i+1 + a)
IF ABS(a-a0)<=EPS(0) THEN EXIT FOR
NEXT i
IF cSGN>0 THEN
LET ATN2=PI/2-x/(1+a)
ELSEIF cSGN<0 THEN
LET ATN2=-PI/2-x/(1+a)
ELSE
LET ATN2=x/(1+a)
END IF
END FUNCTION
MERGE "EXP.LIB" !指数関数 ※級数展開で関数を計算する
MERGE "LOG.LIB" !対数関数
MERGE "TRIGONOM.LIB" !三角関数（正弦，余弦，正接）