package Current.popups.Deg100; import Current.*; import Current.basic.*; import java.applet.Applet; import java.awt.event.*; import java.awt.*; import java.math.*; public class Deg100Compute { public Deg100Compute() {} /*****CERFIFICATE*****/ /*the goal of this routine is to tell if the gradient of a function, as it varies over a square in the parameter space, lies in a quadrant. The values (grad1,grad2) give the gradient at the center of the point. The terms (b1,b2) are bounds on the 2nd derivatives. These bounds guarantee that: grad1-bb1*r < grad1(x) 0)&&(gr1-a1>0)) x[1]=1; if((gr2+a2>0)&&(gr2-a2>0)) x[2]=1; if((gr1+a1<0)&&(gr1-a1<0)) x[3]=1; if((gr2+a2<0)&&(gr2-a2<0)) x[4]=1; if((x[1]>0)&&(x[2]>0)) return(1); if((x[2]>0)&&(x[3]>0)) return(3); if((x[3]>0)&&(x[4]>0)) return(5); if((x[4]>0)&&(x[1]>0)) return(7); //standard quadrant test has failed. //now we look at a square twice as //large and try to show that the gradient //lies in some quadrant. step=Math.pow(2.0,k-1); //twice bound on sq.radius gr1=grad1*step; gr2=grad2*step; a1=Math.max(a[0],a[1]); a2=Math.max(a[1],a[2]); x[1]=gr1+a1; x[2]=gr2+a2; x[3]=gr1-a1; x[4]=gr2-a2; x[5]=Math.abs(x[1]); if(x[5]x[5])&&( x[4]>x[5])) return(2); if((-x[1]>x[6])&&(-x[3]>x[6])) return(4); if((-x[2]>x[5])&&(-x[4]>x[5])) return(6); if(( x[1]>x[6])&&( x[3]>x[6])) return(8); return(0); } //rigorous version public static int certifyAlgebra2(BigInteger grad1,BigInteger grad2,int[] a,int k) { if(k<=1) return(0); int[] x=new int[5]; BigInteger[] y=new BigInteger[7]; BigInteger step=new BigInteger("0"); BigInteger a1=new BigInteger("0"); BigInteger a2=new BigInteger("0"); BigInteger gr1=new BigInteger("0"); BigInteger gr2=new BigInteger("0"); /**convert the a integers to Big integers. We multiply by 2^106 because all the other terms are thus inflated */ Integer A0=new Integer(a[0]); Integer A1=new Integer(a[1]); Integer A2=new Integer(a[2]); BigInteger aa0=new BigInteger(A0.toString()); BigInteger aa1=new BigInteger(A1.toString()); BigInteger aa2=new BigInteger(A2.toString()); BigInteger BIG=new BigInteger("81129638414606681695789005144064"); aa0=aa0.multiply(BIG); aa1=aa1.multiply(BIG); aa2=aa2.multiply(BIG); //test standard quadrant: step=new BigInteger("2"); step=step.pow(k); //anything larger then 2^{k+2}/Pi is legit. gr1=grad1.multiply(step); gr2=grad2.multiply(step); a1=aa0.add(aa1); a2=aa1.add(aa2); x[1]=0;x[2]=0;x[3]=0;x[4]=0; if((gr1.compareTo(a1.negate())==1)&&(gr1.compareTo(a1)==1)) x[1]=1; if((gr2.compareTo(a2.negate())==1)&&(gr2.compareTo(a2)==1)) x[2]=1; if((gr1.compareTo(a1.negate())==-1)&&(gr1.compareTo(a1)==-1)) x[3]=1; if((gr2.compareTo(a2.negate())==-1)&&(gr2.compareTo(a2)==-1)) x[4]=1; if((x[1]>0)&&(x[2]>0)) return(1); if((x[2]>0)&&(x[3]>0)) return(3); if((x[3]>0)&&(x[4]>0)) return(5); if((x[4]>0)&&(x[1]>0)) return(7); //standard quadrant test has failed. //now we look at a square twice as //large and try to show that the gradient //lies in some quadrant. step=new BigInteger("2"); step=step.pow(k-1); gr1=grad1.multiply(step); gr2=grad2.multiply(step); a1=aa0.max(aa1); a2=aa1.max(aa2); y[1]=gr1.add(a1); y[2]=gr2.add(a2); y[3]=gr1.subtract(a1); y[4]=gr2.subtract(a2); y[5]=y[1].abs(); y[5]=y[5].max(y[3].abs()); y[6]=y[2].abs(); y[6]=y[6].max(y[4].abs()); if((y[2].compareTo(y[5])==1)&&(y[4].compareTo(y[5])==1)) return(2); if((y[6].compareTo(y[1].negate())==-1)&&(y[6].compareTo(y[3].negate())==-1)) return(4); if((y[5].compareTo(y[2].negate())==-1)&&(y[5].compareTo(y[4].negate())==-1)) return(6); if((y[1].compareTo(y[6])==1)&&(y[3].compareTo(y[6])==1)) return(8); return(0); } /******FUNCTION EVALUATION*************/ /**this routine selects the appropriate point relative to a dyadic square. The relevant point is chosen by the certificate above. Returns a list of dyadic rationals: Z[0],Z[1] are the coordinates for the first point Z[2],Z[3] are the coordinates for the second point */ public static DyadicRational[] getRelevant(DyadicSquare Q,int direction) { DyadicRational width=new DyadicRational(1,Q.k+1); DyadicRational width2=new DyadicRational(1,Q.k); DyadicRational x=new DyadicRational(Q.x[0],Q.x[1]); int yy=(int)(Math.pow(2,Q.y[1]))-Q.y[0]; DyadicRational y=new DyadicRational(yy,Q.y[1]); DyadicRational[] Z=new DyadicRational[4]; if(direction==1) { Z[0]=DyadicRational.plus(x,width); Z[1]=DyadicRational.plus(y,width); Z[2]=DyadicRational.minus(x,width); Z[3]=DyadicRational.minus(y,width); } if(direction==5) { Z[2]=DyadicRational.plus(x,width); Z[3]=DyadicRational.plus(y,width); Z[0]=DyadicRational.minus(x,width); Z[1]=DyadicRational.minus(y,width); } if(direction==2) { Z[0]=x; Z[1]=DyadicRational.plus(y,width2); Z[2]=x; Z[3]=DyadicRational.minus(y,width2); } if(direction==6) { Z[2]=x; Z[3]=DyadicRational.plus(y,width2); Z[0]=x; Z[1]=DyadicRational.minus(y,width2); } if(direction==3) { Z[0]=DyadicRational.minus(x,width); Z[1]=DyadicRational.plus(y,width); Z[2]=DyadicRational.plus(x,width); Z[3]=DyadicRational.minus(y,width); } if(direction==7) { Z[2]=DyadicRational.minus(x,width); Z[3]=DyadicRational.plus(y,width); Z[0]=DyadicRational.plus(x,width); Z[1]=DyadicRational.minus(y,width); } if(direction==4) { Z[0]=DyadicRational.minus(x,width2); Z[1]=y; Z[2]=DyadicRational.plus(x,width2); Z[3]=y; } if(direction==8) { Z[2]=DyadicRational.minus(x,width2); Z[3]=y; Z[0]=DyadicRational.plus(x,width2); Z[1]=y; } return(Z); } //numerical version public static Complex derivative(int[][] g,int a,int b,DyadicRational X,DyadicRational Y) { Complex total=new Complex(); Complex zz=new Complex(X.n/Math.pow(2,X.d),Y.n/Math.pow(2,Y.d)); Complex z=new Complex(Math.PI*zz.x/2.0,Math.PI*zz.y/2.0); int cong=(a+b)%4; //only use cong=0 or 1. for(int i=1;i<=g[0][0];++i) { double coeff=g[3][i]*Math.pow(g[1][i],a)*Math.pow(g[2][i],b); double arg=g[1][i]*z.x+g[2][i]*z.y; if(cong==0) { total.x=total.x+coeff*Math.cos(arg); total.y=total.y+coeff*Math.sin(arg); } if(cong==1) { total.x=total.x-coeff*Math.sin(arg); total.y=total.y+coeff*Math.cos(arg); } } return(total); } //rigorous version public static BigComplexInterval derivative2(int[][] g,int a,int b,DyadicRational X,DyadicRational Y) { BigInterval TX=new BigInterval(); BigInterval TY=new BigInterval(); TX.L=new BigInteger("0"); TX.R=new BigInteger("0"); TY.L=new BigInteger("0"); TY.R=new BigInteger("0"); int cong=(a+b)%4; //only use cong=0 or 1. DyadicRational ARG1=new DyadicRational(0,0); DyadicRational ARG2=new DyadicRational(0,0); DyadicRational ARG3=new DyadicRational(0,0); BigInterval SIN=new BigInterval(); BigInterval COS=new BigInterval(); for(int i=1;i<=g[0][0];++i) { BigInterval COEFF=new BigInterval(g[3][i]*(int)(Math.pow(g[1][i],a))*(int)(Math.pow(g[2][i],b))); ARG1=new DyadicRational(g[1][i]*X.n,X.d); ARG2=new DyadicRational(g[2][i]*Y.n,Y.d); ARG3=DyadicRational.plus(ARG1,ARG2); COS=Deg100Trig.intervalCos(ARG3.n,ARG3.d); SIN=Deg100Trig.intervalSin(ARG3.n,ARG3.d); if(cong==0) { TX=BigInterval.add(TX,BigInterval.multiply(COEFF,COS)); TY=BigInterval.add(TY,BigInterval.multiply(COEFF,SIN)); } if(cong==1) { TX=BigInterval.subtract(TX,BigInterval.multiply(COEFF,SIN)); TY=BigInterval.add(TY,BigInterval.multiply(COEFF,COS)); } } return(new BigComplexInterval(TX,TY)); } //numerical version public static double evaluate(int[][] n,int[][] d,DyadicRational X,DyadicRational Y) { Complex E=derivative(n,0,0,X,Y); Complex Q=derivative(d,0,0,X,Y); E=Complex.times(E,Complex.conjugate(Q)); return(E.y); } //rigorous version public static BigInterval evaluate2(int[][] n,int[][] d,DyadicRational X,DyadicRational Y) { BigComplexInterval E=derivative2(n,0,0,X,Y); BigComplexInterval Q=derivative2(d,0,0,X,Y); E=BigComplexInterval.times(E,BigComplexInterval.conjugate(Q)); return(E.y); } //numerical version public static int signEvaluate(int[][] n,int[][] d,DyadicSquare Q,int direction) { if(direction==0) return(0); DyadicRational[] DR=getRelevant(Q,direction); double min=evaluate(n,d,DR[2],DR[3]); if(min>=0) return(2); double max=evaluate(n,d,DR[0],DR[1]); if(max<=0) return(1); return(0); } //rigorous version public static int signEvaluate2(int[][] n,int[][] d,DyadicSquare Q,int direction) { if(direction==0) return(0); BigInteger ZILCH=new BigInteger("0"); DyadicRational[] DR=getRelevant(Q,direction); BigInterval min=evaluate2(n,d,DR[2],DR[3]); if(min.L.compareTo(ZILCH)==1) return(2); BigInterval max=evaluate2(n,d,DR[0],DR[1]); if(max.R.compareTo(ZILCH)==-1) return(1); return(0); } //numrical public static double gradient(int q,int[][] n,int[][] d,DyadicRational X,DyadicRational Y) { Complex P=derivative(n,0,0,X,Y); Complex DP=new Complex(); if(q==1) DP=derivative(n,1,0,X,Y); if(q==2) DP=derivative(n,0,1,X,Y); Complex Q=derivative(d,0,0,X,Y); Complex DQ=new Complex(); if(q==1) DQ=derivative(d,1,0,X,Y); if(q==2) DQ=derivative(d,0,1,X,Y); Complex w1=Complex.times(P,Complex.conjugate(DQ)); Complex w2=Complex.times(DP,Complex.conjugate(Q)); Complex w=Complex.plus(w1,w2); return(w.y); } //rigorous public static BigInterval gradient2(int q,int[][] n,int[][] d,DyadicRational X,DyadicRational Y) { BigComplexInterval P=derivative2(n,0,0,X,Y); BigComplexInterval DP=new BigComplexInterval(); if(q==1) DP=derivative2(n,1,0,X,Y); if(q==2) DP=derivative2(n,0,1,X,Y); BigComplexInterval Q=derivative2(d,0,0,X,Y); BigComplexInterval DQ=new BigComplexInterval(); if(q==1) DQ=derivative2(d,1,0,X,Y); if(q==2) DQ=derivative2(d,0,1,X,Y); BigComplexInterval w1=BigComplexInterval.times(P,BigComplexInterval.conjugate(DQ)); BigComplexInterval w2=BigComplexInterval.times(DP,BigComplexInterval.conjugate(Q)); BigComplexInterval w=BigComplexInterval.plus(w1,w2); return(w.y); } public static void inflatePrint(double d) { BigDecimal D=new BigDecimal(d); D=D.multiply(new BigDecimal(new BigInteger("81129638414606681695789005144064"))); BigInteger DI=D.toBigInteger(); System.out.println(DI.toString()); } }