『多項式演算で利用できる一般化された拡張互除法』

多項式演算の基本演算であるユークリッドの一般化拡張互除法を Maple 内で実装する方法を紹介します。一般化拡張互除法とは、与えられた一変数多項式 p(x), q(x) に対して A(x)*p(x) + B(x)*q(x) = x^i を計算するための算法です。ここで i = 0..deg(A)+deg(B) となります。

>    restart:

ある多項式 p(x), q(x) に対して拡張互除法 A(x)*p(x) + B(x)*q(x) = 1 となる A(x), B(x) を求めるための関数 gcdex が Maple には用意されています。

>    p := x^4 - 3*x^3 + x^2 - x + 1:
q :=       2*x^3 - x^2 + x + 3:

>    gcd(p, q);

1

拡張互除法を行うには、gcdex 関数を利用します。3番目・4番目の引数には、A(x), B(x) として割り当てられる変数を指定します。

>    gcdex(p, q, x, 's', 't');

1

変数 s, t には A(x), B(x) となる多項式が割り当てられています。

>    {s, t};

{404/2219+9/317*x+94/2219*x^2, 605/2219-88/2219*x+86/2219*x^2-47/2219*x^3}

実際に、s*p + t*q を計算してみると、

>    expand(s*p + t*q);

1

確かに1となっていることがわかります。

さて、一般化拡張互除法は A(x)*p(x) + B(x)*q(x) = x^i となる A(x), B(x) を計算するものです。

算法自体は難しいものではなく、両辺の 恒等性 を利用して求めることができます。たとえば右辺=1となる場合を例に考えてみます。

まず、A(x), B(x) の多項式を生成します。それぞれの次数に注意してください。

>    Apoly := add(a[i]*x^i, i=0..(degree(q,x)-1));

Apoly := a[0]+a[1]*x+a[2]*x^2

>    Bpoly := add(b[i]*x^i, i=0..degree(p,x));

Bpoly := b[0]+b[1]*x+b[2]*x^2+b[3]*x^3+b[4]*x^4

この2つの多項式 Apoly および Bpoly を、それぞれ p(x), q(x) に掛け合わせた多項式を作ります。

>    euclidPoly := expand(Apoly*p + Bpoly*q);

