#include "pari.h"
#include <math.h>

double G(double t, double w, long upper, long sign)
{
  double x=1;
  double xstep = .0005;
  double total=0;
  int flag=1;
  for (x=1; x<=1.9998; x+=0.0002)	
        {
         if (sign == 1)
          total+=.333333*.0002*flag*(pow(M_E,((-1)*w*x)))*2*(sin(1.0*t*log(1.0*x)));
         else
          total+=.333333*.0002*flag*(pow(M_E,((-1)*w*x)))*2*(cos(1.0*t*log(1.0*x)));        
         if (x != 2)
	   {
	     switch(flag)
	       {
	        case 1: flag = 4; break;
                case 4: flag = 2; break;
                case 2: flag = 4; break;
               }
	   }
	}
      if (sign == 1)
       total+=.333333*.0002*(pow(M_E,((-1)*w*2)))*2*(sin(1.0*t*log(1.0*2)));
      else
       total+=.333333*.0002*(pow(M_E,((-1)*w*2)))*2*(cos(1.0*t*log(1.0*2)));
      for (x=2; x<=40; x+=xstep)
	{
         if (sign == 1)
          total+=0.5*xstep*((pow(M_E,((-1)*w*x)))*2*(sin(1.0*t*log(1.0*x)))+(pow(M_E,((-1)*w*(x+xstep))))*2*(sin(1.0*t*log(1.0*(x+xstep)))));
         else
          total+=0.5*xstep*((pow(M_E,((-1)*w*x)))*2*(cos(1.0*t*log(1.0*x)))+(pow(M_E,((-1)*w*(x+xstep))))*2*(cos(1.0*t*log(1.0*(x+xstep)))));
        }
  return total;
}

int main()
{
 GEN F, E, blah, y;
 long sign, n, N, upper, nbound;
 long a;
 double w, t,finalt;
 int i, signcurve,family;
 double temp;
 double values[40];
 double answer;
 double epsilon=0.00001;
 double tstep=0.1;
 double start;
 double numepsilon = (-1)*(log(epsilon));
 //printf("familyt = %d\n",family);
 pari_init(10000000, 2);

 y = cgeti(32);
 blah = cgeti(32);
 F = cgetg(6, t_VEC);
 for (i=1; i<=5; i++)
 {
  F[i] = lgeti(MEDDEFAULTPREC);
 }
 E = cgetg(14, t_VEC);
 for (i=1; i<=13; i++)
 {
  E[i] = lgeti(MEDDEFAULTPREC);
 }
 for (family = 85; family<=100; family++)
 {
  tstep=0.1;
  gaffsg(0,(GEN) F[1]);
  gaffsg(family+1,(GEN) F[2]);
  gaffsg(0,(GEN) F[3]);
  gaffsg(family,(GEN) F[4]);
  gaffsg(0,(GEN) F[5]);

  gaffsg(1,(GEN) y);
 /*0,t+1,0,t,0*/
 /* 0,-1,1,-10,-20 */
 /* 0,0,0,0,13 */

  E=initell(F,MEDDEFAULTPREC);
  signcurve=0;
  sign = ellrootno(E,y);
  //printf("sign = %ld\n",sign);
  gaffect( (GEN) globalreduction(E)[1], (GEN) y);
  N = gtolong( (GEN) y);
  nbound = (int) (2*log(1.0*N)*pow(N,0.5)+1);
  //printf("\nconductor = %ld\n",N);
  //printf("nbound = %ld\n", nbound);

 //w = (6.2831854)/(pow(11,0.5));
 //printf("G = %f\n", G(1,w,10,-1));

  start = 0.1;
  sign = -1 * sign;
  for (t=start; t<29; t+=tstep)
  {  
   answer=0;

   temp=0;
   for (n=1; n< nbound; n++)
     {
       w = (6.2831854*n)/ sqrt(1.0*N);
       upper = (1/w)*numepsilon;
       gaffsg(n,(GEN) y);
       blah = akell(E,(GEN) y);
       a= gtolong(blah);
       if (a!=0)
	 {
	   temp = 10000*G(t,w,upper,sign);
	   temp = 1.0*a*temp;
          answer+=temp;
  //         if ((n%50) == 0)
    //        printf("n=%ld nbound = %ld\n",n,nbound); 
	 }
     }
   if (t==start)
     {
      if (answer>0)
	signcurve=1;
      if (answer<0)
        signcurve=-1;
      if (answer == 0)
        signcurve=0;
     }
   if ((signcurve<=0)&&(answer>0))
   {
        if (tstep>=0.01)
        {
         t-=tstep;
         tstep=tstep/10;
        }
        else
	{
          finalt=t; 
	  t=100;
        }
        
   }
   if ((signcurve>=0)&&(answer<0))
   {
    if (tstep>=0.01)
        {
         t-=tstep;
         tstep=tstep/10;
        }
    else
      {
	finalt=t;
        t=100;
      }
   }
   //if (t!=100)
   // printf("t = %f value = %f\n",t,answer);
   if (t==100)
    {
     printf("zero at %f for t = %d\n",finalt,family);
    // values[family] = finalt;
    }
  }
 }
 //for (family = 2; family<=40; family++)
 // printf("for familyt = %d, value = %f",family,values[family]);
 return 0;
}