ATPとCASのこと

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$

2017-08-26

実行例(5)

入力は,主に http://www.cs.bath.ac.uk/~djw42/triangular/examplebank.pdf の Examples from [CMXY09](https://arxiv.org/pdf/0903.5221.pdf)から選びました.

環境は,Core i5-3210M 2.50GHz/8GB/Ubuntu 14.04/Maxima 5.40.0/SBCL 1.3.20 です.

コントローラー qvpeds の第 4 引数が h1 は Hong タイプ,l は Lazard タイプの処理を表し,%pp はそれぞれの projection set です.

(%i4) [(qvpeds([ex],[a,b,c,x],0,h1,r11,0),
        qe(bfpcad(ext('(a*x^2+b*x+c = 0))))),%pp]
Evaluation took 0.2600 seconds (0.2620 elapsed) using 30.789 MB.
(%o4) [[[a < root(a,1),true,root(4*a*c-b^2,1) <= c,true],
        [a = root(a,1),b < root(b,1),true,true],
        [a = root(a,1),b = root(b,1),c = root(c,1),true],
        [a = root(a,1),root(b,1) < b,true,true],
        [root(a,1) < a,true,c <= root(4*a*c-b^2,1),true]],
       [[a],[b],[c,4*a*c-b^2],[a*x^2+b*x+c]]]
(%i5) [(qvpeds([ex],[a,b,c,x],0,l,r11,0),qe(bfpcad(ext('(a*x^2+b*x+c = 0))))),%pp]
Evaluation took 0.2110 seconds (0.2090 elapsed) using 30.792 MB.
(%o5) [[[a < root(a,1),true,root(4*a*c-b^2,1) <= c,true],
        [a = root(a,1),b < root(b,1),true,true],
        [a = root(a,1),b = root(b,1),c = root(c,1),true],
        [a = root(a,1),root(b,1) < b,true,true],
        [root(a,1) < a,true,c <= root(4*a*c-b^2,1),true]],
       [[a],[b],[c,4*a*c-b^2],[a*x^2+b*x+c]]]
(%i6) [(qvpeds([all],[c,b,a,x],0,h1,r11,0),
        qe(bfpcad(ext('(x^4+a*x^2+b*x+c >= 0))))),%pp]
Evaluation took 2.7640 seconds (2.7610 elapsed) using 665.230 MB.
(%o6) [[[c = root(c,1),b = root(126976*c^3-135*b^4,1),root(a,1) <= a,true],
        [root(c,1) < c,b < root(b,1),
         root(256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2,2)
           <= a,true],
        [root(c,1) < c,b = root(b,1),
         root(256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2,1)
           <= a,true],
        [root(c,1) < c,root(b,1) < b,
         root(256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2,2)
           <= a,true]],
       [[c],
        [b,256*c^3-27*b^4,1024*c^3-2187*b^4,1024*c^3+3*b^4,4096*c^3+27*b^4,
         126976*c^3-135*b^4],
        [a,8*a*c-9*b^2-2*a^3,
         256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2],
        [x^4+a*x^2+b*x+c]]]
(%i7) [(qvpeds([all],[c,b,a,x],0,l,r11,0),
        qe(bfpcad(ext('(x^4+a*x^2+b*x+c >= 0))))),%pp]
Evaluation took 0.7110 seconds (0.7090 elapsed) using 174.679 MB.
(%o7) [[[c = root(c,1),b = root(4096*c^3+27*b^4,1),root(a,1) <= a,true],
        [root(c,1) < c,b < root(b,1),
         root(256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2,2)
           <= a,true],
        [root(c,1) < c,b = root(b,1),
         root(256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2,1)
           <= a,true],
        [root(c,1) < c,root(b,1) < b,
         root(256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2,2)
           <= a,true]],
       [[c],[b,256*c^3-27*b^4,4096*c^3+27*b^4],
        [256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2],
        [x^4+a*x^2+b*x+c]]]
(%i8) [(qvpeds([all],[d,c,a,b,x],0,h1,r11,0),
        qe(bfpcad(ext('(a*x^4+b*x^2+c*x+d >= 0))))),%pp]
Evaluation took 5.4500 seconds (5.4580 elapsed) using 1195.996 MB.
(%o8) [[[d = root(d,1),c = root(c,1),root(a,1) <= a,root(b,1) <= b,true],
        [root(d,1) < d,c < root(c,1),root(a,1) <= a,
         root(256*a^2*d^3-128*a*b^2*d^2+(144*a*b*c^2+16*b^4)*d-27*a*c^4
                         -4*b^3*c^2,2)
           <= b,true],
        [root(d,1) < d,c = root(c,1),root(126976*a*d^3-135*c^4,1) <= a,
         root(256*a^2*d^3-128*a*b^2*d^2+(144*a*b*c^2+16*b^4)*d-27*a*c^4
                         -4*b^3*c^2,1)
           <= b,true],
        [root(d,1) < d,root(c,1) < c,root(a,1) <= a,
         root(256*a^2*d^3-128*a*b^2*d^2+(144*a*b*c^2+16*b^4)*d-27*a*c^4
                         -4*b^3*c^2,2)
           <= b,true]],
       [[d],[c],
        [a,256*a*d^3-27*c^4,1024*a*d^3-2187*c^4,1024*a*d^3+3*c^4,
         4096*a*d^3+27*c^4,126976*a*d^3-135*c^4],
        [b,8*a*b*d-9*a*c^2-2*b^3,
         256*a^2*d^3-128*a*b^2*d^2+(144*a*b*c^2+16*b^4)*d-27*a*c^4-4*b^3*c^2],
        [a*x^4+b*x^2+c*x+d]]]
(%i9) [(qvpeds([all],[d,c,a,b,x],0,l,r11,0),
        qe(bfpcad(ext('(a*x^4+b*x^2+c*x+d >= 0))))),%pp]
Evaluation took 1.5370 seconds (1.5370 elapsed) using 368.276 MB.
(%o9) [[[d = root(d,1),c = root(c,1),root(a,1) <= a,root(b,1) <= b,true],
        [root(d,1) < d,c < root(c,1),root(a,1) <= a,
         root(256*a^2*d^3-128*a*b^2*d^2+(144*a*b*c^2+16*b^4)*d-27*a*c^4
                         -4*b^3*c^2,2)
           <= b,true],
        [root(d,1) < d,c = root(c,1),root(4096*a*d^3+27*c^4,1) <= a,
         root(256*a^2*d^3-128*a*b^2*d^2+(144*a*b*c^2+16*b^4)*d-27*a*c^4
                         -4*b^3*c^2,1)
           <= b,true],
        [root(d,1) < d,root(c,1) < c,root(a,1) <= a,
         root(256*a^2*d^3-128*a*b^2*d^2+(144*a*b*c^2+16*b^4)*d-27*a*c^4
                         -4*b^3*c^2,2)
           <= b,true]],
       [[d],[c],[a,256*a*d^3-27*c^4,4096*a*d^3+27*c^4],
        [256*a^2*d^3-128*a*b^2*d^2+(144*a*b*c^2+16*b^4)*d-27*a*c^4-4*b^3*c^2],
        [a*x^4+b*x^2+c*x+d]]]
(%i10) [(qvpeds([],[x,z,y],0,h1,r11,0),
         qe(bfpcad(ext('(z^2+y^2+x^2-1 = 0 and z^3+x*z+y = 0))))),%pp]
Evaluation took 0.2330 seconds (0.2330 elapsed) using 44.868 MB.
(%o10) [[[x = root(x+1,1),z = root(z^6+2*x*z^4+(x^2+1)*z^2+x^2-1,1),
          y = root(z^3+x*z+y,1)],
         [root(x+1,1) < x and x < root(x-1,1),
          z = root(z^6+2*x*z^4+(x^2+1)*z^2+x^2-1,1)
            or z = root(z^6+2*x*z^4+(x^2+1)*z^2+x^2-1,2),
          y = root(z^3+x*z+y,1)],
         [x = root(x-1,1),z = root(z^6+2*x*z^4+(x^2+1)*z^2+x^2-1,1),
          y = root(z^3+x*z+y,1)]],
        [[x-1,x,x+1,x^2-3,3*x^3-2*x^2-3*x-2,
          4*x^5-31*x^4+32*x^3+46*x^2-36*x-31],[z^6+2*x*z^4+(x^2+1)*z^2+x^2-1],
         [z^3+x*z+y]]]
(%i11) [(qvpeds([],[x,z,y],0,l,r11,0),
         qe(bfpcad(ext('(z^2+y^2+x^2-1 = 0 and z^3+x*z+y = 0))))),%pp]
Evaluation took 0.1250 seconds (0.1300 elapsed) using 22.110 MB.
(%o11) [[[x = root(x+1,1),z = root(z^6+2*x*z^4+(x^2+1)*z^2+x^2-1,1),
          y = root(z^3+x*z+y,1)],
         [root(x+1,1) < x and x < root(x-1,1),
          z = root(z^6+2*x*z^4+(x^2+1)*z^2+x^2-1,1)
            or z = root(z^6+2*x*z^4+(x^2+1)*z^2+x^2-1,2),
          y = root(z^3+x*z+y,1)],
         [x = root(x-1,1),z = root(z^6+2*x*z^4+(x^2+1)*z^2+x^2-1,1),
          y = root(z^3+x*z+y,1)]],
        [[x-1,x+1,4*x^5-31*x^4+32*x^3+46*x^2-36*x-31],
         [z^6+2*x*z^4+(x^2+1)*z^2+x^2-1],[z^3+x*z+y]]]
(%i12) [(qvpeds([],[y,x],0,h1,r11,0),
         qe(bfpcad(ext('(y^4-2*y^3+y^2+(-3)*x^2*y+2*x^4 = 0))))),%pp]
Evaluation took 0.1530 seconds (0.1530 elapsed) using 20.637 MB.
(%o12) [[[y = root(y,1),x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,1)],
         [root(y,1) < y and y < root(y-1,1),
          x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,1)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,2)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,3)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,4)],
         [y = root(y-1,1),
          x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,1)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,2)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,3)],
         [root(y-1,1) < y and y < root(8*y^2-16*y-1,2),
          x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,1)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,2)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,3)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,4)],
         [y = root(8*y^2-16*y-1,2),
          x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,1)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,2)]],
        [[y-1,y,8*y^2-16*y-1],[y^4-2*y^3+y^2-3*x^2*y+2*x^4]]]
(%i13) [(qvpeds([],[y,x],0,l,r11,0),
         qe(bfpcad(ext('(y^4-2*y^3+y^2+(-3)*x^2*y+2*x^4 = 0))))),%pp]
Evaluation took 0.1170 seconds (0.1170 elapsed) using 20.605 MB.
(%o13) [[[y = root(y,1),x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,1)],
         [root(y,1) < y and y < root(y-1,1),
          x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,1)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,2)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,3)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,4)],
         [y = root(y-1,1),
          x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,1)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,2)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,3)],
         [root(y-1,1) < y and y < root(8*y^2-16*y-1,2),
          x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,1)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,2)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,3)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,4)],
         [y = root(8*y^2-16*y-1,2),
          x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,1)
            or x = root(y^4-2*y^3+y^2-3*x^2*y+2*x^4,2)]],
        [[y-1,y,8*y^2-16*y-1],[y^4-2*y^3+y^2-3*x^2*y+2*x^4]]]
(%i14) [(qvpeds([],[y,x],0,h1,r11,0),
         qe(bfpcad(ext('(144*y^2+96*x^2*y+9*x^4+105*x^2+70*x-98 = 0
                         and x*y^2+6*x*y+x^3+9*x = 0))))),%pp]
Evaluation took 0.0780 seconds (0.0780 elapsed) using 9.584 MB.
(%o14) [[[y = root(72*y^2-49,1) or y = root(72*y^2-49,2),
          x = root(1132535906187264*y^9+2648683072534464*y^8
                                       -11812443576362496*y^7
                                       -144938986561480320*y^6
                                       -169128077155358976*y^5
                                       +764415088047207513*y^4
                                       +4188447484701833676*y^3
                                       +3270406842477283485*y^2
                                       -2768658167053333674*y
                                       +419177363858840230*x
                                       -2534619807825082586,1)]],
        [[72*y^2-49],
         [1132535906187264*y^9+2648683072534464*y^8-11812443576362496*y^7
                              -144938986561480320*y^6-169128077155358976*y^5
                              +764415088047207513*y^4+4188447484701833676*y^3
                              +3270406842477283485*y^2-2768658167053333674*y
                              +419177363858840230*x-2534619807825082586]]]
(%i15) [(qvpeds([],[y,x],0,l,r11,0),
         qe(bfpcad(ext('(144*y^2+96*x^2*y+9*x^4+105*x^2+70*x-98 = 0
                         and x*y^2+6*x*y+x^3+9*x = 0))))),%pp]
Evaluation took 0.0780 seconds (0.0800 elapsed) using 9.506 MB.
(%o15) [[[y = root(72*y^2-49,1) or y = root(72*y^2-49,2),
          x = root(1132535906187264*y^9+2648683072534464*y^8
                                       -11812443576362496*y^7
                                       -144938986561480320*y^6
                                       -169128077155358976*y^5
                                       +764415088047207513*y^4
                                       +4188447484701833676*y^3
                                       +3270406842477283485*y^2
                                       -2768658167053333674*y
                                       +419177363858840230*x
                                       -2534619807825082586,1)]],
        [[72*y^2-49],
         [1132535906187264*y^9+2648683072534464*y^8-11812443576362496*y^7
                              -144938986561480320*y^6-169128077155358976*y^5
                              +764415088047207513*y^4+4188447484701833676*y^3
                              +3270406842477283485*y^2-2768658167053333674*y
                              +419177363858840230*x-2534619807825082586]]]
(%i16) [(qvpeds([ex,ex],[x,y,z,a,b],0,h1,r11,0),
         qe(bfpcad(ext('(x-a*b = 0 and y-a*b^2 = 0 and z-a^2 = 0))))),%pp]
Evaluation took 0.4380 seconds (0.4390 elapsed) using 97.523 MB.
(%o16) [[[x < root(x,1),y < root(y,1) or root(y,1) < y,z = root(y^2*z-x^4,1),
          true,true],[x = root(x,1),y = root(y,1),root(z,1) <= z,true,true],
         [root(x,1) < x,y < root(y,1) or root(y,1) < y,z = root(y^2*z-x^4,1),
          true,true]],
        [[x],[y],[z,z*y^2-x^4],[a,a*y-x^2,(-a^2)+z,z*y-a*x^2],
         [x-a*b,y-b*x,(-a*x)+z*b]]]
(%i17) [(qvpeds([ex,ex],[x,y,z,a,b],0,l,r11,0),
         qe(bfpcad(ext('(x-a*b = 0 and y-a*b^2 = 0 and z-a^2 = 0))))),%pp]
Evaluation took 0.4270 seconds (0.4270 elapsed) using 96.663 MB.
(%o17) [[[x < root(x,1),y < root(y,1) or root(y,1) < y,z = root(y^2*z-x^4,1),
          true,true],[x = root(x,1),y = root(y,1),root(z,1) <= z,true,true],
         [root(x,1) < x,y < root(y,1) or root(y,1) < y,z = root(y^2*z-x^4,1),
          true,true]],
        [[x],[y],[z,(-x^4)+z*y^2],[a,(-x^2)+y*a,(-a^2)+z,(-a*x^2)+z*y],
         [x-a*b,(-b*x)+y,(-a*x)+z*b]]]
(%i18) [(qvpeds([ex,ex,ex],[x,y,z],0,h1,r11,0),
         qe(bfpcad(ext('(x^2+y^2+z^2-1 < 0 and x^2+(y+z-2)^2-1 < 0))))),%pp]
Evaluation took 0.3330 seconds (0.3330 elapsed) using 73.647 MB.
(%o18) [[[true,true,true]],
        [[x-1,x,x+1,x^4+22*x^2-7],
         [y^2+x^2-1,y^4-4*y^3+(x^2+7)*y^2+((-4*x^2)-4)*y+4*x^2],
         [z^2+y^2+x^2-1,z^2+(2*y-4)*z+y^2-4*y+x^2+3]]]
(%i19) [(qvpeds([ex,ex,ex],[x,y,z],0,l,r11,0),
         qe(bfpcad(ext('(x^2+y^2+z^2-1 < 0 and x^2+(y+z-2)^2-1 < 0))))),%pp]
Evaluation took 0.3160 seconds (0.3170 elapsed) using 65.861 MB.
(%o19) [[[true,true,true]],
        [[x-1,x,x+1,x^4+22*x^2-7],
         [y^2+x^2-1,y^4-4*y^3+(x^2+7)*y^2+((-4*x^2)-4)*y+4*x^2],
         [z^2+y^2+x^2-1,z^2+(2*y-4)*z+y^2-4*y+x^2+3]]]
(%i20) [(qvpeds([ex,all,all],[r,x,y],0,h1,r11,0),
         qe(bfpcad(ext('(x-r > 0 and y-r > 0
                         implies x^2*(1+2*y)^2-y^2*(1+2*x^2) > 0))))),%pp]
Evaluation took 0.9550 seconds (0.9550 elapsed) using 189.163 MB.
(%o20) [[[true,true,true]],
        [[r,r+2,2*r+1,4*r+1,2*r^2-1,2*r^2+4*r+1],
         [x,x-r,2*x^2-1,2*x^2+1,(2*r^2+4*r+1)*x^2-r^2],
         [y-r,(2*x^2-1)*y^2+4*x^2*y+x^2]]]
(%i21) [(qvpeds([ex,all,all],[r,x,y],0,l,r11,0),
         qe(bfpcad(ext('(x-r > 0 and y-r > 0
                         implies x^2*(1+2*y)^2-y^2*(1+2*x^2) > 0))))),%pp]
Evaluation took 0.9320 seconds (0.9330 elapsed) using 175.511 MB.
(%o21) [[[true,true,true]],
        [[r,r+2,2*r+1,4*r+1,2*r^2-1,2*r^2+4*r+1],
         [x,x-r,2*x^2-1,2*x^2+1,(2*r^2+4*r+1)*x^2-r^2],
         [y-r,(2*x^2-1)*y^2+4*x^2*y+x^2]]]
(%i22) [(qvpeds([ex],[a,b,r],0,h1r,r11,0),
         qe(bfpcad(ext('(3*a^2*r+3*b^2+(-2)*a*r-a^2-b^2 < 0
                         and 3*a^2*r+3*b^2*r+(-4)*a*r+r+(-2)*a^2+(-2)*b^2+2*a
                          > 0 and a-1/2 >= 0 and b > 0 and r > 0
                         and r-1 < 0))))),%pp]
Evaluation took 2.5190 seconds (2.5200 elapsed) using 659.301 MB.
(%o22) [[[root(2*a-1,1) <= a and a < root(3*a-2,1),
          root(b,1) < b and b < root(b^2+a^2-a,2),true],
         [a = root(3*a-2,1),root(b,1) < b and b <= root(b^2+a^2-a,2),true],
         [root(3*a-2,1) < a and a < root(a-1,1),
          root(b,1) < b and b < root(6*b^4+(9*a^2-12*a+2)*b^2+3*a^4-6*a^3
                                          +3*a^2,4),true],
         [a = root(a-1,1),
          root(b,1) < b and b < root(6*b^4+(9*a^2-12*a+2)*b^2+3*a^4-6*a^3
                                          +3*a^2,3),true],
         [root(a-1,1) < a and a < root(9*a^4-72*a^3+108*a^2-48*a+4,3),
          root(6*b^4+(9*a^2-12*a+2)*b^2+3*a^4-6*a^3+3*a^2,3) < b
            and b < root(6*b^4+(9*a^2-12*a+2)*b^2+3*a^4-6*a^3+3*a^2,4),true]],
        [[9*a^4-72*a^3+108*a^2-48*a+4,9*a^2-12*a+2,3*a-1,3*a-2,2*a-1,a,a-1],
         [6*b^4+(9*a^2-12*a+2)*b^2+3*a^4-6*a^3+3*a^2,3*b^2+3*a^2-4*a+1,
          2*b^2-a^2,b^2+a^2-a,b^2+a^2-2*a+1,b],
         [(3*b^2+3*a^2-4*a+1)*r-2*b^2-2*a^2+2*a,(3*a^2-2*a)*r+2*b^2-a^2,r,
          r-1]]]
(%i23) [(qvpeds([ex],[a,b,r],0,l,r11,0),
         qe(bfpcad(ext('(3*a^2*r+3*b^2+(-2)*a*r-a^2-b^2 < 0
                         and 3*a^2*r+3*b^2*r+(-4)*a*r+r+(-2)*a^2+(-2)*b^2+2*a
                          > 0 and a-1/2 >= 0 and b > 0 and r > 0
                         and r-1 < 0))))),%pp]
Evaluation took 2.1620 seconds (2.1620 elapsed) using 514.650 MB.
(%o23) [[[root(2*a-1,1) <= a and a < root(3*a-2,1),
          root(b,1) < b and b < root(b^2+a^2-a,2),true],
         [a = root(3*a-2,1),root(b,1) < b and b < root(2*b^2-a^2,2),true],
         [root(3*a-2,1) < a and a < root(a-1,1),
          root(b,1) < b and b < root(6*b^4+(9*a^2-12*a+2)*b^2+3*a^4-6*a^3
                                          +3*a^2,4),true],
         [a = root(a-1,1),
          root(6*b^4+(9*a^2-12*a+2)*b^2+3*a^4-6*a^3+3*a^2,2) < b
            and b < root(6*b^4+(9*a^2-12*a+2)*b^2+3*a^4-6*a^3+3*a^2,3),true],
         [root(a-1,1) < a and a < root(9*a^4-72*a^3+108*a^2-48*a+4,3),
          root(6*b^4+(9*a^2-12*a+2)*b^2+3*a^4-6*a^3+3*a^2,3) < b
            and b < root(6*b^4+(9*a^2-12*a+2)*b^2+3*a^4-6*a^3+3*a^2,4),true]],
        [[a-1,a,2*a-1,3*a-2,3*a-1,9*a^4-72*a^3+108*a^2-48*a+4],
         [b,b^2+a^2-2*a+1,b^2+a^2-a,2*b^2-a^2,3*b^2+3*a^2-4*a+1,
          6*b^4+(9*a^2-12*a+2)*b^2+3*a^4-6*a^3+3*a^2],
         [r-1,r,(3*a^2-2*a)*r+2*b^2-a^2,
          (3*b^2+3*a^2-4*a+1)*r-2*b^2-2*a^2+2*a]]]
(%i24) [(qvpeds([all,all,all,ex],[y,x,a,b,c,z],0,h1,r11,0),
         qe(bfpcad(ext('(a > 0 and a*z^2+b*z+c # 0
                         implies a*x^2+b*x+c-y > 0))))),%pp]
Evaluation took 3.0010 seconds (3.0030 elapsed) using 622.586 MB.
(%o24) [[[y <= root(y,1),true,true,true,true,true]],
        [[y],[x],[a,(-x^2*a)+y],
         [b,(-x*b)-x^2*a+y,(-b^2)-4*x*a*b-4*x^2*a^2+4*y*a],
         [c,4*a*c-b^2,(-c)-x*b-x^2*a+y],[a*z^2+b*z+c]]]
(%i25) [(qvpeds([all,all,all,ex],[y,x,a,b,c,z],0,l,r11,0),
         qe(bfpcad(ext('(a > 0 and a*z^2+b*z+c # 0
                         implies a*x^2+b*x+c-y > 0))))),%pp]
Evaluation took 2.9450 seconds (2.9450 elapsed) using 598.900 MB.
(%o25) [[[y <= root(y,1),true,true,true,true,true]],
        [[y],[x],[a,(-x^2*a)+y],
         [b,(-x*b)-x^2*a+y,(-b^2)-4*x*a*b-4*x^2*a^2+4*y*a],
         [c,4*a*c-b^2,(-c)-x*b-x^2*a+y],[a*z^2+b*z+c]]]
(%i26) [(qvpeds([all,all],[b,a,c,x,y],0,h1r,r11,0),
         qe(bfpcad(ext('(0 < a and 0 < b
                               and ((x-c)^2/a^2+y^2/b^2 = 1
                                implies x^2+y^2 <= 1)))))),%pp]
Evaluation took 76.6110 seconds (76.6430 elapsed) using 15531.535 MB.
(%o26) [[[root(b,1) < b and b < root(b-1,1),
          root(a,1) < a and a < root(b^2-a,1),
          root(b^2*c^2+b^4+((-a^2)-1)*b^2+a^2,1) <= c
            and c <= root(b^2*c^2+b^4+((-a^2)-1)*b^2+a^2,2),true,true],
         [root(b,1) < b and b < root(b-1,1),
          root(b^2-a,1) <= a and a < root(a-1,1),
          root(c-a+1,1) <= c and c <= root(c+a-1,1),true,true],
         [root(b,1) < b and b < root(b-1,1),a = root(a-1,1),c = root(c,1),
          true,true],
         [b = root(b-1,1),root(a,1) < a and a <= root(a-1,1),c = root(c,1),
          true,true]],
        [[b+1,b,b-1],
         [(2*a+1)*b^2+a^2,(2*a-1)*b^2-a^2,b^2+a,b^2-a,b+a,b-a,a+1,a,a-1],
         [b^2*c^2+b^4+((-a^2)-1)*b^2+a^2,b^2*c^2-a^2*b^2+a^2,c+a+1,c+a-1,
          c-a+1,c-a-1,c],
         [(b^2-a^2)*x^2-2*b^2*c*x+b^2*c^2-a^2*b^2+a^2,x-c+a,x-c-a,x+1,x-1],
         [a^2*y^2+b^2*x^2-2*b^2*c*x+b^2*c^2-a^2*b^2,y^2+x^2-1]]]
(%i27) [(qvpeds([all,all],[b,a,c,x,y],0,l,r11,0),
         qe(bfpcad(ext('(0 < a and 0 < b
                               and ((x-c)^2/a^2+y^2/b^2 = 1
                                implies x^2+y^2 <= 1)))))),%pp]
Evaluation took 62.6640 seconds (62.6860 elapsed) using 11842.095 MB.
(%o27) [[[root(b,1) < b and b < root(b-1,1),
          root(a,1) < a and a <= root(b^2-a,1),
          root(b^2*c^2+b^4+((-a^2)-1)*b^2+a^2,1) <= c
            and c <= root(b^2*c^2+b^4+((-a^2)-1)*b^2+a^2,2),true,true],
         [root(b,1) < b and b < root(b-1,1),
          root(b^2-a,1) < a and a < root(a-1,1),
          root(c-a+1,1) <= c and c <= root(c+a-1,1),true,true],
         [root(b,1) < b and b < root(b-1,1),a = root(a-1,1),c = root(c+a-1,1),
          true,true],
         [b = root(b-1,1),root(a,1) < a and a <= root((2*a-1)*b^2-a^2,1),
          c = root(b^2*c^2+b^4+((-a^2)-1)*b^2+a^2,1),true,true]],
        [[b-1,b,b+1],
         [a-1,a,a+1,b-a,b+a,b^2-a,b^2+a,(2*a-1)*b^2-a^2,(2*a+1)*b^2+a^2],
         [c-a-1,c-a+1,c+a-1,c+a+1,b^2*c^2-a^2*b^2+a^2,
          b^2*c^2+b^4+((-a^2)-1)*b^2+a^2],
         [x-1,x+1,x-c-a,x-c+a,(b^2-a^2)*x^2-2*b^2*c*x+b^2*c^2-a^2*b^2+a^2],
         [y^2+x^2-1,a^2*y^2+b^2*x^2-2*b^2*c*x+b^2*c^2-a^2*b^2]]]
(%i28) [(qvpeds([ex,all,all],[d,c,a,b],0,h1r,r11,0),
         qe(bfpcad(ext('((a-d = 0 and b-c = 0 or a-c = 0 and b-1 = 0)
                         implies a^2-b = 0))))),%pp]
Evaluation took 9.7020 seconds (9.7060 elapsed) using 1360.997 MB.
(%o28) [[[d = root(d+1,1) or d = root(d-1,1),true,true,true]],
        [[d+1,d,d-1],[(-c)+d^2,(-c)+d,c+1,c,c-1],
         [(-a)+d,(-a^2)+c,(-a)+c,a+1,a-1],[(-b)+c,b-a^2,b-1]]]
(%i29) [(qvpeds([ex,all,all],[d,c,a,b],0,l,r11,0),
         qe(bfpcad(ext('((a-d = 0 and b-c = 0 or a-c = 0 and b-1 = 0)
                         implies a^2-b = 0))))),%pp]
Evaluation took 9.5980 seconds (9.6010 elapsed) using 1325.779 MB.
(%o29) [[[d = root(d+1,1) or d = root(d-1,1),true,true,true]],
        [[d-1,d,d+1],[c-1,c,c+1,(-c)+d,(-c)+d^2],
         [a-1,a+1,(-a)+c,(-a^2)+c,(-a)+d],[b-1,b-a^2,(-b)+c]]]
(%i30) [(qvpeds([ex,ex,ex],[b,a,x,y,z],0,h1,r11,0),
         qe(bfpcad(ext('(x+y+z = 0 and x*y+y*z+z*x-a = 0 and x*y*z-b = 0))))),
        %pp]
Evaluation took 1.1200 seconds (1.1250 elapsed) using 273.509 MB.
(%o30) [[[true,a <= root(27*b^2+4*a^3,1),true,true,true]],
        [[b],[a,4*a^3+27*b^2],[3*x^2+4*a,x^3+a*x-b],[y^2+x*y+x^2+a],[y+x+z]]]
(%i31) [(qvpeds([ex,ex,ex],[b,a,x,y,z],0,l,r11,0),
         qe(bfpcad(ext('(x+y+z = 0 and x*y+y*z+z*x-a = 0 and x*y*z-b = 0))))),
        %pp]
Evaluation took 0.7980 seconds (0.7980 elapsed) using 145.864 MB.
(%o31) [[[true,a <= root(27*b^2+4*a^3,1),true,true,true]],
        [[b],[a,4*a^3+27*b^2],[3*x^2+4*a,x^3+a*x-b],[y^2+x*y+x^2+a],[y+x+z]]]
(%i32) [(qvpeds([ex,ex,ex],[a,b,c,x,y,z],0,h1r,r11,0),
         qe(bfpcad(ext('(x+y+z = a and x*y+y*z+z*x = b and x*y*z = c))))),%pp]
Evaluation took 0.7660 seconds (0.7660 elapsed) using 184.901 MB.
(%o32) [[[true,b < root(3*b-a^2,1),
          root(27*c^2-18*a*b*c+4*a^3*c+4*b^3-a^2*b^2,1) <= c
            and c <= root(27*c^2-18*a*b*c+4*a^3*c+4*b^3-a^2*b^2,2),true,true,
          true],
         [true,b = root(3*b-a^2,1),
          c = root(27*c^2-18*a*b*c+4*a^3*c+4*b^3-a^2*b^2,1),true,true,true]],
        [[],[3*b-a^2],[27*c^2-18*a*b*c+4*a^3*c+4*b^3-a^2*b^2],
         [x^3-a*x^2+b*x-c,3*x^2-2*a*x+4*b-a^2],[y^2+x*y-a*y+x^2-a*x+b],
         [z+y+x-a]]]
(%i33) [(qvpeds([ex,ex,ex],[a,b,c,x,y,z],0,l,r11,0),
         qe(bfpcad(ext('(x+y+z = a and x*y+y*z+z*x = b and x*y*z = c))))),%pp]
Evaluation took 0.6550 seconds (0.6560 elapsed) using 140.720 MB.
(%o33) [[[true,b < root(3*b-a^2,1),
          root(27*c^2-18*a*b*c+4*a^3*c+4*b^3-a^2*b^2,1) <= c
            and c <= root(27*c^2-18*a*b*c+4*a^3*c+4*b^3-a^2*b^2,2),true,true,
          true],
         [true,b = root(3*b-a^2,1),
          c = root(27*c^2-18*a*b*c+4*a^3*c+4*b^3-a^2*b^2,1),true,true,true]],
        [[],[3*b-a^2],[27*c^2-18*a*b*c+4*a^3*c+4*b^3-a^2*b^2],
         [3*x^2-2*a*x+4*b-a^2,x^3-a*x^2+b*x-c],[y^2+x*y-a*y+x^2-a*x+b],
         [z+y+x-a]]]
(%i34) [(qvpeds([ex,ex,ex],[t,x,y],0,h1r,r11,0),
         qe(bfpcad(ext('((17/16)*t-6 >= 0 and (17/16)*t-10 <= 0
                                          and x-(17/16)*t+1 >= 0
                                          and x-(17/16)*t-1 <= 0
                                          and y-(17/16)*t+9 >= 0
                                          and y-(17/16)*t+7 <= 0
                                          and (x-t)^2+y^2-1 <= 0))))),%pp]
Evaluation took 3.1300 seconds (3.1330 elapsed) using 890.252 MB.
(%o34) [[[true,true,true]],
        [[145*t^2-1920*t+6272,145*t^2-2464*t+10368,17*t-96,17*t-112,17*t-128,
          17*t-144,17*t-160,t+32,t,t-32],
         [256*x^2-512*t*x+545*t^2-3808*t+12288,
          256*x^2-512*t*x+545*t^2-4896*t+20480,16*x-17*t+16,16*x-17*t-16,
          x-t+1,x-t-1],[y^2+x^2-2*t*x+t^2-1,16*y-17*t+144,16*y-17*t+112]]]
(%i35) [(qvpeds([ex,ex,ex],[t,x,y],0,l,r11,0),
         qe(bfpcad(ext('((17/16)*t-6 >= 0 and (17/16)*t-10 <= 0
                                          and x-(17/16)*t+1 >= 0
                                          and x-(17/16)*t-1 <= 0
                                          and y-(17/16)*t+9 >= 0
                                          and y-(17/16)*t+7 <= 0
                                          and (x-t)^2+y^2-1 <= 0))))),%pp]
Evaluation took 1.5770 seconds (1.5780 elapsed) using 246.028 MB.
(%o35) [[[true,true,true]],
        [[t-32,t,t+32,17*t-160,17*t-144,17*t-128,17*t-112,17*t-96,
          145*t^2-2464*t+10368,145*t^2-1920*t+6272],
         [x-t-1,x-t+1,16*x-17*t-16,16*x-17*t+16,
          256*x^2-512*t*x+545*t^2-4896*t+20480,
          256*x^2-512*t*x+545*t^2-3808*t+12288],
         [16*y-17*t+112,16*y-17*t+144,y^2+x^2-2*t*x+t^2-1]]]
(%i36) [fpprec,%ez,ratepsilon]:[48,1.0b-12,1.0b-48]
Evaluation took 0.0010 seconds (0.0000 elapsed) using 0 bytes.
(%o36) [48,1.0b-12,1.0b-48]
(%i37) f0:x < 0 and 50000*(x^2+y^2)-49719 < 0
           implies 360000+720000*x+720000*x^2+480000*x^3+240000*x^4+96000*x^5
                         +32200*x^6+9200*x^7+2225*x^8+450*x^9+75*x^10+10*x^11
                         +x^12+(-3000)*x^4*y^2+1200*x^5*y^2+2100*x^6*y^2
                         +1000*x^7*y^2+275*x^8*y^2+50*x^9*y^2+6*x^10*y^2
                         +3000*x^2*y^4+(-6000)*x^3*y^4+(-2250)*x^4*y^4
                         +300*x^5*y^4+350*x^6*y^4+100*x^7*y^4+15*x^8*y^4
                         +(-200)*y^6+2000*x*y^6+(-1900)*x^2*y^6+(-600)*x^3*y^6
                         +150*x^4*y^6+100*x^5*y^6+20*x^6*y^6+225*y^8
                         +(-350)*x*y^8+(-25)*x^2*y^8+50*x^3*y^8+15*x^4*y^8
                         +(-25)*y^10+10*x*y^10+6*x^2*y^10+y^12-360000
            < 0
Evaluation took 0.0290 seconds (0.0300 elapsed) using 3.558 MB.
(%i38) [(qvpeds([all,all],[x,y],0,h1,r11,0),qe(bfpcad(ext(f0)))),%pp]
Evaluation took 10.0710 seconds (10.0740 elapsed) using 2590.940 MB.
(%o38) [[[true,true]],
        [[x,6*x^2+10*x-25,28*x^2-68*x-17,50000*x^2-49719,
          x^5+5*x^4+25*x^3+100*x^2+300*x+600,
          28*x^6-204*x^5-101*x^4-190*x^3+5269*x^2-5710*x+1875,
          1088*x^6-7552*x^5+12736*x^4-3600*x^3+21000*x^2-41240*x+14039,
          200000000000000000000000000000000*x^6
           +599438000000000000000000000000000*x^5
           +1225000789610000000000000000000000*x^4
           +2442297336707065300000000000000000*x^3
           +3754710227977240336626750000000000*x^2
           +3758510381317956468954765766500000*x-5402844825834997735210773,
          1088*x^12-15104*x^11+59904*x^10-3344*x^9+112936*x^8-298200*x^7
                   -6467281*x^6+17310870*x^5-5750825*x^4+6290450*x^3
                   -14463075*x^2+6471500*x-580150,
          120832*x^12-989184*x^11+1965568*x^10+1699840*x^9+3373952*x^8
                     -10134400*x^7-57194464*x^6+73892480*x^5-12046600*x^4
                     -5395400*x^3-126535850*x^2+56676900*x-290075,
          120832*x^19-1648640*x^18+5940736*x^17+8692736*x^16-23252608*x^15
                     -246864000*x^14-1190564192*x^13-3079579840*x^12
                     +12973976600*x^11+1232865000*x^10-84076534050*x^9
                     -209944983400*x^8+39826614875*x^7+1457720531250*x^6
                     +2010504800625*x^5+25224378750*x^4-957912804375*x^3
                     -949924375000*x^2+750565675000*x-43511250000,
          4571136*x^19+2818048*x^18-25075712*x^17+34385920*x^16+930661376*x^15
                      +3433646080*x^14+8198400*x^13-47876089600*x^12
                      -184510347200*x^11-302452144000*x^10+194216259600*x^9
                      +2048648030000*x^8+4813290745000*x^7+5892848620000*x^6
                      +2971088911875*x^5-473950680000*x^4-1572986706250*x^3
                      +258053312500*x^2+440250359375*x+1631671875,
          4571136*x^28+4227072*x^27-51011584*x^26+208793600*x^25
                      +4487152640*x^24+21158707200*x^23-2748128000*x^22
                      -513743699200*x^21-1859278084800*x^20
                      +11836781184000*x^19+177227571163600*x^18
                      +1147497743010000*x^17+5032673524115000*x^16
                      +16422964823940000*x^15+40583034110816875*x^14
                      +73030624007188750*x^13+78801621529646875*x^12
                      -18223494396843750*x^11-288264068603906250*x^10
                      -681410411614046875*x^9-946722681541984375*x^8
                      -792370510167421875*x^7-273080203974921875*x^6
                      +181512821587890625*x^5+197555235933593750*x^4
                      -7722981467968750*x^3-126478640906250000*x^2
                      -6832921593750000*x-1599038437500000,
          855113728*x^28+21377843200*x^27+283256422400*x^26+2608096870400*x^25
                        +18521228902400*x^24+106982744064000*x^23
                        +518383694233600*x^22+2147375910912000*x^21
                        +7692855189760000*x^20+23980542444800000*x^19
                        +65142992800320000*x^18+153700917730560000*x^17
                        +311998085014880000*x^16+534257219535200000*x^15
                        +740482905349000000*x^14+746555405314000000*x^13
                        +321582900808375000*x^12-619086555599625000*x^11
                        -1781201761577531250*x^10-2496929229230000000*x^9
                        -2168728237020312500*x^8-905392388410156250*x^7
                        +397713892669921875*x^6+915322655607031250*x^5
                        +609971037097656250*x^4+165439917011718750*x^3
                        -20690934755859375*x^2+38769456621093750*x
                        -1019794921875],
         [50000*y^2+50000*x^2-49719,
          y^12+(6*x^2+10*x-25)*y^10+(15*x^4+50*x^3-25*x^2-350*x+225)*y^8
              +(20*x^6+100*x^5+150*x^4-600*x^3-1900*x^2+2000*x-200)*y^6
              +(15*x^8+100*x^7+350*x^6+300*x^5-2250*x^4-6000*x^3+3000*x^2)*y^4
              +(6*x^10+50*x^9+275*x^8+1000*x^7+2100*x^6+1200*x^5-3000*x^4)*y^2
              +x^12+10*x^11+75*x^10+450*x^9+2225*x^8+9200*x^7+32200*x^6
              +96000*x^5+240000*x^4+480000*x^3+720000*x^2+720000*x]]]
(%i39) [(qvpeds([all,all],[x,y],0,l,r11,0),qe(bfpcad(ext(f0)))),%pp]
Evaluation took 1.6160 seconds (1.6180 elapsed) using 603.005 MB.
(%o39) [[[true,true]],
        [[x,50000*x^2-49719,x^5+5*x^4+25*x^3+100*x^2+300*x+600,
          200000000000000000000000000000000*x^6
           +599438000000000000000000000000000*x^5
           +1225000789610000000000000000000000*x^4
           +2442297336707065300000000000000000*x^3
           +3754710227977240336626750000000000*x^2
           +3758510381317956468954765766500000*x-5402844825834997735210773,
          855113728*x^28+21377843200*x^27+283256422400*x^26+2608096870400*x^25
                        +18521228902400*x^24+106982744064000*x^23
                        +518383694233600*x^22+2147375910912000*x^21
                        +7692855189760000*x^20+23980542444800000*x^19
                        +65142992800320000*x^18+153700917730560000*x^17
                        +311998085014880000*x^16+534257219535200000*x^15
                        +740482905349000000*x^14+746555405314000000*x^13
                        +321582900808375000*x^12-619086555599625000*x^11
                        -1781201761577531250*x^10-2496929229230000000*x^9
                        -2168728237020312500*x^8-905392388410156250*x^7
                        +397713892669921875*x^6+915322655607031250*x^5
                        +609971037097656250*x^4+165439917011718750*x^3
                        -20690934755859375*x^2+38769456621093750*x
                        -1019794921875],
         [50000*y^2+50000*x^2-49719,
          y^12+(6*x^2+10*x-25)*y^10+(15*x^4+50*x^3-25*x^2-350*x+225)*y^8
              +(20*x^6+100*x^5+150*x^4-600*x^3-1900*x^2+2000*x-200)*y^6
              +(15*x^8+100*x^7+350*x^6+300*x^5-2250*x^4-6000*x^3+3000*x^2)*y^4
              +(6*x^10+50*x^9+275*x^8+1000*x^7+2100*x^6+1200*x^5-3000*x^4)*y^2
              +x^12+10*x^11+75*x^10+450*x^9+2225*x^8+9200*x^7+32200*x^6
              +96000*x^5+240000*x^4+480000*x^3+720000*x^2+720000*x]]]

以下の出力は長い為,省略しました.

(%i40) [fpprec,%ez,ratepsilon]: [32,1.0b-6,1.0b-32];
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.
(%o40) [32,1.0b-6,1.0b-32]
(%i41) [(qvpeds ([],[x,y,z],0,h1,r11,0 ),
         qe( bfpcad(ext( '(  (y-1)*z^4+x*z^3+x*(1-y)*z^2+(y-x-1)*z+y=0   ) ))),%pp)];
Evaluation took 163.0310 seconds (163.1330 elapsed) using 57799.604 MB.
(%o41) [[[x-6,x,x+1,x+2,3*x+8,8*x+3,4*x^2-27,2*x^3-8*x-9,2*x^3-24*x^2+72*x+27,
         8*x^3+7*x^2-42*x-36,x^4+38*x^3+18*x^2-156*x-108,
         2*x^4+47*x^3-124*x-72,4*x^4+x^3-16*x^2-22*x-6,
         4*x^4+21*x^3-4*x^2-36*x-18,5*x^4+26*x^3-9*x^2-54*x-27,
         8*x^4+96*x^3-416*x^2+201*x+282,16*x^4-95*x^3+24*x^2+704*x+512,
         16*x^4+4*x^3-128*x^2-144*x+229,44*x^4+52*x^3-391*x^2-318*x+303,
         4*x^5-88*x^4+645*x^3-1854*x^2+1920*x-1024,
         8*x^5+281*x^4+146*x^3-1806*x^2-1404*x+363,
         15*x^5+278*x^4-30*x^3-1171*x^2-972*x-142,x^6-4*x^5+7*x^3-x^2-3*x+1,
         x^6+46*x^5+142*x^4-74*x^3-426*x^2-342*x-81,
         x^6+78*x^5+556*x^4-208*x^3-1898*x^2-1608*x-405,
         16*x^6+4*x^5-128*x^4-272*x^3+85*x^2+485*x+256,
         192*x^9+1952*x^8-7786*x^7-18888*x^6+41108*x^5+45845*x^4+11523*x^3
                -47238*x^2-78912*x-104832,
         x^10-36*x^9+540*x^8-4212*x^7+16875*x^6-23679*x^5-43659*x^4+114669*x^3
             +96228*x^2-58320*x-46656,
         4*x^11-24*x^10-640*x^9+8377*x^8-38996*x^7+57424*x^6+130396*x^5
               -397672*x^4-88296*x^3+476802*x^2+291240*x+52488,
         2304*x^11-68736*x^10+753424*x^9-2747520*x^8-12830688*x^7
                  +156357720*x^6-532952352*x^5+497201760*x^4+874470897*x^3
                  -1352650752*x^2-339738624*x+905969664,
         192*x^13-1408*x^12-7536*x^11+72256*x^10-85432*x^9-418072*x^8
                 +1080992*x^7+265376*x^6-3323600*x^5+1586496*x^4+3750759*x^3
                 -1186326*x^2-1959552*x-373248,
         16*x^15-640*x^14+11060*x^13-106516*x^12+610680*x^11-1996383*x^10
                +2610948*x^9+4680808*x^8-20830994*x^7+14014728*x^6
                +32414544*x^5-20223465*x^4-40399848*x^3-20947032*x^2-5792544*x
                -944784,
         4096*x^19-76032*x^18-122880*x^17+10690448*x^16-84608248*x^15
                  +268079584*x^14-170213648*x^13-1481638472*x^12
                  +4738394744*x^11-1497816385*x^10-17768059852*x^9
                  +18831296016*x^8+43729432710*x^7-32989362828*x^6
                  -76761910980*x^5+8582283771*x^4+65151311040*x^3
                  +16307847936*x^2-14837686272*x-5792808960,
         4096*x^20-163840*x^19+2852864*x^18-27393408*x^17+149542144*x^16
                  -388286848*x^15-214815103*x^14+4154362488*x^13
                  -7635762592*x^12-10267177106*x^11+44575897180*x^10
                  -943203000*x^9-117037926820*x^8+32371882616*x^7
                  +192201053844*x^6-23391280800*x^5-205468815744*x^4
                  -49971769344*x^3+106333876224*x^2+77778911232*x+16052649984,
         1024*x^25-36864*x^24+258048*x^23+9652320*x^22-275282941*x^21
                  +3486936746*x^20-25967275698*x^19+114334396176*x^18
                  -227290439617*x^17-360243822520*x^16+3107900678073*x^15
                  -4769190186712*x^14-9207670651692*x^13+33891866135830*x^12
                  +3248200641279*x^11-103251848335494*x^10+24822591383424*x^9
                  +196171087421376*x^8-26483105274624*x^7-245526417483264*x^6
                  -55109356535808*x^5+147054424817664*x^4+90730087317504*x^3
                  -5655119265792*x^2-16731447754752*x-3522410053632,
         63488*x^29-3822208*x^28+101436112*x^27-1558519760*x^26
                   +15218355663*x^25-96366411132*x^24+374174103239*x^23
                   -650077067955*x^22-1132725702993*x^21+8060985767681*x^20
                   -9709931306620*x^19-26781342092492*x^18+67248318526264*x^17
                   +73969079197117*x^16-283280225707950*x^15
                   -213264871269207*x^14+938560561534635*x^13
                   +452329212683124*x^12-1991414684849397*x^11
                   -890275611320538*x^10+2384378846393901*x^9
                   +1407423445298100*x^8-1270657944171792*x^7
                   -1145677558363200*x^6+23176788277248*x^5
                   +301610526289920*x^4+117605153636352*x^3-9039935176704*x^2
                   -19813556551680*x-3962711310336,
         7936*x^36-758304*x^35+33419327*x^34-901346428*x^33+16631964287*x^32
                  -222094789346*x^31+2209871734238*x^30-16568887999552*x^29
                  +92997245415767*x^28-377363024540523*x^27
                  +989223612173741*x^26-867722678128324*x^25
                  -5048534162562257*x^24+24309955173059337*x^23
                  -41581904973045947*x^22-28328921944006469*x^21
                  +277043222823124631*x^20-427140634741112450*x^19
                  -394584630361333315*x^18+2108413813884902713*x^17
                  -1146104394578571495*x^16-4964331326135378292*x^15
                  +6113780241822578916*x^14+7684481244922767486*x^13
                  -13180637957355086364*x^12-8973842356231860348*x^11
                  +17034125906519891856*x^10+7340829016726477440*x^9
                  -16471368758915006976*x^8-5867495981300284416*x^7
                  +12273511269086613504*x^6+7406908509570859008*x^5
                  -2449053403743780864*x^4-2741599838700306432*x^3
                  -226468951385702400*x^2+292162779488452608*x
                  +68475651442606080],
        [y-1,y,y-x-1,8*y^2-16*y+3*x+8,
         (4*x^3-16*x-18)*y^4+((-16*x^3)-14*x^2+84*x+72)*y^3
                            +(x^4+38*x^3+18*x^2-156*x-108)*y^2
                            +((-2*x^4)-47*x^3+124*x+72)*y+4*x^4+21*x^3-4*x^2
                            -36*x-18,
         (16*x^4+4*x^3-128*x^2-144*x+229)*y^6
          +((-88*x^4)-104*x^3+782*x^2+636*x-606)*y^5
          +(8*x^5+281*x^4+146*x^3-1806*x^2-1404*x+363)*y^4
          +((-30*x^5)-556*x^4+60*x^3+2342*x^2+1944*x+284)*y^3
          +(x^6+78*x^5+556*x^4-208*x^3-1898*x^2-1608*x-405)*y^2
          +((-2*x^6)-92*x^5-284*x^4+148*x^3+852*x^2+684*x+162)*y+5*x^6+36*x^5
          +48*x^4-46*x^3-144*x^2-108*x-27],
        [(y-1)*z^4+x*z^3+((-x*y)+x)*z^2+(y-x-1)*z+y]]]
(%i42) [(qvpeds ([],[x,y,z],0,l,r11,0 ),
         qe( bfpcad(ext( '(  (y-1)*z^4+x*z^3+x*(1-y)*z^2+(y-x-1)*z+y=0   ) ))),%pp)];
Evaluation took 8.0460 seconds (8.0510 elapsed) using 3347.636 MB.
(%o42) [[[x,x+1,4*x^2-27,5*x^4+26*x^3-9*x^2-54*x-27,
         16*x^4-95*x^3+24*x^2+704*x+512,16*x^4+4*x^3-128*x^2-144*x+229,
         x^6-4*x^5+7*x^3-x^2-3*x+1,
         x^10-36*x^9+540*x^8-4212*x^7+16875*x^6-23679*x^5-43659*x^4+114669*x^3
             +96228*x^2-58320*x-46656],
        [y-1,y,
         (16*x^4+4*x^3-128*x^2-144*x+229)*y^6
          +((-88*x^4)-104*x^3+782*x^2+636*x-606)*y^5
          +(8*x^5+281*x^4+146*x^3-1806*x^2-1404*x+363)*y^4
          +((-30*x^5)-556*x^4+60*x^3+2342*x^2+1944*x+284)*y^3
          +(x^6+78*x^5+556*x^4-208*x^3-1898*x^2-1608*x-405)*y^2
          +((-2*x^6)-92*x^5-284*x^4+148*x^3+852*x^2+684*x+162)*y+5*x^6+36*x^5
          +48*x^4-46*x^3-144*x^2-108*x-27],
        [(y-1)*z^4+x*z^3+((-x*y)+x)*z^2+(y-x-1)*z+y]]]

2017-07-29

Lazard's method(その1)

1994年(Unpublished manuscript は 1990年)に Daniel Lazard が発表した(https://link.springer.com/chapter/10.1007/978-1-4612-2628-4_29CAD の構成方法は,一般化された根の重複度を保つ分割によるもので,その projection set(主係数,定数項,判別式,終結式のみ)は,当時知られていた Collins,McCallum,Hong タイプの何れよりも小型でしたが,正当性の証明に不備がありました.

そのため,2001年に発表された,より小型の McCallum-Brown projection がこれまでの標準でしたが,その正しさが保証されているは well-oriented という条件を満たす場合であるため,他の場合には Hong タイプと使い分ける,或いは,別途多項式を追加するといった制約がありました.

しかし,ここ数年,Scott McCallum らによる検証(https://www.researchgate.net/project/Validity-proof-of-Lazards-method-for-CAD-construction)が進み,遂に 2016年7月,Lazard's method の正当性の証明が発表されたのです(https://arxiv.org/abs/1607.00264).これは非常に大きな進展であり,本年 6 月の MEGA 2017(https://mega2017.inria.fr/files/2017/06/Parusinski.pdf)に続き,9 月の JARCS 2017(http://www.maths.usyd.edu.au/u/laurent/RCSW/)は今回の研究メンバーが主催者に含まれていることもあり,注目しています.

なお,Lazard' method の lifting では,多項式の係数への sample point の代入を拡張し,sample point を根にもつ 1 次式で割り切れる限り割った商に代入したもの,つまり,sample の変数についての multi-index(Lazard valuation)が辞書式順序で最小である 0 でない偏微分係数の 0 でない定数倍(Lazard evaluation)を用います(実装においては,この処理は多項式が消失する場合のみで十分です).

実行例:d = 0, c = 0 のとき,256*a^2*d^3-128*a*b^2*d^2+(144*a*b*c^2+16*b^4)*d-27*a*c^4-4*b^3*c^2 は消失しますが,上記の処理により,Lazard evaluation -4*b^3(実行例の 72*a*b*d-81*a*c^2-2*b^3 に対応)が得られます.

(%i1) (qvpeds ([all],[d,c,a,b,x],0,l,r11,0 ),
       g1:qe( bfpcad(ext( '(  a*x^4+b*x^2+c*x+d>=0 ) )))  );
Evaluation took 3.6700 seconds (3.6800 elapsed) using 506.869 MB.
(%o1) [[d = root(d,1),c = root(c,1),root(a,1) <= a,
        root(72*a*b*d-81*a*c^2-2*b^3,1) <= b,true],
       [root(d,1) < d,c < root(c,1),root(a,1) <= a,
        root(256*a^2*d^3-128*a*b^2*d^2+(144*a*b*c^2+16*b^4)*d-27*a*c^4
                        -4*b^3*c^2,2)
          <= b,true],
       [root(d,1) < d,c = root(c,1),root(4096*a*d^3+27*c^4,1) <= a,
        root(256*a^2*d^3-128*a*b^2*d^2+(144*a*b*c^2+16*b^4)*d-27*a*c^4
                        -4*b^3*c^2,1)
          <= b,true],
       [root(d,1) < d,root(c,1) < c,root(a,1) <= a,
        root(256*a^2*d^3-128*a*b^2*d^2+(144*a*b*c^2+16*b^4)*d-27*a*c^4
                        -4*b^3*c^2,2)
          <= b,true]]

(%i2) %pp;
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.
(%o2) [[d],[c],[a,256*a*d^3-27*c^4,4096*a*d^3+27*c^4],
       [256*a^2*d^3-128*a*b^2*d^2+(144*a*b*c^2+16*b^4)*d-27*a*c^4-4*b^3*c^2],
       [a*x^4+b*x^2+c*x+d]]

(%i3) laz_eval([[],[d=0,c=0,a=-1/2]],%pp[4]);
Evaluation took 0.0000 seconds (0.0100 elapsed) using 368.906 KB.
(%o3) [[b],[4*(72*a*b*d-81*a*c^2-2*b^3)]]

(%i4) qex([],F2G(g1) %eq qex([[A,x]], a*x^4+b*x^2+c*x+d>=0 )  );
Evaluation took 30.2000 seconds (35.1300 elapsed) using 4898.136 MB.
(%o4) true