多項式クラス(13)

n-1行目をもうちょっと詳しく見てみましょう。

    b_n-1 = a_{n-1 0}c₀ + ... + a_{n-1 n}c_n
    a_n-1i = (-1)^n-1-i_n-1C_i/(n-1)! (i < n)

一般的にはこうなりますが、今やってる

    2 + 2x + x¹⁶

で最後に8次式で割れるなら、

    b₇ = (-c₀ + 7c₁ - 21c₂ + 35c₃ - 35c₄ + 21c₅ - 7c₆ + c₇) / 5040 + 12c₈

が整数になる必要があります。
これをちょっと変形してみましょう。

    b₇ = (-c₀ + c₇ + 7(c₁ - 3c₂ + 5c₃ - 5c₄ + 3c₅ - c₆)) / 5040 + 12c₈

これが整数ならば、

    -c₀ + c₇

が7で割れなければなりません。すなわち、b₇ を計算するまでもなく、 c₀ と c₇ だけで整数でないと判定できることがあることになります。しかもこれは結構確率が高いです。
しかし、これはいつも成り立つわけではありません。 n-1が素数のときはうまくいくことがすぐに分かりますが、例えば、7次式で割れるかどうか判定するときは、

    b₆ = (c₀ + c₆ + (-6c₁ + 15c₂ - 20c₃ + 15c₄ - 6c₅)) / 720

が整数であるかチェックするわけですが、 (6, 15, 20)の最大公約数は1なので、8次式のときのようなことはできません。
このように全然括れないことはあまりないようですが、この手はあきらめて反対側から考えましょう。すなわち、

    b₆ = ((c₀ + c₆ - 6c₁ - 6c₅ + 15c₂ + 15c₄) - 20c₃) / 720

としてみると、

    c₀ + c₆ - 6c₁ - 6c₅ + 15c₂ + 15c₄

が20で割れなければならないことが分かります。
これならどんなときも必ず成り立ちます。 b₆ をあと少しで計算できるので効果は小さそうに見えますが、効果は確実にあります。しかも、20で割れる確率は単純に考えれば1/20なので、割と効果があります。
ちなみに、8次式で割る場合は、

    b₇ = ((-c₀ + c₇ + 7c₁ - 7c₆ - 21c₂ + 21c₅) + 35(c₃ - c₄)) / 5040 + 12c₈

となり、

    -c₀ + c₇ + 7c₁ - 7c₆ - 21c₂ + 21c₅

が35で割れる必要があります。
これをコードに組み入れたものが次です。

    factorize1l.cpp

これで、

    2 + 2x + x¹⁶

を因数分解したところ、1.2秒となりました。
もう少し次元を上げてみましょう。

    2 + 2x + x¹⁸

で、48.4秒かかりました。

先ほどの8次式で割る場合の式をもうちょっと変形してみましょう。

    b₇ = (-c₀ + c₇ + 7c₁ - 7c₆ - 21c₂ + 21c₅ + 35(c₃ - c₄)) / 5040 + 12c₈
      = (-c₀ + c₇ + 7c₁ - 7c₆ + 7(-3c₂ + 3c₅ + 5c₃ - 5c₄)) / 5040 + 12c₈
      = (-c₀ + c₇ + 7(c₁ - c₆ - 3c₂ + 3c₅ + 5c₃ - 5c₄)) / 5040 + 12c₈

これは下から順に、

    -c₀ + c₇ は7で割れなければならない
    -c₀ + c₇ + 7c₁ - 7c₆ は7で割れなければならない
    -c₀ + c₇ + 7c₁ - 7c₆ - 21c₂ + 21c₅ は35で割れなければならない

ことを表しています。2番目は実質的に意味はありませんが。
この数は、係数の最大公約数を次々に求めることにより得られます。すなわち、係数は、

    7 21 35

ですが、まず右から35、次にこれと21の最大公約数を取り7、これと7の最大公約数を取り7、となります。
これをコードに組み入れたものが次です。

    factorize1m.cpp

これで、

    2 + 2x + x¹⁸

で、12.0秒かかりました。意外と速くなりますね。

