package Current.popups.Unfold;
import Current.basic.*;
import Current.*;
import Current.plot.Dyadic.*;

import java.applet.Applet;
import java.awt.event.*;
import java.awt.*;
import java.math.*;

public class Function {
    public Function() {}


    /**Section 1:  Function definitions.

      for a generic word, the functions have a numerator and a denominator.
      One tricky part of this business is getting the sign right.  So,
      we have a separate routine for computing the sign.*/


   public static int[][] computeNumerator(VertexPair V,CombinatorialTriangle[] CT) {
      int[] n=Spine.computeSpine(V,CT);
      int spine=n[n[0]+2];
      int[][] m=new int[4][n[0]+2];
      for(int j=1;j<=n[0];++j) {
        m[1][j]=CT[n[j]].d[spine][0];
        m[2][j]=CT[n[j]].d[spine][1];
        m[3][j]=(2*(j%2)-1)*n[n[0]+1];
     }
      m[0][0]=n[0];      //length
      m[0][1]=n[1];      //start
      return(m);
   }


    public static int[][] computeDenominator(int spine,CombinatorialTriangle[] CT) {
      int[] n=Spine.fullSpine(spine,CT);
      int[][] m=new int[4][n[0]+1];
      for(int j=1;j<=n[0];++j) {
        m[0][j]=n[j];
        m[1][j]=CT[n[j]].d[spine][0];
        m[2][j]=CT[n[j]].d[spine][1];
        m[3][j]=2*(j%2)-1;
      }
      m[0][0]=n[0];
      return(m);
    }


    public static int[][] computeDenominator(VertexPair V,CombinatorialTriangle[] CT) {
	int spine=Converter.spineNumber(V,CT);
	return(computeDenominator(spine,CT));
    }


    public static int getSign(int nn,int[][] d) {
       if(d[0][d[0][0]]<nn) return(0);   //a special case

       if(d[0][d[0][0]]==nn) return((int)(Math.pow(-1,nn)));
       int count=1;
       while(d[0][count]<nn) ++count;

       if(d[0][count]!=nn) --count;
       if(count%2==0) return(1);
       return(-1);
    }

    
    public static int[][] getSignedNumerator(int[][] d,VertexPair V,CombinatorialTriangle[] CT) {
       int[][] n=computeNumerator(V,CT); 
       int sign=getSign(n[0][1],d);
       for(int i=0;i<=n[0][0];++i) n[3][i]=sign*n[3][i];
       return(n);
    }





    /**Section 2: Function Computations.   We have set things up so
       that one can take arbitrarily high partial derivatives */


    //McScaled during routine