euclidPoly := a[0]+2*b[1]*x^4+3*b[0]+a[2]*x^6-3*a[2]*x^5+a[2]*x^4-a[2]*x^3+2*b[0]*x^3-b[0]*x^2+b[0]*x-b[1]*x^3+b[1]*x^2+2*b[2]*x^5-b[2]*x^4+b[2]*x^3+2*b[3]*x^6-b[3]*x^5+b[3]*x^4+2*b[4]*x^7-b[4]*x^6+b[4...
euclidPoly := a[0]+2*b[1]*x^4+3*b[0]+a[2]*x^6-3*a[2]*x^5+a[2]*x^4-a[2]*x^3+2*b[0]*x^3-b[0]*x^2+b[0]*x-b[1]*x^3+b[1]*x^2+2*b[2]*x^5-b[2]*x^4+b[2]*x^3+2*b[3]*x^6-b[3]*x^5+b[3]*x^4+2*b[4]*x^7-b[4]*x^6+b[4...

この多項式が右辺の項と一致するように、各係数を決めていきます。

euclidPoly の各係数を取り出してみましょう。

>    coeffList := seq(coeff(euclidPoly, x, i)=d(i), i=0..degree(euclidPoly,x));

coeffList := a[0]+3*b[0] = d(0), b[0]-a[0]+3*b[1]+a[1] = d(1), -b[0]+b[1]+a[0]-a[1]+3*b[2]+a[2] = d(2), -a[2]+2*b[0]-b[1]+b[2]-3*a[0]+a[1]+3*b[3] = d(3), 2*b[1]+a[2]-b[2]+b[3]+a[0]-3*a[1]+3*b[4] = d(4)...
coeffList := a[0]+3*b[0] = d(0), b[0]-a[0]+3*b[1]+a[1] = d(1), -b[0]+b[1]+a[0]-a[1]+3*b[2]+a[2] = d(2), -a[2]+2*b[0]-b[1]+b[2]-3*a[0]+a[1]+3*b[3] = d(3), 2*b[1]+a[2]-b[2]+b[3]+a[0]-3*a[1]+3*b[4] = d(4)...

ここで各右辺は d(i) としました。いま考えているのは、右辺=1の場合なので、d(0)、つまり x^0 の項以外はすべて0となるように d の関数を決めます。

>    d := n -> if n=0 then 1 else 0 end if:

>    coeffList;

a[0]+3*b[0] = 1, b[0]-a[0]+3*b[1]+a[1] = 0, -b[0]+b[1]+a[0]-a[1]+3*b[2]+a[2] = 0, -a[2]+2*b[0]-b[1]+b[2]-3*a[0]+a[1]+3*b[3] = 0, 2*b[1]+a[2]-b[2]+b[3]+a[0]-3*a[1]+3*b[4] = 0, -3*a[2]+2*b[2]-b[3]+b[4]+a...
a[0]+3*b[0] = 1, b[0]-a[0]+3*b[1]+a[1] = 0, -b[0]+b[1]+a[0]-a[1]+3*b[2]+a[2] = 0, -a[2]+2*b[0]-b[1]+b[2]-3*a[0]+a[1]+3*b[3] = 0, 2*b[1]+a[2]-b[2]+b[3]+a[0]-3*a[1]+3*b[4] = 0, -3*a[2]+2*b[2]-b[3]+b[4]+a...

この方程式を各 a_i, b_i に関して解いてみます。

>    solvedCoeff := solve({coeffList}, {coeffs(Apoly,x), coeffs(Bpoly,x)});

solvedCoeff := {b[4] = 0, b[1] = -88/2219, b[0] = 605/2219, b[3] = -47/2219, b[2] = 86/2219, a[1] = 9/317, a[0] = 404/2219, a[2] = 94/2219}

この結果を Apoly, Bpoly にそれぞれ適用(代入)することで、A(x), B(x) に相当する多項式を得ることができます。

>    spoly := subs(solvedCoeff, Apoly);

spoly := 404/2219+9/317*x+94/2219*x^2

>    tpoly := subs(solvedCoeff, Bpoly);

tpoly := 605/2219-88/2219*x+86/2219*x^2-47/2219*x^3

上記の多項式が gcdex 関数で得られた結果と一致していることは自明です。

>    expand(spoly*p + tpoly*q);

1

>    evalb(s=spoly);

true

>    evalb(t=tpoly);

true

ここで行った一般化拡張互除法をひとつの関数としてまとめておきます。

>    monoTerm := (n::integer, d::integer) -> if n=d then 1 else 0 end if:
GeneralExtPolyGCD := proc(p::algebraic, q::algebraic, var::name, deg::integer)
  local Adeg, Bdeg, Apoly, Bpoly, tmpPoly, coefList, solvedCoefs;
  option `Copyright (c) 2003 by Cybernet Systems, Co.,Ltd. All rights reserved.`;
 
  Adeg := degree(q, var)-1;
  Bdeg := degree(p, var);
  Apoly := add(a[i]*x^i, i=0..Adeg);
  Bpoly := add(b[i]*x^i, i=0..Bdeg);

  tmpPoly := expand(Apoly*p+Bpoly*q);

  coefList := seq(
                coeff(tmpPoly, x, i)=monoTerm(i, deg),
                i=0..degree(tmpPoly, x)
              );
                  
  solvedCoefs := solve({coefList}, {coeffs(Apoly,x), coeffs(Bpoly,x)});
  return( subs(solvedCoefs, Apoly), subs(solvedCoefs, Bpoly) );
 
end proc:

以下は、上記関数のテスト用のコードです;

>    p := x^5+2*x^3-x+1:
q := x^4-x+1:

>    gcd(p, q);

1

>    ss, tt := GeneralExtPolyGCD(p, q, x, 2); # this generates A*p+B*q = x^2

ss, tt := 20/29-3/29*x-4/29*x^2-15/29*x^3, -20/29+3/29*x+33/29*x^2+4/29*x^3+15/29*x^4

>    expand(ss*p + tt*q);

x^2

>   

閉じる