/*
Generated by the Cosmo Code Visual Builder, Version 2.3.
*/

import java.awt.*;
import java.applet.*;
import java.lang.*;

import Complex;
//since Complex numbers are not part of the Java 1.1 standard libraries. 
//I used Alexander Anderson's fine class http://www.netlib.no/netlib/java/ (see
//http://www.almide.demon.co.uk/html/Complex/Complex_for_Java.html )
//except the complicated path ORG/... the file is unchanged

//{{ import_classes - Begin Non-Editable Code Block
//}} import_classes - End Non-Editable Code Block

/**
* Put the class description here.
* @author Raymond Manzoni
* @version 0.6 (04 nov 01)
**/
public class critical extends java.awt.Canvas {

    //public DecimalFormat dF;

    /**
    * Initialize this object.
    **/
    public critical() {
        //{{ init_applet - Begin Non-Editable Code Block
        setSize(680, 523);
        //}} init_applet - End Non-Editable Code Block

        //{{ init_objects - Begin Non-Editable Code Block

        //}} init_objects - End Non-Editable Code Block

        //Format init
        //dF = DecimalFormat.getNumberInstance();  //dF.format(
    }

    //replace Double.toString by formatDouble for better control
    public String formatDouble(double d, int N)
    {
       String res= null;
       try {
          Double db=new Double(d);
          res=db.toString();
          if (res.length()>N)
          {
             int ex=res.indexOf('E');
             int pt=res.indexOf('.');
             if (ex>=0) {
                int l=N-res.length();
                int st= N-(res.length()-ex);
                if (st>=pt)
                   res= res.substring(0,st)+res.substring(ex,res.length());
                else
                   res= res.substring(0,pt)+res.substring(ex,res.length());
             }
             if (pt>=N) //cast to integer
                res=Integer.toString(db.intValue());
             else
                res=res.substring(0,N);
          }
       }   
       catch (Exception e) {
          res="1.0";
       } 
       return res;
    }

    public void update(Graphics g) {	
     /* don't erase background
	    g.setColor(getBackground());
	    g.fillRect(0, 0, width, height);
	    g.setColor(getForeground());
 	  */
      paint(g);
    }

    /**
    * Handle an event on an object in this class.
    * @param event The event which occurred.
    * @return true if the event is handled, false if not
    **/
    public boolean handleEvent(Event event) {
        Object arg = event.arg;
        boolean eventHandled = false;

        //{{ wiring - Begin Non-Editable Code Block

        //}} wiring - End Non-Editable Code Block

        if (!eventHandled) {
            super.handleEvent(event);
        }

        // always return false so that others may handle the event as well
        return false;
    }

    //Complete redraw of the canvas rectangle
    public void updateAll() { 
       repaint();
    }   
    
    public void ScaledLine(Graphics gr, int sx, int sy, int dx, int dy)
    {
       gr.drawLine((int)(G_m*sx+G_ox), (int)(G_oy-G_m*sy)
                  ,(int)(G_m*dx+G_ox), (int)(G_oy-G_m*dy));      
    }
    
	//Clear frame and draw x=0, x=1 lines, ...
	public void initScr(Graphics g, Complex z, int choice, double outScaling) {	
      
	   Dimension d = size();
       g.setColor(getBackground());
	   g.fillRect(0, 0, d.width, d.height);
	   g.setColor(getForeground());
	   
      // background lines
       ScaledLine(g, ((int)-outScaling)-1,0, ((int)outScaling)+1,0); //real axis
       ScaledLine(g, 0,((int)-outScaling)-1, 0,((int)outScaling)+1); //imaginary axis
       ScaledLine(g, 1,((int)-outScaling)-1, 1,((int)outScaling)+1); //Re(z)=1 line
     
       g.drawString("X="+ formatDouble(z.re(),10)+ " Y="+formatDouble(z.im(),9),250,460);
 	}

    //restore earlier computed values
    public void cacheRestore(int choice, double inScaling, double outScaling, 
                             double start, double offset) {
       cache_first=0; cache_last=cache_nb-1;
       if (cache_choice==choice && cache_inScaling==inScaling && cache_outScaling==outScaling 
           && cache_offset==offset) {
          int posStart= Math.round((float)((cache_start-start)/offset));
          if (posStart>0) {
             if (posStart<cache_nb) {
                for (int i=cache_nb-1; i>=posStart; --i) {
                   xP[i]= xP[i-posStart];
                   yP[i]= yP[i-posStart];
                }
                cache_last=posStart-1; //0..posStart to fill
             }
          }
          else if (posStart<0) {
             if (posStart+cache_nb>0) {
                for (int i=0; i<cache_nb+posStart; i++) {
                   xP[i]= xP[i-posStart];
                   yP[i]= yP[i-posStart];
                }
                cache_first=cache_nb+posStart; //cache_nb+posStart..nb-1 to fill
             }
          }
          else { //same picture
             cache_last= -1;
          }
       }
    }
    