  public static Complex derivative(int[][] g,int a,int b,Complex zz) {
       Complex total=new Complex();
       Complex z=new Complex(Math.PI*zz.x/2.0,Math.PI*zz.y/2.0);
       int cong=(a+b)%4;
       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 phase=g[1][i]*z.x+g[2][i]*z.y;
	   double C=coeff*Math.cos(phase);
	   double S=coeff*Math.sin(phase);

	   if(cong==0) {
	       total.x=total.x+C;
	       total.y=total.y+S;
	   }
	   if(cong==1) {
	       total.x=total.x-S;
	       total.y=total.y+C;
	   }
	   if(cong==2) {
	       total.x=total.x-C;
	       total.y=total.y-S;
	   }
	   if(cong==3) {
	       total.x=total.x+S;
	       total.y=total.y-C;
	   }
       }
       return(total);
  }


    /**This routine is used in the product formula for partial derivatives.
       If you ever write the chain rule for the product rule in 2
       variables, you will see that it requires a routine like the
       following one.*/

    public static int[] partialContent(int k,int[] a) {
	int b=a[0]+a[1];
	int kk=k;
	int[] c=new int[2];
	c[0]=0;
	c[1]=0;
	for(int i=0;i<b;++i) {
	    int bit=kk&1;
	    if((bit==1)&&(i<a[0])) ++c[0];
	    if((bit==1)&&(i>=a[0])) ++c[1];
	    kk=kk>>1;
	}
	    return(c);
    }


    public static double termInProductRule(int[][] n,int[][] d,int k,int[] a,Complex z) {
	int[] b=partialContent(k,a);
	int[] c=new int[2];
	c[0]=a[0]-b[0];
	c[1]=a[1]-b[1];
	Complex z1=derivative(n,b[0],b[1],z);
	Complex z2=derivative(d,c[0],c[1],z);
	Complex z3=Complex.times(z1,Complex.conjugate(z2));
	return(z3.y);
    }

    /**Here is the routine which evaluates arbitrary partial derivatives of
       the function F=Im(PQ_bar)*/

    public static double evaluate(int[][] n,int[][] d,int a0,int a1,Complex z) {
	double w=0;
	int[] a=new int[2];
	a[0]=a0;
	a[1]=a1;
	for(int k=0;k<Math.pow(2,a0+a1);++k) {
	    w=w+termInProductRule(n,d,k,a,z);
	}
	return(w);
    }






    /**Section 3:  Foiling */

    /*In one place in our plotting algorithm we need to take the
      product of the numerator and the (conjugate of the) denominator.
      We do this using the foil method, which is slow but reliable.
      We don't know a better way to get the bounds which use this
      method*/
    

