多倍長LOG(π)を求める

 投稿者：しばっち  投稿日：2010年 1月23日(土)19時22分46秒
  多倍長 LOG(π)を求める (2000桁)

計算式に

    LOG(1 - X) = -(X + X^2/2 + X^3/3 + X^4/4 +...)

を使用した。

X は試行錯誤した結果

     X = 1 - 544482330679994391053312457583 / 1710541690073718870111737129379 * π

とした。計算後に

  LOG(17105416...) - LOG(54448233...)を加算して、LOG(π)を求める。

また、計算時間短縮のため、π、LOG(17105416...)、LOG(54448233...)値は予め用意した。

　共通ルーチンが必要。

PUBLIC NUMERIC BIAS, KETA, SIGN, EPS
LET  EPS=12
LET  N=2^9-EPS
LET  BIAS = 0
LET  KETA =-BIAS+N+EPS
LET  SIGN = -BIAS - 1
DIM X(-BIAS - 1 TO KETA)
CALL LLOGPI(X)
CALL DISPLAY2(X)
END

EXTERNAL  SUB LLOGPI(X())
DIM L(-BIAS - 1 TO KETA)
DIM A(-BIAS - 1 TO KETA), B(-BIAS - 1 TO KETA)
DIM S(-BIAS - 1 TO KETA), C(-BIAS - 1 TO KETA)
CALL PI(L)
CALL SDIV(L,9*31*14699) !'1710541690073718870111737129379 で割る
CALL SDIV(L,304439)
CALL SDIV(L,49368533)
CALL SDIV2(L,27751800277)
CALL SMUL(L,101*1657) !' 544482330679994391053312457583 倍する
CALL SMUL(L,3580259)
CALL SMUL(L,5392007)
CALL SMUL2(L,168529144463)
CALL LCLR(X)
LET X(0)=1
CALL LSUB2(X, L)
CALL LCOPY(B, X)
CALL LCOPY(S, X)
LET  NN = 2
DO
   CALL LMUL2(B, X)
   CALL LCOPY(C, B)
   CALL SDIV(C, NN)
   LET  NN = NN + 1
   CALL LADD2(S, C)
LOOP UNTIL ZERO(B)<>0
LET S(SIGN)=-1
CALL LN1710541690073718870111737129379(L)
CALL LADD2(S,L)
CALL LN544482330679994391053312457583(L)
CALL LSUB(S,L,X)
END SUB

EXTERNAL  SUB SMUL2(X(),XA) !'多倍長数X = 多倍長数X * 整数XA (1E+11～1E+14程度まで)
DIM A(-BIAS*4 TO KETA*4+3)  !'(基数10000 → 基数10 に変換して計算) 精度落ち対策
LET  SIGNA = X(SIGN)
IF XA >= 0 THEN LET SG=1 ELSE LET SG=-1
LET  XX = ABS(XA)
FOR I=-BIAS TO KETA
   LET A(I*4)=MOD(INT(X(I)/1000),10)
   LET A(I*4+1)=MOD(INT(X(I)/100),10)
   LET A(I*4+2)=MOD(INT(X(I)/10),10)
   LET A(I*4+3)=MOD(X(I),10)
NEXT I
MAT A=(XX)*A
FOR I=KETA*4 TO -BIAS*4+1 STEP -1
   IF A(I) >= 10 THEN
      LET  R = INT(A(I) / 10)
      LET  A(I) = MOD(A(I),10)
      LET  A(I - 1) = A(I - 1) + R
   END IF
NEXT I
IF A(-BIAS*4) >= 10 THEN
   PRINT "OVER FLOW in SMUL2"
   STOP
END IF
FOR I=-BIAS TO KETA-1
   LET X(I)=A(I*4)*1000+A(I*4+1)*100+A(I*4+2)*10+A(I*4+3)
NEXT I
LET X(SIGN)= SIGNA * SG
END SUB