   /*============================================================*
    *                 Mathematical   Part                        *
    *============================================================*/
   
    
/* Eratosthenes' algorithm for identifying prime numbers

    Written by:   Peter Holcroft
                  University of Hull, United Kingdom
    Contact:      peter@holcroft.cjb.net
    Date written: 27/3/00
    Updated:      12/9/00 (it now works! Sorry)
    Modified for zeta:20/10/01 R.M.
    
    This isn't meant to be an all singing all dancing prime number finder. It shows you one
    way of implemeting Eratosthenes' algorithm in Java and has been deliberately kept simple.
    Feel free to do whatever you want with this source code.

    Just change "max" to set the size of N, and "startat" for setting where to start 
    displaying the prime numbers from.
 */
 
    static final double ROOT_1=14.13472514d;
    
    static final int max=20000;
    static byte[] PrimeOffs=null;
    static int primeCount=0;
 
    public void Sieve(int n) {

       // Create array and set all values to true
       boolean numberpool [] = new boolean [n];
       for (int i=2; i<n; i++)
          numberpool[i]=true;
       int limit = (int) (Math.sqrt(n)+1);         // Find square root of n
       for (int j=3; j<=limit; j+=2)               // do up to sqrt of n
          if (numberpool[j])              // only do if j is a prime
             for (int i=j+j; i<n; i+=j)        // start with 2j as j is prime
                numberpool[i]=false;             // set all multiples to false
       //convert it in prime offset table
       primeCount=1;
       for (int j=3; j<n; j+=2)              // do up to sqrt of n
          if (numberpool[j])              // only do if j is a prime
             primeCount++;
       PrimeOffs= new byte[primeCount];
       PrimeOffs[0]=2;
       int i=0;
       int p=2;
       for (int j=3; j<n; j+=2)              // test every odd number
          if (numberpool[j]) {            // only do if j is a prime
             PrimeOffs[++i]= (byte)(j-p);         // keep the prime offset
             p=j;
          }
    }
    
    // Gamma and lnGamma functions
    // 'Numerical recipes in C' method (10 digits accurate in critical strip)
    // http://lib-www.lanl.gov/numerical/bookcpdf.html
    static final double sqrt2Pi = 2.50662827463100050242d;  //sqrt(2*Pi)
    static final double ln2Pi   = 1.83787706640934548356d;  //ln(2*Pi)
       
    static final double coeff[]= { 1.000000000190015d,
                         76.18009172947146d,
                        -86.50532032941677d,
                         24.01409824083091d,
                         -1.231739572450155d,
                          1.208650973866179E-3d,
                         -5.395239384953E-6d};
    
     public Complex Gamma(Complex zz) {
    
       Complex sum= new Complex(coeff[0]);
       Complex unit= new Complex(1.0);
       Complex z= zz;
       for (int i=1; i<7; i++) {
          sum= sum.add((new Complex(coeff[i])).div(z));
          z = z.add(unit);
       }
       sum= sum.mul(new Complex(sqrt2Pi));
       Complex z45= zz.add(new Complex(4.5d));
       Complex zm05= zz.sub(new Complex(0.5d));
       sum= sum.mul(z45.pow(zm05));
       sum= sum.mul(z45.neg().exp());
       return sum;
    }

    //same reference (log to avoid over and underflow)
    public Complex lnGamma(Complex zz) {
    
       Complex sum= new Complex(coeff[0]);
       Complex unit= new Complex(1.0);
       Complex z= zz;
       for (int i=1; i<7; i++) {
          sum= sum.add((new Complex(coeff[i])).div(z));
          z = z.add(unit);
       }
       sum= sum.mul(new Complex(sqrt2Pi));      
       Complex z45= zz.add(new Complex(4.5d));
       Complex zm05= zz.sub(new Complex(0.5d));
       sum= sum.log();
       sum= sum.add(z45.log().mul(zm05));
       sum= sum.sub(z45);
       return sum;
    }
     
    //============== Riemann-Siegel formula ===============================
    // Ref : H.M. Edwards fine 'Riemann's zeta function' Dover (53f)
    // for C code and details see Glen Pugh's master thesis : 'The Riemann-Siegel 
    // formula and large scale computation of the Riemann zeta function' at :
    // http://www.math.ubc.ca/~pugh (he has an applet too for this function)
    // for more advanced topics see Andrew Odlyzko's papers at 
    // http://www.dtc.umn.edu/~odlyzko/
    
    //set A0  An:= series(exp(2*PI*I*y^2),y);  //to Order 40 for example
    //iterate An:= series(-1/2*y*expr(An)-1/(32*PI^2)*diff(expr(An)/y,y $2),y); //A1,A2..
    //compute series(An*sum(Psi^(m)*y^m/m!,m=0..16),y); 
    //keeping the real part of the constant term should give the n case below
    // (coeff obtained by expanding the taylor series around (z+1)/2)
    
