正規部分群の構成

今回のオリジナルプログラムでは各元の一点集合の conjugate closure から極大正規部分群を選んでいますが,それが必ずしも正しい結果を与えないことも,例えば

G:[[1,2,3,4,5,6],[1,2,3,4,6,5],[1,2,4,3,5,6],[1,2,4,3,6,5],
     [2,1,3,4,5,6],[2,1,3,4,6,5],[2,1,4,3,5,6],[2,1,4,3,6,5]];

正規部分群

[[1,2,3,4,5,6],[1,2,3,4,6,5],[1,2,4,3,5,6],[1,2,4,3,6,5]]

のように指摘されています.

一方,共役類に分割された群  G=\bigcup_{i\in A}C_i とその部分群 S について,次の2つは等価でした.
(1)SG正規部分群である.
(2) S=\bigcup_{i\in B}C_i となる空でない B \subseteq A がある.

従って,有限群 G を共役類に分割すれば,G の位数 d正規部分群の全体は,元の個数の合計が d,かつ,和が積閉である共役類の和の全体に一致します.

以下,この性質による組成列計算のコードを記します.なお,mnsg では powerset を用いましたが,入力の群のサイズが程々(笑)なら共役類も多くはないので,爆発しないと思われます(母関数を用いると,元の個数の合計が d になる組のみを抽出できますが,…).