EXTERNAL  SUB SDIV2(X(),XA) !'多倍長数X = 多倍長数X / 整数XA (1E+11～1E+14程度まで)
DIM A(-BIAS*4 TO KETA*4+3)  !'(基数10000 → 基数10 に変換して計算) 精度落ち対策
IF XA=0 THEN
   PRINT "ERROR in SDIV2"
   STOP
END IF
LET  SIGNA = X(SIGN)
LET  SG=SGN(XA)
LET  XX = ABS(XA)
FOR I=-BIAS TO KETA
   LET A(I*4)=MOD(INT(X(I)/1000),10)
   LET A(I*4+1)=MOD(INT(X(I)/100),10)
   LET A(I*4+2)=MOD(INT(X(I)/10),10)
   LET A(I*4+3)=MOD(X(I),10)
NEXT I
FOR I=-BIAS*4 TO KETA*4-1
   LET  R = A(I) - INT(A(I) / XX) * XX
   LET  A(I) = INT(A(I) / XX)
   LET  A(I + 1) = A(I + 1) + R * 10
NEXT I
LET A(KETA*4)=INT(A(KETA*4)/XX)
FOR I=-BIAS TO KETA-1
   LET X(I)=A(I*4)*1000+A(I*4+1)*100+A(I*4+2)*10+A(I*4+3)
NEXT I
LET X(SIGN)= SIGNA * SG
END SUB

EXTERNAL  SUB LMUL(A(),B(),C())!'多倍長同士の乗算 C=A*B
LET N=(KETA+BIAS)*2
IF INT(LOG2(N))<>LOG2(N) THEN
   PRINT "ERROR in LMUL"
   STOP
END IF
OPTION BASE 0
DIM AA(N*2),BB(N*2),CC(N*2)
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, COS(PI/N), SIN(PI/N), AA)
CALL CDFT(2*N, COS(PI/N), SIN(PI/N), BB)
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, COS(PI/N), -SIN(PI/N), CC)
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
   ELSEIF C(I)<0 THEN
      LET C(I)=C(I)+10000
      LET C(I-1)=C(I-1)-1
   END IF
NEXT I
LET C(SIGN)=A(SIGN)*B(SIGN)
END SUB

EXTERNAL  SUB LMUL2 (A(), B()) !'多倍長同士の乗算 A=A*B
DIM C(-BIAS-1 TO KETA)
CALL LMUL(A,B,C)
CALL LCOPY(A,C)
END SUB

EXTERNAL  SUB CDFT(N, WR, WI, A()) !'ネット上から入手した(※原版はFORTRAN)
LET  WMR = WR
LET  WMI = WI
LET  M = N
DO WHILE M > 4
   LET  L = M / 2
   LET  WKR = 1
   LET  WKI = 0
   LET  WDR = 1 - 2 * WMI * WMI
   LET  WDI = 2 * WMI * WMR
   LET  SS = 2 * WDI
   LET  WMR = WDR
   LET  WMI = WDI
   FOR J = 0 TO N - M STEP M
      LET  I = J + L
      LET  XR = A(J) - A(I)
      LET  XI = A(J + 1) - A(I + 1)
      LET  A(J) = A(J) + A(I)
      LET  A(J + 1) = A(J + 1) + A(I + 1)
      LET  A(I) = XR
      LET  A(I + 1) = XI
      LET  XR = A(J + 2) - A(I + 2)
      LET  XI = A(J + 3) - A(I + 3)
      LET  A(J + 2) = A(J + 2) + A(I + 2)
      LET  A(J + 3) = A(J + 3) + A(I + 3)
      LET  A(I + 2) = WDR * XR - WDI * XI
      LET  A(I + 3) = WDR * XI + WDI * XR
   NEXT J
   FOR K = 4 TO L - 4 STEP 4
      LET  WKR = WKR - SS * WDI
      LET  WKI = WKI + SS * WDR
      LET  WDR = WDR - SS * WKI
      LET  WDI = WDI + SS * WKR
      FOR J = K TO N - M + K STEP M
         LET  I = J + L
         LET  XR = A(J) - A(I)
         LET  XI = A(J + 1) - A(I + 1)
         LET  A(J) = A(J) + A(I)
         LET  A(J + 1) = A(J + 1) + A(I + 1)
         LET  A(I) = WKR * XR - WKI * XI
         LET  A(I + 1) = WKR * XI + WKI * XR
         LET  XR = A(J + 2) - A(I + 2)
         LET  XI = A(J + 3) - A(I + 3)
         LET  A(J + 2) = A(J + 2) + A(I + 2)
         LET  A(J + 3) = A(J + 3) + A(I + 3)
         LET  A(I + 2) = WDR * XR - WDI * XI
         LET  A(I + 3) = WDR * XI + WDI * XR
      NEXT J
   NEXT  K
   LET  M = L
