ATPとCASのこと

2019-01-23

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

2019-01-07

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

2018-03-25

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 関数です(ので,速度は関心の外ということに).

2018-03-24

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<a<0 を表します.そして,第2成分 [-oo,oo] は -oo<b<+oo,第3成分 [[ 4*c*a-b^2 1 ],oo] のベクトル [ 4*c*a-b^2 1 ] は第1,2成分の a,b に対して,c についての多項式 4*c*a-b^2 の小さい方から数えて1個目の実根(つまり,(b^2)/(4*a) )の閉端点を表し,[ [ 4*c*a-b^2 1 ],oo] は (b^2)/(4*a)<=c<+oo を表します.同様に,2個目のリストの第1成分 [ [a 1],[a 1] ] は 0<=a<=0(つまり,a=0)を表します.

なお,引数の個数が3以上の真理関数は定義済みの真理関数を組み合わせて定義します.また,真理関数が連言,原子論理式が等式の場合,オプション grob=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)$
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 を超えます).

CAD-QE on Risa/Asir(ソースコード)

Risa/Asir 上に CAD ベースの RCF-QE を実装しましたので,ソースコードと使用例を公開します.

環境は,Core i5-3210M 2.50GHz/8GB/Ubuntu 14.04,Risa/Asir Version 20170917 (Kobe Distribution),OpenXM/Risa/Asir-Contrib $Revision: 1.143 $ (20170210) です.

/*%% 18.03.24.Sat. 20:33:27*/
load("sp")$load("gr")$load("sturm")$load("os_muldif.rr")$
def deld(L) " delete dupl. " {if(L==[]){return [];}
 else {M=[C=car(L=qsort(L))];L=cdr(L);
 while(L!=[]) {
  if(C!=(C=car(L))) {L=cdr(L);M=cons(C,M);} else L=cdr(L);}
 return M;}
}
def delc(L) " without os_md " {A=deld(L);B=[];
 while(A!=[]) {if(type(F=car(A))==2) B=cons(F,B);A=cdr(A);}
 return B;}
def factors(L) { return L==[] ? [] : delc(taka_base_flatten(map(fctr,delc(L)))); }
def proj(VV,L) " without os_md "
{W=reverse(cdr(VV));A=L;PP=[];
 while(W!=[]) {V=car(W);W=cdr(W);C=factors(A);A=B=[]; /*print(V);*/
  while(C!=[]) {F=car(C);C=cdr(C);
   if((DG=deg(F,V))>0)
    {if(vars(F)!=[V] || count_real_roots(F)>0) {A=cons(F,A);}
     B=cons(coef(F,DG,V),B);
     B=cons(res(V,F,diff(F,V)),B);
     CL=[];for(K=0;K<=DG;K++) {if((CF=coef(F,K,V))!=0) CL=cons(CF,CL);}
     if((GR=nd_f4(CL,vars(CL),0,0))!=[1] && GR!=[-1])
      B=cons(coef(F,0,V),B);}
   else B=cons(F,B);}
  PP=cons(A,PP);
  for(I=0;I<(LA=length(A))-1;I++) for(J=I+1;J<LA;J++) 
   B=cons(res(V,A[I],A[J]),B);
  A=B;}
 A=factors(A);B=[];V=car(VV);while(A!=[]) {F=car(A);A=cdr(A);
  if(count_real_roots(F)>0) B=cons(F,B);}
 return cons(B,PP);
}
def prod(L) {P=1;while(L!=[]){P*=car(L);L=cdr(L);} return P;}
def adp(L)
{R=car(L)[1];/*R=rdp(L); maybe long */L=cdr(L);
 while(L!=[]){R=res((C=car(L))[0],C[1],R);L=cdr(L);}
 return polsqf1(R);}