    //Psi(p)= cos(2*pi*(p^2-p-1/16))/cos(2*pi*p)
    // Psi:=func(cos(2*PI*((p^2-p)-1/16))/cos(2*PI*p),p);
    static final int SIZE=23;
    static final double tabC[]= {
    //n=0 case  corresponding to Taylor series of Psi((z+1)/2)
    // float(series(Psi(z/2+1/2),z));
    +.38268343236508977173d, //0 (z^0 coefficient)
    +.43724046807752044936d, //2 (z^2 coefficient)
	+.13237657548034352332d, //4
	-.01360502604767418865d, //6
	-.01356762197010358089d, //8 
	-.00162372532314446528d, //10
	+.00029705353733379691d, //12
	+.00007943300879521470d, //14
	+.00000046556124614505d, //16
	-.00000143272516309551d, //18 
	-.00000010354847112313d, //20
	+.00000001235792708386d, //22
	+.00000000178810838580d, //24 
	-.00000000003391414390d, //26 
	-.00000000001632663390d, //28
	-.000000000000378510932d,//30
	+.0000000000000932742326d,     //32
	+.0000000000000052218430d,     //34
	-.00000000000000033506731d,    //36
	-.000000000000000034124265d,   //38
	+.00000000000000000057512033d, //40
	+.00000000000000000014895301d, //42
	+.000000000000000000001256537d,//44
	//n=1 case  -Psi^(3)(p)/(2^5*3*pi^2)
	// float(-1/(2^5*3*PI^2)*series(subs(diff(Psi(x),x $ 3),x=(z+1)/2),z));
	-.02682510262837534703d, //1
	+.01378477342635185305d, //3 
	+.03849125048223508223d, //5
	+.00987106629906207647d, //7
	-.00331075976085840433d, //9 
	-.00146478085779541508d, //11
	-.00001320794062487696d, //13 
	+.00005922748701847141d, //15
	+.00000598024258537345d, //17
	-.00000096413224561698d, //19
	-.00000018334733722714d, //21
	+.00000000446708756272d, //23 
	+.00000000270963508218d, //25
	+.00000000007785288654d, //27
	-.00000000002343762601d, //29
	-.00000000000158301728d, //31
	+.00000000000012119942d, //33
	+.00000000000001458378d, //35 
	-.00000000000000028786d, //37 28786.305
	-.00000000000000008663d, //39 8662.8629
	-.00000000000000000084d, //41 84.307227
	+.00000000000000000036d, //43 36.308072
	+.00000000000000000001d, //45 1.1626698
	//n=2 case  Psi^(6)(p)/(2^11*3^2*pi^4)+Psi^(2)(p)/(2^6*pi^2)
	// float(1/(2^11*3^2*PI^4)*series(subs(diff(Psi(x),x $ 6),x=(z+1)/2),z)+\
	//       1/(2^6*PI^2)*series(subs(diff(Psi(x),x $ 2),x=(z+1)/2),z));
	+.00518854283029316849d, //0
	+.00030946583880634746d, //2 
	-.01133594107822937338d, //4
	+.00223304574195814477d, //6
	+.00519663740886233021d, //8
	+.00034399144076208337d, //10
	-.00059106484274705828d, //12
	-.00010229972547935857d, //14
	+.00002088839221699276d, //16
	+.00000592766549309654d, //18
	-.00000016423838362436d, //20
	-.00000015161199700941d, //22
	-.00000000590780369821d, //24
	+.00000000209115148595d, //26
	+.00000000017815649583d, //28
	-.00000000001616407246d, //30
	-.00000000000238069625d, //32 
	+.00000000000005398265d, //34 
	+.00000000000001975014d, //36 
	+.00000000000000023333d, //38 
	-.00000000000000011188d, //40
	-.00000000000000000416d, //42 
	+.00000000000000000044d, //44 
	//n=3 case  
	// float(-1/(2^16*3^4*PI^6)*series(subs(diff(Psi(x),x $ 9),x=(z+1)/2),z)\
	//       -1/(2^8*3*5*PI^4)*series(subs(diff(Psi(x),x $ 5),x=(z+1)/2),z)\
	//       -1/(2^6*PI^2)*series(subs(diff(Psi(x),x $ 1),x=(z+1)/2),z));
	-.00133971609071945690d, //1
	+.00374421513637939370d, //3
	-.00133031789193214681d, //5
	-.00226546607654717871d, //7
	+.00095484999985067304d, //9
	+.00060100384589636039d, //11
	-.00010128858286776622d, //13
	-.00006865733449299826d, //15 
	+.00000059853667915386d, //17
	+.00000333165985123995d, //19
	+.00000021919289102435d, //21
	-.00000007890884245681d, //23
	-.00000000941468508130d, //25
	+.00000000095701162109d, //27
	+.00000000018763137453d, //29
	-.00000000000443783768d, //31
	-.00000000000224267385d, //33
	-.00000000000003627687d, //35
	+.00000000000001763981d, //37
	+.00000000000000079608d, //39
	-.00000000000000009420d, //41
	-.00000000000000000713d, //43
	+.00000000000000000033d, //45
	//n=4 case
	// float( 1/(2^23*3^5*PI^8)*series(subs(diff(Psi(x),x $ 12),x=(z+1)/2),z)\
	//       +11/(2^17*3^2*5*PI^6)*series(subs(diff(Psi(x),x $ 8),x=(z+1)/2),z)\
	//       +19/(2^13*3*PI^4)*series(subs(diff(Psi(x),x $ 4),x=(z+1)/2),z)\
	//       +1/(2^7*PI^2)*series(subs(diff(Psi(x),x $ 0),x=(z+1)/2),z));
    +.00046483389361763382d, //0 
	-.00100566073653404708d, //2
	+.00024044856573725793d, //4 
	+.00102830861497023219d, //6 
	-.00076578610717556442d, //8
	-.00020365286803084818d, //10
	+.00023212290491068728d, //12 
	+.00003260214424386520d, //14 
	-.00002557906251794953d, //16
	-.00000410746443891574d, //18
	+.00000117811136403713d, //20 
	+.00000024456561422485d, //22
	-.00000002391582476734d, //24
	-.00000000750521420704d, //26 
	+.00000000013312279416d, //28
	+.00000000013440626754d, //30
	+.00000000000351377004d, //32
	-.00000000000151915445d, //34
	-.00000000000008915418d, //36
	+.00000000000001119589d, //38
	+.00000000000000105160d, //40
	-.00000000000000005179d, //42 
	-.00000000000000000807d};//44 
   // n=5 case
   // -2879/(1769472*PI^4)*Psi^(3) -901/(82575360*PI^6)*Psi^(7) 
   // -7/(849346560*PI^8)*Psi^(11) -1/(978447237120*PI^10)*Psi^(15)
   // will be added for better precision (im(z)<100)
   
   
   
