|
|
OPTION ARITHMETIC NATIVE
PUBLIC NUMERIC BIAS, KETA, SIGN, EPS
LET EPS = 15
LET NN=2^13-EPS
LET BIAS = 0
LET KETA =-BIAS+NN+EPS
LET SIGN = -BIAS - 1
DIM A(-BIAS - 1 TO KETA), B(-BIAS - 1 TO KETA), T(-BIAS - 1 TO KETA), D(-BIAS - 1 TO KETA)
DIM Y(-BIAS - 1 TO KETA)
PRINT (KETA-EPS)*4;"桁の計算"
PRINT "初期値計算開始 ";TIME$
LET TT=TIME
LET B(1)=5000
LET B(SIGN)=1
CALL LLET(A,"1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623091229702492483605585073721264412149709993583141322266592750559275579995050115278206057147010955997160597027453459686201472851741864088919860955232923048430871432145083976260362799525140798968725339654633180882964062061525835239505474575028775996172983557522033753185701135437460340849884716038689997069900481503054402779031645424782306849293691862158057846311159666871301301561856898723723528850926486124949771542183342042856860601468247207714358548741556570696776537202264854470158588016207584749226572260020855844665214583988939443709265918003113882464681570826301005948587040031864803421948972782906410450726368813137398552561173220402450912277002269411275736272804957381089675040183698683684507257993647290607629969413804756548237289971803268024744206292691248590521810044598421505911202494413417285314781058036033710773091828693147101711116839165817268894197587165821521282295184884720896946")
CALL LROOT(B,2,A)
CALL LCLR(A)
LET A(0)=1
LET A(SIGN)=1
LET T(0)=1
LET T(SIGN)=1
PRINT "本計算開始 ";TIME$
LET N=INT(LOG2((KETA-EPS)*4))
PRINT N;"回のループ"
FOR I=1 TO N
PRINT I; "回目ループ開始 ";TIME$
CALL LCOPY(Y, A) !'Y=A
CALL LADD2(A, B) !'A=A+B
CALL SDIV(A, 2) !'A=(A+B)/2
CALL LMUL2(B, Y) !'B=B*Y
CALL LLET(D,"1.1803405990160962260453379405584885872337166348814472995158643994043041807207157949784586161958079542094501011740292391558983114363307772537179901431278907198180022631848551446228568968727290914938763567035131526302916966294160360929620596994708114931694510735455482288368551127932729308542681350952930150721252050562501929985875965308823634425612598226755198069995211240731261603162039175693541218683276929413095721235428538964841422109753898357049667178621838238977087340954876250984589313456259797012783650441595738187089403712449669273254118014112496048274590183868547091977755666227180131335726631660131266131004593333344606488043269036157940955769174454232077882420546288372747664608765565451867337191134341713055280848497823085519339914098700300563588272756017194524101464447383612115436206799071109638039056449850030469593376362703990715387837080990494695017569187571119930919000206494006335605678681262111333568785739621375941779126327769125705948757162148053434214663157870501672695252082809896437")
CALL LROOT(B,2,D) !'B=SQR(B)
CALL LSUB(Y, A, D) !'D=Y-A
CALL LMUL2(D,D) !'D=(Y-A)^2
CALL SMUL(D,2^(I+1)) !'D=(Y-A)^2*2^(I+1)
CALL LSUB2(T, D) !'T=T-D
NEXT I
PRINT "本計算終了 ";TIME$
CALL LADD2(A, B) !'A=A+B
CALL LMUL2(A,A) !'A=(A+B)^2
CALL LLET(D,"1.09421980761323831941838497035223227162968636141278336105964310528972596922296647388551657483878043159010324264134932996553974460016759121846133050698996332288555287910236240242455177433717867102689074846117652736155807179216436398690608159709925338697120024700336549551674314126512619334904422167530426298919613374763698493815324524156559140599889382645126987244457105603196228753398324115087573478288902030415549530952126504651807549640134536466732889257628998153939987431623998315595746787310366295077371569110455025751324445372136919929577952495837146837832273056504149196082794537931421338179675190631479791945279819574381912180875059411423540455151542190868866660269903546477431500319054432488573938743984765851668467635521405993224813414120248734750833760132048057438717569869091694761936201474148039365220149010688572137927945794255839641679283999821834698109484233041245570229569281484483555944122238475703510331242049350614542146739175784519109295251742429195053833668410032388383123379673265018793")
PRINT "割り算計算開始 ";TIME$
CALL REC(T,D) !'逆数
CALL LMUL2(A,T)
PRINT "計算終了 ";TIME$;TIME-TT;"秒"
CALL DISPLAY(A,"PI.TXT")
END
EXTERNAL SUB REC(Y(),X()) !'逆数
OPTION ARITHMETIC NATIVE
DIM YY(-BIAS - 1 TO KETA)
DIM XX(-BIAS - 1 TO KETA),A(-BIAS - 1 TO KETA)
LET A(0)=2
LET A(SIGN)=1
CALL LCOPY(YY,Y)
DO
CALL LCOPY(XX,X)
CALL LCOPY(Y,YY)
LET Y(SIGN)=-1
CALL LMUL2(Y,X) !'Y=Y*X
CALL LADD2(Y,A) !'Y=Y+2
CALL LMUL2(Y,X) !'Y=Y*X
CALL LCOPY(X,Y)
LOOP UNTIL EQUAL(X,XX)=-1
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
IF A$<>"" THEN PRINT #1:A$ !' <--追加訂正
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
LET EQUAL=0
EXIT FUNCTION
END IF
NEXT I
LET EQUAL=-1
END FUNCTION
以下省略
※「n乗根の計算」から必要なサブルーチンをコピペしてください
|
|