def polrootsreal(F,P) {R=pari(roots,F,P);
 if(type(R)!=5) return [R];
 else
  {R=vtol(R);RR=[];
   while(R!=[]){ if(imag(C=car(R))==0) {RR=cons(C,RR);} R=cdr(R); }
 return reverse(RR);}
}
def setG(N) { pari(allocatemem,10^9*N); }$
setG(1)$
/* 1G */
ctrl("bigfloat",1)$
extern REALPRECISION$
REALPRECISION=50$setprec(REALPRECISION)$
def polrootsreal4(P) {A=polrootsreal(P,LP=REALPRECISION);
 if((NA=length(A))!=0)
  {if(NA==1) R=+oo;
   else while((R=(os_md.lmin(ltov(cdr(A))-ltov(reverse(cdr(reverse(A)))))*6/25))*10^LP<=os_md.lmax(map(number_abs,A))) {A=polrootsreal(P,LP=2*LP);/*print(A);*/}
   return([A,R]);}
 else return [];
}
def toint(L) {R=L[1];return base_makelist([s-R,s+R],s,L[0],0);}
def toint2(VF,L) {V=VF[0];F=VF[1];R=L[1];
 return base_makelist([subst(F,V,s-R),subst(F,V,s+R)],s,L[0],0);}
def toint4(VF,L) {V=VF[0];F=VF[1];M=L[0];R=(L[1]==oo)?1:L[1];Z=[];
 while(M!=[]){S=car(M);M=cdr(M);
  if(subst(F,V,S-R)*subst(F,V,S+R)<0) Z=cons(S,Z);}
 return Z;
}
extern CAD,CAD0;
def ans() { if(CAD==[]){ print("false");}
 else {for(K=0;K<length(CAD);K++){printf(rtostr(1+K));print(CAD[K]);}}
}
def lambda1(X) { return car(X); }
extern PP0,PPj,Vj,FF,AFS,US,UN,UA,QQj,EXS,ALLS,TC,FC,TCS,TFUNC;
def id(F) {F;}
def not(F){ 1-F; }
def or(F,G){ 1-(1-F)*(1-G); } /* F || G */
def and(F,G){ F * G; } /* F && G */
def imp(F,G){ 1+F*(G-1); }
def eq(F,G){ (1+F*(G-1))*(1+G*(F-1)); }
def tfunc(F,G){ F * G; }
def tuf(N,S,A) {L=AFS;M=[];SUS=cons([Vj,S],US);NUN=cons(Vj,cons(N,UN));
 while(L!=[]){F=car(car(L));G=car(L)[1];L=cdr(L);  /*print([F,SUS]);*/
  F=call(subst,cons(F2=simpalg_noalg(F,SUS),NUN));  /*print("->"+rtostr(F2));print("->"+rtostr(F));*/
  if(type(F)==2) M=cons(1/2,M);
  else M=cons(eval_str(rtostr(F)+G),M);
}
 FC=TC=0;
 V=call(TFUNC,reverse(M));
 if(V!=0 && V!=1) {return [NUN,SUS,cons(A,UA)];}
 else {if(V==0) {FC=1;return 0;} else {TC=1;TCS++;
  return [cons(A,UA)];}}
}
def addinf(U) {
 return (length(U)==1)?[cons([-oo,oo],U[0])]
        :[cons(Vj,cons(0,U[0])),cons([Vj,Vj],U[1]),cons([-oo,oo],U[2])];
}
def lambda3(L) { return reverse(car(L)); }
def makeANS1(ANS0) {
if(ANS0==[[]]){
TUF=tuf(0,Vj,[-oo,oo]);
    if(TUF) {ANS1=[TUF];}}
else {
   ANS1=[];LAST=[car(car(ANS0))-2,0,-oo];
   while(ANS0!=[])
   {Last=car(LAST);TOP=car(ANS0);Top=car(TOP);DDD=2;
    do {DDD*=2;Mid=floor((Last+Top)/2*DDD)/DDD; /*print(DDD);*/ }
     while(Mid<=Last || Top<=Mid);
TUF=tuf(Mid,ptozp(Vj-Mid),[LAST[2],TOP[2]]);
if(EXS*TC) {return [TUF];}
if(ALLS*FC) {return [];}
    if(TUF) {ANS1=cons(TUF,ANS1);}
TUF=tuf(Top,TOP[1],map(vect,[TOP[2],TOP[2]]));
if(EXS*TC) {return [TUF];}
if(ALLS*FC) {return [];}
    if(TUF) {ANS1=cons(TUF,ANS1);}
    LAST=TOP;ANS0=cdr(ANS0);}
TUF=tuf(Etr=ceil(Top+1),ptozp(Vj-Etr),[LAST[2],oo]);
if(EXS*TC) {return [TUF];}
if(ALLS*FC) {return [];}
   if(TUF) {ANS1=cons(TUF,ANS1);}}
return ANS1;  }
def deltopn(N,L) { for(K=1;K<=N;K++){L=cdr(L);} return L; }
def consinfn(N,L) { for(K=1;K<=N;K++){L=cons([-oo,oo],L);} return L; }
def delnth(L,N)
{ M=[];
for(K=0;K<N;K++){M=cons(car(L),M);L=cdr(L);}
L=cdr(L);
for(K=0;K<N;K++){L=cons(car(M),L);M=cdr(M);}
return L; }
def repnth(L,N,A)
{ M=[];
for(K=0;K<N;K++){M=cons(car(L),M);L=cdr(L);}
L=cdr(L);L=cons(A,L);
for(K=0;K<N;K++){L=cons(car(M),L);M=cdr(M);}
return L; }
def cnn(LenVV,ANS2)
{
LAST0=base_makelist([[0,0],[0,0]],k,1,LenVV);
for(K=0;K<LenVV;K++){
ANS3=[];LAST=LAST0;
while(ANS2!=[]){TOP=car(ANS2);ANS2=cdr(ANS2);
if(delnth(LAST,K)!=delnth(TOP,K)){ANS3=cons(TOP,ANS3);LAST=TOP;}
else
 {LASTK=LAST[K];TOPK=TOP[K];
  if(ltov(LASTK[1])==TOPK[0] || LASTK[1]==ltov(TOPK[0]))
   {LAST=repnth(LAST,K,[LASTK[0],TOPK[1]]);
    ANS3=cdr(ANS3);ANS3=cons(LAST,ANS3);}
  else{ANS3=cons(TOP,ANS3);LAST=TOP;}
 }}
ANS2=reverse(ANS3);
}
return ANS2;
}
def lambda4(L) { C=car(L);
if(C==8){C="==0";}else{
if(C==14){C="!=0";}else{
if(C==9){C="<=0";}else{
if(C==11){C="<0";}else{
if(C==13){C=">=0";}else{
if(C==15){C=">0";}}}}}}
return [L[1],C];
}
def at2afs(AT) { return map(lambda4,map(fopargs,AT)); }
extern RVV,TOP;
def lambda5(F) { TpF=type(F);
 if(TpF==2){ return F; }
 if(TpF==4){ return [prod(flatmf(fctr(gr([F[0],TOP],RVV,2)[1]))),F[1]]; }
 if(TpF==5){ return ltov([prod(flatmf(fctr(gr([F[0],TOP],RVV,2)[1]))),F[1]]); }
}
extern LenQQ;
def lambda6(F) { for(K=1;K<=LenQQ;K++){F=cdr(F);} return F; }
def ppzeros(V,PP/*must be !=[]*/,L1,L2)
 {
 M=[];while(PP!=[]){M=append(M,zeros(V,car(PP),L1,L2));PP=cdr(PP);}
 M=qsort(M,sort2);
 N=[LAST=car(M)];M=cdr(M);
 while(M!=[]){
  if(LAST[0]!=((TOP=car(M))[0])) {N=cons(TOP,N);LAST=TOP;}
/*  else {if(lambda2(LAST[1])!=lambda2(TOP[1])) {*/
/*  else {if(subst(LAST[1],U0)!=subst(TOP[1],U0)) { */
  else {if(ptozp(simpalg_noalg(ptozp( LAST[1] ),L1))
         !=ptozp(simpalg_noalg(ptozp( TOP[1] ),L1))) {
/*  else {if(LAST[1]!=TOP[1]) {*/
   print("  *** common-roots fault in ppzeros:");print([LAST,TOP,US,UN]);
   print(" *** Increase the value of REALPRECISION.");}}
  M=cdr(M);
 }
/* print(N); */
 return N;
}
def cadqe(QQ,VV,Tfunc,AF)
"
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)$
"
{
JJ0=(LenVV=length(VV))-(LenQQ=length(QQ));EXS=ALLS=0;
if(type(car(AF))==14){AF=at2afs(AF);}
TFUNC=Tfunc;
/*if(length(AF)==1){TFUNC=id;}else{TFUNC=tfunc;}*/
FF=map(lambda1,AF);
if(type(getopt(grob))!=-1){FF=gr(FF,reverse(VV),2);}
PP=proj(VV,FF);ANS2=[[[],[],[]]];PP0=PP;AFS=AF;
for(J=0;J<LenVV;J++){Vj=VV[J];PPj=PP[J];
if(J==LenVV-1){PPZ=ppzeros;}else{PPZ=ppzeros;}
if((JJ=J-JJ0)>=0){QQj=QQ[JJ];
 if(QQj==ex){EXS++;ALLS=0;}else{EXS=0;ALLS++;}; /*print([QQj,EXS,ALLS]);*/ }
if(PPj==[]){ANS3=[];while(ANS2!=[])
 {/*print(ANS2);*/U=car(ANS2);ANS2=cdr(ANS2);
  if(length(U)==1){TUF=[cons([-oo,oo],UA=car(U))]; /*print(TUF);*/}
   else{US=U[1];UN=U[0];UA=U[2];TUF=tuf(0,Vj,[-oo,oo]);}
  if(TUF) {ANS3=cons(TUF,ANS3);}}
 ANS2=reverse(ANS3);}
else{ANS3=[];while(ANS2!=[])
 {U=car(ANS2);ANS2=cdr(ANS2);ETC=TCA=TCS=0;
  if(length(U)==1){ANS1=[[cons([-oo,oo],UA=car(U))]]; /*print(ANS1);*/}
  else {ANS0=( *PPZ)(Vj,PPj,US=U[1],UN=U[0]);UA=U[2]; /*if(Vj==x){print(ANS0);};*/
        ANS1FULL=2*length(ANS0)+1;
        ANS1=makeANS1(ANS0); /*print([CNT,ANS1]);*/
/*print([EXS,TC,TCS,ANS1FULL,ALLS,FC,J]);*/
  if((ETC=EXS*TC) || ALLS*FC) {
   UA=deltopn(EAS=(EXS+ALLS-1),UA); /*print(UA);*/
/*print([ETC,TCA,TCS,UA,EAS]);print(ANS2);*/
   while(deltopn(EAS,car(reverse(car(ANS2))))==UA){ANS2=cdr(ANS2);}
/*print(ANS2);print(" ");*/
   while(deltopn(EAS+1,car(reverse(car(car(ANS3)))))==UA){ANS3=cdr(ANS3);}}
  if(ETC /*|| TCS==ANS1FULL*/){ANS1=[[consinfn(EAS+1,UA)]];}
   else{if(TCS==ANS1FULL){ANS1=[[cons([-oo,oo],UA)]];}} }
  if(ANS1){ANS3=cons(reverse(ANS1),ANS3);}; /*print(["ANS3",ANS3]);print(["ANS1",ANS1]);*/ }
 ANS2=os_md.m2l(reverse(ANS3)|flat=1); /*print([J,ANS2]);*/
    }
}
if(ANS2==[]){ return CAD=[];}
CAD0=ANS2=map(lambda1,ANS2);
do {ANS3=ANS2;ANS2=cnn(LenVV,ANS2);} while (ANS3!=ANS2);
for(K=JJ;0<=K;K--){
 if(QQ[K]==all){ANS3=[];
  while(ANS2!=[]){TOP=car(ANS2);ANS2=cdr(ANS2);
   if(TOP[JJ-K]==[-oo,oo]){ANS3=cons(TOP,ANS3);}}
  ANS2=reverse(ANS3);
  do {ANS3=ANS2;ANS2=cnn(LenVV,ANS2);} while (ANS3!=ANS2);
 }
 if(QQ[K]==ex){ANS3=[LAST=repnth(car(ANS2),JJ-K,[-oo,oo])];
  while(ANS2!=[]){TOP=repnth(car(ANS2),JJ-K,[-oo,oo]);ANS2=cdr(ANS2);
   if(LAST!=TOP){ANS3=cons(TOP,ANS3);LAST=TOP;}}
  ANS2=reverse(ANS3);
  do {ANS3=ANS2;ANS2=cnn(LenVV,ANS2);} while (ANS3!=ANS2);
 }
}
LenVV=JJ0;
ANS2=map(reverse,map(lambda6,ANS2));
ANS2=call(mat,ANS2);RVV=reverse(VV);
for(I=0;I<size(ANS2)[0];I++){
 for(J=0;J<LenVV-1;J++){Aij=ANS2[I][J];
  if((Aij0=Aij[0])==Aij[1]){TOP=Aij0[0];
   for(K=J+1;K<LenVV;K++){ANS2[I][K]=map(lambda5,ANS2[I][K]);}
  }
 }
}
CAD=ANS2=os_md.m2ll(ANS2);
ans();
return CAD;
}
def ansm() { print(call(mat,CAD)); }
def af_noalg(F,DL)
{
	DL = reverse(DL);
	N = length(DL);
	Tab = newvect(N);
	/* Tab = [[a1,r1],...]; ri is a root of di(t,r(i-1),...,r1). */
	AL = [];
	for ( I = 0; I < N; I++ ) {
		T = DL[I];
		for ( J = 0, DP = T[1]; J < I; J++ )
			DP = subst(DP,Tab[J][0],Tab[J][1]);
		B = newalg(DP);
		Tab[I] = [T[0],B];
		F = subst(F,T[0],B);
		AL = cons(B,AL);
	}
	FL = af(F,deld(AL)); /* modified */
	for ( T = FL, R = []; T != []; T = cdr(T) )
		R = cons([conv_noalg(T[0][0],Tab),T[0][1]],R);
	return reverse(R);
}
def sort2(X,Y) { return (X[0]<Y[0])?1:0; }
def lambda2(F) { return ptozp(simpalg_noalg(ptozp(F),US)); }
def sort1(X,Y) { return (X[0]<Y[0])?1:0; }
def zeros(V,F,L1,L2)
 {
 G=simpalg_noalg(F,L1); /*print(G);*/
 if(G==0) F=G=subst3(F,L1); else 0; /*print([V,F,L1,L2,G]);*/
 L=flatmf(af_noalg(G,L1));
 if(L==[1]) return [];
 W=[];
 while(L!=[]){G=car(L);L=cdr(L);
 if((S=polrootsreal4(adp(cons([V,G],L1))))!=[])
  { if((S=toint4([V,call(subst,cons(G,L2))],S))!=[])
    W=append(W,base_makelist([s,G],s,S,0));}
 }
 W=qsort(W,sort1);U=[];C=length(W);
 while(W!=[]){U=cons(append(car(W),[[F,C]]),U);W=cdr(W);C--;}
 return U;
}
def zeros2(V,F,L1,L2)
 {
 G=simpalg_noalg(F,L1); /*print(G);*/
 if(G==0) G=subst3(F,L1); else 0; /*print([V,F,L1,L2]);*/
 G=prod(factors([prod(flatmf(asq_noalg(cons([V,G],L1))))]));
 if(G==1) return [];
 W=[];L=factors([adp(cons([V,G],L1))]);
 while(L!=[]){A=car(L);L=cdr(L);
 if((S=polrootsreal4(A))!=[])
  { if((S=toint4([V,call(subst,cons(G,L2))],S))!=[])
    W=append(W,base_makelist([s,A],s,S,0));}
 }
 W=qsort(W,sort1);U=[];C=length(W);
 while(W!=[]){U=cons(append(car(W),[[G,C],[F,G]]),U);W=cdr(W);C--;}
 return U;
}
def asq_noalg(DL)
{   F = car(DL)[1]; V = car(DL)[0]; DL = cdr(DL);
	DL = reverse(DL);
	N = length(DL);
	Tab = newvect(N);
	/* Tab = [[a1,r1],...]; ri is a root of di(t,r(i-1),...,r1). */
	AL = [];
	for ( I = 0; I < N; I++ ) {
		T = DL[I];
		for ( J = 0, DP = T[1]; J < I; J++ )
			DP = subst(DP,Tab[J][0],Tab[J][1]);
		B = newalg(DP);
		Tab[I] = [T[0],B];
		F = subst(F,T[0],B);
		AL = cons(B,AL);
	}
		FL = asq(red(monica(simpalg(F),V)));
	for ( T = FL, R = []; T != []; T = cdr(T) )
		R = cons([conv_noalg(T[0][0],Tab),T[0][1]],R);
	return reverse(R);
}
def simpalg_noalg(F,DL)
{
	DL = reverse(DL);
	N = length(DL);
	Tab = newvect(N);
	/* Tab = [[a1,r1],...]; ri is a root of di(t,r(i-1),...,r1). */
	AL = [];
	for ( I = 0; I < N; I++ ) {
		T = DL[I];
		for ( J = 0, DP = T[1]; J < I; J++ )
			DP = subst(DP,Tab[J][0],Tab[J][1]);
		B = newalg(DP);
		Tab[I] = [T[0],B];
		F = subst(F,T[0],B);
		AL = cons(B,AL);
	}
	return conv_noalg(simpalg(F),Tab);
}
def subst3(F/* must be banish */,L)
/*
 subst3((a+1)*x^2+(a+b)*x+b^2-1,[[b,b+a],[a,a+1]]);
*/
 {L=reverse(L);M=[];
 while(L!=[]) {M=cons(car(L),M);L=cdr(L);
  while((G=simpalg_noalg(F,M))==0) {F=diff(F,M[0][0]);}
  F=G;
 }
 return F;
}
def rdp0(L)
/* relative def. poly. using af_noalg.
   rdp0([[x,a*x^2+b*x+1],[b,b^2-4*a],[a,a^2-2]]); */
{ return prod(flatmf(af_noalg(car(L)[1],cdr(L)))); }
def rdp(L)
/* relative def. poly. using asq_noalg. faster than by af_noalg.
   rdp([[x,a*x^2+b*x+1],[b,b^2-4*a],[a,a^2-2]]);
*/
{ return prod(flatmf(asq_noalg(L))); }
def monica(F,V) {return F==0?0:F/coef(F,deg(F,V),V);}
def polsqf(F,V) {return F==0?0:sdiv(F,gcd(F,diff(F,V)),V);}
def polsqf1(F) {return prod(flatmf(sqfr(F)));}
def polsqfa(A,P,B,X) " polsqfa(p^2-2,p,(x-p)^4*(p*x^2+4*x+2*p)^5,x); "
{ 
 if(B==0) return 0;
 else {C=simpcoef(simpalg(monica(subst(B,P,P1=newalg(A)),X)));
       return subst(algptorat(simpalg(
        sdiv(C,cr_gcda(C,diff(C,X)),X))),algptorat(P1),P);}
}
cputime(1)$
end$

