再帰呼出しでの不具合

 投稿者：山中和義  投稿日：2011年 9月18日(日)16時53分11秒
  十進モードでは、不正な値（答え）になる。
1000桁モードでは、INT関数が桁あふれする。


!多項式に対する拡張ユークリッド互除法

!多項式f(x)=A[m]x^m+ … +A[1]x+A[0]、g(x)=B[n]x^n+ … +B[1]x+B[0]、m≧nとして、
!f(x)S(x)+g(x)T(x)=gcd(f(x),g(x))=C(x)となるS(x),T(x),C(x)を求める。

PUBLIC NUMERIC MAX_DEGREE !次数　※m^2以上
LET MAX_DEGREE=10
!------------------------------


DATA 3 !次数
DATA 1,0,-3,1 !g(x)=x^3-3x+1
DATA 3
DATA 1,1,0,6 !f(x)=x^3+x^2+6


READ M !次数
DIM A(0 TO M) !係数
FOR i=M TO 0 STEP -1 !f(x)=A[m]x^m+ … +A[1]x+A[0]
   READ A(i)
NEXT i

READ N !次数
DIM B(0 TO N) !係数
FOR i=N TO 0 STEP -1 !g(x)=B[n]x^n+ … +B[1]x+B[0]
   READ B(i)
NEXT i


DIM P(0 TO MAX_DEGREE),Q(0 TO MAX_DEGREE),C(0 TO MAX_DEGREE)
MAT P=ZER
MAT Q=ZER
MAT C=ZER
CALL ExGCD(M,A,N,B, pp,P,qq,Q,cc,C) !拡張ユークリッド互除法


MAT PRINT P; !結果を表示する
MAT PRINT Q;
MAT PRINT C;

END


!拡張ユークリッド互除法
!　f(x)=A[m]x^m+ … +A[1]x+A[0]、g(x)=B[n]x^n+ … +B[1]x+B[0]、m≧nとして、
!　f(x)S(x)+g(x)T(x)=gcd(f(x),g(x))=C(x)となるS(x),T(x),C(x)を求める。

EXTERNAL SUB ExGCD(aa,A(),bb,B(), ss,S(),tt,T(),cc,C()) !拡張ユークリッド互除法
IF bb=0 AND B(0)=0 THEN !IF b=0 THEN
   LET S(0)=1 !s=1　※f(x)*1+0*0=f(x)とする
   LET ss=0
   LET T(0)=0 !t=0
   LET tt=0
   MAT C=A !c=a
   LET cc=aa
ELSE
   DIM Q(0 TO MAX_DEGREE),R(0 TO MAX_DEGREE)
   IF aa=0 AND bb=0 THEN !定数項のみ
      LET Q(0)=INT(A(0)/B(0))
      LET qq=0
      LET R(0)=MOD(A(0),B(0))
      LET rr=0
   ELSE
      CALL poly_div(aa,A,bb,B, qq,Q,rr,R) !q=INT(a/b), r=MOD(a,b)
   END IF

   CALL ExGCD(bb,B,rr,R, tt,T,ss,S,cc,C) !k=n-1,…,3,2 まで続ける

   DIM W(0 TO MAX_DEGREE)
   CALL poly_mul(ss,S,qq,Q, ww,W) ! t=u-v*q
   MAT T=T-W
   LET tt=ww
END IF
END SUB


!補助ルーチン

!演算関連

EXTERNAL SUB poly_mul(aa,A(),bb,B(), ss,S()) !乗算 S=A*B　※S≠A、S≠B
FOR i=aa TO 0 STEP -1
   FOR j=bb TO 0 STEP -1
      LET S(i+j)=S(i+j)+A(i)*B(j) !すべての係数をかける
   NEXT j
NEXT i
LET ss=aa+bb !次数
END SUB


EXTERNAL SUB poly_div(aa,A(),bb,B(), qq,Q(),rr,R()) !除算　※被除数=商*除数+余り
IF bb=0 AND B(0)=0 THEN !0なら
   PRINT "0で割ることはできません。"
   STOP
ELSE
   MAT Q=ZER !商
   MAT R=ZER !余り
   MAT R=A
   FOR t=aa TO bb STEP -1 !被除数の次数が除数のより大きいなら
      IF R(t)<>0 THEN !係数が０以外なら
         LET k=R(t)/B(bb) !商の係数、その次数
         LET w=t-bb
         LET Q(w)=k !商

         FOR i=bb TO 0 STEP -1 !余り　※R=A-k*B
            LET R(w+i)=R(w+i)-k*B(i)
         NEXT i
      END IF
   NEXT t
   LET qq=MAX(aa-bb,0) !次数
   IF aa>=bb THEN LET rr=MAX(bb-1,0) ELSE LET rr=aa
END IF
END SUB
 		 

Re: 再帰呼出しでの不具合

 投稿者：白石和夫  投稿日：2011年 9月18日(日)17時56分27秒
  > No.1663[元記事へ]

2005行と2065行を追加しました。
有理数モードで正しい答えになるようであれば，2000行と2050行の計算誤差が原因ではないでしょうか。


