2023年12月23日土曜日

231223

SINGULAR


n次多項式の判別式の項数について

この記事は
日曜数学 Advent Calendar 2023
の12/23 分として書いております。

2013年にn=17について計算されて以降、新しい情報はないようですが、
現段階で、普通のPCでどの程度計算可能か確認してみた。

一般にn次多項式f(x)(=a_0*x^n + a_1*x^(n-1) + ... +a_n)の判別式は、
(-1)^(n*(n-1)/2) * 1/a_0 * resultant(f(x), f'(x))
と表される。
判別式の項数を考えているので、a_0=1としても良い。

様々な計算方法があるが、
http://syskiso.fuee.u-fukui.ac.jp/~kkimur/FINAL.pdf
のP.20-21のアルゴリズムを使うことにする。

このアルゴリズムを使うと、
x^3 + a_1*x^2 + a_2*x + a_3 = 0 の判別式は次のように計算できる。
① x^2 + a_1*x + a_2 = 0 の判別式はa_1^2 - 4*a_2.
② D0 = a_2^2*(a_1^2 - 4*a_2) = a_1^2*a_2^2 - 4*a_2^3.
③ D1 = -4*a_1^3 + 18*a_1*a_2, D2 = -27.
④ D = D0 + D1*a_3 + D_2*a_3^2.

SINGULARを使うと、n=14までは現実的な時間で算出することができる。

ring r=(0),(a14,a13,a12,a11,a10,a9,a8,a7,a6,a5,a4,a3,a2,a1),lp;
list v=1,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14;
poly s=v[2]^2-4*v[3]; // 2次の場合
poly t0;
poly u0;
int n,j,k,w;
for(n=3; n<15; n++){
  t0=v[n]^2*s; // D0=a{n-1}^2*disc.(f{n-1}(x))
  s=t0;
  for(k=1; k<n; k++){
    u0=n*diff(t0,v[2]);
    for(j=2; j<n; j++){
      u0=u0+(n+1-j)*v[j]*diff(t0,v[j+1]);
    }
    t0=-u0/(k*v[n]);
    s=s+t0*v[n+1]^k; // D=D0+D1*an+D2*an^2+...+D{n-1}*an^(n-1)
  }
  print(size(s));
}

quit;

出力結果
5
16
59
246
1103
5247
26059
133881
706799
3815311
20979619
117178725
Auf Wiedersehen.

0 件のコメント:

コメントを投稿

注: コメントを投稿できるのは、このブログのメンバーだけです。