Groebner基底の利用(その1)

jurupapaさん( http://maxima.hatenablog.jp/ )が公開しておられる「可解な方程式を冪根で解く」プログラムの処理の一部に Groebner 基底を用いてみました.例は以下のものです.
http://maxima.hatenablog.jp/entry/2018/10/09/005917
http://maxima.hatenablog.jp/entry/2018/10/11/012019
http://maxima.hatenablog.jp/entry/2018/10/13/134605

load("SolveSolvable2.mac")$
load("grobner")$

gal_phi(xlist):=ev(sum(xlist[i]*cl[i],i,1,length(xlist)),cl:[1,-1,2,0,3,2,5])$
p1:x^4+2*x^3+3*x^2+4*x+5$

ldc(f,s):=coeff(f,s,hipow(f,s))$
tst():=(
PS:gal_init_polynomial_info(p1),
gal_minimal_polynomial_V2(PS),
gal_sol_V4(PS),
gal_galois_group2(PS),
get_matrix_permutation_group(PS),
FG:new(FiniteGroup),
gr_gen_tables(FG,PS@solperm,PS@galois_group_index),
NSGS:gr_subnormal_series(FG),
g0v:PS@gal_minpoly,
C:[],

gr1:listify(NSGS[2]@index_elements),
gr12:makelist(gr_mult(2,elem,NSGS[2]),elem,gr1),
h0:product(x-remainder(ev(V[gr1[i]],PS@vncond),g0v),i,1,length(gr1)),
h1:product(x-ev(V[gr12[i]],PS@vncond),i,1,length(gr12)),
R:remainder(rat(h0-h1)/2,g0v,V),
R:ldc(R,x),

[g1v,Da]:poly_reduced_grobner([g0v,R-a],[V,a]),
push([a,Da],C),

gr2:listify(NSGS[3]@index_elements),
gr22:makelist(gr_mult(4,elem,NSGS[3]),elem,gr2),
gr23:makelist(gr_mult(4,elem,NSGS[3]),elem,gr22),
h0:product(x-remainder(ev(V[gr2[i]],PS@vncond),g1v),i,1,length(gr2)),
h1:product(x-remainder(ev(V[gr22[i]],PS@vncond),g1v),i,1,length(gr22)),
h2:product(x-remainder(ev(V[gr23[i]],PS@vncond),g1v),i,1,length(gr23)),
R:remainder(rat(h0+w*h1+w^2*h2)/3,g1v,V),
push([w,Dw:w^2+w+1],C),
R:ef_polynomial_reduction(R,C),
R:ldc(R,x),

[Da,Dw,Db,g2v]:poly_reduced_grobner([Da,Dw,g1v,R-b],[V,b,w,a]),
push([b,Db],C),

gr3:listify(NSGS[4]@index_elements),
gr32:makelist(gr_mult(17,elem,NSGS[4]),elem,gr3),
h0:product(x-ev(V[gr3[i]],PS@vncond),i,1,length(gr3)),
h1:product(x-ev(V[gr32[i]],PS@vncond),i,1,length(gr32)),
R:remainder(rat(h0-h1)/2,g2v,V),
R:ef_polynomial_reduction(R,C),
R:ldc(R,x),

[Da,Dw,Db,Dc,g3v]:poly_reduced_grobner([Da,Dw,Db,g2v,R-c],[V,c,b,w,a]),
push([c,Dc],C),

gr4:listify(NSGS[5]@index_elements),
gr42:makelist(gr_mult(8,elem,NSGS[5]),elem,gr4),
h0:product(x-ev(V[gr4[i]],PS@vncond),i,1,length(gr4)),
h1:product(x-ev(V[gr42[i]],PS@vncond),i,1,length(gr42)),
R:remainder(rat(h0-h1)/2,g3v,V),
R:ef_polynomial_reduction(R,C),
R:ldc(R,x),

[Da,Dw,Db,Dc,g4v,Dd]:poly_reduced_grobner([Da,Dw,Db,Dc,g3v,R-d],[V,d,c,b,w,a]),
push([d,Dd],C),

RROOTS:[],
for i:1 thru 4 do (
x[i]:remainder(ev(solV[i](V),PS@solcond),g1v,V),
for gv in [g2v,g3v,g4v] do (
x[i]:remainder(x[i],gv,V)),
push((x[i]:ratexpand(ef_polynomial_reduction(x[i],C))),RROOTS)))$

showtime:all$
SolveSolvable(p1)$
tst()$[RROOTS,C];

実行結果です.

Minimal polynomial of V 
V^24+24*V^23+336*V^22+3344*V^21+25740*V^20+159984*V^19+820856*V^18
    +3519504*V^17+12721926*V^16+39075680*V^15+104485896*V^14+257189424*V^13
    +603068156*V^12+1264487184*V^11+1484791560*V^10-3707413456*V^9
    -23515353279*V^8-53513746296*V^7-7075256024*V^6+299352120960*V^5
    +770653544880*V^4+869309952000*V^3+1145273500800*V^2+1451723788800*V
    +1818528595200
  
Galois Group of x^4+2*x^3+3*x^2+4*x+5 
matrix([a,b,c,d],[a,b,d,c],[a,c,b,d],[a,c,d,b],[a,d,b,c],[a,d,c,b],[b,a,c,d],
       [b,a,d,c],[b,c,a,d],[b,c,d,a],[b,d,a,c],[b,d,c,a],[c,a,b,d],[c,a,d,b],
       [c,b,a,d],[c,b,d,a],[c,d,a,b],[c,d,b,a],[d,a,b,c],[d,a,c,b],[d,b,a,c],
       [d,b,c,a],[d,c,a,b],[d,c,b,a])
  
Subnormal series and quotients of orders 
FiniteGroup[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24] 
FiniteGroup[1,4,5,8,9,12,13,16,17,20,21,24] 
FiniteGroup[1,8,17,24] 
FiniteGroup[1,8] 
FiniteGroup[1] 
x^4+2*x^3+3*x^2+4*x+5 is solvable. 
Solutions 
'x[1] = (alpha[1]*alpha[2]^2*Zeta[3]*alpha[3]*alpha[4])/292032000
      -(3*alpha[2]^2*Zeta[3]*alpha[3]*alpha[4])/10816
      +(3*alpha[2]*Zeta[3]*alpha[3]*alpha[4])/416
      +(23*alpha[1]*alpha[2]^2*alpha[3]*alpha[4])/292032000
      -(7*alpha[2]^2*alpha[3]*alpha[4])/54080+(alpha[2]*alpha[3]*alpha[4])/104
      -(alpha[3]*alpha[4])/16-alpha[4]/4+alpha[3]/4-1/2
  
'x[2] = (-(alpha[1]*alpha[2]^2*Zeta[3]*alpha[3]*alpha[4])/292032000)
      +(3*alpha[2]^2*Zeta[3]*alpha[3]*alpha[4])/10816
      -(3*alpha[2]*Zeta[3]*alpha[3]*alpha[4])/416
      -(23*alpha[1]*alpha[2]^2*alpha[3]*alpha[4])/292032000
      +(7*alpha[2]^2*alpha[3]*alpha[4])/54080-(alpha[2]*alpha[3]*alpha[4])/104
      +(alpha[3]*alpha[4])/16+alpha[4]/4+alpha[3]/4-1/2
  
'x[3] = (-(alpha[1]*alpha[2]^2*Zeta[3]*alpha[3]*alpha[4])/292032000)
      +(3*alpha[2]^2*Zeta[3]*alpha[3]*alpha[4])/10816
      -(3*alpha[2]*Zeta[3]*alpha[3]*alpha[4])/416
      -(23*alpha[1]*alpha[2]^2*alpha[3]*alpha[4])/292032000
      +(7*alpha[2]^2*alpha[3]*alpha[4])/54080-(alpha[2]*alpha[3]*alpha[4])/104
      +(alpha[3]*alpha[4])/16-alpha[4]/4-alpha[3]/4-1/2
  
'x[4] = (alpha[1]*alpha[2]^2*Zeta[3]*alpha[3]*alpha[4])/292032000
      -(3*alpha[2]^2*Zeta[3]*alpha[3]*alpha[4])/10816
      +(3*alpha[2]*Zeta[3]*alpha[3]*alpha[4])/416
      +(23*alpha[1]*alpha[2]^2*alpha[3]*alpha[4])/292032000
      -(7*alpha[2]^2*alpha[3]*alpha[4])/54080+(alpha[2]*alpha[3]*alpha[4])/104
      -(alpha[3]*alpha[4])/16+alpha[4]/4-alpha[3]/4-1/2
  
with 
[[alpha[4],
  alpha[4]^2+(23*alpha[1]*alpha[2]^2*Zeta[3])/4563000
            -(7*alpha[2]^2*Zeta[3])/845-(2*alpha[2]*Zeta[3])/13
            +(11*alpha[1]*alpha[2]^2)/2281500+(8*alpha[2]^2)/845
            +(6*alpha[2])/13+4],
 [alpha[3],
  alpha[3]^2-(11*alpha[1]*alpha[2]^2*Zeta[3])/2281500
            -(8*alpha[2]^2*Zeta[3])/845+(8*alpha[2]*Zeta[3])/13
            +(alpha[1]*alpha[2]^2)/4563000-(3*alpha[2]^2)/169+(2*alpha[2])/13
            +4],
 [alpha[2],
  (14*alpha[1]*Zeta[3])/27-1440*Zeta[3]+alpha[2]^3-(19*alpha[1])/135-2120],
 [Zeta[3],Zeta[3]^2+Zeta[3]+1],[alpha[1],alpha[1]^2-38880000]]
  
Verification 
Numeric calcuration of the above solutions 
[(-0.8578967583284902*%i)-1.287815479557649,
 1.416093080171908*%i+0.2878154795576489,
 0.2878154795576489-1.416093080171908*%i,
 0.85789675832849*%i-1.287815479557649]
  
Numeric solutions with allroots( x^4+2*x^3+3*x^2+4*x+5 ) 
[x = 1.416093080171908*%i+0.287815479557648,
 x = 0.287815479557648-1.416093080171908*%i,
 x = 0.8578967583284904*%i-1.287815479557648,
 x = (-0.8578967583284904*%i)-1.287815479557648]
  
Evaluation took 16.7930 seconds (16.8840 elapsed) using 6984.509 MB.
(%i8) Evaluation took 6.1570 seconds (6.1670 elapsed) using 2897.350 MB.
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.
(%o9) [[(a*b^2*c*d*w)/292032000-(3*b^2*c*d*w)/10816+(3*b*c*d*w)/416
                               +(23*a*b^2*c*d)/292032000-(7*b^2*c*d)/54080
                               +(b*c*d)/104-(c*d)/16+d/4-c/4-1/2,
        (-(a*b^2*c*d*w)/292032000)+(3*b^2*c*d*w)/10816-(3*b*c*d*w)/416
                                  -(23*a*b^2*c*d)/292032000+(7*b^2*c*d)/54080
                                  -(b*c*d)/104+(c*d)/16-d/4-c/4-1/2,
        (-(a*b^2*c*d*w)/292032000)+(3*b^2*c*d*w)/10816-(3*b*c*d*w)/416
                                  -(23*a*b^2*c*d)/292032000+(7*b^2*c*d)/54080
                                  -(b*c*d)/104+(c*d)/16+d/4+c/4-1/2,
        (a*b^2*c*d*w)/292032000-(3*b^2*c*d*w)/10816+(3*b*c*d*w)/416
                               +(23*a*b^2*c*d)/292032000-(7*b^2*c*d)/54080
                               +(b*c*d)/104-(c*d)/16-d/4+c/4-1/2],
      [[d,
          (-23*a*b^2*w)+37800*b^2*w+702000*b*w-4563000*d^2-22*a*b^2-43200*b^2
                       -2106000*b-18252000],
         [c,
          22*a*b^2*w+43200*b^2*w-2808000*b*w-4563000*c^2-a*b^2+81000*b^2
                    -702000*b-18252000],
         [b,(-70*a*w)+194400*w-135*b^3+19*a+286200],[w,w^2+w+1],
         [a,a^2-38880000]]]