	//C coeff for Riemann-Siegel formula
    public double Coeff(int n, double z) {
    
 	   double p;
	   double s=0.0d;
	   double z2=z*z;
       if (n<5) {
          int st=n*SIZE;
          if ((n&1)!=0)
             p= z;
          else
             p= 1.0;
          int i=SIZE; 
          while (--i>=0) {
             s += tabC[st]*p;
             p *= z2;
             st++;
          }
       }
       return s;
    }

    static final double pi      = 3.1415926535897932385d; // java.lang.Math.PI
    static final double piD2    = 1.5707963267948966192d; // pi/2
    static final double invTwoPi= 0.1591549430918953358d; // 1/(2*Pi)
     
    // theta(t)=Im{log[Pi(it/2-3/4)]}-t/2*log(pi)
    // using the approximate
    // t/2*(log(t/(2*Pi))-1)-Pi/8+1/(48*t)+7/(5760*t^3)+31/(80640*t^5)+127/(430080*t^7)+511/(1216512*t^9)
    // +1414477/(1476034560*t^11)+1898482392450447133975223/(2642676384606731475419894*t^13)
    // should not be called for t<=10 (absolute error of order 4*10^(-9))
    public double theta(double t) {
   		return t*0.5d*(Math.log(t*invTwoPi)-1.0) - pi/8.0d + 1.0d/(48.0d*t) + 7.0d/(5760.0d*t*t*t);
    }

	//Riemann-Siegel Z function (0<=n<5)
	//should not be called for t<100 or so even if n=4 (if precision is vital)
	public double Z(double t, int n) { 
	    if (t<10.0d) {
    	   if (t<-10.0d)
    	      return Z(-t, n);
    	   return -EulerML(new Complex(0.5d, t)).abs();
	    }
		double w= t*invTwoPi;   	//t/(2*pi)
		double p= Math.sqrt(w); 	//sqrt(t/(2*pi))
		double ip= 1.0d/p;          //sqrt(2*pi/t)
		double sip= Math.sqrt(ip); 	//(2*pi/t)^(1/4)
		int N = (int) p;
		p -= N;
		double tt = theta(t);
		int j=1;
		double ZZ= Math.cos(tt); // % (2.0d*pi));	
		while (++j <= N) { 
			ZZ += Math.cos((tt-t*Math.log((double) j))) // % (2.0d*pi))
   			  	 /Math.sqrt((double) j);
		}
		ZZ *= 2.0d;
		double po= 1.0d;
		double R= 0.0d;	
		int k=-1; 
		while (++k <= n) { 
			R += Coeff(k, 2.0d*p-1.0d)*po;
			po *= ip;
		}	
		if ((N & 1)!=0)
		   return ZZ+R*sip;
		else
		   return ZZ-R*sip;
	} 
	