LOOP
IF M > 2 THEN
   FOR J = 0 TO N - 4 STEP 4
      LET  XR = A(J) - A(J + 2)
      LET  XI = A(J + 1) - A(J + 3)
      LET  A(J) = A(J) + A(J + 2)
      LET  A(J + 1) = A(J + 1) + A(J + 3)
      LET  A(J + 2) = XR
      LET  A(J + 3) = XI
   NEXT J
END IF
IF N > 4  THEN CALL BITRV2(N, A)
END SUB 		 

Re: 多倍長LOG(π)を求める

 投稿者：しばっち  投稿日：2010年 1月23日(土)19時23分47秒
  > No.984[元記事へ]

続き

EXTERNAL  SUB BITRV2(N, A())
LET  M = N / 4
LET  M2 = 2 * M
LET  N2 = N - 2
LET  K = 0
FOR J = 0 TO M2 - 4 STEP 4
   IF J < K THEN
      LET  XR = A(J)
      LET  XI = A(J + 1)
      LET  A(J) = A(K)
      LET  A(J + 1) = A(K + 1)
      LET  A(K) = XR
      LET  A(K + 1) = XI
   ELSEIF J > K THEN
      LET  J1 = N2 - J
      LET  K1 = N2 - K
      LET  XR = A(J1)
      LET  XI = A(J1 + 1)
      LET  A(J1) = A(K1)
      LET  A(J1 + 1) = A(K1 + 1)
      LET  A(K1) = XR
      LET  A(K1 + 1) = XI
   END IF
   LET  K1 = M2 + K
   LET  XR = A(J + 2)
   LET  XI = A(J + 3)
   LET  A(J + 2) = A(K1)
   LET  A(J + 3) = A(K1 + 1)
   LET  A(K1) = XR
   LET  A(K1 + 1) = XI
   LET  L = M
   DO WHILE K >= L
      LET  K = K - L
      LET  L = L / 2
   LOOP
   LET  K = K + L
NEXT J
END SUB

EXTERNAL  SUB PI(X()) !'π値
CALL LCLR (X)
LET X(SIGN)=1
FOR I=0 TO KETA
   READ IF MISSING THEN EXIT FOR:X(I)
