プログラミング雑感  新JavaScript入門  JavaScript,Neo-Generation  掲示板  表紙
16. 多項式クラス(9)  18. 多項式クラス(11) 
プログラミング雑感
Written 1/17/03
17. 多項式クラス(10)
前回いそいでプログラムを作ったような気がしますが、 もう少しゆっくり考えていきたいと思います。
■基本的な考え方
例えば次のような多項式を因数分解することを考えましょう。 
    x2 + 2xy + y2 - z2 - 3y - z + 2
3変数の多項式ですが、これを x の多項式として見てみましょう。
1変数多項式のときの考え方を参考にすると、 因数分解できるとしたら x の1次多項式同士の積になるはずです。 上の多項式を f とすると、 
  f(x,y,z) = a0 + a1x + a2x2
  g(x,y,z) = b0 + b1x
  h(x,y,z) = c0 + c1x
ai, bi, ci は y, z の多項式です。 
  f(0,y,z) = g(0,y,z)h(0,y,z) = y2 - z2 - 3y - z + 2
  f(1,y,z) = g(1,y,z)h(1,y,z) = y2 - z2 - y - z + 3
g(0, y, z) は y2 - z2 - 3y - z + 2 を割り切らなければならないので、 
    g(0,y,z) = ±1, ±(y2 - z2 - 3y - z + 2)
あるいは、y2 - z2 - 3y - z + 2 を因数分解した式になるはずです。g(1, y, z) も同様です。
こうして3変数多項式の因数分解は2変数多項式の因数分解に置き換えられました。 結局多変数多項式は1変数多項式に帰着することができます。
次からは1変数多項式をスピードアップすることを考えます。
■とりあえず書いてみる
せっかくのwebなのでJavaScriptで書けばいいようなものですが、 やっぱりちょっと書きにくいということで、C++で書きます。 最後にはJavaScriptで書きたいと思います。
1変数多項式のときの考え方に沿って、 というかこのときのプログラムと基本的に同じように書きます。
多項式クラスを次のように vector<int> を継承して、 四則演算などを定義します。 
    class Polynomial : public vector<int> {
        vector<Polynomial>    factors;    //因数分解の結果を格納する配列
    public:
        Polynomial() { }
        Polynomial(const int *a, int size);
        Polynomial(const vector<int>& a);
        bool operator ==(const Polynomial& a);
        bool operator ==(int a);
        bool operator !=(const Polynomial& a);
        bool operator !=(int a);
        Polynomial& operator +=(const Polynomial& a);
        Polynomial& operator -=(const Polynomial& a);
        Polynomial& operator *=(const Polynomial& a);
        Polynomial& operator /=(const Polynomial& a);
        Polynomial& operator %=(const Polynomial& a);
        void factorize();
        int value(int x) const;           //多項式関数の値
        void print() const;               //因数分解の形での出力
    };
STLなんか使ったら遅いという話もありますが、 そんなことは問題にならず何桁もスピードアップする予定です。
とりあえず作ってみたものが次です。 
    factorize1a.cpp
このプログラムは、 
    f(x) = 1 + 2x + x11
を因数分解します。 Athlon1.2GHzで9秒かかります。
書いていませんでしたが、私のとりあえずの目標は、 16次の多項式を数秒で因数分解するというものなので、かなり遅いです。
■代入値を選ぶ
今の方法は、 
    f(x0), f(x1), ..., f(xn)
のそれぞれの約数の組 
    p0, p1, ..., pn
全てについて行列の計算を行っているのでした（因数分解できない限り）。 だから、 
    f(x0), f(x1), ..., f(xn)
の約数の数が少なければそれだけ速く演算できます。
実際、どれくらいの数の約数があるのかを調べてみると次のようになりました。 
    f(-5) = -48828134 32
    f(-4) =  -4194311  8
    f(-3) =   -177152 44
    f(-2) =     -2051  8
    f(-1) =        -2  4
    f(0)  =         1  2
    f(1)  =         4  6
    f(2)  =      2053  4
    f(3)  =    177154 16
    f(4)  =   4194313  8
    f(5)  =  48828136 32
右端に書いてある数字が約数の個数です。
上ではこのうち下の f(0) 〜 f(5) 使っていますが、 これを f(-2) 〜 f(3) とすれば、 計算の回数はそれぞれの約数の個数の積を2で割ったものだから、 
    98,304 -> 12,288
と激減します（別に連番でなくてもよいのですが、 逆行列がきれいにならない気がするので）。 これを実際にプログラムにしたものが次です。 
    factorize1b.cpp
これはプログラムもちょっと改良して速くなってます。 上の例で、0.6秒と15倍速くなりました。
■オーバーフロー
この勢いで12次の因数分解もしたいと思います。 例えば、次のような多項式です。 
    1 + 3x + x12
ところがこれを最初の方法で計算しようとするとオーバーフローするんですね。 
    612 = 2176782336 > 2147483647 = 231 - 1
その次の方法ではもちろんそんなことはないのですが、 多項式関数の値を求めるときのみあらかじめオーバーフローに備えておきます。 
    try {
        f[i] = value(i - nhalf);
    }
    catch(int e) {
        f[i] = 0;
        continue;
    }
メンバ関数 value の中でオーバーフローが発生したら例外を投げるようにしています。
これが、 
    factorize1c.cpp
です。これで、 
    1 + 3x + x12
を因数分解したところ、4.6秒かかりました。
first, prev, next, exit