|
|
OPTION ARITHMETIC NATIVE !'二進モード
PUBLIC NUMERIC BIAS, KETA, SIGN, EPS
LET EPS = 16 !'打ち切りによる収束誤差分(末尾 EPS*4桁分)
LET NN=2^15-EPS !' (※ NN=2^18-EPSで100万桁オーバー)
LET BIAS = 2 !'10000^(BIAS+1)まで
LET KETA =-BIAS+NN+EPS !'10000^(-KETA) (KETA*4桁)
LET SIGN = -BIAS - 1 !'多倍長数の符号
DIM Y(-BIAS-1 TO KETA),X(-BIAS-1 TO KETA)
INPUT PROMPT "K^(1/N) K=":K
DO
INPUT PROMPT "K^(1/N) N=":N
LET N=INT(ABS(N))
LOOP UNTIL N>=2
CALL LLET(Y,STR$(K))
CALL LLET(X,USING$("#.################",K^(-(N-1)/N))) !'ニュートン法の初期値
PRINT (KETA-EPS)*4;"桁の計算開始 ";TIME$
LET T=TIME
CALL LROOT(Y,N,X)
PRINT "計算終了 ";TIME$;TIME-T;"秒"
PRINT "ファイルへ書き出し中"
CALL DISPLAY(Y,"OUT.txt") !'ファイルへ書き出し
END
EXTERNAL SUB LROOT(T(),N,X()) !'ニュートン法
OPTION ARITHMETIC NATIVE
DIM TT(-BIAS-1 TO KETA),XX(-BIAS-1 TO KETA),Y(-BIAS-1 TO KETA)
CALL LCOPY(TT,T)
FOR I=1 TO N-2
CALL LMUL2(TT,T)
NEXT I
DO
CALL LCOPY(Y,X)
CALL LCOPY(XX,X)
CALL SMUL(XX,N+1)
FOR I=1 TO N
CALL LMUL2(X,Y)
NEXT I
CALL LMUL2(X,TT)
CALL LSUB2(XX,X)
CALL SDIV(XX,N)
CALL LCOPY(X,XX)
LOOP UNTIL EQUAL(Y,X)<>0
CALL LMUL2(T,X)
END SUB
EXTERNAL SUB DISPLAY(X(),N$)
OPTION ARITHMETIC NATIVE
IF N$="" THEN
OPEN #1: TEXTWINDOW1
ELSE
OPEN #1:NAME N$
END IF
ERASE #1
FOR K=-BIAS TO 0
IF X(K)<>0 THEN EXIT FOR
NEXT K
IF X(SIGN) = -1 THEN PRINT #1:"- ";
IF K>=0 THEN
LET K=0
PRINT #1:STR$(X(0));"."
ELSE
PRINT #1:STR$(X(K));
FOR I=K+1 TO 0
LET A$=A$ & RIGHT$("000" & STR$(X(I)),4)
IF LEN(A$)=100 THEN
PRINT #1:A$
LET A$=""
END IF
NEXT I
IF LEN(A$)>0 THEN
PRINT #1:A$;"."
LET A$=""
END IF
END IF
LET S=0
FOR I=1 TO KETA-EPS
LET A$=A$ & RIGHT$("000" & STR$(X(I)),4)
IF LEN(A$)=100 THEN
LET S=S+100
FOR J=1 TO 10
PRINT #1:LEFT$(A$,10);" ";
IF J=5 THEN PRINT #1:" ";
LET A$=RIGHT$(A$,LEN(A$)-10)
NEXT J
PRINT #1:":";S
LET A$=""
IF MOD(S,1000)=0 THEN PRINT #1
END IF
NEXT I
IF LEN(A$)>0 THEN
LET S=S+LEN(A$)
LET A$=A$ & REPEAT$(" ",10)
FOR J=1 TO 9
PRINT #1:RTRIM$(LEFT$(A$,10));" ";
IF J=5 THEN PRINT #1:" ";
LET A$=RIGHT$(A$,LEN(A$)-10)
IF RTRIM$(A$)="" THEN EXIT FOR
NEXT J
!' PRINT #1:":";S
END IF
CLOSE #1
END SUB
EXTERNAL FUNCTION EQUAL(A(), B())
OPTION ARITHMETIC NATIVE
FOR I=-BIAS-1 TO KETA-EPS
IF A(I)<>B(I) THEN
FOR J=0 TO 3
IF INT(A(I)/10^(3-J))<>INT(B(I)/10^(3-J)) THEN EXIT FOR
NEXT J
PRINT (I-1)*4+J;"桁まで計算終了 ";TIME$
LET EQUAL = 0
EXIT FUNCTION
END IF
NEXT I
PRINT (I-1)*4;"桁まで計算終了 ";TIME$
LET EQUAL=-1
END FUNCTION
EXTERNAL FUNCTION GREAT(A(), B())
OPTION ARITHMETIC NATIVE
LET SIGNA = A(SIGN)
LET SIGNB = B(SIGN)
IF SIGNA = -1 AND SIGNB = 1 THEN
LET GREAT = 0
EXIT FUNCTION
END IF
IF SIGNA = 1 AND SIGNB = -1 THEN
LET GREAT = -1
EXIT FUNCTION
END IF
FOR I = -BIAS TO KETA
IF SIGNA = -1 AND SIGNB = -1 THEN
IF A(I) < B(I) THEN
LET GREAT = -1
EXIT FUNCTION
END IF
IF A(I) > B(I) THEN
LET GREAT = 0
EXIT FUNCTION
END IF
ELSE
IF A(I) > B(I) THEN
LET GREAT = -1
EXIT FUNCTION
END IF
IF A(I) < B(I) THEN
LET GREAT = 0
EXIT FUNCTION
END IF
END IF
NEXT I
LET GREAT = 0
END FUNCTION
EXTERNAL SUB LCLR(A())
OPTION ARITHMETIC NATIVE
MAT A=ZER
LET A(SIGN) = 1
END SUB
EXTERNAL SUB LCOPY(A(), B())
OPTION ARITHMETIC NATIVE
MAT A=B
END SUB
EXTERNAL SUB LLET(A(), X$)
OPTION ARITHMETIC NATIVE
CALL LCLR(A)
LET X$ = LTRIM$(RTRIM$(X$))
FOR I = 1 TO LEN(X$)
IF POS("0123456789.-", MID$(X$, I, 1)) = 0 THEN
PRINT "ERROR in LLET"
STOP
END IF
NEXT I
IF LEFT$(X$, 1) = "-" THEN
LET A(SIGN) = -1
LET X$ = RIGHT$(X$, LEN(X$)-1)
ELSE
LET A(SIGN) = 1
END IF
LET K1 = POS(X$, ".")
LET X$ = "000" & X$
IF K1 = 0 THEN LET X$ = X$ & "."
LET KK = POS(X$, ".")
IF KK - 4 > BIAS * 4 + 4 THEN
PRINT "OVER FLOW in LLET"
STOP
END IF
FOR J = 0 TO INT((KK - 1) / 4) - 1
LET A(-J) = VAL(MID$(X$, KK - 4 * J - 4, 4))
NEXT J
LET X$ = X$ & "000"
FOR J = 1 TO(LEN(X$) - KK) / 4
LET A(J) = VAL(MID$(X$, KK + 4 * J - 3, 4))
NEXT J
END SUB
EXTERNAL SUB SDIV(A(), XA)
OPTION ARITHMETIC NATIVE
LET SIGNA = A(SIGN)
LET SG = SGN(XA)
LET XA = ABS(XA)
FOR I = -BIAS TO KETA - 1
LET R = A(I) - INT(A(I) / XA) * XA
LET A(I) = INT(A(I) / XA)
LET A(I + 1) = A(I + 1) + R * 10000
NEXT I
LET A(KETA) = INT(A(KETA) / XA)
LET A(SIGN)=SIGNA * SG
END SUB
EXTERNAL SUB SMUL(A(), XA)
OPTION ARITHMETIC NATIVE
LET SIGNA = A(SIGN)
LET SG = SGN(XA)
LET XA = ABS(XA)
MAT A = XA * A
FOR I = KETA TO -BIAS + 1 STEP -1
IF A(I) >= 10000 THEN
LET R = INT(A(I) / 10000)
LET A(I) = MOD(A(I),10000)
LET A(I - 1) = A(I - 1) + R
END IF
NEXT I
IF A(-BIAS) >= 10000 THEN
PRINT "OVER FLOW in SMUL"
STOP
END IF
LET A(SIGN)=SIGNA * SG
END SUB
EXTERNAL SUB LADD2(A(), B())
OPTION ARITHMETIC NATIVE
DIM C(-BIAS - 1 TO KETA)
CALL LADD(A,B,C)
CALL LCOPY(A,C)
END SUB
EXTERNAL SUB LSUB2(A(), B())
OPTION ARITHMETIC NATIVE
DIM C(-BIAS - 1 TO KETA)
CALL LSUB(A,B,C)
CALL LCOPY(A,C)
END SUB
EXTERNAL SUB LADD(A(), B(), C())
OPTION ARITHMETIC NATIVE
LET SIGNA = A(SIGN)
LET SIGNB = B(SIGN)
IF SIGNA = 1 AND SIGNB = -1 THEN
LET B(SIGN) = 1
CALL LSUB(A, B, C)
LET B(SIGN) = -1
EXIT SUB
ELSEIF SIGNA = -1 AND SIGNB = 1 THEN
LET A(SIGN) = 1
CALL LSUB(B, A, C)
LET A(SIGN) = -1
EXIT SUB
END IF
MAT C=A+B
FOR I = KETA TO -BIAS + 1 STEP -1
IF C(I) >= 10000 THEN
LET C(I) = C(I) - 10000
LET C(I - 1) = C(I - 1) + 1
END IF
NEXT I
IF C(-BIAS) >= 10000 THEN
PRINT "OVER FLOW in LADD"
STOP
END IF
IF SIGNA = -1 AND SIGNB = -1 THEN LET C(SIGN) = -1 ELSE LET C(SIGN) = 1
END SUB
EXTERNAL SUB LSUB(A(), B(), C())
OPTION ARITHMETIC NATIVE
LET SIGNA = A(SIGN)
LET SIGNB = B(SIGN)
LET A(SIGN) = 1
LET B(SIGN) = 1
IF SIGNA * SIGNB = -1 THEN
CALL LADD(A, B, C)
LET C(SIGN) = SIGNA
LET A(SIGN) = SIGNA
LET B(SIGN) = SIGNB
EXIT SUB
END IF
LET GR = GREAT(A, B)
IF SIGNA = 1 AND SIGNB = 1 THEN
IF GR<>0 THEN
MAT C=A-B
LET C(SIGN) = 1
ELSE
MAT C=B-A
LET C(SIGN) = -1
END IF
ELSE
IF GR<>0 THEN
MAT C=B-A
LET C(SIGN) = 1
ELSE
MAT C=A-B
LET C(SIGN) = -1
END IF
END IF
FOR I = KETA TO -BIAS + 1 STEP -1
IF C(I) < 0 THEN
LET C(I) = C(I) + 10000
LET C(I - 1) = C(I - 1) - 1
END IF
NEXT I
LET A(SIGN) = SIGNA
LET B(SIGN) = SIGNB
END SUB
EXTERNAL SUB LMUL(A(),B(),C())
OPTION ARITHMETIC NATIVE
LET N=(KETA+BIAS)*2
IF INT(LOG2(N))<>LOG2(N) THEN
PRINT "ERROR in LMUL"
STOP
END IF
OPTION BASE 0
DIM AA(2*N),BB(2*N),CC(2*N),IP(2+2^(INT(LOG(N+0.5)/LOG(2))/2)),W(N/2-1)
FOR I = 0 TO N/2-1
LET AA(2*I) = A(-BIAS+I)
LET BB(2*I) = B(-BIAS+I)
NEXT I
CALL CDFT(2*N, 1, AA, IP, W)
CALL CDFT(2*N, 1, BB, IP, W)
FOR I = 0 TO N-1
LET CC(2*I) = AA(2*I) * BB(2*I) - AA(2*I+1) * BB(2*I+1)
LET CC(2*I+1) = AA(2*I) * BB(2*I+1) + BB(2*I) * AA(2*I+1)
NEXT I
CALL CDFT(2*N, -1,CC, IP, W)
FOR I=0 TO N/2-1
IF -2*BIAS+I>=-BIAS AND -2*BIAS+I<=KETA THEN
LET C(-2*BIAS+I)=INT(CC(2*I)/N+.5)
END IF
NEXT I
FOR I=KETA TO -BIAS STEP -1
IF C(I) >=10000 THEN
LET R = INT(C(I) / 10000)
LET C(I) = C(I) - R * 10000
LET C(I - 1) = C(I - 1) + R
ELSE
DO WHILE C(I)<0
LET C(I)=C(I)+10000
LET C(I-1)=C(I-1)-1
LOOP
END IF
NEXT I
LET C(SIGN)=A(SIGN)*B(SIGN)
END SUB
EXTERNAL SUB LMUL2(A(),B())
OPTION ARITHMETIC NATIVE
DIM C(-BIAS - 1 TO KETA)
CALL LMUL(A,B,C)
CALL LCOPY(A,C)
END SUB
! これより以下のFFTルーチンは下記サイトより入手したものを移植した
!
! -------- Complex DFT (Discrete Fourier Transform) --------
! [definition]
! <case1>
! X(k) = sum_j=0^n-1 x(j)*exp(2*pi*i*j*k/n), 0<=k<n
! <case2>
! X(k) = sum_j=0^n-1 x(j)*exp(-2*pi*i*j*k/n), 0<=k<n
! (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
! [usage]
! <case1>
! ip(0) = 0 ! first time only
! call cdft(2*n, 1, a, ip, w)
! <case2>
! ip(0) = 0 ! first time only
! call cdft(2*n, -1, a, ip, w)
! [parameters]
! 2*n :data length (integer)
! n >= 1, n = power of 2
! a(0:2*n-1) :input/output data (real*8)
! input data
! a(2*j) = Re(x(j)),
! a(2*j+1) = Im(x(j)), 0<=j<n
! output data
! a(2*k) = Re(X(k)),
! a(2*k+1) = Im(X(k)), 0<=k<n
! ip(0:*) :work area for bit reversal (integer)
! length of ip >= 2+sqrt(n)
! strictly,
! length of ip >=
! 2+2**(int(log(n+0.5)/log(2.0))/2).
! ip(0),ip(1) are pointers of the cos/sin table.
! w(0:n/2-1) :cos/sin table (real*8)
! w(),ip() are initialized if ip(0) = 0.
! [remark]
! Inverse of
! call cdft(2*n, -1, a, ip, w)
! is
! call cdft(2*n, 1, a, ip, w)
! do j = 0, 2 * n - 1
! a(j) = a(j) / n
! end do
EXTERNAL SUB CDFT(N, ISGN, A(), IP(), W())
OPTION ARITHMETIC NATIVE
IF N > 4 * IP(0) THEN
CALL MAKEWT(N / 4, IP, W)
END IF
IF N > 4 THEN
IF ISGN >= 0 THEN
CALL BITRV2(N, IP, A)
CALL CFTFSUB(N, A, W)
ELSE
CALL BITRV2CONJ(N, IP, A)
CALL CFTBSUB(N, A, W)
END IF
ELSEIF N = 4 THEN
CALL CFTFSUB(N, A, W)
END IF
END SUB
EXTERNAL SUB MAKEWT(NW, IP(), W())
OPTION ARITHMETIC NATIVE
LET IP(0) = NW
LET IP(1) = 1
IF NW > 2 THEN
LET NWH = NW / 2
LET DELTA = ATN(1.0) / NWH
LET W(0) = 1
LET W(1) = 0
LET W(NWH) = COS(DELTA * NWH)
LET W(NWH + 1) = W(NWH)
IF NWH > 2 THEN
FOR J = 2 TO NWH - 2 STEP 2
LET X = COS(DELTA * J)
LET Y = SIN(DELTA * J)
LET W(J) = X
LET W(J + 1) = Y
LET W(NW - J) = Y
LET W(NW - J + 1) = X
NEXT J
CALL BITRV2(NW, IP, W)
END IF
END IF
END SUB
EXTERNAL SUB CFTFSUB(N, A(), W())
OPTION ARITHMETIC NATIVE
LET L = 2
IF N > 8 THEN
CALL CFT1ST(N, A, W)
LET L = 8
DO WHILE 4 * L < N
CALL CFTMDL(N, L, A, W)
LET L = 4 * L
LOOP
END IF
IF 4 * L = N THEN
FOR J = 0 TO L - 2 STEP 2
LET J1 = J + L
LET J2 = J1 + L
LET J3 = J2 + L
LET X0R = A(J) + A(J1)
LET X0I = A(J + 1) + A(J1 + 1)
LET X1R = A(J) - A(J1)
LET X1I = A(J + 1) - A(J1 + 1)
LET X2R = A(J2) + A(J3)
LET X2I = A(J2 + 1) + A(J3 + 1)
LET X3R = A(J2) - A(J3)
LET X3I = A(J2 + 1) - A(J3 + 1)
LET A(J) = X0R + X2R
LET A(J + 1) = X0I + X2I
LET A(J2) = X0R - X2R
LET A(J2 + 1) = X0I - X2I
LET A(J1) = X1R - X3I
LET A(J1 + 1) = X1I + X3R
LET A(J3) = X1R + X3I
LET A(J3 + 1) = X1I - X3R
NEXT J
ELSE
FOR J = 0 TO L - 2 STEP 2
LET J1 = J + L
LET X0R = A(J) - A(J1)
LET X0I = A(J + 1) - A(J1 + 1)
LET A(J) = A(J) + A(J1)
LET A(J + 1) = A(J + 1) + A(J1 + 1)
LET A(J1) = X0R
LET A(J1 + 1) = X0I
NEXT J
END IF
END SUB
|
|