十進BASIC 第2掲示板 過去ログ

!定数
DECLARE EXTERNAL NUMERIC MultiPrecision.RADIX, MultiPrecision.L2
DECLARE EXTERNAL NUMERIC MultiPrecision.c0(), MultiPrecision.c1()

!演算
DECLARE EXTERNAL FUNCTION MultiPrecision.lNum2Str$, MultiPrecision.lcomp
DECLARE EXTERNAL SUB MultiPrecision.lStr2Num, MultiPrecision.lcopy
DECLARE EXTERNAL SUB MultiPrecision.ladd, MultiPrecision.lsub
DECLARE EXTERNAL SUB MultiPrecision.lmul, MultiPrecision.ldiv

DECLARE EXTERNAL SUB MultiPrecision.llmul, MultiPrecision.ldivqr
!----------------------------------------


!Σ[k=1,n]1/k の計算

LET t0=TIME


DIM P(0 TO L2),Q(0 TO L2)
CALL lStr2Num("0",P) !t=P/Q=1
CALL lStr2Num("1",Q)

DIM T1(0 TO L2),T2(0 TO L2)
FOR k=1 TO 1129
   CALL lmul(P,k, T1) !t=t+1/k=(P*k+Q)/(Q*k)　通分
   CALL ladd(T1,Q, P)
   CALL lmul(Q,k, T2)
   CALL lcopy(T2,Q)
NEXT k
PRINT lNum2Str$(P)
PRINT lNum2Str$(Q)


DIM A(0 TO L2),B(0 TO L2)
CALL lcopy(P,A) !LET A=P !copy it
CALL lcopy(Q,B) !LET B=Q

LET c=0 !繰り返し回数 debug
DO UNTIL lcomp(B,c0)=0 !最大公約数を求める
   LET c=c+1 !debug
   CALL ldivqr(A,B, T1,T2) !LET R=MOD(A,B) !MOD(x,y)=x-y*INT(x/y)
   CALL lcopy(B,A) !LET A=B
   CALL lcopy(T2,B) !LET B=R
LOOP
PRINT "c=";c !debug
PRINT "gcd=";lNum2Str$(A) !PRINT "gcd=";A !debug

CALL ldivqr(P,A, T1,T2) !LET P=P/A !約分
CALL lcopy(T1,P)
CALL ldivqr(Q,A, T1,T2) !LET Q=Q/A
CALL lcopy(T1,Q)

PRINT "P=";lNum2Str$(P) !PRINT "P=";P !結果の表示
PRINT "Q=";lNum2Str$(Q) !PRINT "Q=";Q


PRINT TIME-t0

END


MODULE MultiPrecision !正の多桁（多倍長）整数の計算

!桁の数字の列を格納する配列a()を考える。
!上位桁　　　　　　下位桁
!a(L2),a(L2-1),…,a(1),a(0)　の構造で表すことができる。
!言い換えると、n進数表記の各桁がa()となる。
!たとえば、100進数とすると、a(k)は正の整数２桁（0～99）となる。
!例. 12345は1*100*100+23*100+45だから、a(2)=1、a(1)=23、a(0)=45。


SHARE NUMERIC KETA
PUBLIC NUMERIC RADIX,L2

LET KETA=4 !桁数 1,2,3,4　※32bitか64bit整数の範囲
LET RADIX=10^KETA !基数

LET L=3000 !求める桁数　※
LET L2=INT(L/KETA)+1 !配列の大きさ


PUBLIC NUMERIC c0(0 TO 1000) !※
CALL lStr2Num("0",c0) !定数０
PUBLIC NUMERIC c1(0 TO 1000) !※
CALL lStr2Num("1",c1) !定数１

SHARE NUMERIC aa(0 TO 1000),bb(0 TO 1000),m(0 TO 1000) !作業用　※


!----- 下位の演算ルーチン

PUBLIC SUB ladd
EXTERNAL SUB ladd(a(),b(), c()) !多桁＋多桁　※C=A+B
   LET cy=0 !桁上がり
   FOR i=0 TO L2 !下の桁から
      LET d=a(i)+b(i)+cy !同じ桁で
      IF d<RADIX THEN
         LET c(i)=d
         LET cy=0
      ELSE
         LET c(i)=d-RADIX
         LET cy=1 !上の桁へ
      END IF
   NEXT i
   IF cy>0 THEN
      PRINT "加算オーバーフロー"
      STOP
   END IF
