|
|
サブルーチン mPOLY.LIB
!mPOLY.LIB
!多変数多項式の演算（加減乗算） p(a,b,c, … ,y,z)=ΣN*a^A*b^B*c^C* … *y^Y*z^Z
! 式 p(項,変数)
! 係数 変数
! 項＼ 0 v1 v2 v3 …
! 0: k - - - …
! 1: c1 e11 e12 e13 … 項 c1*(v1^e11)*(v2^e12)*(v3^e13)* … の意
! 2: c2 e21 e22 e23 …
! 3: c3 e31 e32 e33 …
! :
! :
! k: ck ek1 ek2 ek3 …
! 例 p=a^2*b+3*a*b+2*b
! p(0 TO 3, 0 TO 2) として、
! 係数 変数
! 項＼ 0 a b
! 0: 3 - -
! 1: 1 2 1
! 2: 3 1 1
! 3: 2 0 1
!表示関連
EXTERNAL SUB PolyPrint(p(,)) !多変数多項式を表示する
OPTION ARITHMETIC RATIONAL
LET K=p(0,0) !項の数
IF K>0 THEN
FOR i=1 TO K !各項について
LET t=p(i,0) !係数の部分
IF t<>0 THEN !０は非表示
LET FLG=0
IF t>0 THEN !符号を演算子として表示する
IF i>1 THEN PRINT "+"; !第1項目は表示しない
ELSE
PRINT "-";
END IF
LET t=ABS(t) !定数の部分
IF t<>1 THEN !１は非表示
CALL DispNum(t)
LET FLG=1 !分子を表示した
END IF
LET s=0
FOR e=1 TO N !次数の部分
IF p(i,e)>0 THEN !べき乗数が正の場合
CALL DispVar(v$(e:e),p(i,e))
LET FLG=1 !分子を表示した
ELSEIF p(i,e)<0 THEN
LET s=s+1
END IF
NEXT e
IF s>0 THEN !べき乗数が負の場合、分数で表示する
IF FLG=0 THEN PRINT "1"; !分子が定数±１のとき
IF s>1 THEN PRINT "/("; ELSE PRINT "/";
FOR e=1 TO N !次数の部分
IF p(i,e)<0 THEN CALL DispVar(v$(e:e),ABS(p(i,e)))
NEXT e
IF s>1 THEN PRINT ")";
LET FLG=1 !分母を表示した
END IF
IF FLG=0 THEN PRINT "1"; !定数±１のみのとき
END IF
NEXT i
IF K=1 AND p(1,0)=0 THEN PRINT "0"; !定数項０のみのとき
ELSE
PRINT "NUL"; !未定義である
END IF
END SUB
EXTERNAL SUB DispNum(t) !数値を表示する ※1,(-2),(3/4),(-5/6)など
OPTION ARITHMETIC RATIONAL
IF t=INT(t) AND t>=0 THEN !非負の整数なら
PRINT STR$(t);
ELSE !負の整数、有理数なら
PRINT "(";STR$(t);")";
END IF
END SUB
EXTERNAL SUB DispVar(x$,t) !x^t形式で変数を表示する
OPTION ARITHMETIC RATIONAL
IF t<>0 THEN !x^0は非表示
PRINT x$; !変数
IF t<>1 THEN !x^1の1は非表示
PRINT "^"; !次数
CALL DispNum(t)
END IF
END IF
END SUB
EXTERNAL SUB PolyPrintCollect(p(,),x$) !変数xの多項式とみなして表示する
OPTION ARITHMETIC RATIONAL
DIM t(0 TO MAX_TERM, 0 TO N)
CALL PolyCopy(p, t) !次数の昇順に並べる
LET e=idx(x$)
LET K=t(0,0)
FOR i=1 TO K-1 !バブルソート O(n^2)
LET A=t(i,e)
FOR J=i+1 TO K
LET B=t(J,e)
IF A>B THEN !iとjを交換する
LET A=B
FOR m=0 TO N !swap it
LET t(K+1,m)=t(i,m) !iを作業領域へ
NEXT m
FOR m=0 TO N
LET t(i,m)=t(J,m) !j→i
NEXT m
FOR m=0 TO N
LET t(J,m)=t(K+1,m) !i（作業領域）→j
NEXT m
END IF
NEXT J
NEXT i
DIM s(0 TO MAX_TERM, 0 TO N)
LET u=t(1,e)
LET i=1
DO
PRINT "+("; !同類項（次数が同じ値のもの）で括る
LET J=0
DO
LET J=J+1
FOR m=0 TO N !copy it
LET s(J,m)=t(i,m)
NEXT m
LET s(J,e)=0 !この次数は除く
LET i=i+1 !次へ
LET w=t(i,e)
LOOP WHILE i<=K AND w=u !異なる次数が現れたら
LET s(0,0)=J
CALL PolyPrint(s) !表示する
PRINT ")";
CALL DispVar(x$,u) !この次数の変数を表示する
LET u=w !次へ
LOOP UNTIL i>K
END SUB
!演算関連
EXTERNAL SUB PolyAdd(p(,),q(,), r(,)) !加算 r=p+q ※r≠p、r≠qの配列を指定する
OPTION ARITHMETIC RATIONAL
CALL PolyCopy(p, r) !そのままコピーする
LET M=p(0,0) !項の数
FOR J=1 TO q(0,0) !項を末尾に加える
FOR e=0 TO N !次数、係数
LET r(M+J,e)=q(J,e)
NEXT e
NEXT J
LET r(0,0)=M+q(0,0) !項の数
CALL PolySimplify(r) !同類項をまとめる、0サプレス
END SUB
EXTERNAL SUB PolySubtract(p(,),q(,), r(,)) !減算 r=p-q ※r≠p、r≠qの配列を指定する
OPTION ARITHMETIC RATIONAL
DIM t(0 TO MAX_TERM, 0 TO N)
CALL PolyMultiplyC(q,-1, t) !(-1)倍
CALL PolyAdd(p,t, r) !p-q
END SUB
EXTERNAL SUB PolyMultiply(p(,),q(,), r(,)) !乗算 r=p*q ※r≠p、r≠qの配列を指定する
OPTION ARITHMETIC RATIONAL
LET M=0
FOR i=1 TO p(0,0) !項と項をかける
FOR J=1 TO q(0,0)
LET M=M+1
FOR e=1 TO N !次数
LET r(M,e)=p(i,e)+q(J,e)
NEXT e
LET r(M,0)=p(i,0)*q(J,0) !係数
NEXT J
NEXT i
LET r(0,0)=M !項の数
CALL PolySimplify(r) !同類項をまとめる、0サプレス
END SUB
EXTERNAL SUB PolyMultiplyC(p(,),k, r(,)) !乗算（定数倍） r=k*p
OPTION ARITHMETIC RATIONAL
LET M=p(0,0)
FOR i=1 TO M
FOR e=1 TO N !次数
LET r(i,e)=p(i,e)
NEXT e
LET r(i,0)=k*p(i,0) !係数
NEXT i
LET r(0,0)=M !項の数
END SUB
EXTERNAL SUB PolyDivide(p(,), f1(,), r(,)) !除算（単項式による） r=p/f1
OPTION ARITHMETIC RATIONAL
LET M=p(0,0)
FOR i=1 TO M
FOR e=1 TO N !次数
LET r(i,e)=p(i,e)-f1(1,e)
NEXT e
LET r(i,0)=p(i,0)/f1(1,0) !係数
NEXT i
LET r(0,0)=M !項の数
END SUB
EXTERNAL SUB PolyCommonFactor(p(,), f1(,)) !共通因数を得る（共通因数で括る）
OPTION ARITHMETIC RATIONAL
LET G=p(1,0) !copy it
FOR e=1 TO N
LET f1(1,e)=p(1,e)
NEXT e
FOR i=2 TO p(0,0)
LET G=gcd(G,p(i,0)) !係数
FOR e=1 TO N !次数
LET f1(1,e)=MIN(f1(1,e),p(i,e))
NEXT e
NEXT i
LET f1(1,0)=G
LET f1(0,0)=1 !項の数
END SUB
EXTERNAL FUNCTION gcd(a,b) !最大公約数
OPTION ARITHMETIC RATIONAL
DO UNTIL b=0
LET t=b
LET b=MOD(a,b)
LET a=t
LOOP
LET gcd=ABS(a)
END FUNCTION
!演算関連2
EXTERNAL SUB PolyComposition(p(,),x$,q(,), r(,)) !合成関数 P(x=Q) ※r≠p、r≠qの配列を指定する
OPTION ARITHMETIC RATIONAL
MAT r=ZER
LET r(0,0)=1 !定数０
LET r(1,0)=0
CALL PolyDegree(p,x$, mx,mn) !次数の最大値、最小値を得る
DIM T2(0 TO MAX_TERM, 0 TO N), T3(0 TO MAX_TERM, 0 TO N)
FOR ii=mx TO mn STEP -1 !ホーナー法による
CALL PolyMultiply(r,q, T3) !T[n]*X+T[n-1]
CALL PolyCoefficient(p,x$,ii, T2) !係数 ※多項式
CALL PolyAdd(T3,T2, r)
NEXT ii
FOR ii=0 TO mn-1 !残り部分
CALL PolyMultiply(r,q, T3)
CALL PolyCopy(T3,r)
NEXT ii
END SUB
EXTERNAL SUB PolyPowN(f(,),K, fk(,)) !べき乗 f^k
OPTION ARITHMETIC RATIONAL
DIM t(0 TO MAX_TERM, 0 TO N)
MAT fk=ZER !定数 1
LET fk(0,0)=1
LET fk(1,0)=1
FOR i=1 TO K
CALL PolyMultiply(fk,f, t) !fk=fk*f
CALL PolyCopy(t, fk)
NEXT i
END SUB
!deg(f)≧deg(g)≧0、gの最高次の係数は単項式 として、f÷g=商q余りr を求める
EXTERNAL SUB PolyQuotientRemainder(f(,),g(,),x$, q(,),r(,)) !f÷g=商q余りr
OPTION ARITHMETIC RATIONAL
DIM t1(0 TO MAX_TERM, 0 TO N), t2(0 TO MAX_TERM, 0 TO N)
DIM p(0 TO MAX_TERM, 0 TO N), u(0 TO MAX_TERM, 0 TO N)
CALL PolyDegree(f,x$, mx1,mn1) !次数を得る
CALL PolyDegree(g,x$, mx2,mn2)
MAT q=ZER !q,rの初期値
LET q(0,0)=1
CALL PolyCopy(f, r)
CALL PolyCoefficient(g,x$,mx2, p) !P
FOR i=mx1 TO mx2 STEP -1 !筆算による
CALL PolyCoefficient(r,x$,i, u) !A
LET e=idx(x$)
FOR J=0 TO u(0,0) !x^k倍
LET u(J,e)=u(J,e)+(i-mx2)
NEXT J
CALL PolyDivide(u,p, u) !Q=A/P
CALL PolyAdd(q,u, t2) !q=ΣQ
CALL PolyCopy(t2, q)
CALL PolyMultiply(u,g, t1) !r=r-Qg
CALL PolySubtract(r,t1, t2)
CALL PolyCopy(t2, r)
NEXT i
CALL PolySimplify(r) !0に注意する
END SUB
!共通ルーチン
EXTERNAL SUB PolyCopy(p(,), q(,)) !コピーする LET q=p
OPTION ARITHMETIC RATIONAL
LET K=p(0,0) !項の数
FOR i=1 TO K
FOR e=0 TO N !係数、次数
LET q(i,e)=p(i,e)
NEXT e
NEXT i
LET q(0,0)=K !項の数
END SUB
EXTERNAL SUB PolySimplify(r(,)) !同類項をまとめるなどの簡約化を行う
OPTION ARITHMETIC RATIONAL
LET K=r(0,0) !項の数
LET W=0
FOR i=1 TO K
LET T=r(i,0) !係数が０のものを篩う
IF T<>0 THEN
FOR j=i+1 TO K !同類項を探す
CALL IsEqualDegree(r,i, r,j, rc) !次数が同じなら、同類項！
IF rc<>0 THEN
LET T=T+r(j,0) !前方へ吸収する
LET r(j,0)=0 !後方は削除する
END IF
NEXT j
IF T<>0 THEN !係数が０なら、無効！
LET W=W+1
LET r(W,0)=T !係数
IF i>W THEN !前方へ吸収する ※ガーベジ・コレクション
FOR e=1 TO N !次数
LET r(W,e)=r(i,e)
NEXT e
END IF
END IF
END IF
NEXT i
IF W=0 THEN
LET r(0,0)=1 !定数０
FOR e=0 TO N
LET r(1,e)=0
NEXT e
ELSE
LET r(0,0)=W !項の数
END IF
END SUB
EXTERNAL SUB PolySimplify2(r(,)) !同類項をまとめるなどの簡約化を行う ※i^2=-1、ω^2+ω+1=0 を適用する
OPTION ARITHMETIC RATIONAL
LET K=r(0,0) !項の数
LET W=K
FOR i=1 TO K !虚数i、ωの条件式を適用する ※x^k=Aのとき、x^m≡A^int(m/k) x^mod(m,k)
LET e=idx("i") !虚数i
IF e>0 THEN
LET T=INT(r(i,e)/2) !i^2=-1を適用する ※i^3=-iなど
LET r(i,0)=(-1)^MOD(T,2)*r(i,0) !係数
LET r(i,e)=MOD(r(i,e),2) !次数
END IF
LET e=idx("w") !ω（x^3-1=0の3つの解の内、虚数解の1つをωとする）
IF e>0 THEN
LET r(i,e)=MOD(r(i,e),3) !ω^3=1を適用する ※ω^5=ω^2など
!ω^2+ω+1=0を適用する。すなわち、Aω^2=-Aω-Aとする。
IF r(i,e)=2 THEN !Aω^2なら
LET r(i,0)=-r(i,0) !係数 -Aω
LET r(i,e)=1 !次数
LET W=W+1 !-Aを後方に追加する
FOR J=0 TO N !copy it
LET r(W,J)=r(i,J)
NEXT J
LET r(W,e)=0
END IF
END IF
NEXT i
LET r(0,0)=W !項の数
CALL PolySimplify(r) !a+bi、a+bω（a+bω+cω^2）の形へ
END SUB
EXTERNAL SUB PolyCoefficient(p(,),x$,m, t(,)) !変数xにおけるm次の係数（多項式）を得る
OPTION ARITHMETIC RATIONAL
LET K=0
LET J=idx(x$) !桁位置を確認する
FOR i=1 TO p(0,0) !すべての項を検索する
IF p(i,J)=m THEN !同じ次数なら
LET K=K+1
FOR e=0 TO N !copy it
LET t(K,e)=p(i,e)
NEXT e
LET t(K,J)=0 !該当する変数は除く
END IF
NEXT i
IF K=0 THEN
LET t(0,0)=1 !定数０
FOR e=0 TO N
LET t(1,e)=0
NEXT e
ELSE
LET t(0,0)=K !項の数
END IF
END SUB
EXTERNAL SUB PolyDegree(p(,),x$, mx,mn) !変数xにおける次数の最大値と最小値を得る
OPTION ARITHMETIC RATIONAL
LET e=idx(x$) !桁位置を確認する
LET mx=p(1,e) !最大
LET mn=p(1,e) !最小
FOR i=2 TO p(0,0) !すべての項を検索する
IF p(i,e)>mx THEN LET mx=p(i,e)
IF p(i,e)<mn THEN LET mn=p(i,e)
NEXT i
END SUB
つづく
|
|