2017-09-23

REDUCE 更新

https://sourceforge.net/projects/reduce-algebra/files/snapshot_2017-09-22

インストールスクリプト

svn co http://svn.code.sf.net/p/reduce-algebra/code/trunk reduce-algebra
cd ./reduce-algebra
./configure --with-psl
make
sudo apt-get install libffi-dev
./configure --with-csl
make
cd ./generic/redfront
sudo make install

実行例.

$ rfpsl -b

Redfront 3.2, built 23-Sep-2017 ...
Reduce (Free PSL version, revision 4218), 23-Sep-2017 ...

1: load redlog;rlset ofsf;on rlverbose;off nat;

5: rlcad all({x},a*x^2+b*x+c>=0);

+++ Preparation Phase
+ Kernel order set to (x c b a)
+++ Projection Phase
+ projection order: (x c b a)
+ number of all projection factors: 5
+ number of projection factors of level r,...,1: 1,2,1,1
+++ Extension Phase
+ Building partial CAD tree...
+ number of partial CAD tree nodes of level 0,...,r: 1,3,6,16,0
+++ Simple Solution Formula Construction Phase
+ levels to be considered: (1 2 3)
+ Generating signatures for 4 polynomials and 19 cells.
+ Number of cells on level 3 is 16.
+ SSFC succeded.
+++ Finish Phase
+ Kernel order was restored.