    // Euler-MacLaurin expansion of zeta
    // should not be called for im(z) >> 1 (slow!)
	public Complex EulerML(Complex z) {
	
	    // works fine if  im(z) < 2*Pi*N
	    int N= (int) (z.im()*invTwoPi)+10;
	    
	    // sum(1/k^x, k=1..N)
	    Complex nz= z.neg();
	    double xx=nz.re();
	    double yy=nz.im();
	    double lk;
	    double ek;
	    double py;
	    double sx=1.0d;
	    double sy=0.0d;
	    //faster than s= s.add(Complex.pow((double) k, nz)); and heap friendly
	    for (int k=2; k<=N; k++) {
	       lk= Math.log((double) k);
	       ek= Math.exp(xx*lk);
	       sx += Math.cos(py=lk*yy)*ek;
	       sy += Math.sin(py)*ek;
	    }
	    Complex s= new Complex(sx, sy);
	    	
	    // + 1/((x-1)*N^(x-1))
	    Complex unit= new Complex(1.0d);
	    Complex xd= z.sub(unit);
	    Complex p= unit.div(xd.mul(Complex.pow((double)N, xd)));
	    s= s.add(p);
	    
	    //Bernoulli terms : 
	    
        // -1/(2*N^x)
        xd= z;
        p= unit.div(Complex.pow((double)N, xd).mul(new Complex(2.0d)));
        s= s.sub(p);
        
        // +x/(12*N^(x+1))
        p= p.mul(xd).div(new Complex(6.0d*N));
        s= s.add(p);
        
        //  -x*(x+1)*(x+2)/(720*N^(x+3))
        double NN= N*N;
        xd= xd.add(unit);
        p= p.mul(xd);
        xd= xd.add(unit);
        p= p.mul(xd).div(new Complex(120.0d*NN));
        s= s.sub(p);
        
        //  +x*(x+1)*(x+2)*(x+3)*(x+4)/(30240*N^(x+5))  		  
        xd= xd.add(unit);
        p= p.mul(xd);
        xd= xd.add(unit);
        p= p.mul(xd).div(new Complex(84.0d*NN));
        s= s.add(p);
        
        // result
        return s; 
    }

    //functional equation : zeta(z)= fE(z)*zeta(1-z)
    // with fE(z)= Gamma(1-z)*(2*pi)^(z-1)*2*sin(z*pi/2)   
    public Complex functionalEquation(Complex z) {
 
       Complex unit= new Complex(1.0d); // 1
       Complex oneMz= unit.sub(z); // 1-z
       Complex zMone= oneMz.neg(); // z-1
       /* Complex coef= Gamma(oneMz); //=> underflow for im(z)> 451.86
          coef= coef.mul((new Complex(pi*2)).pow(zMone));
          coef= coef.mul(z.mul(new Complex(piD2)).sin()).mul(new Complex(2.0));*/
       Complex coef= lnGamma(oneMz);
       coef= coef.add(zMone.mul(new Complex(ln2Pi)));
       if (z.im()>400.0) //sinus overflow
          coef= coef.add(new Complex(piD2*z.im(),piD2*(1.0d-z.re())));
       else if (z.im()<-400.0)
          coef= coef.add(new Complex(-piD2*z.im(),piD2*(z.re()-1.0d)));
       else //
          coef= coef.add(z.mul(new Complex(piD2)).sin().mul(new Complex(2.0d)).log());
       return coef.exp();
	}
	
	// use the appropriate method to evaluate the zeta function 
	// on the complex plane (for Re(z) > -6)
	public Complex zeta(Complex z) {
	   if (z.im()<0)
	      return zeta(z.conj()).conj();
	   else {
	      if ((z.re()==0.5d) && (z.im()>20.0d)) { //far critical line
	         // use Riemann-Siegel formula  ~ sqrt(y)
	         double ZZ= Z(z.im(),4);
	         double tt= theta(z.im());
	         return new Complex(ZZ*Math.cos(tt), -ZZ*Math.sin(tt));
	      }
	      else {
	         // use Euler-MacLaurin expansion  ~ y
	         return EulerML(z);
 	      }
	   }
	}
	

