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

/*this is an improvement from t_1.c as it doesn't input rational points for the two fixed points since they will always be integers. */


int main() 
{

  /*declare pari vars*/
  /*pr and qr points, b1 – b3 are entries in bilinear height matrix*/
  GEN E, F, fp1, fp2, b1, b2, b3, pt, pt2, pt3, pt4, y, hp1, hp2;

  /*A is  coeffs of the elliptic curve (both based on t)*/
  long  A, ctr, t, t_squared, sign, ind;
  double det, a11, a12, a22, ht1, ht2;
  FILE *fp;

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

  b1 = cgetr(MEDDEFAULTPREC);
  b2 = cgetr(MEDDEFAULTPREC);
  b3 = cgetr(MEDDEFAULTPREC);

/*initialize the heights of the two forced points*/

  hp1 = cgetr(MEDDEFAULTPREC);
  hp2 = cgetr(MEDDEFAULTPREC);

  y = cgeti(32);

/*pari versions of the powers of t*/

  pt  = cgeti(32);
  pt2 = cgeti(32);
  pt3 = cgeti(32);
  pt4 = cgeti(32);

/*initialize two forced points as vectors*/

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

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


/*inintialize curves before assigning values*/

  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);
  
  gaffsg(1, (GEN) y);

  fp = fopen("bhtdat.c", "w"); /*opens the file for writing*/  

  /*  printf("This program returns the values of t for which the points (t, t^2) and (t^2, t^3) are linearly independent on an elliptic curve of the form y^2 = x^3 – t^2x + t^4.  Please enter an INTEGER for t.\n"); */
  
  for(t = 1; t<=1000; t++){
    
    t_squared = t*t;
    A = -t_squared;
    
    /*assign powers of t in pari*/
    gaffsg(t_squared, (GEN) pt2);
    gaffsg(t, (GEN) pt);
    gaffect(gmul( (GEN) pt2, (GEN) pt2), (GEN) pt4);
    gaffect(gmul( (GEN) pt2, (GEN) pt), (GEN) pt3);
    
    /* 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(A, (GEN) F[4]);
    gaffect((GEN) pt4, (GEN) F[5]);
    
    E=initell(F, MEDDEFAULTPREC);
    
    /*Assign the vectors (t, t^2) and (t^2, t^3) to our points*/
    
    gaffsg(t, (GEN) fp1[1]);
    gaffsg(t_squared, (GEN) fp1[2]);
    
    gaffsg(t_squared, (GEN) fp2[1]);
    gaffect( (GEN) pt3, (GEN) fp2[2]);
    
    /*calculate the heights of these points and turn them to C vars*/
    
    gaffect( (GEN) ghell2(E, fp1, MEDDEFAULTPREC), (GEN) hp1);
    gaffect( (GEN) ghell2(E, fp2, MEDDEFAULTPREC), (GEN) hp2);
    
    ht1 = gtodouble( (GEN) hp1);
    ht2 = gtodouble( (GEN) hp2);
    
    gaffect( (GEN) bilhell(E,fp1, fp1,MEDDEFAULTPREC), (GEN) b1);
    gaffect( (GEN) bilhell(E,fp1, fp2,MEDDEFAULTPREC), (GEN) b2);
    gaffect( (GEN) bilhell(E,fp2, fp2,MEDDEFAULTPREC), (GEN) b3);
    
    a11 = gtodouble( (GEN) b1);
    a12 = gtodouble( (GEN) b2);
    a22 = gtodouble( (GEN) b3);
    
    /*calculate determinant in C, take its absolute value*/

    det = a11*a22 - a12*a12;
    if(det<0)
      det = a12*a12 - a11*a22;

    /*determine whether or not points are linearly dependent*/

    if(det<0.001)
      ind = 0;
    else ind = 1;

    /*calculate sign of functional equation*/
    
    sign = ellrootno(E,y);

    fprintf(fp, "%ld %f %f %f %ld %ld\n", t, ht1, ht2, det, ind, sign);
  }
  fclose(fp);
  return 0;
}