END SUB
PUBLIC SUB lsub
EXTERNAL SUB lsub(a(),b(), c()) !多桁－多桁　※A>B、C=A-B
   LET brrw=0 !借り
   FOR i=0 TO L2 !下の桁から
      LET d=a(i)-b(i)-brrw !同じ桁で
      IF d>=0 THEN
         LET c(i)=d
         LET brrw=0
      ELSE
         LET c(i)=d+RADIX
         LET brrw=1 !上の桁から借りる
      END IF
   NEXT i
   IF brrw>0 THEN
      PRINT "減算A-Bで、A>Bではありません。"
      STOP
   END IF
END SUB
PUBLIC SUB lmul
EXTERNAL SUB lmul(a(),b, c())  !多桁×正の整数（0～255）※10000進数なら32767
   LET cy=0
   FOR i=0 TO L2 !下の桁から
      LET d=a(i)*b+cy
      LET cy=INT(d/RADIX) !桁上がり
      LET c(i)=MOD(d,RADIX) !この桁
   NEXT i
   IF cy>0 THEN
      PRINT "乗算オーバーフロー"
      STOP
   END IF
END SUB
PUBLIC SUB ldiv
EXTERNAL SUB ldiv(a(),b, c()) !多桁÷正の整数（0～255）※10000進数なら32767
   IF b=0 THEN
      PRINT "０では割れません。"
      STOP
   END IF

   LET r=0 !余り
   FOR i=L2 TO 0 STEP -1 !上の桁から
      LET d=a(i)+r
      LET c(i)=INT(d/b) !商はこの桁
      LET r=MOD(d,b)*RADIX !余りを下の桁へ
   NEXT i
END SUB

PUBLIC SUB lcopy
EXTERNAL SUB lcopy(a(), b()) !コピー B=A
   FOR i=0 TO L2 !mat b=a
      LET b(i)=a(i) !同じ桁で
   NEXT i
END SUB
PUBLIC FUNCTION lcomp
EXTERNAL FUNCTION lcomp(a(),b()) !比較　A>Bなら1、A=Bなら0、A<Bなら-1
   FOR i=L2 TO 0 STEP -1 !上の桁から
      IF a(i)>b(i) THEN
         LET lcomp=1
         EXIT FUNCTION
      ELSEIF a(i)<b(i) THEN
         LET lcomp=-1
         EXIT FUNCTION
      END IF
   NEXT i
   LET lcomp=0
END FUNCTION

!----- 入力、出力ルーチン

PUBLIC SUB lStr2Num
EXTERNAL SUB lStr2Num(a$, a()) !数字列を多桁整数に変換する
   FOR i=0 TO L2
      LET A(i)=0
   NEXT i
   FOR i=LEN(a$) TO 1 STEP -KETA
      LET d=INT((LEN(a$)-i)/KETA)
      IF d>L2 THEN
         PRINT "桁数が足りません。"
         STOP
      END IF
      LET a(d)=VAL(a$(i-KETA+1:i))
   NEXT i
END SUB
PUBLIC FUNCTION lNum2Str$
EXTERNAL FUNCTION lNum2Str$(a()) !多桁整数を数字列に変換する
   LET aMSD=GetMSD(a)
   IF aMSD<0 THEN !０なら
      LET a$=" 0"
   ELSE
      LET a$=" "&STR$(a(aMSD)) !１桁目
      FOR i=aMSD-1 TO 0 STEP -1 !２桁目以降
         LET a$=a$&right$("000"&STR$(a(i)),4)
      NEXT i
   END IF
   LET lNum2Str$=a$
END FUNCTION


!----- 上位の演算ルーチン

PUBLIC SUB llmul
EXTERNAL SUB llmul(a(),b(), c())  !多桁×多桁　※C=A*B
   FOR i=0 TO L2*2 !mat c=zer
      LET c(i)=0
   NEXT i

   LET aMSD=GetMSD(a)
   LET bMSD=GetMSD(b)
   IF aMSD<0 OR bMSD<0 THEN EXIT SUB !乗数、被乗数が０なら

   FOR j=0 TO bMSD !乗数：下の桁から
      LET cy=0
      IF b(j)>0 THEN !０は計算しない
         FOR i=0 TO aMSD !被乗数：下の桁から
            LET d=a(i)*b(j)+cy + c(i+j) !累積　※筆算参照 O(n^2)
            LET c(i+j)=MOD(d,RADIX) !この桁
            LET cy=INT(d/RADIX) !桁上がり
         NEXT i
         IF cy>0 THEN LET c(j+aMSD+1)=cy !上位桁へ
      END IF
   NEXT j
