/*
******tries to find an extra point(s) independent on y^2 = x^3 - t^2*x + t^4******


v 1.1 Precursor to a simple program to determine rational solutions to elliptic curves:
This program returns rational x-values and floating point squares of y-values of points on an elliptic curve with user-specified integer coefficients.  The x-values of the points returned all have naive height less than a user-specified value.
  v 1.2 Computes y value instead of y^2; uses pow function instead of power 
  v 1.3 Added RAT.h functions, which compute x and y_squared as rational numbers; 
computes y_squared as a rational number instead of a floating point; 
program must be run with RATbetter.c functions; 
says whether this is a rational point on curve or not
  v 1.4 Returns only those points that are rational, have low height, and are on the curve
  v 1.5 Ooops!  The previous versions only work for curves y^2=Ax^3+B.  This problem has been fixed, and the new version returns points on curves y^2=x^3+Ax+B.
  v 1.5-3 Checks to see whether number is a square mod p
  v 1.5-4 Includes sqrt2 function which is alternative method for calculating square roots... based on Rudin, ch. 3, exercise 16... This should help when dealing with large numbers ????? (I don't really see how it will help at the moment, but it made sense when Steve discussed this with me.).  Please note that the accuracy can be changed around line 6 where it says #define Accuracy; this sayshow accurate the square root should be.  Also, #define MaxIterations says max number of iterations to be done by sqrt2.  Both of these variables would probably be useful in other functions, too.*/

/* This is now the program to determine low heights of the family of curves: */
/* Y^3 + X^3 = A. It asks for what rational point search method to use, determines */
/* the number of points using this method and the rank of the curve, and prints all */
/* the information to a file */

/* To Do:
   Check what's wrong with methods
   Fix method 2
   - Make the switch from t to (9t+1)^2
   Check that 9t+1 is square-free (?)
   - progress bar (sort of done)
   - Check canonical height is bounded
   - Track S(N) = Si(N) + St(N) where Si and St are 
     the number solutions that are independent and torsion, respectively
   N -> log N [x_1, c^2 -> e^(x_1), e^(c^2)] (?)
   - eliminate methods 1 and 2
*/

/* Note: Y^3 + X^3 = A  <==> Y^2 = X^3 + 2^4*(-3)^3*(9t+1)^2 [see Miller] */

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

/*accuracy for calculating square root with sqrt2 function*/
#define Accuracy 0.00001
/* accuracy for checking a number is 0 - 10^-9 */
#define Accuracy0 0.0000000001
/*maximum number of iterations to be done by sqrt2 function*/
#define MaxIterations 1000
/* number of primes checked in method 3
 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);

/*Functions from RATbetter.c to simplify numerator and denominator*/
int SIMPLIFYnum (long int, long int);	/*fcn simplifies numerator*/
int SIMPLIFYdenom (long int, long int); /*fcn simplifies denominator*/
int isSquare(int); /* says whether an integer is a square */
double sqrt2(double); /*algorithm for calculating square of large number*/

/* calculates if a number is an integer */
int isInt(double x)
{
  int x2=0;
  x2=floor(x);
  if ( (x/x2) <= (1+pow(10, -12)) )
    return 1;
  else
    return 0;
}