/* 置換の積,逆元 */
mul(A,B):=map(lambda([s],A[s]),B)$
inv(A):=map(second,sort(map("[",A,sort(A)),lambda([s,t],s[1]<t[1])))$
/* 共役類分割 */
CC(GG):=block([g,C:[],G:setify(GG),H],
unless G={} do (g:listify(G)[1],
H:setify(map(lambda([s],mul(s,mul(g,inv(s)))),GG)),
G:setdifference(G,H),push(H,C)),map(listify,C))$
/* 積閉判定 */
isG(GG):=is(setify(apply('append,makelist(makelist(mul(i,j),i,GG),j,GG)))=setify(GG))$
/* 極大正規部分群 */
mnsg(GG):=block([S:[],PS:map(lambda([s],apply('append,listify(s))),listify(powerset(setify(CC(GG)))))],
for d in rest(reverse(listify(divisors(length(GG))))) while S=[] do
S:sublist(PS,lambda([s],length(s)=d and isG(s))),S[1])$
/* Galois 群の組成列 */
cs(GG):=block([S:[GG],m:GG],unless length(m)=1 do push(m:mnsg(m),S),
[S:reverse(S),S:map(length,S),makelist(S[i-1]/S[i],i,2,length(S))])$

例えば,冒頭のGに対して,cs(G); と入力すれば

[[[[1,2,3,4,5,6],[1,2,3,4,6,5],[1,2,4,3,5,6],[1,2,4,3,6,5],
          [2,1,3,4,5,6],[2,1,3,4,6,5],[2,1,4,3,5,6],[2,1,4,3,6,5]],
         [[1,2,3,4,5,6],[1,2,3,4,6,5],[1,2,4,3,5,6],[1,2,4,3,6,5]],
         [[1,2,3,4,5,6],[1,2,3,4,6,5]],[[1,2,3,4,5,6]]],[8,4,2,1],[2,2,2]]

と出力(最後の2つは各部分群の位数,隣接するそれらの比の値)されます.

次の位数9の部分群は,冒頭の例と同じく,singleton の conjugate closure にはならないものです.

[x^6+x^5+x^4+x^3-4*x^2+5,
 [[[[1,2,3,4,5,6],[1,3,2,5,4,6],[3,4,1,6,2,5],[3,1,4,2,6,5],[2,4,1,3,6,5],
    [2,1,4,6,3,5],[4,2,3,5,1,6],[4,3,2,1,5,6],[3,5,1,2,6,4],[3,1,5,6,2,4],
    [2,5,1,6,3,4],[2,1,5,3,6,4],[5,2,3,1,4,6],[5,3,2,4,1,6],[3,5,4,6,2,1],
    [3,4,5,2,6,1],[2,5,4,3,6,1],[2,4,5,6,3,1],[1,6,3,5,4,2],[1,3,6,4,5,2],
    [1,6,2,4,5,3],[1,2,6,5,4,3],[6,4,1,2,3,5],[6,1,4,3,2,5],[4,6,3,1,5,2],
    [4,3,6,5,1,2],[4,6,2,5,1,3],[4,2,6,1,5,3],[6,5,1,3,2,4],[6,1,5,2,3,4],
    [5,6,3,4,1,2],[5,3,6,1,4,2],[5,6,2,1,4,3],[5,2,6,4,1,3],[6,5,4,2,3,1],
    [6,4,5,3,2,1]],
   [[1,2,3,4,5,6],[1,2,6,5,4,3],[1,3,2,5,4,6],[1,6,3,5,4,2],[4,2,6,1,5,3],
    [4,3,2,1,5,6],[4,6,3,1,5,2],[5,2,6,4,1,3],[5,3,2,4,1,6],[5,6,3,4,1,2],
    [1,3,6,4,5,2],[1,6,2,4,5,3],[4,2,3,5,1,6],[5,2,3,1,4,6],[4,3,6,5,1,2],
    [4,6,2,5,1,3],[5,3,6,1,4,2],[5,6,2,1,4,3]],
   [[1,2,3,4,5,6],[1,3,6,4,5,2],[1,6,2,4,5,3],[4,2,3,5,1,6],[5,2,3,1,4,6],
    [4,3,6,5,1,2],[5,6,2,1,4,3],[4,6,2,5,1,3],[5,3,6,1,4,2]],
   [[1,2,3,4,5,6],[1,3,6,4,5,2],[1,6,2,4,5,3]],[[1,2,3,4,5,6]]],[36,18,9,3,1],
  [2,2,3,3]]]

冪乗に分解できるタイプ

今回のオリジナルプログラムでは,例えば,SolveSolvable(x^2-2)$ の根の表示には [[alpha[1],alpha[1]^2-8]] が現れるので,SolveSolvable(x^2-8)$ と問うと,[[alpha[1],alpha[1]^2-32]] が現れるので,...となってしまいます.

プログラムの趣意に反する可能性もありますが,例えば,x^2-2,x^6+x^3+1,x^20+x^15+3*x^5+4 といった冪乗に分解できるタイプについては,それぞれ x-2,x^2+x+1,x^4+x^3+3*x+4 の根の表示を得た上で,その各根の2,3,5乗根のすべてを出力するコードを書いてみました.

まず,1の原始N乗根の最小多項式を出力する関数を用意します.

ROU(x,N):=block([D:reverse(listify(divisors(N)))],
 num(factor((x^pop(D)-1)/apply("*",map(lambda([s],x^s-1),D)))))$
/* テスト */
makelist(ROU(x,i),i,1,10);

次に SolveSolvable の定義の局所変数の宣言をコメントアウトした SolveSolvable2.mac をリロードして,以下を実行すると,上記の分解を挟んだ結果が得られます.

/* サンプル選択 */
p:[x^2-2,x^6+x^3+1,x^20+x^15+3*x^5+4][2];
/* 分解(polydecompの出力は気まぐれなので...) */
for i:(N:hipow(p,x)) step -1 unless polynomialp(q:rat(subst(x=x^(1/i),p)),[x]) do N:i-1$[q,x^N];
/* ベースの方程式の求根 */
if hipow(q,x)=1 then [C,SolN,x[1]]:[[],1,rhs(solve(q,x)[1])] else (SolveSolvable(q),C:SI[StageN]@ExtensionList)$
/* 各根x[i]のN乗根y[i]に1の原始N乗根の冪乗を掛ける(得られるもの全体はy[i]の選び方に拠らない) */
for i:1 thru SolN do C:push([y[i],y[i]^N-x[i]],C)$push([w,ROU(w,N)],C);
RROOTS:flatten(makelist(makelist(y[i]*w^j,j,0,N-1),i,1,SolN));
/* 検算 */
ef_polynomial_reduction(P:apply("*",map(lambda([s],x-s),RROOTS)),C);

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]]]

CAD-QE on PARI+Singular(実行例)

同じ例を開発中の PARI+Singular 上での CAD-QE で処理すると...

? tst11Sv2s01([ex,ex],[a,b,x,y],and,[sgn np (x^2+y^2-1),sgn 0 (a*x+b*y-1)],2*7);Ans();
  *** connecting adjacent 16/796 cells.
1[a <= [a+1,1],true,true,true]
2[[a+1,1] < a < [a,1],b <= [a^2+(b^2-1),1] or [a^2+(b^2-1),2] <= b,true,true]
3[a = [a,1],b <= [b^2-1,1] or [b^2-1,2] <= b,true,true]
4[[a,1] < a < [a-1,1],b <= [a^2+(b^2-1),1] or [a^2+(b^2-1),2] <= b,true,true]
5[[a-1,1] <= a,true,true,true]
time = 107 ms.
? tst11Sv2s01([ex,ex],[a,b,x,y],and,[sgn np (x^4+y^3-1),sgn 0 (a*x+b*y-1)],2*7);Ans();
12 no real. 9 12 10 12 9 12 12 9 12 10 12 9 12 no real. 
  *** polrootsreal: Warning: increasing stack size to 16000000.
  *** connecting adjacent 71/1817 cells.
1[a <= [a^8+6*a^4-3,1],true,true,true]
2[[a^8+6*a^4-3,1] < a < [a,1],b <= [27*a^12+(-432*b^3-54)*a^8+(1728*b^6+432*b^3+27)*a^4+(256*b^12-768*b^9+768*b^6-256*b^3),1] or [27*a^12+(-432*b^3-54)*a^8+(1728*b^6+432*b^3+27)*a^4+(256*b^12-768*b^9+768*b^6-256*b^3),2] <= b,true,true]
3[a = [a,1],b < [b^4-b,1] or [b^4-b,2] <= b,true,true]
4[[a,1] < a < [a^8+6*a^4-3,2],b <= [27*a^12+(-432*b^3-54)*a^8+(1728*b^6+432*b^3+27)*a^4+(256*b^12-768*b^9+768*b^6-256*b^3),1] or [27*a^12+(-432*b^3-54)*a^8+(1728*b^6+432*b^3+27)*a^4+(256*b^12-768*b^9+768*b^6-256*b^3),2] <= b,true,true]
5[[a^8+6*a^4-3,2] <= a,true,true,true]
time = 9,464 ms.

QEPCAD Version B 1.72

https://www.usna.edu/Users/cs/wcbrown/qepcad/B/QEPCAD.html
上記の通り,更新されていましたので,インストールスクリプトをUPします.

export CAS=~/CAS
mkdir $CAS
export qe=$CAS/qesource
export saclib=$CAS/saclib2.2.7
cd $CAS
wget https://www.usna.edu/Users/cs/wcbrown/qepcad/INSTALL/saclib2.2.7.tar.gz
tar zxvf ./saclib2.2.7.tar.gz -C ./
cd $saclib/bin
./sconf
./mkproto
./mkmake
./mklib all
cd $CAS
wget https://www.usna.edu/Users/cs/wcbrown/qepcad/INSTALL/qepcad-B.1.72.tar.gz
tar zxvf ./qepcad-B.1.72.tar.gz -C ./
cd $qe
sed -i "s/csh/sh/g" $qe/Makefile
make
sed -i "s/\#SINGULAR/SINGULAR/g" $qe/default.qepcadrc
cd $CAS
wget https://www.usna.edu/Users/cs/wcbrown/qepcad/SLFQ/simplify-1.20.tar.gz
tar zxvf ./simplify-1.20.tar.gz -C ./
cd ./simplify
make
sudo ln -s $CAS/qesource/bin/qepcad /usr/local/bin/qepcad
sudo ln -s $CAS/simplify/slfq /usr/local/bin/slfq

実行例(1)

echo -e "[] \n (a,b,x,y) \n 2 \n (Ex)(Ey)[x^2+y^2<=1 /\ a x+b y=1]. \n finish" | qepcad
=======================================================
                Quantifier Elimination                 
                          in                           
            Elementary Algebra and Geometry            
                          by                           
      Partial Cylindrical Algebraic Decomposition      
                                                       
      Version B 1.72, Sun Mar 18 22:49:08 EDT 2018
                                                       
                          by                           
                       Hoon Hong                       
                  (hhong@math.ncsu.edu)                
                                                       
With contributions by: Christopher W. Brown, George E. 
Collins, Mark J. Encarnacion, Jeremy R. Johnson        
Werner Krandick, Richard Liska, Scott McCallum,        
Nicolas Robidoux, and Stanly Steinberg                 
=======================================================
Enter an informal description  between '[' and ']':
[]Enter a variable list:
 (a,b,x,y)Enter the number of free variables:
 2 Enter a prenex formula:
 (Ex)(Ey)[x^2+y^2<=1 /\ a x+b y=1].

=======================================================

Before Normalization >
 finish

An equivalent quantifier-free formula:

b^2 + a^2 - 1 >= 0


=====================  The End  =======================

-----------------------------------------------------------------------------
0 Garbage collections, 0 Cells and 0 Arrays reclaimed, in 0 milliseconds.
334662 Cells in AVAIL, 500000 Cells in SPACE.

System time: 31 milliseconds.
System time after the initialization: 17 milliseconds.
-----------------------------------------------------------------------------

実行例(2)

echo -e "[] \n (a,b,x,y) \n 2 \n (Ex)(Ey)[x^4+y^3<=1 /\ a x+b y=1]. \n finish" | qepcad +L100000 +N10000000
=======================================================
                Quantifier Elimination                 
                          in                           
            Elementary Algebra and Geometry            
                          by                           
      Partial Cylindrical Algebraic Decomposition      
                                                       
      Version B 1.72, Sun Mar 18 22:49:08 EDT 2018
                                                       
                          by                           
                       Hoon Hong                       
                  (hhong@math.ncsu.edu)                
                                                       
With contributions by: Christopher W. Brown, George E. 
Collins, Mark J. Encarnacion, Jeremy R. Johnson        
Werner Krandick, Richard Liska, Scott McCallum,        
Nicolas Robidoux, and Stanly Steinberg                 
=======================================================
Enter an informal description  between '[' and ']':
[]Enter a variable list:
 (a,b,x,y)Enter the number of free variables:
 2 Enter a prenex formula:
 (Ex)(Ey)[x^4+y^3<=1 /\ a x+b y=1].

=======================================================

Before Normalization >
 finish

An equivalent quantifier-free formula:

a^8 + 6 a^4 - 3 > 0 \/ 256 b^12 - 768 b^9 + 1728 a^4 b^6 + 768 b^6 - 432 a^8 b^3 + 432 a^4 b^3 - 256 b^3 + 27 a^12 - 54 a^8 + 27 a^4 > 0 \/ [ b > 0 /\ 256 b^12 - 768 b^9 + 1728 a^4 b^6 + 768 b^6 - 432 a^8 b^3 + 432 a^4 b^3 - 256 b^3 + 27 a^12 - 54 a^8 + 27 a^4 = 0 ]


=====================  The End  =======================

-----------------------------------------------------------------------------
638 Garbage collections, -1442980436 Cells and 12 Arrays reclaimed, in 25548 milliseconds.
881113 Cells in AVAIL, 5000000 Cells in SPACE.

System time: 205804 milliseconds.
System time after the initialization: 205697 milliseconds.
-----------------------------------------------------------------------------

CAD-QE on Risa/Asir(特徴)

CAD の要は,projection と代数体上での根の分離です.

昨年の maxima における実装は,Lazard method と複素数値根のクラスタリングのための pscs の計算(Mathematica 似)との併用というスタイルであり,それは primitive elements の計算(squarefree norm や代数体上の多項式の GCD の計算)が maxima では非常に重いことによる選択でしたが,肝心の多倍長浮動小数点数の計算も maxima では高速とは言えず,求根には外部の MPSolve を用いていました.

そうした経緯から,Lazard method の軽さを損なわない,primitive elements 生成の実装を探した結果,PARI+Singular での実装(後日公開予定,かなり高速)に至る訳ですが,当初,多変数多項式因数分解(PARI での実装は実験的かつ遅い)などを含むルーチンには,Risa を利用しており,今回公開したコードは,それと Risa の特長の一つである逐次拡大体上の各種関数とを組み合わせた,Risa/Asir(PARI ライブラリも含む)プロパーな CAD-QE 関数です(ので,速度は関心の外ということに).

CAD-QE on Risa/Asir(使用例)

先のソースコードを cadqe.rr として保存,Asir を起動,それをロードします.

asir -quiet
load("cadqe.rr")$

幾つかのメッセージが表示された後,入力待ちになるので

help("cadqe")$

と入力すると,次のような使用例が表示されます.

cadqe([ex],[a,b,c,x],id,[a*x^2+b*x+c@==0])$
cadqe([all],[a,b,x],imp,[x^2-1@<=0,x^3+a*x+b@<=0])$
cadqe([all,all],[a,b,c,y,x],eq,[a*x^2+b*y^2@<=c,x^2+y^2@<=1])$
def _(A,B,C){and(A,and(B,C));}$cadqe([ex,ex,ex],[a,b,c,x,y,z],_,[x+y+z@==a,x*y+y*z+z*x@==b,x*y*z@==c]|grob=1)$

書式は

cadqe(量化子のリスト,自由変数と束縛変数とのリスト,真理関数または恒等関数id,原子論理式のリスト,オプション)$

ここで,量化子∈{ex,all},定義済みの真理関数∈{not,and,or,imp,eq},述語記号∈{@<,@<=,@>,@>=,@==,(@=,)@!=} であり,量化子の個数が n(>0) の場合,変数リストの最後の n 個が各量化子に対する束縛変数です.

例えば

cadqe([ex],[a,b,c,x],id,[a*x^2+b*x+c@==0])$

と入力すると

1[[-oo,[a,1]],[-oo,oo],[[ 4*c*a-b^2 1 ],oo]]
2[[[ a 1 ],[ a 1 ]],[-oo,[b,1]],[-oo,oo]]
3[[[ a 1 ],[ a 1 ]],[[ b 1 ],[ b 1 ]],[[ c 1 ],[ c 1 ]]]
4[[[ a 1 ],[ a 1 ]],[[b,1],oo],[-oo,oo]]
5[[[a,1],oo],[-oo,oo],[-oo,[ 4*c*a-b^2 1 ]]]
0.06575sec + gc : 0.04272sec(0.1256sec)

のように,5個のリストが表示され,各リストはその成分が表す論理式の連言を表し,それらの選言が結果となります.ここで,1個目のリストの第1成分 [-oo,[a,1] ] のリスト [a,1] は a についての多項式 a の小さい方から数えて1個目の実根(つまり,0)の開端点を表し,[-oo,[a,1] ] は -oo

def _(A,B,C){and(A,and(B,C));}$
cadqe([ex,ex,ex],[a,b,c,x,y,z],_,[x+y+z@==a,x*y+y*z+z*x@==b,x*y*z@==c]|grob=1)$
1[[-oo,oo],[-oo,[a^2-3*b,1]],[[ 4*c*a^3-b^2*a^2-18*c*b*a+4*b^3+27*c^2 1 ],[ 4*c*a^3-b^2*a^2-18*c*b*a+4*b^3+27*c^2 2 ]]]
2[[-oo,oo],[[ a^2-3*b 1 ],[ a^2-3*b 1 ]],[[ a^3-27*c 1 ],[ a^3-27*c 1 ]]]
0.2741sec + gc : 0.2713sec(0.5929sec)

直前の projection set は

PP0;

と入力すると

[[],[a^2-3*b],[4*c*a^3-b^2*a^2-18*c*b*a+4*b^3+27*c^2],[x^3-a*x^2+b*x-c,3*x^2-2*a*x-a^2+4*b],[x^2+(y-a)*x+y^2-a*y+b],[x+y+z-a]]
0sec(3.099e-06sec)

のように確認できます(grob=1 オプションなしの場合,この例の projection factors の個数は 200 を超えます).