#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. */

/*uses the curve y^2 = x^3 - t^2x + t^4. checks heights of (t, t^2) and (t^2, t^3). computes the determinant of the bilinear height matrix of these points, and determines whether or not they are linearly independent. finally, checks the sign of the functional equation. loops through many values of t*/


int main() 
{

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

  /*A and B are coeffs of the elliptic curve (both based on t)*/
  long  A, B, ctr, t, t_squared, sign, ind1,ind2, ind3, ind4, max_t;
  double det1, a11, a12, a22, a13, a23, a33, ht1, ht2, ht3, det2, det3, det4, pair;
  FILE *fp;
  long lTop;

  max_t = 1000; /*sets maximum value of t to calculate*/

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

  b1 = cgetr(MEDDEFAULTPREC);
  b2 = cgetr(MEDDEFAULTPREC);
  b3 = cgetr(MEDDEFAULTPREC);
  b4 = cgetr(MEDDEFAULTPREC);
  b5 = cgetr(MEDDEFAULTPREC);
  b6 = cgetr(MEDDEFAULTPREC);

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

  hp1 = cgetr(MEDDEFAULTPREC);
  hp2 = cgetr(MEDDEFAULTPREC);
  hp3 = 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);

  fp3=cgetg(3, t_VEC);
  fp3[1] = cgeti(MEDDEFAULTPREC);
  fp3[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);

  fp = fopen("r3dat10.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<=max_t; t++){
    
    t_squared = t*t;
    A = -16*t_squared;
    B = 64*t_squared;
    
    lTop = avma;

    /*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(5, (GEN) F[2]);
    gaffsg(0, (GEN) F[3]);
    gaffsg(A, (GEN) F[4]);
    gaffsg(B, (GEN) F[5]);
    
    E=initell(F, MEDDEFAULTPREC);
    
    /*Assign the vectors (t, t^2) and (t^2, t^3) to our points*/
    
    gaffsg(4*t, (GEN) fp1[1]);
    gaffsg(12*t, (GEN) fp1[2]);
    
    gaffsg((-4)*t, (GEN) fp2[1]);
    gaffsg(12*t, (GEN) fp2[2]);

    gaffsg(16*t*t, (GEN) fp3[1]);
    gaffsg(8*t*(8*t*t+1), (GEN) fp3[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);
    gaffect( (GEN) ghell2(E, fp3, MEDDEFAULTPREC), (GEN) hp3);

    ht1 = gtodouble( (GEN) hp1);
    ht2 = gtodouble( (GEN) hp2);
    ht3 = gtodouble( (GEN) hp3);
    
    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);
    gaffect( (GEN) bilhell(E,fp1, fp3,MEDDEFAULTPREC), (GEN) b4);
    gaffect( (GEN) bilhell(E,fp2, fp3,MEDDEFAULTPREC), (GEN) b5);
    gaffect( (GEN) bilhell(E,fp3, fp3,MEDDEFAULTPREC), (GEN) b6);

    a11 = gtodouble( (GEN) b1);
    a12 = gtodouble( (GEN) b2);
    a22 = gtodouble( (GEN) b3);
    a13 = gtodouble( (GEN) b4);
    a23 = gtodouble( (GEN) b5);
    a33 = gtodouble( (GEN) b6);

    /*calculate sign of functional equation*/

    gaffsg(1, (GEN) y);
    sign = ellrootno(E,y);

    avma = lTop;
    
    /*calculate determinants in C, take absolute values*/

    det1 = a11*a22*a33 + 2*a12*a23*a13 - a13*a13*a22 - a23*a23*a11 - a12*a12*a33;
    det2 = a11*a22 - a12*a12;
    det3 = a22*a33 - a23*a23;
    det4 = a11*a33 - a13*a13;

    if(det1 < 0)
      det1 = -det1;

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

    if(det1<0.001)
      ind1 = 0;
    else ind1 = 1;

    if(det2 < 0)
      det2 = -det2;

    /*determine whether or not first 2 points are linearly dependent*/
    
    if(det2<0.001)
      ind2 = 0;
    else ind2 = 1;


    if(det3 < 0)
      det3 = -det3;

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

    if(det3<0.001)
      ind3 = 0;
    else ind3 = 1;

    if(det4 < 0)
      det4 = -det4;

    /*determine whether or not first + last points are linearly dependent*/

    if(det4<0.001)
      ind4 = 0;
    else ind4 = 1;
    
    fprintf(fp, "%ld %f %f %f %f %f %f %f %ld %ld %ld %ld %ld\n", t, ht1, ht2, ht3, det1, det2, det3, det4, ind1, ind2, ind3, ind4, sign);
  }
  fclose(fp);
  return 0;
}