/* 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;
  GEN sp1, sp2, b11, b12, b13, b22, b23, b33;
  double a11, a12, a13, a22, a23, a33, det;
  long ind;

  /*N=maximum naive height*/
  /*A and B are coeffs of the elliptic curve*/
  long  N, A, B, ctr, GCD=0, c_squared=0, c_cubed=0;
  long templ=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, x_2, c, lower_bound;
  Rational x, x_squared, x_cubed, y, y_squared;
  long tempx3p, tempc3p, tempc6p, tempk3p, tempt2p, p, tempt1p, tempxp, tempt1xp, tempt1xkp;
  double doubleD;
  double canonical_height; //canonical height as computed by ghell2
  double d=0.0, temp=0.0;
  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 ismp[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},
		     {7, 0, 1, 4, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
		     {11, 0, 1, 4, 3, 9, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1}, 
		     {13, 0, 1, 9, 10, 3, 4, 12, -1, -1, -1, -1, -1, -1, -1, -1, -1},
		     {17, 0, 1, 9, 13, 15, 2, 4, 8, 16, -1, -1, -1, -1, -1, -1, -1},
		     {19, 0, 1, 5, 4, 16, 11, 17, 7, 6, 9, -1, -1, -1, -1, -1, -1},
		     {23, 0, 1, 12, 8, 6, 4, 3, 18, 2, 16, 13, 9, -1, -1, -1, -1},
		     {29, 0, 1, 22, 6, 5, 25, 13, 9, 20, 16, 4, 24, 23, 7, 28, -1},
                     {31, 0, 1, 16, 8, 25, 9, 4, 7, 28, 20, 2, 19, 18, 14, 5, 10} };
  long ctr1=0, ctr2=0, temp2=0, k=0, D=0;
  /* pass verifies that a point passed the first 10 primes */
  int pass=0;

  double progress=0, increment=0;

  /* Si and St are the number of solutions in the independent and torsion groups, respectively */
  long Si=0, St=0;

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

  pari_init(10000000, 2);
  printf("Initializing pari variables\n");
  tempg=cgeti(32);
  a=cgeti(32);

  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);

  /*new code end*/

  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("Checking new mod function:\n");
  printf("mod(1039,19)=%ld mod(-1894,19)=%ld mod(-242,31)=%ld mod(242,31)=%ld mod(-423,8)=%ld mod(423,8)=%ld\n",
	 mod(1039, 19), mod(-1894, 19), mod(-242, 31), mod(242, 31), mod(-423, 8), mod(423, 8));
  */
 
  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 - t^2*x + t^4.  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\n", N, t_lower, t_upper);
  
  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 = 1;
      // t2=16*t*t;

      t1=-t*t;
      t2=t*t*t*t;
      fprintf(fp, "Curve: y^2 = x^3 - %ld*x + %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]);
      
      /*
	printf("\tChecking F\n");
	for (ctr=1; ctr<6; ctr++)
	{
	printf("F[%ld]=", ctr);
	outbrute((GEN) F[ctr]); printf("\n");
	}
      */
      
      E=initell(F, MEDDEFAULTPREC);
      
      /*
	printf("\tChecking E\n");
	for (ctr=1; ctr<14; ctr++)
	{
	printf("E[%ld]=", ctr);
	outbrute((GEN) E[ctr]); printf("\n");
	}
      */

      /*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]);
      
      /*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*/
      
      /*lower_bound = pt of order 2 in E(R) w/ min x coord...HOW DO YOU DO THIS??? For now, let lower_bound = 1*/
      lower_bound = 1;
      
      /* Debugging statement */
      // printf("N=%ld, t2=%ld\n\n",N, t2);
      
      fprintf(fp, "x\ty\tcanonical height\tdeterminant\tindependent?\n");
      
      /* 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. */
      
      // printf("Third method results:\n");
      for (x_1 = (-1 * N); x_1 <= N;  x_1++) {
	for (c = lower_bound; c*c <=N; c++)  
	  {
	    // if (x_1 == -1071 && c == 10) printf("\nHELLO WORLD");
	    c_squared=c*c;
	    c_cubed=c_squared*c;
	    pass=0;
	    if ( gcd(x_1*x_1,c) ==1)
	      {
		// if (x_1 == -1071 && c == 10) printf("\nHELLO WORLD2");
		for (ctr1=0; ctr1<num_primes; ctr1++)            /* cycle through primes in smp[][] */
		  {
		    
		    // printf("\n\t\tx_1=%ld c=%ld c^2=%ld c^3=%ld p=%ld ", 
		    //           x_1, c, c_squared, c_cubed, smp[ctr1][0]); 
		    
		    if ( gcd(c, smp[ctr1][0]) == 1)
		      {
			// if (x_1 == -1071 && c == 10) printf("\nHELLO WORLD2a"); 
			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;
			  }
			
			// printf("c^2 = %ld inv = %ld p = %ld\n", temp2, k,smp[ctr1][0]);
			
			/* 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 = tempc6p * (tempx3p * tempk3p + tempt2p);
			// D = mod(tempx3p * tempk3p, p) + tempt1xkp + tempt2p;
			D = mod(D, p);
			
			/* D=c_cubed*c_cubed*( (x_1*x_1*x_1*k*k*k) + t2); */
			
			// printf("k=%ld D=%ld ", k, D);
			
			/* 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 (x_1==-1071 && c==10) printf("\n\t\tsmp[%ld][%ld]=%ld D=%ld", 
			    //ctr1, ctr2, smp[ctr1][ctr2], (D%smp[ctr1][0]));
			    if ( mod(D, smp[ctr1][0]) == smp[ctr1][ctr2])
			      pass++; 
				//printf(" smp[%ld][%ld]=%ld ", ctr1, ctr2, smp[ctr1][ctr2]); 
			  }
			  
			// printf("pass=%d\n", pass);
			  
			// if (x_1 == -1071 && c == 10) printf("\n Here = %ld, pass = %d\n", ctr1, pass);
			
			/* if D is a square mod p, for first 10 primes, pass should be 10 */
			if (pass==num_primes)
			  {
			    // if (x_1 == -1071 && c == 10) printf("\nHELLO WORLD3\n");
			    /* calculates the numerator of Y = D^(1/2)/c^3 */
			    /* so D=x_1^3 + t2*c^6 */
			    // doubleD=(1.0*x_1*x_1*x_1)+(1.0*t1*x_1*c*c_cubed)+(1.0*t2*c_cubed*c_cubed);
			    doubleD=(1.0*x_1*x_1*x_1)+(1.0*t2*c_cubed*c_cubed);
			    // printf("doubleD = %f\n", doubleD);
			    
			    /* templ is the integer closest to D which is less than D */
			    // templ=floor(sqrt(1.0*doubleD));
			    /*templ = 100;*/
			    doublesqrtD = sqrt(doubleD);
			    longD1 = floor(doublesqrtD);
			    sqrtdifference = doublesqrtD - 1.0*longD1; //This is a double
			    
			    /* trouble here - templ was exceeding long */
			    
			    // if (x_1==-1071 && c==10) printf("\n\t\t\tnew D=%ld templ=%ld\n", D, templ);
			    
			    /* 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)
			      {
				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)
				      {
					if ( gtodouble((GEN) ghell2(E,pr,MEDDEFAULTPREC)) < Accuracy0 )
					  St++; //solutions in torsion gp
					else
					  Si++; //solutions in free gp
					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);
					/*
					  outbrute((GEN) fp1[1]); printf("/");
					  outbrute((GEN) fp1[2]); printf("\t\t");
					  outbrute((GEN) fp2[1]); printf("/");
					  outbrute((GEN) fp2[2]); printf("\n\n");
					*/
					
					/*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*/
				      }
				  }
			      }
			  }
		      }
		    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");
      /* end of third method */
      
      /*
	for (ctr1=0; ctr1<10; ctr1++)
	for (ctr2=0; ctr2<17; ctr2++)
	{
	if (smp[ctr1][ctr2] < 0) printf("%ld ", smp[ctr1][ctr2]);
	if (smp[ctr1][ctr2] >=0) printf("%ld  ", smp[ctr1][ctr2]);
	if (ctr2==16)
	printf("\n");
	}
      */
      fprintf(fp, "S(%ld)=%ld\tSi(%ld)=%ld\tSt(%ld)=%ld\n", N, (Si+St), N, Si, N, St);
      fprintf(fp, "-------------------------------------------------\n");
      Si=0; St=0;
    }
  fprintf(fp, "======================================================\n");
  fprintf(fp, "======================================================\n");
  fclose(fp);
  // printf("\n\n %d\n", gcd(-1071,10));
  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);
} 



