暗号や符号などで既約多項式がよく使われることは誰しも知っていると思う。
でも、じゃあその既約多項式ってどうやって見つけてくるの？ということはよくわからない。
そこで私が調べた結果わかったことを記事にしてみようと思った。

今回紹介するのは、Ben-Orの（多項式の）既約性判定テストの方法だ。
完全なソースは以下のURLです。

私が知っていたのは手元に以下の本を持っていたからだ。
計算機代数の辞書的な本に相当するが、数学系のプログラミングをしている人にはぜひともおすすめしたい。ネットで調べて出てくるのはアイゼンシュタインの判定法とかだが、計算機で使う判定法とはあまり関係ない。

（値段がすごいことになってますが、自分が買ったときは1万円行かなかったと思います。ノイキルヒの代数的整数論より高い？）

この本には以下のような疑似コードが書かれている。

１．ランダムに有限体上の多項式を取り、それを$f(x) \in F_q[x]$とする。
２．for i in 1 to floor(n/2) do
$g_i←gcd(x^{q^i}-x,f)$,if $g_i≠1$ then goto 1.
3.else retturn f.

たったこれだけで既約性を判定できるのだ。
コードで示そう。

//GF(2^m) then set m in this function.
int ben_or(OP f)
{
  int i, n, flg = 0;
  OP s = {0}, u = {0}, r = {0};
  vec x = {0};
  //if GF(8192) is 2^m and m==13 or if GF(4096) and m==12 if GF(16384) is testing
  int m = 13;

  x.x[1] = 1;
  s = v2o(x); //s(x)=x
  r = s;
  n = deg(o2v(f));

  if (n == 0)
    return -1;

  i = 0;

  //r(x)=x^{q^i} square pow mod
  while (i < n / 2 + 1)
  {
    flg = 1;

    //over GH(8192) 2^13
    r = opowmod(r, f, m); //r(x)=r(x)^m % f

    //over GF2
    //r=omod(opow(r,2),f);

    u = oadd(r, s);
    if (LT(u).n == 0 && LT(u).a == 0)
      return -1;
    if (LT(u).n == 0 && LT(u).a == 1)
    {
      i++;
      flg = 0;
    }
    if (LT(u).n > 0)
      u = gcd(f, u);

    if (LT(u).n > 0)
      return -1;

    if (flg == 1)
      i++;
  }

  return 0;
}

ここでちょっと問題なのは、$x^{q^i}$が２重指数なので、多項式の次数$i=n/2$が上がると簡単にオーバーフローしてしまう。これを正直に計算すると、$x^{{8192}^{64}}$なんてことになる。
仮に巨大整数型を使うにしても、xを$q^i$回計算するのは効率が悪すぎる。
しかし安心してほしい。この計算には秘訣があるのだ。

まずコードで見てみよう。

//多項式のべき乗余
OP opowmod(OP f, OP mod, int n) {
    vec v={0};
    OP ret;

     v.x[0]= 1;
    ret=v2o(v);
    while (n > 0) {
        if (n & 1) ret = omod(omul(ret , f),mod) ;  // n の最下位bitが 1 ならば x^(2^i) をかける
        f = omod(omul(f , f),mod);
        n >>= 1;  // n を1bit 左にずらす
    }
    return ret;
}

悪い例

OP opowmod(OP f, OP mod, int n)
{
  int i, j = 0;

  //繰り返し２乗法
  for (i = 1; i < n + 1; i++)
  {
    f = omul(f, f);
    if (odeg(f) > odeg(mod))
      f = omod(f, mod);
  }

  return f;
}

ここで出てくる繰り返し２乗法は楕円曲線やRSAの実装をしたことがあれば知っていると思うので割愛する。q=2ならば、f同士を掛け算することを繰り返していけば、$x^{2^i}$は$(x*x)^i$回だから、i回２乗すれば計算できる。例えば有限体がGF(8192)の場合でも、8192回ではなく$q=8192=2^{13}$なので、$m=13$回かければ終わる。

更に多項式の次数nが128次のとき、$n/2=64$回のループを繰り返す。この計算は、もしfが可約の場合、ｆを構成する既約因子の次数を上げていく事に相当する。
このように、このアルゴリズムの計算量は多項式の次数nと体の拡大次数mに関係している。
計算の途中でgcdが出てくるが、これは$x^{q^i}-x$が有限体上既約な、次数$q^{i-1}$次以下のすべての多項式の積として表されるからで、この式に含まれる既約因子とのgcdをとり、fがdeg(f)以下の共通因子を持てば可約であることがわかるからだ。この時、$i=n/2$回のループ全てにおいて$gcd(x^{q^i}-x,f)==1$（互いに素）であれば、既約と判定できる。

$x^{q^i}$の次数が巨大整数にならないのは、fの次数が計算の途中で法となる多項式modの次数がdeg(f)>deg(mod)になった時点で剰余を計算するからだ。$x^{q^i}$%modのように剰余に対して掛け算を繰り返し実行すればいいだけなのでfの次数はdeg(mod)の大きさにとどまる（のでオーバーフローしないですむ）ことがわかる。つまり、$f_i=x^{q^i}$　÷ mod、$f_0=x,q=2$とおくと、$f_i=(f_{i-1}*f_{i-1})÷ mod$ を$m=13$回、逐次計算すればいい。
（蛇足だが、q=3の場合は$f_{i}=(f_{i-1}*f_{i-1}*f_{i-1})$ ÷　modとなり、ちょっと複雑になる。やってないのでわからないが、基礎体$q$が大きなときは巨大整数のときのように$f_i$に関する演算テーブルを用意すれば効率化できるかもしれない。）

今回の実装ではランダムな多項式を生成して、その中から$GF(2^{13})$上既約になる128次多項式を40秒位で決定することが出来た。有限体を生成する際に必要とされるのは、主にGF(2)上で既約な多項式だと思うが、このアルゴリズムは$GF(2^m)$などの任意の拡大体における既約多項式の判定にも使えることに注意しよう。

以上簡単ではあるが、有限体上定義された既約多項式の判定法の紹介である。
プログラムを動かしてみたい方はコメントください。

詳しい理屈が知りたい方は、以下の参考文献をお読みいただきたい。

20210324（追記）
なんで探している間は出てこないんだろう・・・
本に乗ってることがネットで見つからないなんてありえないと思っていたけど、検索内容が変わるというのはよくあるらしい。