#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "pari.h"
#include "RAT.h"

/* number of primes checked
 should be bounded by number of primes in smp[][]*/
#define num_primes 10

/*Euclid's algorithm to return gcd*/
int gcd(long int, long int);

/* new mod function that handles negative values */
long mod(long x, long y)
{
  long x2=1, tempmody=0;
  if (x<0)
    {
      x2=0-x;
      tempmody=y-(x2%y);
      tempmody=tempmody%y;
      return tempmody;
    }
  else if (x>0)
    return (x%y);
  else
    return 0;
}


int main(int argc, char *argv[]) 
{
  /*declare pari vars*/
  GEN E, F, pr, fp1, fp2, a, tempg, pp1, pp2, pp11, pp12, pp21, pp22;

  GEN sp1, sp2, b11, b12, b13, b22, b23, b33;
  double a11, a12, a13, a22, a23, a33, det, tempdet, lowestdet;
  long ind;


  /*N=maximum naive height*/
  /*A and B are coeffs of the elliptic curve*/
  long  N, ctr, GCD=0, c_squared=0, c_cubed=0;
  long t=0, t1=0, t2=0, t_upper=0, t_lower=0;
  /*x = x_1/x_2*/
  /* x = a/c^2  y = b/c^3 */
  /*lower_bound = point of order 2 in E(R) w/ min x-coord */
  long x_1, c, lower_bound;
  long tempx3p, tempc3p, tempc6p, tempk3p, tempt2p, p, tempt1p, tempxp, tempt1xp, tempt1xkp;
  double doubleD;
  double canonical_height; //canonical height as computed by ghell2
  double doublesqrtD, sqrtdifference;
  long longD1;

  /* variables for third method - square mod p. smp is a 2x2 array of the first 
     10 primes and their squares (including 0). ismp is a 2x2 array of the inverse 
     of these squares mod p. ismp is not currently used, but i included it anyway, 
     since it may be helpful */
  long smp[10][17]={ {3, 0, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
		     {5, 0, 1, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
		     {7, 0, 1, 2, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
		     {11, 0, 1, 3, 4, 5, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
		     {13, 0, 1, 3, 4, 9, 10, 12, -1, -1, -1, -1, -1, -1, -1, -1, -1},
		     {17, 0, 1, 2, 4, 8, 9, 13, 15, 16, -1, -1, -1, -1, -1, -1, -1},
		     {19, 0, 1, 4, 5, 6, 7, 9, 11, 16, 17, -1, -1, -1, -1, -1, -1},
		     {23, 0, 1, 2, 3, 4, 6, 8, 9, 12, 13, 16, 18, -1, -1, -1, -1},
		     {29, 0, 1, 4, 5, 6, 7, 9, 13, 16, 20, 22, 23, 24, 25, 28, -1},
		     {31, 0, 1, 2, 4, 5, 7, 8, 9, 10, 14, 16, 18, 19, 20, 25, 28} };

  long ctr1=0, ctr2=0, temp2=0, k=0, D=0;
  /* pass verifies that a point passed the first 10 primes */
  int pass=0;

  int method=1;

  double progress=0, increment=0;

  FILE *fp;
  fp=fopen("output", "w");

  pari_init(10000000, 2);
  printf("Initializing pari variables\n");

  b11 = cgetr(MEDDEFAULTPREC);   /*new code I added*/
  b12 = cgetr(MEDDEFAULTPREC);
  b13 = cgetr(MEDDEFAULTPREC);
  b22 = cgetr(MEDDEFAULTPREC);
  b23 = cgetr(MEDDEFAULTPREC);
  b33 = cgetr(MEDDEFAULTPREC);

  sp1=cgetg(3, t_VEC);
  sp1[1] = cgeti(MEDDEFAULTPREC);
  sp1[2] = cgeti(MEDDEFAULTPREC);

  sp2=cgetg(3, t_VEC);
  sp2[1] = cgeti(MEDDEFAULTPREC);
  sp2[2] = cgeti(MEDDEFAULTPREC);

  pp1=cgetg(3, t_VEC);
  (GEN) pp1[1]=cgetg(3, t_FRAC);
  (GEN) pp1[2]=cgetg(3, t_FRAC);
  mael2(pp1,1,1) = lgeti(MEDDEFAULTPREC);
  mael2(pp1,1,2) = lgeti(MEDDEFAULTPREC);
  mael2(pp1,2,1) = lgeti(MEDDEFAULTPREC);
  mael2(pp1,2,2) = lgeti(MEDDEFAULTPREC);
  pp11 = cgetg(3, t_FRAC);          /* makes fp1 a vector */
  pp11[1] = lgeti(MEDDEFAULTPREC); /* initializes components */
  pp11[2] = lgeti(MEDDEFAULTPREC);
  pp12 = cgetg(3, t_FRAC);          /* makes fp2 a vector */
  pp12[1] = lgeti(MEDDEFAULTPREC); /* initializes components */
  pp12[2] = lgeti(MEDDEFAULTPREC);

  pp2=cgetg(3, t_VEC);
  (GEN) pp2[1]=cgetg(3, t_FRAC);
  (GEN) pp2[2]=cgetg(3, t_FRAC);
  mael2(pp2,1,1) = lgeti(MEDDEFAULTPREC);
  mael2(pp2,1,2) = lgeti(MEDDEFAULTPREC);
  mael2(pp2,2,1) = lgeti(MEDDEFAULTPREC);
  mael2(pp2,2,2) = lgeti(MEDDEFAULTPREC);
  pp21 = cgetg(3, t_FRAC);          /* makes fp1 a vector */
  pp21[1] = lgeti(MEDDEFAULTPREC); /* initializes components */
  pp21[2] = lgeti(MEDDEFAULTPREC);
  pp22 = cgetg(3, t_FRAC);          /* makes fp2 a vector */
  pp22[1] = lgeti(MEDDEFAULTPREC); /* initializes components */
  pp22[2] = lgeti(MEDDEFAULTPREC);

  /*new code end*/

  tempg=cgeti(32);
  a=cgeti(32);
  pr=cgetg(3, t_VEC);
  (GEN) pr[1]=cgetg(3, t_FRAC);
  (GEN) pr[2]=cgetg(3, t_FRAC);
  mael2(pr,1,1) = lgeti(MEDDEFAULTPREC);
  mael2(pr,1,2) = lgeti(MEDDEFAULTPREC);
  mael2(pr,2,1) = lgeti(MEDDEFAULTPREC);
  mael2(pr,2,2) = lgeti(MEDDEFAULTPREC);
  fp1 = cgetg(3, t_FRAC);          /* makes fp1 a vector */
  fp1[1] = lgeti(MEDDEFAULTPREC); /* initializes components */
  fp1[2] = lgeti(MEDDEFAULTPREC);
  fp2 = cgetg(3, t_FRAC);          /* makes fp2 a vector */
  fp2[1] = lgeti(MEDDEFAULTPREC); /* initializes components */
  fp2[2] = lgeti(MEDDEFAULTPREC);
  F=cgetg(6, t_VEC);
  for (ctr=1; ctr<6; ctr++)
    F[ctr]=lgeti(MEDDEFAULTPREC);
  E=cgetg(14, t_VEC);
  for (ctr=1; ctr<14; ctr++)
    E[ctr]=lgeti(MEDDEFAULTPREC);
  printf("Finished initializing pari variables\n");

 
  printf("Choose an upper bound for the height:\n");
  scanf("%ld", &N);
  printf("This program returns the rational points with height <= %ld on an elliptic curve of the form y^2=x^3+Ax+B.  Please enter a lower bound for t:\n", N);
  scanf("%ld", &t_lower);
  printf("Now enter an an upper bound for t:\n");
  scanf("%ld", &t_upper);

  increment=1.0/(N*(t_upper - t_lower+1));
 
  printf("        ");
  fprintf(fp, "N=%ld; t from %ld to %ld; method %d\n", N, t_lower, t_upper, method);
  
  for (t=t_lower; t<=t_upper; t++)
    {
     /*HERE IS WHERE THE CURVE IS HARD-WIRED IN*/
      /*y^2=x^3 + t1*x + t2*/
	t1 = -1*t*t;
	t2 = t*t*t*t;
	fprintf(fp, "Curve: Y^2 = X^3 + %ldX + %ld\n", t1, t2);
	/* storing these values in the F and E vector */
	gaffsg(0, (GEN) F[1]);
	gaffsg(0, (GEN) F[2]);
	gaffsg(0, (GEN) F[3]);
	gaffsg(t1, (GEN) F[4]);
	gaffsg(t2, (GEN) F[5]);


      /*new code added by me*/
      /* assign forced points*/

      gaffsg(t, (GEN) sp1[1]);
      gaffsg(t*t, (GEN) sp1[2]);
      
      gaffsg(t*t, (GEN) sp2[1]);
      gaffsg(t*t*t, (GEN) sp2[2]);

      gaffsg(t, (GEN) pp11[1]); gaffsg(1, (GEN) pp11[2]);
      gaffsg(t*t, (GEN) pp12[1]); gaffsg(1, (GEN) pp12[2]);
      
      gaffsg(t*t, (GEN) pp21[1]); gaffsg(1, pp21[2]);
      gaffsg(t*t*t, (GEN) pp22[1]); gaffsg(1, (GEN) pp22[2]);

      gaffect((GEN) pp11[1], (GEN) mael2(pp1,1,1));
      gaffect((GEN) pp11[2], (GEN) mael2(pp1,1,2));
      gaffect((GEN) pp12[1], (GEN) mael2(pp1,2,1));
      gaffect((GEN) pp12[2], (GEN) mael2(pp1,2,2));

      gaffect((GEN) pp21[1], (GEN) mael2(pp2,1,1));
      gaffect((GEN) pp21[2], (GEN) mael2(pp2,1,2));
      gaffect((GEN) pp22[1], (GEN) mael2(pp2,2,1));
      gaffect((GEN) pp22[2], (GEN) mael2(pp2,2,2));


      /*assign values to some entries of height matrix*/
      
      gaffect( (GEN) bilhell(E,sp1, sp1,MEDDEFAULTPREC), (GEN) b11);
      gaffect( (GEN) bilhell(E,sp1, sp2,MEDDEFAULTPREC), (GEN) b12);
      gaffect( (GEN) bilhell(E,sp2, sp2,MEDDEFAULTPREC), (GEN) b22);
      
      a11 = gtodouble( (GEN) b11);
      a12 = gtodouble( (GEN) b12);
      a22 = gtodouble( (GEN) b22);
      
      /*end of new code*/

	E=initell(F, MEDDEFAULTPREC);
	lower_bound = 1;
	
	
	
	fprintf(fp, "x\t\ty\t\tcanonical height\tdeterminant\tindependent?\n");
	/*Debugging line: Making sure all the expected points are there*/
	fprintf(fp, "(t, t^2)=(%ld, %ld); (t^2, t^3)=(%ld, %ld); 
(-t, t^2)=(%ld, %ld); (1/4, t^2-1/8)=(1/4, %ld - 1/8)\n\n", t, t*t, t*t, t*t*t, -1*t, t*t, t*t);
	

	    /* third method for determining y's. note: x_1/c^2 ~ xk (mod p) where ~ means congruent to.
	       this method determines whether x^3+Ax+B ~ (xk)^3 + t2 (modp) is a square mod p,
	       where k=(c^2)^(-1). i will test for the first num_primes. */
  
	    for (x_1 = (-1 * N); x_1 <= N;  x_1++) {
	      for (c = lower_bound; c*c <=N; c++)  
		{
		  c_squared=c*c;
		  c_cubed=c_squared*c;
		  pass=0;
		  if ( gcd(x_1*x_1,c) ==1)
		    {
		     
		      for (ctr1=0; ctr1<num_primes; ctr1++)            /* cycle through primes in smp[][] */
			{
			  			 
			  if ( gcd(c, smp[ctr1][0]) == 1)
			    {
			      temp2=mod(c_squared, smp[ctr1][0]);      /* used to determine (c^2)^-1 */
			      
			      /* Determines the inverse of temp2 and assigns it to k */
			      /* cycles through integetes < p and if one of them multiplied by */
			      /* temp2 (mod p) is 1, it is assigned to k */
			      for (ctr=0; ctr<smp[ctr1][0]; ctr++)
				{
				  if (temp2==0)
				    k=0;
				  else if ( (mod((temp2*ctr), smp[ctr1][0])) == 1)
				    k=ctr;
				}
			      
			      /* Determines c^6 * ((x_1*k)^3 + t2) */

			      p = smp[ctr1][0];
			      tempc3p = mod(c, p) * mod(c, p) * mod(c, p);
			      tempc3p = mod(tempc3p, p);
			      tempc6p = mod((tempc3p * tempc3p), p);
			      tempx3p = mod(x_1, p) * mod(x_1, p) * mod(x_1, p);
			      tempx3p = mod(tempx3p, p);

			      tempk3p = mod(k, p) * mod(k, p) * mod(k, p);
			      tempk3p = mod(tempk3p, p);
			      tempt2p = mod(t2, p);

			      tempt1p = mod(t1, p);
			      tempxp = mod(x_1, p);
			      tempt1xp = mod(tempt1p*tempxp, p);
			      tempt1xkp = mod(tempt1xp*k, p);


			      
			      D = mod(tempx3p * tempk3p, p) + tempt1xkp + tempt2p;
			      
			      D = mod(D, p);
			     
			      
			      /* Determines if D is a square mod p for the first 10 primes */
			      /* cycles through all the points stored in the row of a prime */
			      /* and checks if D(mod p) is any of them (-1 is ignored) */
			      for (ctr2=1; ctr2<17; ctr2++)
				{
				 
				  if ( mod(D, smp[ctr1][0]) == smp[ctr1][ctr2])
				    pass++; 
			
				}
			    			      
			      /* if D is a square mod p, for first 10 primes, pass should be 10 */
			      if (pass==num_primes)
				{
				  
				  doubleD=(1.0*x_1*x_1*x_1)+(1.0*t1*x_1*c*c_cubed)+(1.0*t2*c_cubed*c_cubed);
				 
				 if (doubleD >= 0)
				 {
				  doublesqrtD = sqrt(doubleD);
				  longD1 = floor(doublesqrtD);
				  sqrtdifference = doublesqrtD - 1.0*longD1; //This is a double


				  
				  
				  /* if the difference between (templ)^2 and D is very small,
				     print out the points.
				     i think this can be changed, though. since floor() will
				     give you the integer closest to the given number. thus
				     the difference betwen (templ)^2 and D must be an integer.
				     so, instead, we could just check that the two are equal, since
				     if there difference is less than some decimal, they would
				     have to be equal. i left it the way it is, following Steve's
				     suggestion, though */
				  if (sqrtdifference<0.0001)
				    {
				      /*debugging statement*/
				      /*printf("sqrtdifference<0.0001");*/
				      gaffsg( x_1, (GEN) fp1[1]);
				      gaffsg( c_squared, (GEN) fp1[2]);
				      GCD=gcd(longD1, c_cubed);
				      gaffsg( (longD1/GCD), (GEN) fp2[1]);
				      gaffsg( (c_cubed/GCD), (GEN) fp2[2]);

				      gaffect((GEN) fp1[1], (GEN) mael2(pr,1,1));
				      gaffect((GEN) fp1[2], (GEN) mael2(pr,1,2));
				      gaffect((GEN) fp2[1], (GEN) mael2(pr,2,1));
				      gaffect((GEN) fp2[2], (GEN) mael2(pr,2,2));

				      if (oncurve(E, pr))
					{
					  /* check that canonical height is bounded by N */
					  canonical_height = gtodouble((GEN) ghell2(E, pr, MEDDEFAULTPREC));
					  if ( canonical_height <= N)
					    {
					      fprintf(fp, "%ld", gtolong((GEN) fp1[1])); fprintf(fp, "/");
					      fprintf(fp, "%ld", gtolong((GEN) fp1[2])); fprintf(fp, "\t");
					      fprintf(fp, "%ld", gtolong((GEN) fp2[1])); fprintf(fp, "/");
					      fprintf(fp, "%ld", gtolong((GEN) fp2[2])); fprintf(fp, "\t");
					      fprintf(fp, "%f\t\t", canonical_height);

					      /*new code added*/
					      gaffect( (GEN) bilhell(E,sp1, pr,MEDDEFAULTPREC), (GEN) b13);
					      gaffect( (GEN) bilhell(E,sp2, pr,MEDDEFAULTPREC), (GEN) b23);
					      gaffect( (GEN) bilhell(E,pr, pr,MEDDEFAULTPREC), (GEN) b33);
					      
					      gaffect( (GEN) bilhell(E,sp1, pr,MEDDEFAULTPREC), (GEN) b13);
					      gaffect( (GEN) bilhell(E,sp2, pr,MEDDEFAULTPREC), (GEN) b23);
					      gaffect( (GEN) bilhell(E,pr, pr,MEDDEFAULTPREC), (GEN) b33);
					      
					      a13 = gtodouble( (GEN) b13);
					      a23 = gtodouble( (GEN) b23);
					      a33 = gtodouble( (GEN) b33);
					      
					      det = a11*a22*a33 + 2*a12*a23*a13 - a13*a13*a22 - a23*a23*a11 - a12*a12*a33;
					      if(det < 0)
						det = -det;
					      
					      if(det < 0.001)
						ind = 0;
					      else ind = 1;
					
					      fprintf(fp, "%f\t%ld\n", det, ind);
					      /*new code end*/
					    }
					}
				    }
}
				} /*end of "if pass == number" of primes statement*/
			    }
			  else
			    {
			      pass++;
			    }
			} //end of for statement
		    } //end of if statement
		}
	      progress=progress+increment;
	      printf("\b\b\b\b\b\b\b\b%8.6f ", progress);
	    }
	    printf("\n");
	fprintf(fp, "-------------------------------------------------\n");
    }
  
  fprintf(fp, "======================================================\n");
  fprintf(fp, "======================================================\n");
  fclose(fp);
  
  return 0;
  
}

/*Euclid's algorithm*/
int gcd(long int m, long int n)
{
  if (n==0) return m;  /*return gcd*/
  return gcd(n, m%n);
}