|
|
十進モードでは、不正な値（答え）になる。
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
|
|