    public void paint(Graphics gr) {
      
      if (G_POINTSNUMBER<=0)
         return;
      //use local variables (global may be modified)           
      Complex z= G_z;      // z
      Complex mz= z.neg(); // -z
      int choice=G_CHOICE;
      int nb= G_POINTSNUMBER;
      double inScaling=G_INSCALING;
      double outScaling=G_OUTSCALING;
      int count= nb;  //impose count <= nb
       
      try {
       
       //array and number of points
       if (cache_nb!=nb) { //new arrays only if size change
          xP= new int[nb];
          yP= new int[nb];
          cache_nb= nb; //allocated size
          cache_choice= -1; 
       }
       
       //points generation corresponding to choice              
       switch (choice)  
       {
       case 0: //finite sum (convergent for Re(z)>1)
          {
             double d=2.0;
             Complex sum= new Complex(1.0d);
             xP[0]=(int)(G_m*sum.re()+G_ox);
             yP[0]=(int)(G_oy-G_m*sum.im());
             for (int i= 1; i< nb; i++) {                              
                sum= sum.add(Complex.pow(d, mz));
                xP[i]=(int)(G_m*sum.re()+G_ox);
                yP[i]=(int)(G_oy-G_m*sum.im());
                d += 1.0;
             }
          }
          break;
       case 1: //mirror finite sum (using functional equation) Re(z)<1 
          {
             //use the functional equation :
             Complex coef=functionalEquation(z);
             double d=1.0;
             Complex res= null;
             Complex zMone= z.sub(new Complex(1.0d)); // z-1
             Complex sum= new Complex(0.0d);
             for (int i= 0; i< nb; i++) {                              
                sum= sum.add(Complex.pow(d, zMone));
                res= sum.mul(coef);
                xP[i]=(int)(G_m*res.re()+G_ox);
                yP[i]=(int)(G_oy-G_m*res.im());
                d += 1.0;
             }
          }
          break;
       case 2: //alternate sum (convergent for Re(z)>0)
          {
             double d=2.0;
             Complex cSum;
             boolean pos = false;
             //compute  coef= 1/(2^(1-x)-1)
             Complex unit= new Complex(1.0d);
             Complex coef=unit.div(Complex.pow(2.0,unit.add(mz)).sub(unit));
             Complex sum= new Complex(-1.0d);
             cSum= sum.mul(coef);
             xP[0]=(int)(G_m*cSum.re()+G_ox);
             yP[0]=(int)(G_oy-G_m*cSum.im());
             for (int i= 1; i< nb; i++) { 
                pos = !pos;                             
                if (pos)
                   sum= sum.add(Complex.pow(d, mz));  //mz = -z
                else
                   sum= sum.sub(Complex.pow(d, mz));
                cSum= sum.mul(coef);
                xP[i]=(int)(G_m*cSum.re()+G_ox);
                yP[i]=(int)(G_oy-G_m*cSum.im());
                d += 1.0;
             }
          }
          break;
       case 3: //mirror alternate sum
          {
              //use the functional equation :
             Complex coef=functionalEquation(z);
             double d=2.0;
             Complex cSum;
             boolean pos = false;
             Complex unit= new Complex(1.0d);
             //multiply by 1/(2^(1-x)-1) with z=1-x here
             coef= coef.div(Complex.pow(2.0,z).sub(unit));
             Complex sum= new Complex(-1.0d);
             Complex zMone= z.sub(unit); // z-1
             cSum= sum.mul(coef);
             xP[0]=(int)(G_m*cSum.re()+G_ox);
             yP[0]=(int)(G_oy-G_m*cSum.im());
             for (int i= 1; i< nb; i++) { 
                pos = !pos;                             
                if (pos)
                   sum= sum.add(Complex.pow(d, zMone));  //zMone = z-1
                else
                   sum= sum.sub(Complex.pow(d, zMone));
                cSum= sum.mul(coef);
                xP[i]=(int)(G_m*cSum.re()+G_ox);
                yP[i]=(int)(G_oy-G_m*cSum.im());
                d += 1.0;
             }
          }
          break;
       case 4: //finite Euler product
          {
             if (PrimeOffs==null)  //use the Sieve to compute primes offset table
                Sieve(max);
             if (PrimeOffs!=null) {
                           
                double d=2.0;
                Complex prd= new Complex(1.0d);
                Complex unit= new Complex(1.0d);
                xP[0]=(int)(G_m*prd.re()+G_ox);
                yP[0]=(int)(G_oy-G_m*prd.im());
                if (nb>=primeCount)
                   count= primeCount-1;
                else
                  count= nb;
                for (int i= 1; i< count; i++) { 
                   prd= prd.div(unit.sub(Complex.pow(d, mz)));                   
                   xP[i]=(int)(G_m*prd.re()+G_ox);
                   yP[i]=(int)(G_oy-G_m*prd.im());
                   d += PrimeOffs[i];
                }
             }
          }
          break;
       case 5: //Riemann-Siegel function
          {
             //Output y corresponds to y input (output scaling on x only)
             //On critical line use Riemann-Siegel (and zeta for abs(im(z))<10)
             if (G_xRef==0.5d) {
                double xx;
                double yy;
                double maxT= G_yRef+inScaling/2.0; //logical (at top of frame)
                double offset=inScaling/nb;
                double max= Math.round((float)(maxT/offset))*offset; //for cache
                cacheRestore(choice, inScaling, outScaling, max, -offset); //set cache_first and cache_last
             	for (int i= cache_first; i<= cache_last; i++) {
             		yy= max-i*offset;            	    
                	xx= Z(yy, 4);
	                xP[i]=(int)(G_m*xx+G_WIDTH/2); //output scaling acts on X only
                }
                double oy=(maxT-max)/offset;
             	for (int i= 0; i< nb; i++)
    	            yP[i]=(int) ((oy+i*G_HEIGHT)/nb);
                cache_start= max;
                cache_offset= -offset;
             }
          }
          break;
       case 6: //Riemann zeta function
          {
             //along y axis only yet (any constraint later...)
             double max= G_yRef+inScaling/2.0; //logical value at top of frame
             double yy;
             Complex zz;
             double offset= inScaling/nb;
             max= Math.round((float)(max/offset))*offset; //for cache
             cacheRestore(choice, inScaling, outScaling, max, -offset); //set cache_first and cache_last
             for (int i= cache_first; i<= cache_last; i++) {
             	 yy= max-i*offset;				 
             	 zz= zeta(new Complex(G_xRef, yy));
             	 xP[i]=(int)(G_m*zz.re()+G_WIDTH/2);
             	 yP[i]=(int)(G_oy-G_m*zz.im());
             }
             cache_start= max;
             cache_offset= -offset;
          }
          break;
       case 7: //Riemann zeta 3D
          {
              //along y axis only yet (any constraint later...)
             double max= G_yRef+inScaling/2.0; //logical value at top of frame
             double min= max-inScaling; //logical value corresponding to foreground
                                          // (distance 1 from observer)
             double yy;
             Complex zz;
             double rb= 2000.0d/outScaling;
             double offset= inScaling/nb;
             min= Math.round((float)(min/offset))*offset; //for cache
             cacheRestore(choice, inScaling, outScaling, min, offset); //set cache_first and cache_last
             for (int i= cache_first; i<= cache_last; i++) {
             	 yy= min+(i*inScaling)/nb;				 
             	 zz= zeta(new Complex(G_xRef, yy));
         	     xP[i]=Math.round((float) (zz.re()*rb));
         	     yP[i]=Math.round((float) (zz.im()*rb));
         	 }
         	 //drawing
         	 int x,y;
         	 rb= G_HEIGHT/2000.0d;
        	 int ox=(int)(rb*xP[0]+G_WIDTH/2);
             int oy=(int)(G_oy-rb*yP[0]);
             double r;
             initScr(gr,z,choice,outScaling);            	 
             for (int i= 1; i< nb; i++) {
                 r= rb/(1.0d+(i*inScaling)/nb);
         	     x=(int)(r*xP[i]+G_WIDTH/2);
             	 y=(int)(G_oy-r*yP[i]);
        	 	 gr.drawLine(ox, oy, x, y);
             	 ox= x;
             	 oy= y;
             }
             count=0; //drawing already done (xP,yP for buffering only)
            /*for (int i= 0; i< nb; i++) {
             	 yy= min+(i*inScaling)/nb;				 
             	 zz= zeta(new Complex(G_xRef, yy));
                 r= G_HEIGHT/(outScaling+2.0d*(yy-min));
         	     xP[i]=(int)(r*zz.re()+G_WIDTH/2);
             	 yP[i]=(int)(G_oy-r*zz.im());
             }*/
          }
          break;
       case 8: //goodies
          {
             //none yet...
             count=0;
          }
          break;
       }
       //effective drawing
       if (count!=0) {       
          initScr(gr,z,choice,outScaling);   
          gr.drawPolyline(xP, yP, count);
       }
     } 
     catch (Exception e) {
       //error trapping
       nb= -1; //empty cache (in 'finally')
       xP= null;
       yP= null;
     } 
     finally {
       cache_nb= nb;
       cache_choice= choice;
       cache_inScaling= inScaling;
       cache_outScaling= outScaling;
     }     
    }
    
