多変数多項式の計算 mPOLY.LIB

 投稿者:山中和義  投稿日:2014年11月15日(土)23時16分42秒
  サブルーチン 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 !0は非表示
         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 !1は非表示
            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"; !分子が定数±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"; !定数±1のみのとき
      END IF
   NEXT i

   IF K=1 AND p(1,0)=0 THEN PRINT "0"; !定数項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 !定数0
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) !係数が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 !係数が0なら、無効!
         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 !定数0
   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 !定数0
   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


つづく

 

Re: 多変数多項式の計算 mPOLY.LIB

 投稿者:山中和義  投稿日:2014年11月15日(土)23時18分41秒
  > No.3547[元記事へ]

つづき



!マクロ

EXTERNAL FUNCTION idx(x$) !変数a,b,c,…を1,2,3,…へ
OPTION ARITHMETIC RATIONAL
LET idx=POS(v$,x$)
END FUNCTION


EXTERNAL SUB IsEqualDegree(p(,),i, q(,),j, rc) !多項式pの第i項の次数と多項式qの第j項の次数が等しいかどうか確認する
OPTION ARITHMETIC RATIONAL
LET rc=0
FOR e=1 TO N !次数が等しくないなら
   IF p(i,e)<>q(j,e) THEN EXIT SUB
NEXT e
LET rc=1 !等しい
END SUB


EXTERNAL SUB IsConst(p(,),i, rc) !多項式pの第i項が定数かどうか確認する
OPTION ARITHMETIC RATIONAL
LET rc=0
FOR e=1 TO N !次数が0以外なら
   IF p(i,e)<>0 THEN EXIT SUB
NEXT e
LET rc=1 !定数
END SUB



!----------------------------------------------

!文字列で表された多項式を符号化する Σ(-M/N)a^(-P/Q)  例 xy^5+(4/3)x^2y+z^(-1)

EXTERNAL SUB PolySet(s$, a(,)) !多項式を符号化する
OPTION ARITHMETIC RATIONAL
MAT a=ZER

LET p=1 !文字位置

LET K=1 !項数

CALL GetToken(p,s$, t$) !1番目の単項式の±は符号とみなす
LET s=1 !符号
IF t$="+" THEN
   LET p=p+1 !eat it
ELSEIF t$="-" THEN
   LET p=p+1 !eat it
   LET s=-1
END IF

CALL GetToken(p,s$, t$) !トークン
IF t$="" THEN
   PRINT "符号のみで項がありません。"; p
   STOP
END IF
DO WHILE t$<>"" !終端以外なら
   IF s=0 THEN
      IF t$="+" THEN !+演算子
         LET s=1
      ELSEIF t$="-" THEN !-演算子
         LET s=-1
      END IF
      LET p=p+1 !eat it

      LET K=K+1
   ELSE
      LET a(K,0)=s
      DO
         IF t$<>"" AND idx(t$)>0THEN !変数なら
            LET p=p+1 !eat it

            LET t=idx(t$)

            CALL GetToken(p,s$, t$)
            LET v=1
            IF t$="^" THEN !べき乗
               LET p=p+1 !eat it
               CALL GetRational(p,s$, v)
            END IF
            LET a(K,t)=a(K,t)+v
         ELSE
            CALL GetRational(p,s$, v) !係数
            LET a(K,0)=a(K,0)*v
         END IF

         CALL GetToken(p,s$, t$)
      LOOP UNTIL t$="+" OR t$="-" OR t$="" !±演算子、終端が現れるまで

      LET s=0 !次は±演算子である
   END IF

   CALL GetToken(p,s$, t$) !トークン
LOOP
IF s<>0 THEN
   PRINT "項がありません。"; p
   STOP
END IF
LET a(0,0)=K !項数
END SUB


EXTERNAL SUB GetToken(p,s$, t$) !1文字を得る
OPTION ARITHMETIC RATIONAL
LET t$=s$(p:p)
DO WHILE t$=" " !空白文字はスキップする
   LET p=p+1
   LET t$=s$(p:p)
LOOP
END SUB


EXTERNAL SUB GetRational(p,s$, v) !n,(±m/n)形式の数を得る
OPTION ARITHMETIC RATIONAL
CALL GetToken(p,s$, t$)
IF t$="(" THEN !負の数、有理数
   LET p=p+1 !eat it

   CALL GetToken(p,s$, t$)
   LET s=1 !符号を得る
   IF t$="+" THEN
      LET p=p+1 !eat it
   END IF
   IF t$="-" THEN
      LET p=p+1 !eat it
      LET s=-1
   END IF
   CALL GetInteger(p,s$, v1)
   !!!PRINT s; v1

   LET v2=1
   CALL GetToken(p,s$, t$)
   IF t$="/" THEN !有理数なら、分母を得る
      LET p=p+1 !eat it
      CALL GetInteger(p,s$, v2)
      !!!PRINT v2
   END IF
   LET v=s*v1/v2

   CALL GetToken(p,s$, t$)
   IF t$<>")" THEN
      PRINT ")がありません。"; p
      STOP
   END IF
   LET p=p+1 !eat it

ELSE !正の整数
   CALL GetInteger(p,s$, v)

END IF
END SUB


EXTERNAL SUB GetInteger(p,s$, v) !正の整数を得る
OPTION ARITHMETIC RATIONAL
CALL GetToken(p,s$, t$)
IF t$<"0" OR t$>"9" THEN
   PRINT "数字ではありません。"; p
   STOP
END IF
LET v=0
DO WHILE t$>="0" AND t$<="9" !連続する数字列を10進法の数とみなす
   LET v=v*10+VAL(t$)
   LET p=p+1
   LET t$=s$(p:p)
LOOP
END SUB



終わり

 

戻る