END SUB

EXTERNAL FUNCTION GetMSD(a()) !最上位桁の位置を得る
   FOR i=L2 TO 0 STEP -1
      IF a(i)<>0 THEN EXIT FOR
   NEXT i
   LET GetMSD=i !桁位置　※０なら-1
END FUNCTION

PUBLIC SUB ldivqr
EXTERNAL SUB ldivqr(a(),b(), q(),r()) !多桁÷多桁　※商q 余りr
   LET bMSD=GetMSD(b)
   IF bMSD<0 THEN !除数が０かどうか確認する
      PRINT "０では割れません。"
      STOP
   END IF

   FOR i=0 TO L2 !商を０とする
      LET q(i)=0
   NEXT i

   LET v=b(bMSD) !最上位桁の値を得る
   LET nk=INT(RADIX/(v+1)) !b(MSD)をRADIX/2以上にするための最小係数
   IF nk>1 THEN
      CALL lmul(b,nk, bb) !除数の最上位桁をRADIX/2以上、RADIX未満にする
      CALL lmul(a,nk, aa) !aもnk倍する
   ELSE
      CALL lcopy(b, bb) !×１
      CALL lcopy(a, aa)
   END IF

   LET t=GetMSD(bb)
   DO UNTIL lcomp(aa,bb)<0 !a<bなら終了
      LET s=GetMSD(aa)
      IF aa(s)>=bb(t) THEN
         IF CompareOffset(aa,bb,s-t)>=0 THEN !a>=b*RADIX^(s-t)を検査する
            LET u=s-t !商の最上位桁位置uおよび値q(u)の候補を求める
            LET q(u)=1
         ELSE
            LET u=s-t-1
            LET q(u)=RADIX-1
         END IF
      ELSE
         LET u=s-t-1
         LET q(u)=INT( (aa(s)*RADIX+aa(s-1))/bb(t) )
      END IF

      CALL lmul(bb,q(u), m)
      DO WHILE CompareOffset(aa,m,u)<0 !a>=m=b*q(0)を満たす最大のq(u)を求める
         LET q(u)=q(u)-1
         CALL lmul(bb,q(u), m)
      LOOP

      CALL SubOffset(aa,m, u) !a=a-b*q(u)*RADIX^u
   LOOP

   IF nk>1 THEN
      CALL ldiv(aa,nk, r) !余り 1/nk倍
   ELSE
      CALL lcopy(aa, r)
   END IF
END SUB
EXTERNAL FUNCTION CompareOffset(a(),b(),n) !桁をずらして比較する
   LET aMSD=GetMSD(a)
   LET bMSD=GetMSD(b)+n
   IF aMSD>bMSD THEN
      LET CompareOffset=1
   ELSEIF aMSD<bMSD THEN
      LET CompareOffset=-1
   ELSE
      FOR i=aMSD TO n STEP -1
         IF a(i)>b(i-n) THEN
            LET CompareOffset=1
            EXIT FUNCTION
         END IF
         IF a(i)<b(i-n) THEN
            LET CompareOffset=-1
            EXIT FUNCTION
         END IF
      NEXT i
      DO !a(aMSD)～a(n)とb(bMSD)～b(0)が一致する場合、
         IF NOT i>=0 THEN EXIT DO !a(n-1)～a(0)に非零のものがあれば、aの方が大きい
         IF a(i)<>0 THEN
            LET CompareOffset=1
            EXIT FUNCTION
         END IF
         LET i=i-1
      LOOP
      LET CompareOffset=0
   END IF
END FUNCTION
EXTERNAL SUB SubOffset(a(),b(),n) !桁をずらして減算　※a=a-b*RADIX^n
   LET brrw=0 !借り
   LET bMSD=GetMSD(b) !最上位桁までbをひく
   FOR i=0 TO bMSD
      LET d=a(i+n)-b(i)-brrw !同じ桁で
      IF d>=0 THEN
         LET a(i+n)=d
         LET brrw=0
      ELSE
         LET a(i+n)=d+RADIX
         LET brrw=1 !上の桁から借りる
      END IF
   NEXT i
   DO !上位桁への桁借り
      IF NOT brrw<>0 THEN EXIT DO
      IF a(i+n)<>0 THEN
         LET a(i+n)=a(i+n)-1
         LET brrw=0
      ELSE
         LET a(i+n)=RADIX-1
         LET brrw=1
      END IF
      LET i=i+1
   LOOP
END SUB

END MODULE