    // new 'number of points'
    public void newPointsNumber(int pointsNumber) {
       G_POINTSNUMBER = pointsNumber;
       updateAll();      
    }
    
    //Reset Default 
    public void resetDefault() {     
       G_POINTSNUMBER = 200;
       G_INSCALING  = 1;
       G_OUTSCALING = 1;
       G_m= G_HEIGHT; //unit height in pixels
       G_ix = G_WIDTH/2; //-G_HEIGHT/2;
       G_iy = G_HEIGHT/2;
       G_ox = G_ix;
       G_oy = G_iy;
       G_xRef= 0.5d;
       G_yRef= 14.13472514d;
       G_z= new Complex(G_xRef, G_yRef);
       updateAll();      
    }
    
    //new choice
    public void update(int choice) {     
       G_CHOICE= choice;
       updateAll();      
    }

    public void restrainInput(int restrain) {
       G_RESTRAIN= restrain;
    }

    //clic in the top panel
    
    public void newPoint(int x, int y) {
       switch (G_RESTRAIN) {
       case 0 : //none
          G_z= new Complex(G_xRef+(x-G_ix)*G_INSCALING/G_HEIGHT,
          G_yRef-(y-G_iy)*G_INSCALING/G_HEIGHT);
          break;
       case 1 : //x =G_ox only
          G_z= new Complex(G_xRef, G_yRef-(y-G_iy)*G_INSCALING/G_HEIGHT);
          break;
       case 2 : //y =G_oy only
          G_z= new Complex(G_xRef+(x-G_ix)*G_INSCALING/G_HEIGHT, G_yRef);
          break;
       case 3 : { //radius                   
             double dist=Math.sqrt((x-G_ix)*(x-G_ix)+(y-G_iy)*(y-G_iy));
             if (dist==0.0)
                G_z= new Complex(G_xRef, G_yRef);
             else
                G_z= new Complex(G_xRef+(x-G_ix)*G_Radius/dist, 
                                 G_yRef-(y-G_iy)*G_Radius/dist);
          }
          break;      
       case 4 : { //arg
             double prod=((x-G_ix)*Math.cos(G_Arg)-(y-G_iy)*Math.sin(G_Arg))/G_HEIGHT;
             G_z= new Complex(G_xRef+prod*Math.cos(G_Arg), 
                              G_yRef+prod*Math.sin(G_Arg));
          }
          break;      
       }
       updateAll();      
    }
    
   public void newOutputScaling(double outputScaling) {
       G_OUTSCALING = outputScaling;
       G_m= G_HEIGHT/outputScaling; //unit height in pixels
       updateAll();      
    }
    