1000 !多項式に対する拡張ユークリッド互除法
1010
1020 !多項式f(x)=A[m]x^m+ … +A[1]x+A[0]、g(x)=B[n]x^n+ … +B[1]x+B[0]、m≧nとして、
1030 !f(x)S(x)+g(x)T(x)=gcd(f(x),g(x))=C(x)となるS(x),T(x),C(x)を求める。
1040
1050 PUBLIC NUMERIC MAX_DEGREE !次数　※m^2以上
1060 LET MAX_DEGREE=10
1070 !------------------------------
1080
1090
1100 DATA 3 !次数
1110 DATA 1,0,-3,1 !g(x)=x^3-3x+1
1120 DATA 3
1130 DATA 1,1,0,6 !f(x)=x^3+x^2+6
1140
1150
1160 READ M !次数
1170 DIM A(0 TO M) !係数
1180 FOR i=M TO 0 STEP -1 !f(x)=A[m]x^m+ … +A[1]x+A[0]
1190    READ A(i)
1200 NEXT i
1210
1220 READ N !次数
1230 DIM B(0 TO N) !係数
1240 FOR i=N TO 0 STEP -1 !g(x)=B[n]x^n+ … +B[1]x+B[0]
1250    READ B(i)
1260 NEXT i
1270
1280
1290 DIM P(0 TO MAX_DEGREE),Q(0 TO MAX_DEGREE),C(0 TO MAX_DEGREE)
1300 MAT P=ZER
1310 MAT Q=ZER
1320 MAT C=ZER
1330 CALL ExGCD(M,A,N,B, pp,P,qq,Q,cc,C) !拡張ユークリッド互除法
1340
1350
1360 MAT PRINT P; !結果を表示する
1370 MAT PRINT Q;
1380 MAT PRINT C;
1390
1400 END
1410
1420
1430 !拡張ユークリッド互除法
1440 !　f(x)=A[m]x^m+ … +A[1]x+A[0]、g(x)=B[n]x^n+ … +B[1]x+B[0]、m≧nとして、
1450 !　f(x)S(x)+g(x)T(x)=gcd(f(x),g(x))=C(x)となるS(x),T(x),C(x)を求める。
1460
1470 EXTERNAL SUB ExGCD(aa,A(),bb,B(), ss,S(),tt,T(),cc,C()) !拡張ユークリッド互除法
1480 IF bb=0 AND B(0)=0 THEN !IF b=0 THEN
1490    LET S(0)=1 !s=1　※f(x)*1+0*0=f(x)とする
1500    LET ss=0
1510    LET T(0)=0 !t=0
1520    LET tt=0
1530    MAT C=A !c=a
1540    LET cc=aa
1550 ELSE
1560    DIM Q(0 TO MAX_DEGREE),R(0 TO MAX_DEGREE)
1570    IF aa=0 AND bb=0 THEN !定数項のみ
1580       LET Q(0)=INT(A(0)/B(0))
1590       LET qq=0
1600       LET R(0)=MOD(A(0),B(0))
1610       LET rr=0
1620    ELSE
1630       CALL poly_div(aa,A,bb,B, qq,Q,rr,R) !q=INT(a/b), r=MOD(a,b)
1640    END IF
1650
1660    CALL ExGCD(bb,B,rr,R, tt,T,ss,S,cc,C) !k=n-1,…,3,2 まで続ける
1670
1680    DIM W(0 TO MAX_DEGREE)
1690    CALL poly_mul(ss,S,qq,Q, ww,W) ! t=u-v*q
1700    MAT T=T-W
1710    LET tt=ww
1720 END IF
1730 END SUB
1740
1750
1760 !補助ルーチン
1770
1780 !演算関連
1790
1800 EXTERNAL SUB poly_mul(aa,A(),bb,B(), ss,S()) !乗算 S=A*B　※S≠A、S≠B
1810 FOR i=aa TO 0 STEP -1
1820    FOR j=bb TO 0 STEP -1
1830       LET S(i+j)=S(i+j)+A(i)*B(j) !すべての係数をかける
1840    NEXT j
1850 NEXT i
1860 LET ss=aa+bb !次数
1870 END SUB
1880
1890
1900 EXTERNAL SUB poly_div(aa,A(),bb,B(), qq,Q(),rr,R()) !除算　※被除数=商*除数+余り
1910 IF bb=0 AND B(0)=0 THEN !0なら
1920    PRINT "0で割ることはできません。"
1930    STOP
1940 ELSE
1950    MAT Q=ZER !商
1960    MAT R=ZER !余り
1970    MAT R=A
1980    FOR t=aa TO bb STEP -1 !被除数の次数が除数のより大きいなら
1990       IF R(t)<>0 THEN !係数が０以外なら
2000          LET k=R(t)/B(bb) !商の係数、その次数
2005          PRINT k
2010          LET w=t-bb
2020          LET Q(w)=k !商
2030
2040          FOR i=bb TO 0 STEP -1 !余り　※R=A-k*B
2050             LET R(w+i)=R(w+i)-k*B(i)
2060          NEXT i
2065          MAT PRINT R
2070       END IF
2080    NEXT t
2090    LET qq=MAX(aa-bb,0) !次数
2100    IF aa>=bb THEN LET rr=MAX(bb-1,0) ELSE LET rr=aa
2110 END IF
2120 END SUB

 		 

Re: 再帰呼出しでの不具合

 投稿者：山中和義  投稿日：2011年 9月18日(日)19時28分52秒
  > No.1664[元記事へ]

白石和夫さんへのお返事です。

> 2005行と2065行を追加しました。
> 有理数モードで正しい答えになるようであれば，2000行と2050行の計算誤差が原因ではないでしょうか。

> 1480 IF bb=0 AND B(0)=0 THEN !IF b=0 THEN

1480 IF bb=0 AND ABS(B(0))<=1E-13 THEN !IF b=0 THEN

として対応します。
通常の整数の場合は、INT関数で０になるのでよいのですが、
多項式の場合は、「これに相当する処理がないため、有理数モードでの実行が前提となる」ということですね。

お騒がせしました。
 		 

戻る