NEXT I
DATA 0003
DATA 1415,9265,3589,7932,3846,2643,3832,7950,2884,1971,6939,9375,1058,2097,4944,5923,0781,6406,2862,0899,8628,0348,2534,2117,0679
DATA 8214,8086,5132,8230,6647,0938,4460,9550,5822,3172,5359,4081,2848,1117,4502,8410,2701,9385,2110,5559,6446,2294,8954,9303,8196
DATA 4428,8109,7566,5933,4461,2847,5648,2337,8678,3165,2712,0190,9145,6485,6692,3460,3486,1045,4326,6482,1339,3607,2602,4914,1273
DATA 7245,8700,6606,3155,8817,4881,5209,2096,2829,2540,9171,5364,3678,9259,0360,0113,3053,0548,8204,6652,1384,1469,5194,1511,6094
DATA 3305,7270,3657,5959,1953,0921,8611,7381,9326,1179,3105,1185,4807,4462,3799,6274,9567,3518,8575,2724,8912,2793,8183,0119,4912
DATA 9833,6733,6244,0656,6430,8602,1394,9463,9522,4737,1907,0217,9860,9437,0277,0539,2171,7629,3176,7523,8467,4818,4676,6940,5132
DATA 0005,6812,7145,2635,6082,7785,7713,4275,7789,6091,7363,7178,7214,6844,0901,2249,5343,0146,5495,8537,1050,7922,7968,9258,9235
DATA 4201,9956,1121,2902,1960,8640,3441,8159,8136,2977,4771,3099,6051,8707,2113,4999,9998,3729,7804,9951,0597,3173,2816,0963,1859
DATA 5024,4594,5534,6908,3026,4252,2308,2533,4468,5035,2619,3118,8171,0100,0313,7838,7528,8658,7533,2083,8142,0617,1776,6914,7303
DATA 5982,5349,0428,7554,6873,1159,5628,6388,2353,7875,9375,1957,7818,5778,0532,1712,2680,6613,0019,2787,6611,1959,0921,6420,1989
DATA 3809,5257,2010,6548,5863,2788,6593,6153,3818,2796,8230,3019,5203,5301,8529,6899,5773,6225,9941,3891,2497,2177,5283,4791,3151
DATA 5574,8572,4245,4150,6959,5082,9533,1168,6172,7855,8890,7509,8381,7546,3746,4939,3192,5506,0400,9277,0167,1139,0098,4882,4012
DATA 8583,6160,3563,7076,6010,4710,1819,4295,5596,1989,4676,7837,4494,4825,5379,7747,2684,7104,0475,3464,6208,0466,8425,9069,4912
DATA 9331,3677,0289,8915,2104,7521,6205,6966,0240,5803,8150,1935,1125,3382,4300,3558,7640,2474,9647,3263,9141,9927,2604,2699,2279
DATA 6782,3547,8163,6009,3417,2164,1219,9245,8631,5030,2861,8297,4555,7067,4983,8505,4945,8858,6926,9956,9092,7210,7975,0930,2955
DATA 3211,6534,4987,2027,5596,0236,4806,6549,9119,8818,3479,7753,5663,6980,7426,5425,2786,2551,8184,1757,4672,8909,7777,2793,8000
DATA 8164,7060,0161,4524,9192,1732,1721,4772,3501,4144,1973,5685,4816,1361,1573,5255,2133,4757,4184,9468,4385,2332,3907,3941,4333
DATA 4547,7624,1686,2518,9835,6948,5562,0992,1922,2184,2725,5025,4256,8876,7179,0494,6016,5346,6804,9886,2723,2791,7860,8578,4383
DATA 8279,6797,6681,4541,0095,3883,7863,6095,0680,0642,2512,5205,1173,9298,4896,0841,2848,8626,9456,0424,1965,2850,2221,0661,1863
DATA 0674,4278,6220,3919,4945,0471,2371,3786,9609,5636,4371,9172,8746,7764,6575,7396,2413,8908,6583,2645,9958,1339,0478,0275,9009
DATA 9465,7640,7895,1269,4683,9835,2595,7098,2582,2620,5224,8940
END SUB

EXTERNAL  SUB LN1710541690073718870111737129379(X()) !'LOG(1710..)値
CALL LCLR(X)
LET X(SIGN)=1
FOR I=0 TO KETA
   READ IF MISSING THEN EXIT FOR:X(I)