public static int[][] foil(int[][] f1,int[][] f2) {
    int i,j,imax,jmax,sign,v1,v2,count;

    /**computing the range of frequencies.  This determines the array size*/

  int x1min,x1max,y1min,y1max,x2min,x2max,y2min,y2max;

  x1min=10000000;
  x1max=-10000000;
  y1min=10000000;
  y1max=-10000000;
  x2min=10000000;
  x2max=-10000000;
  y2min=10000000;
  y2max=-10000000;

  for(int ii=1;ii<=f1[0][0];++ii) {
      if(x1min>f1[1][ii]) x1min=f1[1][ii];
      if(x1max<f1[1][ii]) x1max=f1[1][ii];
      if(y1min>f1[2][ii]) y1min=f1[2][ii];
      if(y1max<f1[2][ii]) y1max=f1[2][ii];
  }

  for(int ii=2;ii<=f2[0][0];++ii) {
      if(x2min>f2[1][ii]) x2min=f2[1][ii];
      if(x2max<f2[1][ii]) x2max=f2[1][ii];
      if(y2min>f2[2][ii]) y2min=f2[2][ii];
      if(y2max<f2[2][ii]) y2max=f2[2][ii];
  }

  int xoffset1=Math.max(Math.abs(x1min),Math.abs(x1max));
  int xoffset2=Math.max(Math.abs(x2min),Math.abs(x2max));
  int xoffset=xoffset1+xoffset2+1;

  int yoffset1=Math.max(Math.abs(y1min),Math.abs(y1max));
  int yoffset2=Math.max(Math.abs(y2min),Math.abs(y2max));
  int yoffset=yoffset1+yoffset2+1;

  int[][] a=new int[xoffset+1][2*yoffset+1];

  /**done computing the range of frequencies*/


  for(i=0;i<=xoffset;++i) {
    for(j=0;j<=2*yoffset;++j) {
      a[i][j]=0;
    }}

  imax=-2*xoffset;
  jmax=0;
  int freq=0;

  for(i=1;i<=f1[0][0];++i) {
    for(j=1;j<=f2[0][0];++j) {
      sign=f1[3][i]*f2[3][j];
      v1=f1[1][i]-f2[1][j];
      v2=f1[2][i]-f2[2][j];
      if((v1<=0)&&(imax<-v1)) imax=-v1;
      if((v1>=0)&&(imax<v1)) imax=v1;

      if(v1>=0) {
	freq=yoffset+v2;
	a[v1][freq]=a[v1][freq]+sign;
	if(jmax<freq) jmax=freq;


      } 
      if(v1<0) {
	freq=yoffset-v2;
	a[-v1][freq]=a[-v1][freq]-sign;
	if(jmax<freq) jmax=freq;
      }
    }
  }

  a[0][yoffset]=0;
    count=1;
    int big=imax*jmax+1;
    int[][] f3=new int[4][big];

    for(i=0;i<=imax;++i) {
      for(j=0;j<=jmax;++j) {

	if(a[i][j]>0) {
	  f3[0][count]=i;
	  f3[1][count]=j-yoffset;
	  f3[2][count]=a[i][j];
	  ++count;
	}
	if(a[i][j]<0) {
          f3[0][count]=-i;
	  f3[1][count]=-j+yoffset;
	  f3[2][count]=-a[i][j];
	  ++count;
	}
      }
    }
    f3[0][0]=count-1;
    return(f3);
}


    /**This is the bound which requires the foil method*/

  public static int absoluteSecondDerivative(int[][] f,int n1,int n2) {
    int a,b,c,d;
    int total=0;
    d=0;
    for(int i=1;i<=f[0][0];++i) {
      a=f[0][i];
      b=f[1][i];
      c=f[2][i];
      if((n1==2)&&(n2==0)) d=c*a*a;
      if((n1==1)&&(n2==1)) d=c*a*b;
      if((n1==0)&&(n2==2)) d=c*b*b;
      if(d<0) d=-d;
      total=total+d;
    }
    return(total);
  }


    /**Here is a foil-based routine for computing the partial derivatives*/

    public static double foilDerivative(int[][] g,int a,int b,Complex zz) {
	double total=0;
	double X=zz.x;
	double Y=zz.y;
        int cong=(a+b)%4;
        for(int i=1;i<=g[0][0];++i) {  
	   int c1=g[2][i];
	   int c2=g[0][i];
           int c3=g[1][i];
	   double c4=Math.pow(c2,a);
	   double c5=Math.pow(c3,b);
           double coeff=c1*c4*c5;
	   double phase=(c2*zz.x+c3*zz.y)*Math.PI/2;
	   double C=coeff*Math.cos(phase);
	   double S=coeff*Math.sin(phase);

           if(cong==0) {
	       total=total+S;
	   }
	   if(cong==1) {
	       total=total+C;
	   }
	   if(cong==2) {
	       total=total-S;
	   }
	   if(cong==3) {
	       total=total-C;
	   }
	}
       return(total);
    }







    /**Section 4: precomputing for the plotter.
       when we do our plotting algorithm we will have the denominator
       functions precomputed.  This routine evaluates these functions
       at some commonly sampled points.*/

    public static Complex[] assignData(int[][]f,DyadicSquare Q) {
    Complex z=Q.getCenter();
    int diam=Q.k;
    Complex[] zz=new Complex[10];
    int pos1,pos2;
    double max,min,width;
    for(int i=1;i<=8;++i) zz[i]=new Complex();
    width=Math.pow(0.5,diam+1);
    zz[1].x=z.x+width;
    zz[1].y=z.y+width;
    zz[2].x=z.x;
    zz[2].y=z.y+2.0*width;
    zz[3].x=z.x-width;
    zz[3].y=z.y+width;
    zz[4].x=z.x-2.0*width;
    zz[4].y=z.y;
    zz[5].x=z.x-width;
    zz[5].y=z.y-width;
    zz[6].x=z.x;
    zz[6].y=z.y-2.0*width;
    zz[7].x=z.x+width;
    zz[7].y=z.y-width;
    zz[8].x=z.x+2.0*width;
    zz[8].y=z.y;
    Complex[] Z=new Complex[11];
    Z[0]=derivative(f,0,0,z);
    for(int i=1;i<=8;++i) Z[i]=derivative(f,0,0,zz[i]);
    Z[9]=derivative(f,1,0,z);
    Z[10]=derivative(f,0,1,z);
    return(Z);
    }


}