/* Functions from RATbetter.c */

long int simplifynum (long int num, long int denom) /*simplify numerat.*/
{
if (num==denom)
	{
	num=1;
	}	/*if num=denom, simplifies num is 1*/
else 
	{
	num = ((num)/(gcd (num, denom))); /*make num relatively prime to
					denom*/
	}
return num;
}

long int simplifydenom (long int num, long int denom) /*simpl. denom*/
{
denom=(simplifynum (denom, num)); /*make denom prime relative to numer.*/
return denom;
}

Rational RATinit (int p, int q)	/*put p &q in Rational structure*/
{
Rational init;
init.p=simplifynum (p, q);
init.q=simplifydenom (p, q);
return init;	/*initialize simplified Rational*/
}

void RATshow (Rational RAT1) /*print Rational "p/q"*/
{
printf ("%ld/%ld", simplifynum (RAT1.p, RAT1.q), simplifydenom (RAT1.p,
RAT1.q));	/*print simplified Rational*/
}

Rational RATadd (Rational RAT1, Rational RAT2)	/*perform Rat addition*/
{
Rational sum;
long int common, a, b;
a=RAT1.p; b=RAT2.p;	/*set a and b equal to RAT1.p and RAT2.p, resp.*/
RAT1.p=simplifynum (RAT1.p, RAT1.q); /*simplify RAT1*/
RAT1.q=simplifydenom (a, RAT1.q);   /*simplify RAT1*/
RAT2.p=simplifynum (RAT2.p, RAT2.q);/*simplify RAT2*/
RAT2.q=simplifydenom (b, RAT2.q);  /*simplify RAT2*/
common=gcd (RAT1.q, RAT2.q);	/*find gcd of denoms of RAT1 and RAT2*/
RAT1.p=(RAT2.q/common)*RAT1.p;	/*compute numerator of RAT1*/
RAT2.p=(RAT1.q/common)*RAT2.p;	/*compute numerator of RAT2*/
RAT1.q=(RAT2.q/common)*RAT1.q;	/*compute denominator of RAT1*/

sum.p=((RAT1.p)+(RAT2.p));	/*compute numerator of sum*/
sum.q=RAT1.q;			/*denominator of sum=denom of RAT1*/
return sum;
}

