複素数の計算

 投稿者:山中和義  投稿日:2009年11月20日(金)09時17分31秒
  複素数モードをサポートする十進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 !第1象限、第4象限
   LET r=ATN2(y/x)
ELSE
   IF x<0 THEN
      IF y>=0 THEN !第2象限
         LET r=ATN2(y/x)+PI
      ELSE !第3象限
         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 "0では割れません。"
   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" !三角関数(正弦,余弦,正接)
 

Re: 複素数の計算

 投稿者:山中和義  投稿日:2009年11月20日(金)09時21分31秒
  > No.747[元記事へ]

続き  ※必要に応じて追加してください。
!演算関連 べき乗、関数など

EXTERNAL SUB CompPowN(x1,y1, n, x,y) !べき乗 (x1+y1*i) ^ n → x+y*i
DECLARE EXTERNAL FUNCTION SIN,COS !関数のオーバーロード
LET r=CompABS(x1,y1) !r,θ
LET th=CompARG(x1,y1)
LET t=r^n !z^n=r^n*(COS(n*θ)+i*SIN(n*θ))、r=SQR(x1*x1+y1*y1)
LET x=t*COS(n*th)
LET y=t*SIN(n*th)
END SUB

EXTERNAL SUB CompNPow(n, x1,y1, x,y) !べき乗 n ^ (x1+y1*i) → x+y*i
DECLARE EXTERNAL FUNCTION LOG !関数のオーバーロード
CALL CompMul(LOG(n),0, x1,y1, a,b) !n^z1=exp(log(n)*z1)より
CALL CompEXP(a,b, x,y)
END SUB


EXTERNAL SUB CompSIN(x1,y1, x,y) !三角関数 正弦
DECLARE EXTERNAL FUNCTION EXP,SIN,COS !関数のオーバーロード
!sin(x+yi)
!={exp(x*i-y)-exp(-x*i+y)}/(2*i)  ※sin(z)=(exp(i*z)-exp(-i*z))/(2*i)より
!={exp(-y)*(cos(x)+i*sin(x))-exp(y)*(cos(x)-i*sin(x))}/(2*i) ※exp(i*x)=cos(x)+i*sin(x)より
!={(exp(y)-exp(-y))*(-cos(x)) + (exp(y)+exp(-y))*sin(x)*i}/(2*i)
!={(exp(y)-exp(-y))*cos(x)*i + (exp(y)+exp(-y))*sin(x)}/2

LET e=EXP(y1)
LET f=1/e
LET yy=0.5*(e-f)*COS(x1)
LET xx=0.5*(e+f)*SIN(x1)
LET x=xx
LET y=yy
END SUB

EXTERNAL SUB CompCOS(x1,y1, x,y) !三角関数 余弦
DECLARE EXTERNAL FUNCTION EXP,SIN,COS !関数のオーバーロード
LET e=EXP(y1) !(exp(i*z)+exp(-i*z))/2
LET f=1/e
LET yy=0.5*(f-e)*SIN(x1)
LET xx=0.5*(f+e)*COS(x1)
LET x=xx
LET y=yy
END SUB

EXTERNAL SUB CompTAN(x1,y1, x,y) !三角関数 正接
DECLARE EXTERNAL FUNCTION EXP,SIN,COS !関数のオーバーロード
LET e=EXP(2*y1)
LET f=1/e
LET d=COS(2*x1)+0.5*(e+f)
LET x=SIN(2*x1)/d
LET y=(e-f)/d
END SUB

!逆三角関数
! ArcSin(z)=-i*LOG(SQR(1-z^2)+z*i)
! ArcCos(z)=-i*LOG(z+i*SQR(1-z^2))
! ArcTan(z)=i/2*LOG((i+z)/(i-z))


EXTERNAL SUB CompSinh(x1,y1, x,y) !双曲線正弦
DECLARE EXTERNAL FUNCTION EXP,SIN,COS !関数のオーバーロード
LET e=EXP(x1)
LET f=1/e
LET xx=0.5*(e-f)*COS(y1)
LET yy=0.5*(e+f)*SIN(y1)
LET x=xx
LET y=yy
END SUB

EXTERNAL SUB CompCosh(x1,y1, x,y) !双曲線余弦
DECLARE EXTERNAL FUNCTION EXP,SIN,COS !関数のオーバーロード
LET e=EXP(x1)
LET f=1/e
LET xx=0.5*(e+f)*COS(y1)
LET yy=0.5*(e-f)*SIN(y1)
LET x=xx
LET y=yy
END SUB

EXTERNAL SUB ComppTanh(x1,y1, x,y) !双曲線正接
DECLARE EXTERNAL FUNCTION EXP,SIN,COS !関数のオーバーロード
LET e=EXP(2*x1)
LET f=1/e
LET d=0.5*(e+f)+COS(2*y1)
LET y=SIN(2*y1)/d
LET x=0.5*(e-f)/d
END SUB
 

戻る