NEXT I
DATA 0069
DATA 6143,6288,7993,3268,3004,4228,2256,8485,2026,6785,3596,6394,4004,2459,4994,5470,1366,5214,8581,6415,4697,3145,7866,1523,7272
DATA 6551,2059,0127,0534,7351,7820,4920,0008,7989,9541,7800,2163,7623,0895,4887,3922,8357,3520,5004,0375,5683,9209,4332,0242,9540
DATA 1748,4435,3241,8232,0917,6048,6351,8267,6655,4103,1920,3713,6534,2808,6848,0068,3923,2800,1886,4684,2394,0286,5595,7822,8870
DATA 6824,4083,5434,5674,3002,0302,6430,7887,0940,1081,9799,5193,2835,3565,5690,4021,5682,1679,9587,2726,5741,7738,9966,8303,9374
DATA 4065,5626,7370,2317,3186,2152,1129,2984,3621,1141,6888,7106,2708,7294,1105,4881,3419,1360,3002,9792,2275,7776,7209,3100,7189
DATA 1333,2926,5906,5310,6574,2333,0024,4444,6043,7221,5935,4199,0820,3015,5200,1473,4002,4862,2284,7807,2731,8261,8953,2953,7309
DATA 4162,4615,8900,1280,1186,1445,3416,9135,4404,4277,0805,7623,1006,4866,6480,4750,1673,5797,3590,6210,4828,8466,8064,0621,5369
DATA 6822,3112,8752,9922,5756,4306,5244,5723,7782,5053,2215,5848,6850,0720,0010,0839,6060,6790,9162,7960,6519,9616,7812,3819,6481
DATA 2850,4009,5322,6067,3979,0496,8839,2466,5255,0930,4737,8771,8972,6183,6781,6485,6742,6847,7569,0628,3272,4958,9055,0595,1667
DATA 0318,1311,7655,2636,5449,8892,8994,6549,8185,5605,1410,8418,3957,3941,9585,7859,8354,6284,4894,3324,7101,4235,7355,6297,7705
DATA 6284,7948,3668,3827,3890,3916,7911,8682,2761,6156,8577,8088,8066,9387,2851,5278,4084,2809,0879,0882,6939,4694,4014,9549,7450
DATA 2911,5763,0278,3692,1766,1804,8124,9781,4074,8058,6307,1025,2785,0760,5599,0083,4081,9099,0502,6793,1794,6341,6502,0176,2868
DATA 3796,2139,9878,6835,7170,1164,2205,6497,0107,8877,3771,7488,8657,5263,9761,2416,6278,3016,1690,0953,1383,8181,7952,3541,1635
DATA 6891,5224,4409,2332,4874,2319,4524,2884,0889,5022,8222,5650,3677,3971,7593,6214,1500,2361,0686,4206,4017,9737,8034,3987,7461
DATA 5150,7437,9961,5876,6881,2221,9103,9820,1354,6600,2760,3220,1344,0667,5700,0559,0702,6840,8532,2596,9193,2134,9548,5155,8877
DATA 1115,9327,4668,8319,5471,2802,0874,9640,1131,7970,8168,2428,1598,3983,6017,5960,3281,6942,4686,9574,9764,9433,1216,4959,2387
DATA 0873,8907,8785,5093,6238,0367,4603,0384,3624,6755,5657,2459,9473,2847,7735,8517,5154,0299,6744,3440,0577,4508,6886,3276,1136
DATA 3026,8720,1951,3001,5676,9739,1394,0643,7786,6988,0521,2185,7078,9937,0618,4466,8235,0108,8657,3294,3409,0176,9410,4535,3487
DATA 5923,4566,3163,5408,2042,0034,7213,0857,1778,0402,9447,8452,8732,6232,4880,1913,2182,3249,5532,6094,1509,8743,6779,9084,9870
DATA 0148,8381,3012,8627,2066,3865,1370,4664,5664,9339,4367,0190,7057,6215,5163,9757,6009,6007,7540,5813,9957,0761,4465,5132,0577
DATA 5503,8417,6557,8706,8939,7814,1500,5247,4757,4949,9209,4162
END SUB

EXTERNAL  SUB LN544482330679994391053312457583(X()) !'LOG(5444..)値
CALL LCLR(X)
LET X(SIGN)=1
FOR I=0 TO KETA
   READ IF MISSING THEN EXIT FOR:X(I)
