|
|
多倍長 LOG(5)を求める
(旧掲示板、投稿ネタからの掘り起しです)
T=1-X/2^N (0<T<=1) として
LOG(1-T)=-(T+T^2/2+T^3/3+T^4/4+...)
LOG(X)=LOG(T)+N*LOG(2)
又は
T=X/2^(N-1) (1<T<=2) として
LOG((1+U)/(1-U))=2*(U+U^3/3+U^5/5...) U=(T-1)/(T+1)
LOG(X)=LOG(T)+(N-1)*LOG(2)
として求める。
なお、LOG(2)は予め、用意しておく
実行には、2進モードをご使用ください。
PUBLIC NUMERIC BIAS, KETA, SIGN, EPS
INPUT PROMPT "桁数=":KETA !' 2000桁まで
LET KETA = INT(KETA/4)
LET EPS = 5 !'計算誤差分(4*EPS 桁)
LET BIAS = 0 !'10000 ^ (BIAS+1) まで
LET KETA = KETA + EPS !'10000 ^ (-KETA) まで
LET SIGN = -BIAS - 1 !'多倍長数符号 A(SIGN)=1... 正 A(SIGN)=-1...負
DIM A(-BIAS-1 TO KETA)
LET X=5
!'INPUT PROMPT "LOG(X) X=":X
CALL SLOG(X,A)
CALL DISPLAY2(A)
END
EXTERNAL SUB SLOG(X,S())
DIM A(-BIAS-1 TO KETA)
FOR I=1 TO 32
IF X<=2^I THEN EXIT FOR
NEXT I
IF (2^I-X)/2^I<(X-2^(I-1))/2^(I-1) THEN
CALL SLOG2(X,S)
ELSE
CALL SLOG3(X,S)
END IF
END SUB
EXTERNAL SUB SLOG2(AX, X()) !'LOG(1-T) T=1-X/2^N (0<T<=1) → LOG((2^N-X)/2^N)
IF AX < 0 THEN
PRINT "ERROR in SLOG2"
STOP
END IF
DIM U(-BIAS - 1 TO KETA), V(-BIAS - 1 TO KETA), S(-BIAS - 1 TO KETA), N(-BIAS - 1 TO KETA)
LET NN= 1
LET U(0) = 1
LET U(SIGN) = 1
LET S(SIGN) = 1
LET F = 1
DO WHILE 2 ^ F < AX
LET F = F + 1
LOOP
LET XX = 2 ^ F - AX
LET XA = 2 ^ F
CALL LN2(N)
CALL SMUL(N,F)
IF XX=0 THEN
CALL LCOPY(X,N)
EXIT SUB
END IF
DO
CALL SMUL(U,XX)
CALL SDIV(U,XA)
CALL LCOPY(V, U)
CALL SDIV(V, NN)
CALL LADD2(S, V)
LET NN = NN+ 1
LOOP UNTIL ZERO(V)<>0
CALL LSUB (N, S, X)
END SUB
EXTERNAL SUB SLOG3(XA, X()) !'LOG((T-1)/(T+1)) T=X/2^N (1<T<=2) → LOG((X-2^N)/(X+2^N))
IF XA < 0 THEN
PRINT "ERROR in SLOG3"
STOP
END IF
DIM D(-BIAS - 1 TO KETA)
DIM M(-BIAS - 1 TO KETA), S(-BIAS - 1 TO KETA)
DO
LET K=K+1
LOOP WHILE 2 ^ (K+1) < XA
IF XA - 2 ^ (K+1) = 0 THEN
CALL LN2(X)
CALL SMUL(X,K+1)
EXIT SUB
END IF
LET N = 1
CALL LCLR (X)
LET X(0) = XA - 2^K
LET X(SIGN) = 1
CALL SDIV(X,XA+2^K)
CALL LCOPY(S, X)
DO
CALL LCOPY(M, S)
CALL SMUL(X,(XA-2^K)^2)
CALL SDIV(X,(XA+2^K)^2)
CALL LCOPY(D, X)
LET N = N + 2
CALL SDIV(D, N)
CALL LADD2(S, D)
LOOP UNTIL EQUAL(M, S)<>0
CALL SMUL(S,2)
CALL LN2(D)
CALL SMUL(D,K)
CALL LADD(D,S,X)
END SUB
EXTERNAL SUB LN2(X()) !' LOG(2) (2000桁分)
CALL LCLR (X)
LET X(SIGN)=1
FOR I = 1 TO KETA
READ IF MISSING THEN EXIT FOR:X(I)
NEXT I
DATA 6931,4718,0559,9453,0941,7232,1214,5817,6568,0755,0013,4360,2552,5412,0680,0094,9339,3621,9696,9471,5605,8633,2699,6418,6875
DATA 4200,1481,0205,7068,5733,6855,2023,5758,1305,5703,2670,7516,3507,5961,9307,2757,0828,3714,3519,0307,0386,2389,1673,4711,2335
DATA 0115,3644,9795,5239,1204,7517,2681,5749,3206,5155,5247,3413,9525,8829,5045,3007,0953,2636,6642,6541,0423,9157,8149,5204,3740
DATA 4303,8550,0801,9441,7064,1671,5186,4471,2839,9681,7178,4546,9570,2627,1631,0645,4615,0257,2074,0248,1637,7733,8963,8550,6952
DATA 6066,8341,1372,7387,3722,9289,5649,3547,0257,6265,2098,8596,9320,1965,0585,5476,4703,3067,9365,4432,5476,3274,4951,2504,0606
DATA 9438,1471,0468,9946,5062,2016,7720,4245,2452,9612,6879,4654,6193,1651,7468,1392,6725,0410,3802,5462,5965,6869,1441,9287,1608
DATA 2938,0317,2714,3677,8265,4877,5664,8508,5674,0776,4845,1464,4399,4046,1422,6031,9309,6735,4025,7444,6070,3080,9608,5047,4866
DATA 3852,3138,1816,7675,1438,6674,7664,7890,8814,3714,1985,4942,3151,9973,5488,0375,1658,6127,5352,9166,1000,7105,3558,2498,7941
DATA 4729,5092,9311,3897,1559,9820,5654,3928,7170,0072,1808,5761,0252,3688,9213,2449,7138,9320,3784,3935,3088,7748,2597,0171,5591
DATA 0708,8236,8362,7589,8425,8918,5353,0243,6342,1436,7061,1892,3678,9192,3723,1467,2321,7205,3401,6492,5687,2747,7823,4453,5347
DATA 6481,1494,1864,2386,7767,7440,6069,5626,5737,9600,8670,7625,7199,1847,3402,2651,4628,3790,4883,0620,3306,1144,6300,7371,9489
DATA 0027,4364,3965,0025,8093,6519,4430,4119,1150,6080,9487,9306,7865,1588,7090,0605,2034,6842,9736,1938,4128,9652,5565,3968,6022
DATA 1941,2292,4207,5743,2175,7489,0977,0675,2687,1158,1705,1137,0091,5894,2665,4785,9596,4890,6530,5846,0258,6683,8294,0022,8330
DATA 0538,2074,0056,7705,3046,7870,0184,1624,0441,8833,2327,9838,6349,0015,6312,1889,5606,5055,3151,2721,9939,8332,0307,5140,8426
DATA 0914,7900,1265,1682,4344,3893,5724,7278,8205,4862,7155,2741,8772,4300,2489,7945,4019,6187,2339,8086,0831,6648,1149,0930,6675
DATA 1933,9312,8904,3164,1370,6813,9777,6498,1769,7486,8903,8877,8999,1296,5036,1927,0710,8892,6410,5230,9247,8391,7373,5012,2984
DATA 2420,4995,6893,5992,2066,0220,4654,9415,1061,3918,7885,7442,4557,7510,2068,3703,0866,6194,8089,6412,1868,0779,0208,1815,8858
DATA 0001,6881,1597,3056,1866,7619,9187,3952,0076,6719,2145,9223,6720,6025,3959,5436,5416,5531,1295,1759,8994,0056,0003,6651,3567
DATA 5690,5124,5926,8257,4394,6483,1683,3262,4901,8038,2424,0824,2314,5230,6140,9638,0570,0702,5513,8770,2681,7851,6306,9025,5137
DATA 0323,4053,8021,4501,9015,3740,2950,9942,2629,9577,9647,4271,3815,7363,8017,2987,3940,7042,4217,9972,2669,6297,9939,3127,0693
DATA 5747,2404,9338,6530,8797
END SUB
!'以下、多倍長計算共通ルーチン (多倍長LOG(π)を求める、でも使用)
EXTERNAL SUB DISPLAY2(X()) !'表示
FOR K=-BIAS TO 0
IF X(K)<>0 THEN EXIT FOR
NEXT K
IF X(SIGN) = -1 THEN PRINT "- ";
IF K>=0 THEN
LET K=0
PRINT STR$(X(0));"."
ELSE
PRINT STR$(X(K));
FOR I=K+1 TO 0
LET A$=A$ & RIGHT$("000"&STR$(X(I)),4)
IF LEN(A$) = 100 THEN
PRINT A$
LET A$=""
END IF
NEXT I
IF LEN(A$) > 0 THEN
PRINT 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 LEFT$(A$,10);" ";
IF J=5 THEN PRINT " ";
LET A$=RIGHT$(A$,LEN(A$)-10)
NEXT J
PRINT ":";S
LET A$=""
IF MOD(S,1000)=0 THEN PRINT
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 RTRIM$(LEFT$(A$,10));" ";
IF J=5 THEN PRINT " ";
LET A$=RIGHT$(A$,LEN(A$)-10)
IF RTRIM$(A$)="" THEN EXIT FOR
NEXT J
!' PRINT ":";S
END IF
END SUB
EXTERNAL FUNCTION EQUAL(A(), B()) !'等しいかどうか
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
EXTERNAL FUNCTION ZERO(A()) !'0値かどうか
FOR I=-BIAS TO KETA-EPS
IF A(I)<>0 THEN
LET ZERO=0
EXIT FUNCTION
END IF
NEXT I
LET ZERO=-1
END FUNCTION
EXTERNAL FUNCTION GREAT(A(), B()) !'多倍長数A > 多倍長数B なら真
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
END FUNCTION
|
|