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$