NEXT I
DATA 0068
DATA 4696,3300,2143,9266,5590,0800,8743,3179,3315,0312,4115,3479,0888,5308,1371,0786,7772,6582,8354,8087,7911,8425,1658,1624,6863
DATA 3534,2720,9432,7548,7948,6284,7447,1357,2875,4696,1898,3726,1025,5558,4562,0907,4477,6791,4979,6390,9355,8684,7181,2383,6067
DATA 5125,9038,3164,5872,4149,1463,9290,8100,8221,3858,4661,7368,8894,4392,1782,5529,1379,0706,7263,6717,7094,2864,3556,3900,7575
DATA 8220,5281,5426,2279,7293,8843,9725,1183,2375,0876,2708,1484,0473,9670,8817,3223,4819,0068,2413,2727,7063,0621,1175,4750,7107
DATA 8671,3441,3378,9146,9346,1923,8237,0526,7108,5413,6817,1946,4759,6110,2330,5046,3222,6584,6590,6347,9478,9162,1229,8258,0639
DATA 8629,0999,7506,3639,3766,9619,4748,6069,6562,5828,2302,4778,3079,2466,1935,4340,5305,9867,9470,7594,5310,3809,3799,2927,1647
DATA 2356,9062,6748,9596,1245,0147,5975,5100,6274,8333,6949,6146,2909,4481,6231,2366,8112,4857,0275,2529,1010,3037,7042,6637,0392
DATA 1670,0491,9722,4984,9411,3815,8869,2731,5703,6430,4132,8603,8424,8368,8011,9123,4089,8308,1853,8980,8362,2095,2457,1928,3638
DATA 0053,3991,4538,1624,6079,4415,8152,2109,6785,7592,9725,2195,0650,0526,5590,0422,3146,3046,8647,5876,9130,4732,3389,1245,8068
DATA 7812,3293,5439,3929,1890,9399,6703,8878,7077,5335,6941,7431,1782,2323,8457,0834,4553,9586,3123,6397,5711,0712,1356,5910,7625
DATA 2193,7131,0681,5285,6609,1095,4149,6754,7763,2836,8796,8730,4951,2823,6551,9299,6014,9348,6341,1187,0646,7366,9790,3242,0076
DATA 3159,2614,7305,1479,3803,0883,1726,8138,1385,3562,1576,6493,8245,0199,7851,4583,9672,9911,8496,3168,7333,7765,0503,4451,0920
DATA 5090,5214,6787,1314,6772,6862,3736,8572,0251,7508,8096,1615,0470,5139,3677,8073,5561,0597,8884,4253,9075,0348,7834,3251,6913
DATA 6549,9352,6478,6587,3687,6446,5803,7343,1607,6931,6148,2364,8107,6286,4138,6620,4422,4124,9848,5744,8327,3790,8854,9012,5488
DATA 6930,5446,1559,2103,7235,1105,4337,1632,8662,1702,8001,9460,1977,8871,6694,9328,1262,7307,0409,3747,7022,2234,8243,0766,1057
DATA 1794,9925,2660,2292,2678,0591,1757,8530,6030,0003,8275,0990,6180,5366,9308,0990,1612,3368,9265,3857,0729,5854,4739,4278,6728
DATA 2111,2446,8161,9364,0415,6814,3644,1779,3609,4265,9090,3322,2336,7923,4096,4595,1538,1568,3121,3374,7747,7340,2591,1985,8297
DATA 1375,4564,3882,8472,3644,4193,1371,2433,7491,7818,1991,7233,1335,2034,1296,1177,3339,7135,6549,2440,4227,1453,0204,3370,3408
DATA 9449,0546,0281,3840,9872,7548,2470,7503,6138,4063,2497,2268,8601,9381,1855,0714,6065,2571,9613,8541,6472,6325,5532,3632,6764
DATA 0684,0804,6469,2615,3651,3647,2317,6321,2200,3895,4491,6134,5365,5560,5641,2651,0168,7531,0092,4478,8313,3640,3413,7462,9541
DATA 1403,9443,4183,7689,8013,2877,3469,9226,6078,9523,0380,9531
END SUB

以下に共通ルーチンをコピペする(多倍長LOG(5)を求める、より) 		 

戻る