上の例では結構素因数分解に時間がかかっています。
以下ではこの素因数分解を避けられるときは避けるようにしてみます。
まず、上の例でどういう整数を素因数分解しているか調べます。

    f(-8) = 18014398509481970   64
    f(-7) = 1628413597910437    16
    f(-6) = 101559956668406     16
    f(-5) = 3814697265617       64
    f(-4) = 68719476730         16
    f(-3) = 387420485           32
    f(-2) = 262142               8
    f(-1) = 1                    2
    f(0) =  2                    2
    f(1) =  5                    4
    f(2) =  262150              72
    f(3) =  387420497           16
    f(4) =  68719476746         16
    f(5) =  3814697265637        4
    f(6) =  101559956668430    384
    f(7) =  1628413597910465    64
    f(8) =  18014398509482002   32

右端はそれぞれ割り切れる整数の数です。ただし、f(0) だけは正の値だけ取ることにしているので半分になっています。
このうち9つを選んで基本的な考え方の p の組み合わせが最小になるようしていて、この場合は、

    f(-3), f(-2), ... , f(5)

を選んでいて、組み合わせの数は、

    301,989,888

となります。
ところで、この場合 f(-8) の素因数分解は計算する必要がないことが分かります。 f(0) が0,1,-1のときは素因数分解するまでもないし、そうでないときは、割り切れる整数の数は4以上です。仮に4として、f(-8) を選んだとき、すなわち

    f(-8), f(-7), ... , f(0)

を選んだときの組み合わせの数は、

    1,073,741,824

となり、素因数分解するまでもなくこの組み合わせは取り得ないことになります。同様の理由で f(7), f(8) も素因数分解する必要がありません。
これを組み込んだものが、次です。

    factorize1n.cpp

これで、

    2 + 2x + x¹⁸

が11.6秒になりました。
この例はあまり速くなりません。もっといい例は次です。

    f(x) = 1 - x¹⁷

これは確か正17角形の作図と関係あったと記憶していますが、

    f(-7) = 29,078,814,248,401

が素数で、一つ前のプログラムだと、これの素因数分解にほとんどの時間を費やしていて、1.4秒かかってしまいます。これが今回のプログラムならこの素因数分解を避けられて、0.04秒となりました。

素因数分解が速ければ上のようなことはする必要はないのですが、これがなかなか見つかりません。
因子を見つけるアルゴリズムはたくさんあるようですが、１つρ法というのを試してみたところかなり遅くなってしまいました。上の例では効果が出てこないというのと、プログラムが悪いってことでしょう。
素数かどうかの判定は、今まで確率的な判定方法が使われてきたようですが、いくらその確率が1に極めて近いから実用上問題無いと言われても、だったらこんなプログラムなんか書いてないでどっかからいいソースを拾って来い、という話になってしまいます。
そんな中、去年の8月に多項式時間で判定できる決定的な素数判定法が発表されたそうです。論文を読んでみると、意外と簡単そうで、ちょっと組んでみようかな、という感じです。ただしすごく遅いようで、全く実用的ではないようです。けれど、こういうものは一つ突破口が見つかればどんどん速くなっていくものなのでしょう。

さて、そんな夢のある話はおいといて、お手軽に素因数分解を速くしましょう。
今までは、n の因数を求めるとき、ひたすら n^1/2 までの奇数で割っていました。もうちょっと工夫して、2と3で割り切れないなら、 5からは6k+5と6k+1で表される整数で割っていけばよいです。これだと、割とコードが増えないではっきり効果が現れます。

    factorize1o.cpp

これで11.4秒になりました。
以上を組み込んで、かつ、入力と出力をもうちょっとましにしたものを作りました。

    polynomial1.h
    ---------------------------------------------------------------------
    class Polynomial : public vector<int> {
        vector<Polynomial>    factors;
    public:
        Polynomial() { }
        Polynomial(const int *a, int size);
        Polynomial(const vector<int>& a);
        Polynomial(const string& s);
        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);
        friend bool operator <(const Polynomial& a, const Polynomial& b);
        void factorize();
        myint value(int x) const;
        void print() const;
    private:
        void _add(int a, unsigned int dim);
    };

これは次のように使います。

    polynomial1.h
    -----------------------------------------------------------------
    #include "polynomial1.h"
    
    void main(void) {
        Polynomial  p("2+2x+x^18");
        try {
            p.factorize();
        }
        catch(int e) {
            if(e == 1)
                cerr << "オーバーフローのため因数分解できず" << endl;
            else
                cerr << "メモリエラー" << endl;
        }
        p.print(); cout << endl;
    }