    public void newInputScaling(double inputScaling) {
       if (inputScaling!=0.0)
          G_INSCALING = inputScaling;
       if (G_CHOICE>4)
          updateAll();      
    }   

    
    public void newXOffset(double xOffset) {
       G_xRef= xOffset; 
       G_z= new Complex(xOffset, G_z.im());
       G_ix= G_WIDTH/2-xOffset/G_INSCALING; //new center (pixels)
       updateAll();      
    }

    
     public void newYOffset(double yOffset) {
       G_yRef= yOffset;
       G_z= new Complex(G_z.re(), yOffset);
       G_iy= G_HEIGHT/2+yOffset/G_INSCALING; //new center (pixels)     
       updateAll();      
    }
    
    public void newRadius(double radius) {
       G_Radius= radius;
       updateAll();      
    }

    //input in degrees
    public void newArg(double arg) {
       G_Arg= arg*Math.PI/180.0;
       updateAll();      
    }

    //if next=1 next one (->right or ^ top)
    public void findZero(int next) {
       //start near the selected point
       G_xRef= G_z.re();
       G_yRef= G_z.im();
       double inc=0.0;
       if (next>0) {
          if (G_xRef<-2.0d) {
             double nx= Math.ceil(G_xRef/2.0d)*2.0d;
             if (G_xRef==nx)
                nx += 2.0d;
             G_xRef= nx;
             G_yRef= 0.0d;
          }
          else {
             G_xRef=0.5d;  //OK R.H supposed but...
             if (Math.abs(G_yRef)< ROOT_1)
                G_yRef=ROOT_1;
             else { //search needed
                inc=next;
             }
          }                
       }       
       else {
         if (G_xRef<0.0d) {
             double nx= Math.floor(G_xRef/2.0d)*2.0d;
             if (G_xRef==nx)
                nx -= 2.0d;
             G_xRef= nx;
             G_yRef= 0.0d;
          }
          else {
             G_xRef=0.5d;  //OK R.H supposed but...
             if (Math.abs(G_yRef)< ROOT_1)
                G_yRef=-ROOT_1;
             else { //search needed
                inc=next;
             }
          }                
       }
       if (inc!=0.0d) { //real work
          //perhaps use Gram points later...
          double st=G_yRef; //Math.abs(G_yRef);
          inc /= (2.0d*Math.log(Math.abs(st)));
          double zz= Z(st,4); 
          while (Math.abs(zz)<2.5e-6d) { //initial root
             st += inc;
             zz= Z(st,4); 
          }
          double z=st;
          double nz;
          do {         
             z += inc;         
             nz= Z(z,4);
          }
          while (nz*zz>0);
          //use interpolation
          double z1 = z-inc;
          double zz1= Z(z1,4);
          double z2 = z;
          double zz2= nz;
          int cpt=100; //security limit
          do {
             z= z1+zz1*(z2-z1)/(zz1-zz2);
             nz= Z(z,4);
             if (nz*zz1<0.0d) {
                z2=z;
                zz2=nz;
             }
             else if (nz*zz2<0.0d) {
                z1=z;
                zz1=nz;
             }
             else
                break;               
             --cpt;
          }
          while (Math.abs(nz)>1e-7d && zz1!=zz2 &&cpt>=0);
          G_yRef= z;
       }
       G_z= new Complex(G_xRef, G_yRef);
       G_ix= G_WIDTH/2-G_xRef/G_INSCALING; //new center (pixels)
       G_iy= G_HEIGHT/2+G_yRef/G_INSCALING;
       updateAll();      
    }    
    
    //{{ declarations - Begin Non-Editable Code Block
    //}} declarations - End Non-Editable Code Block
    
    public static double G_HEIGHT= 495; //Back height
    public static double G_WIDTH = 700; //Back width
    public static int    G_POINTSNUMBER= 200;
    public static int    G_CHOICE= 2; 
    public static int    G_RESTRAIN= 0; //none
    public static double G_m= G_HEIGHT; //corresponding to G_OUTSCALING = 1
    public static double G_OUTSCALING = 1; //scaling on output
    public static double G_INSCALING = 1;  //scaling on input
    public static double G_ix= G_WIDTH/2;  //input offset
    public static double G_iy= G_HEIGHT/2;
    public static double G_ox= G_ix;    //output offset
    public static double G_oy= G_iy;
    public static double G_Radius= 0.0; //constraint
    public static double G_Arg= 0.0;
    public static double G_xRef= 0.5;
    public static double G_yRef= ROOT_1;
    public static Complex G_z= new Complex(G_xRef, G_yRef);
    public static int[] xP= null;
    public static int[] yP= null;
    public static int cache_nb=-1; //arrays empty
    public static int cache_choice=-1; //none
    public static double cache_inScaling=1.0d; //for comparison purpose
    public static double cache_outScaling=1.0d;
    public static double cache_start=0.0d;
    public static double cache_offset=1.0d;
    public static int cache_first=0; //first not covered by cache
    public static int cache_last=0; //last not covered by cache
 
    //{{ extra_code
    //}} extra_code
}