c >= 0 and b = 0 and a >= 0 or c > 0 and b <> 0 and a > 0 and 4*a*c - b**2 >= 0
$

6: rlslfq ws;

+++ creating /tmp/rlslfq-k-18455.slfq ... done

+++ calling slfq -N 100 -P 200000 < /tmp/rlslfq-k-18455.slfq 2> /dev/null | awk -v rf=/tmp/rlslfq-k-18455.red -v verb=t -v time=nil -v slfqvb=nil -v name=SLFQ -f qepcad.awk
+++ SLFQ raw output: [ c >= 0 /\ a >= 0 /\ 4 c a - b^2 >= 0 /\ [ b = 0 \/ c > 0 ] ]

c >= 0 and a >= 0 and 4*a*c - b**2 >= 0 and (b = 0 or c > 0)$

7: rlqepcad all({x},a*x^2+b*x+c>=0);

+++ creating /tmp/rlqepcad-k-18455.qepcad ... done

+++ calling qepcad +N100000000 +L200000 < /tmp/rlqepcad-k-18455.qepcad | awk -v rf=/tmp/rlqepcad-k-18455.red -v verb=t -v time=nil -v slfqvb=nil -v name=QEPCAD -f qepcad.awk
+++ QEPCAD raw output: c >= 0 /\ a >= 0 /\ 4 c a - b^2 >= 0

c >= 0 and a >= 0 and 4*a*c - b**2 >= 0$