Rational RATmul (Rational RAT1, Rational RAT2) /*multiply 2 RAT nums*/
{
long int a, b, c, d;
Rational product;
a=RAT1.p; b=RAT1.q; c=RAT2.p; d=RAT2.q;

RAT1.p=(simplifynum ((simplifynum (a, b)), d));	/*smplfy num of RAT1*/
RAT1.q=(simplifydenom (c, simplifydenom (a, b))); /*smplfy denom of RAT1*/

RAT2.p=(simplifynum (simplifynum (c, simplifydenom (a, b)), 
	simplifydenom (simplifynum (a,b), d))); /*smplfy num of RAT2*/
RAT2.q=(simplifydenom (simplifynum (c, simplifydenom (a,b)),
	simplifydenom (simplifynum (a,b), d))); /*smplfy denom of RAT2*/
product.p=(RAT1.p)*(RAT2.p);	/*compute num of product*/
product.q=(RAT1.q)*(RAT2.q);	/*compute denom of product*/
return product;
}

int isSquare(int p)
{
  if(sqrt(p)==(int)sqrt(p))
    return 1;
  else
    return 0;
}

/*computes square root; decreasing monotonically... So gives UPPER BOUND for square root; Algorithm from ch.3, exercise 16, Rudin.*/
double sqrt2( double a )
{
  double x_prev=a, x_current=0.0;
  int n = 0;
  while( ((x_prev-x_current)>Accuracy) && (++n<=MaxIterations))
    {
      x_current=0.5*(x_prev + a/x_prev);
    }
  return x_current;
}