// This is slopes2.
// This version was created on 29 June 2011.

// As a consequence of slopes1 completing successfully, we can
// restrict to the following parameter space:
// 1 \leq e_2 \leq 1.4
// e_2 \leq m \leq 2.5
// 1.7 \leq h \leq 4.0
// 0 \leq t \leq 0.5
// In this routine, we assume that e_3 = 51/32.
// We use Tool 1 either to reach a contradiction, or
// to prove that the number of exceptional slopes is at most 10 and 
// that the maximal intersection number between exceptional slopes is at most 8.
// If Tool 1 works for a particular value of e_3, then it works for
// larger values of e_3, and so the successful completion of this routine
// will imply that we may restrict, in later routines, to e_2 \leq e_3 \leq 51/32.

// We include standard libraries. We also include myDouble, which is used to
// define the Double variable type.

#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <math.h>
#include "myDouble.h"
#include "myDouble.cpp"
#include <fenv.h>

// The following variable records underflow and overflow exceptions.

int fetestexcept(int excepts);

// We now declare the variables that are used later in the routine.

Double overlapApprox( Double a, Double b, Double c );
Double psArcSin( Double x );
int crudeSlopeBound( Double L, Double hhh, Double m);
int fancySlopeBound( Double L, Double hhh, Double m, Double ttt);
int crudeIntBound( Double L, Double A);
bool q1slopes( Double L, Double hhh, Double m, Double ttt);
bool q2slopes( Double L, Double hhh, Double m, Double ttt);
bool q3slopes( Double L, Double hhh, Double m, Double ttt);
int minpatq1( Double L, Double hhh, Double m, Double ttt);
int maxpatq1( Double L, Double hhh, Double m, Double ttt);
int minpatq2( Double L, Double hhh, Double m, Double ttt);
int maxpatq2( Double L, Double hhh, Double m, Double ttt);
int minpatq3( Double L, Double hhh, Double m, Double ttt);
int maxpatq3( Double L, Double hhh, Double m, Double ttt);
int intnumber1( Double L, Double hhh, Double m, Double ttt);
int intnumber2( Double L, Double hhh, Double m, Double ttt);
int intnumber3( Double L, Double hhh, Double m, Double ttt);
int intnumber4( Double L, Double hhh, Double m, Double ttt);
int intnumber5( Double L, Double hhh, Double m, Double ttt);
int intnumber6( Double L, Double hhh, Double m, Double ttt);
int intnumber7( Double L, Double hhh, Double m, Double ttt);
int intnumber8( Double L, Double hhh, Double m, Double ttt);
int intnumber9( Double L, Double hhh, Double m, Double ttt);
int maxintersection( Double L, Double hhh, Double m, Double ttt);

// In Section 5, the function overlap(a,b,c) is defined to be
// the area of the overlap of two discs of radius a and b whose
// centres are separated by c. The following function is an
// approximation to this. It is a strict overestimate. This is
// because its only use is to find a lower bound for the area of
// the cusp torus, and in this lower bound, it always occurs with a
// minus sign in front. (See Theorems 5.2, 5.4 and 5.5.)
// For a justification of the formula used in this function, see
// Lemma 3.12 of [17].

Double overlapApprox( Double a, Double b, Double c)
{
	Double result, Pi, oneThird;
	
	Pi.value = 3.14159265358979323846;
    Pi.error = UU;
	oneThird.value = 0.33333333333333333333;
	oneThird.error = UU;
	
	result = a*a*((5.0*oneThird - Pi/2.0)*Pow((a*a-b*b+c*c)/(2.0*a*c),5)
		+Pow((a*a-b*b+c*c)/(2.0*a*c),3)/3.0 - 2.0*((a*a-b*b+c*c)/(2.0*a*c)) + Pi/2.0) +
			 b*b*((5.0*oneThird - Pi/2.0)*Pow((-a*a+b*b+c*c)/(2.0*b*c),5)
		+Pow((-a*a+b*b+c*c)/(2.0*b*c),3)/3.0 - 2.0*((-a*a+b*b+c*c)/(2.0*b*c)) + Pi/2.0);
		
	if (result.value + result.error < 0){
		result.value = 0.0;
		result.error = 0.0;
		}
		
	return result;
}

	
Double psArcSin( Double x )
// This function is an approximation to arcsin.
// It is a strict underestimate, and the upper side is likely to be wrong, but
// it's only use is to compute a lower bound on the critical slope length
// in Theorem 4.1.
{
	Double result;
	result = x + Pow(x,3)/6.0 + 3.0*Pow(x,5)/40.0 + 5.0*Pow(x,7)/112.0 + 35.0*Pow(x,9)/1152.0;
	return result;
}

int crudeSlopeBound( Double L, Double hhh, Double m)
// This function either returns an upper bound on the number of
// slopes of length at most L on a torus, or, if this upper bound is
// more than 10, it returns the value 11. The reason for the value 11
// is that the only significance of this function is whether it
// is more than 10.
// We assume that the torus has a parallelogram as fundamental domain, 
// with base length m, and height hhh. The base of the parallelogram 
// specifies a slope mu on the torus. The routine considers slopes 
// with intersection number 1 with mu, then with intersection number 2, and 
// then with intersection number 3. We assume that
// mu is a slope of shortest length on the cusp torus, and
// that L is at most 6. It is argued in Section 7 that, in this case,
// all curves with intersection number at least 4 with mu have length
// more than 6.
{
	int crudeslopes;
	Double rhs;
	double bigrhs;
	
	crudeslopes = 1;
	rhs = 4.0*(L*L - hhh*hhh)/(m*m);
	bigrhs = bigDsize(rhs);
	if (bigrhs < 16) crudeslopes = crudeslopes + 4;
	else if (bigrhs < 25) crudeslopes = crudeslopes + 5;
	else if (bigrhs < 36) crudeslopes = crudeslopes + 6;
	else if (bigrhs < 49) crudeslopes = crudeslopes + 7;
	else if (bigrhs < 64) crudeslopes = crudeslopes + 8;
	else if (bigrhs < 81) crudeslopes = crudeslopes + 9;
	else 
	{
		crudeslopes = 11;
		return crudeslopes;
	}

	rhs = 4.0*(L*L - 4.0*hhh*hhh)/(m*m);
	bigrhs = bigDsize(rhs);
	if (bigrhs < 0) crudeslopes = crudeslopes;
	else if (bigrhs < 4) crudeslopes = crudeslopes + 1;
	else if (bigrhs < 16) crudeslopes = crudeslopes + 2;
	else if (bigrhs < 36) crudeslopes = crudeslopes + 3;
	else if (bigrhs < 64) crudeslopes = crudeslopes + 4;
	else crudeslopes = crudeslopes + 5;

	rhs = 4.0*(L*L - 9.0*hhh*hhh)/(m*m);
	bigrhs = bigDsize(rhs);
	if (bigrhs < 0) crudeslopes = crudeslopes;
	else if (bigrhs < 1) crudeslopes = crudeslopes + 1;
	else if (bigrhs < 9) crudeslopes = crudeslopes + 2;
	else if (bigrhs < 16) crudeslopes = crudeslopes + 3;
	else crudeslopes = crudeslopes + 4;

	return crudeslopes;
}


int fancySlopeBound( Double L, Double hhh, Double m, Double ttt)
// This function is a slightly more sophisticated version of
// the previous function, crudeSlopeBound. In addition to specifying the
// base length and height of the parallelogram, it also takes as input
// the parameter ttt, which measures how sheared the parallelogram is.
// More specifically, we assume that three vertices of the parallelogram
// have co-ordinates (0,0), (m,0) and (m * ttt, hhh), where ttt is between
// 0 and 1/2. We also make the same assumptions as for the crudeSlopeBound function.
{
	int fancyslopes;
	double p;
	Double rhs;
	
	fancyslopes = 1;
	rhs = (L*L - hhh*hhh)/(m*m);	
	for ( p = -6; p < 6.5; p = p + 1){
		if ( (p + ttt) * (p + ttt) > rhs) fancyslopes = fancyslopes;
		else fancyslopes = fancyslopes + 1;
		}
		
	rhs = (L*L - 4.0*hhh*hhh)/(m*m);
	for ( p = -3; p < 3.5; p = p + 2){
		if ( (p + 2.0 * ttt) * (p + 2.0 * ttt) > rhs) fancyslopes = fancyslopes;
		else fancyslopes = fancyslopes + 1;
		}
	
	rhs = (L*L - 9.0*hhh*hhh)/(m*m);
	for ( p = -2; p < 2.5; p = p + 1){
		if ( p == 0 ) continue;
		else if ( (p + 3.0 * ttt) * (p + 3.0 * ttt) > rhs) fancyslopes = fancyslopes;
		else fancyslopes = fancyslopes + 1;
		}
		
	return fancyslopes;
}


int crudeIntBound( Double L, Double A)
// This gives a crude upper bound on the intersection number between slopes
// of length at least L on a torus with area A.
// We assume that A is at least 3.35 by Cao-Meyerhoff, and that
// L is at most 6. Thus, the intersection number is at most 10.
{
	int result;
	Double rhs;
	double bigrhs;
	
	rhs = (L*L)/A;
	bigrhs = bigDsize(rhs);
	
	result = 10;
	
	if (bigrhs < 10.0) result = 9;
	if (bigrhs < 9.0) result = 8;
	if (bigrhs < 8.0) result = 7;
	if (bigrhs < 7.0) result = 6;
	if (bigrhs < 6.0) result = 5;
	if (bigrhs < 5.0) result = 4;
	if (bigrhs < 4.0) result = 3;
	if (bigrhs < 3.0) result = 2;
	if (bigrhs < 2.0) result = 1;
	if (bigrhs < 1.0) result = 0;
	
	return result;
	
}

// The following functions determine whether there are slopes (p,q)
// with q = 1,2,3 and length at most L. If there are definitely no such slopes, the functions
// return the value FALSE. Otherwise, they return the value TRUE.

bool q1slopes( Double L, Double hhh, Double m, Double ttt)
{
	bool result;
	result = !((L*L - hhh*hhh)/(m*m) < (ttt*ttt));
	return result;
}

bool q2slopes( Double L, Double hhh, Double m, Double ttt)
{
	bool result;
	result = !((L*L - 4.0*hhh*hhh)/(m*m) < (2.0*ttt-1.0)*(2.0*ttt-1.0));
	return result;
}

bool q3slopes( Double L, Double hhh, Double m, Double ttt)
{
	bool result;
	result = !((L*L - 9.0*hhh*hhh)/(m*m) < (3.0*ttt-1.0)*(3.0*ttt-1.0));
	return result;
}
	
/* The following functions give the minimum and maximum values of p
such that (p,q) is a slope of length at most L, for q = 1, 2, 3.
They assume that there is at least one such slope. Otherwise, an
error message is given. They also assume the area of the torus is 
at least 3.35, and that L is at most 6. These assumptions imply, 
for example, that the minimal value of p, when q=1, is at least -6. */

int minpatq1( Double L, Double hhh, Double m, Double ttt)
{
	int result, p;
	Double rhs;
	
	result = 0;
	
	if (! q1slopes(L,hhh,m,ttt)) {
		printf("Error");
		exit(1);
	}

	rhs = (L*L - hhh*hhh)/(m*m);
	
	for ( p = 0; p > -6.5; p = p - 1){
		if ( !((p + ttt) * (p + ttt) > rhs)) result = p;
	}
	
	return result;
}

int maxpatq1( Double L, Double hhh, Double m, Double ttt)
{
	int result, p;
	Double rhs;
	
	result = 0;
	
	if (! q1slopes(L,hhh,m,ttt)) {
		printf("Error");
		exit(1);
	}

	rhs = (L*L - hhh*hhh)/(m*m);
	for ( p = 0; p < 6.5; p = p + 1){
		if ( !((p + ttt) * (p + ttt) > rhs)) result = p;
	}
	
	return result;
}


int minpatq2( Double L, Double hhh, Double m, Double ttt)

{
	int result, p;
	Double rhs;
	
	result = 0;

	if (! q2slopes(L,hhh,m,ttt)) {
		printf("Error");
		exit(1);
	}
	
	rhs = (L*L - 4.0*hhh*hhh)/(m*m);
	for ( p = -1; p > -3.5; p = p - 2){
		if ( !((p + 2.0 * ttt) * (p + 2.0 * ttt) > rhs)) result = p;
	}
	
	return result;
}



int maxpatq2( Double L, Double hhh, Double m, Double ttt)
{
	int result, p;
	Double rhs;
	
	result = 0;

	if (! q2slopes(L,hhh,m,ttt)) {
		printf("Error");
		exit(1);
	}
	
	rhs = (L*L - 4.0*hhh*hhh)/(m*m);
	for ( p = -1; p < 3.5; p = p + 2){
		if ( !((p + 2.0 * ttt) * (p + 2.0 * ttt) > rhs)) result = p;
	}
	
	return result;
}

int maxpatq3( Double L, Double hhh, Double m, Double ttt)
{
	int result, p;
	Double rhs;
	
	result = 0;

	if (! q3slopes(L,hhh,m,ttt)) {
		printf("Error");
		exit(1);
	}
	
	rhs = (L*L - 9.0*hhh*hhh)/(m*m);
	for ( p = -1; p < 1.5; p = p + 1){
		if ( p == 0 ) continue;
		if ( !((p + 3.0 * ttt) * (p + 3.0 * ttt) > rhs)) result = p;
	}
	
	return result;
}

int minpatq3( Double L, Double hhh, Double m, Double ttt)
{
	int result, p;
	Double rhs;
	
	result = 0;

	if (! q3slopes(L,hhh,m,ttt)) {
		printf("Error");
		exit(1);
	}
		
	
	rhs = (L*L - 9.0*hhh*hhh)/(m*m);
	for ( p = -1; p > -2.5; p = p - 1){
		if ( p == 0 ) continue;
		if ( !((p + 3.0 * ttt) * (p + 3.0 * ttt) > rhs)) result = p;
	}
	
	return result;
}

/* The following functions compute the intersection number between certain slopes.
For example, intnumber1 computes the intersection number between (p1,1)
and (p2,1), where p1 and p2 are the minimal and maximal values of p such that
(p,1) has length at most L. If there are no slopes (p,q) for
one of the relevant values of q, then the function returns 0.
*/

int intnumber1( Double L, Double hhh, Double m, Double ttt)
{
	int result;
	result = 0;
	
	if  (q1slopes(L,hhh,m,ttt))
		result = maxpatq1(L,hhh,m,ttt) - minpatq1(L,hhh,m,ttt);
	
	return result;
}
	
int intnumber2( Double L, Double hhh, Double m, Double ttt)
{
	int result;
	result = 0;
	
	if  (q1slopes(L,hhh,m,ttt) && q2slopes(L,hhh,m,ttt))
		result = abs(2*maxpatq1(L,hhh,m,ttt) - minpatq2(L,hhh,m,ttt));
	
	
	return result;
}

int intnumber3( Double L, Double hhh, Double m, Double ttt)
{
	int result;
	result = 0;
	
	if  (q1slopes(L,hhh,m,ttt) && q2slopes(L,hhh,m,ttt))
		result = abs(2*minpatq1(L,hhh,m,ttt) - maxpatq2(L,hhh,m,ttt));
	
	return result;
}

int intnumber4( Double L, Double hhh, Double m, Double ttt)
{
	int result;
	result = 0;
	
	if  (q2slopes(L,hhh,m,ttt))
		result = 2*(maxpatq2(L,hhh,m,ttt) - minpatq2(L,hhh,m,ttt));
	
	return result;
}


int intnumber5( Double L, Double hhh, Double m, Double ttt)
{
	int result;
	result = 0;
	
	if  (q2slopes(L,hhh,m,ttt) && q3slopes(L,hhh,m,ttt))
		result = abs(2*maxpatq3(L,hhh,m,ttt) - 3*minpatq2(L,hhh,m,ttt));
	
	return result;
}

int intnumber6( Double L, Double hhh, Double m, Double ttt)
{
	int result;
	result = 0;
	
	if  (q2slopes(L,hhh,m,ttt) && q3slopes(L,hhh,m,ttt))
		result = abs(2*minpatq3(L,hhh,m,ttt) - 3*maxpatq2(L,hhh,m,ttt));
	
	return result;
}

int intnumber7( Double L, Double hhh, Double m, Double ttt)
{
	int result;
	result = 0;
	
	if  (q1slopes(L,hhh,m,ttt) && q3slopes(L,hhh,m,ttt))
		result = abs(maxpatq3(L,hhh,m,ttt) - 3*minpatq1(L,hhh,m,ttt));
	
	return result;
}

int intnumber8( Double L, Double hhh, Double m, Double ttt)
{
	int result;
	result = 0;
	
	if  (q1slopes(L,hhh,m,ttt) && q3slopes(L,hhh,m,ttt))
		result = abs(minpatq3(L,hhh,m,ttt) - 3*maxpatq1(L,hhh,m,ttt));
	
	return result;
}

int intnumber9( Double L, Double hhh, Double m, Double ttt)
{
	int result;
	result = 0;
	
	if  (q3slopes(L,hhh,m,ttt))
		result = 3*(maxpatq3(L,hhh,m,ttt) - minpatq3(L,hhh,m,ttt));
	
	return result;
}

/* One of the above 9 functions is the maximum intersection number
between two slopes of length at most L. The following function computes
this maximum.
*/

int maxintersection( Double L, Double hhh, Double m, Double ttt)
{
	int result;
	result = 0;
	
	if (intnumber1(L,hhh,m,ttt) > result) result = intnumber1(L,hhh,m,ttt);
	if (intnumber2(L,hhh,m,ttt) > result) result = intnumber2(L,hhh,m,ttt);
	if (intnumber3(L,hhh,m,ttt) > result) result = intnumber3(L,hhh,m,ttt);
	if (intnumber4(L,hhh,m,ttt) > result) result = intnumber4(L,hhh,m,ttt);
	if (intnumber5(L,hhh,m,ttt) > result) result = intnumber5(L,hhh,m,ttt);
	if (intnumber6(L,hhh,m,ttt) > result) result = intnumber6(L,hhh,m,ttt);
	if (intnumber7(L,hhh,m,ttt) > result) result = intnumber7(L,hhh,m,ttt);
	if (intnumber8(L,hhh,m,ttt) > result) result = intnumber8(L,hhh,m,ttt);
	if (intnumber9(L,hhh,m,ttt) > result) result = intnumber9(L,hhh,m,ttt);
			
	return result;
}

// The following is the main routine.

int main()
{
  double ee2, mmtl, tt, hh;
  Double e2, e3, mtl, Area, t, h, critSlopeLength;
  double widthe2, widthe3, widthm, widtht, widthh;
  Double Pi;
  int badslopes, maxint;
   	
  printf("We analyze Exceptional slopes over a 4-dimensional parameter space,\n");
		
  /* mtl is the minimum translation length (first generator), en is the Euclideanized nth
	 ortholength, t is the second-generator offset, h is the height of
	 the second generator.  We compute upper bounds on the number of exceptional
	 surgeries and on the maximal intersection number between exceptional slopes. 
	 We achieve rigor by using Double's to carry out the calculations.  */

  widthe2 = 1.0/1024.0;
  widthe3 = 1.0/1024.0;
  widthm = 1.0/1024.0;
  widtht = 1.0/512.0;
  widthh = 1.0/512.0;
  Pi.value = 3.14159265358979323846;
  Pi.error = UU;
  
  for (ee2 = 1 + 0.0/512.0 + widthe2; ee2 < 1 + 13.0/32.0; ee2 = ee2 + 2.0 * widthe2){
	
	e2.value = ee2;
	e2.error = widthe2;
	
	// The following variable critSlopeLength is taken from Theorem 4.1.
	
	critSlopeLength = (Pi * e2)/psArcSin(e2/2.0);

	printf("e2 = %6.5f  critSlopeLength = %6.5f  \n",ee2,critSlopeLength.value);
	
	e3.value = 1 + 19.0/32.0;
	e3.error = widthe3;
			
	// The following is the lower bound on the area of the cusp torus given by Theorem 5.2.
	  
	  Area = Pi*e3*e3/2.0 + 2.0 *Pi*(e3/e2 - e3/2.0)*(e3/e2 - e3/2.0) -
				overlapApprox(e3/2.0,e3/2.0,e2) - 
				2.0 * overlapApprox(e3/2.0,e3/e2 - e3/2.0,1.0/e2);
	
      for (mmtl = ee2 + widthm; mmtl < 2 + 1.0/2.0 + widthm; mmtl = mmtl + 2.0*widthm){
	
	  mtl.value = mmtl;
	  mtl.error = widthm;

	// The following test determines whether we can eliminate this parameter point using
	// the functions crudeSlopeBound and crudeIntBound, which do not depend on the parameter t.

	  	  if ((crudeSlopeBound( critSlopeLength, Area/mtl, mtl) < 11) && 
		  (crudeIntBound( critSlopeLength, Area) < 9))
	    continue;

		for (tt = 0.0 + widtht; tt < 0.5 + widtht; tt = tt + 2.0*widtht){
		
		t.value = tt;
		t.error = widtht;
		
//			Let's see if we can get out of this without an h loop.

			if ((fancySlopeBound( critSlopeLength, Area/mtl, mtl, t) < 11)
				&& (maxintersection ( critSlopeLength, Area/mtl, mtl, t) < 9))
				continue;

			for (hh =  1 + 21.0/32.0 + widthh; hh < 4.0 + widthh; hh = hh + 2.0*widthh){

			h.value = hh;
			h.error = widthh;

			// If h is too small, h*mtl will be strictly lower than the known lower bound
			// on cusp area, which would be a contradiction.
			
			if (h*mtl < Area) continue;
			
			if ((fancySlopeBound( critSlopeLength, h, mtl, t) < 11)
				&& (maxintersection ( critSlopeLength, h, mtl, t) < 9))
				continue;
				
			badslopes = fancySlopeBound( critSlopeLength, h, mtl, t);
			maxint = maxintersection ( critSlopeLength, h, mtl, t);
		
	  printf("\nOoops.  It looks like we have a bit of a problem when \n");
	  printf("e2 = %6.5f, mtl = %6.5f\n",
		 ee2,mmtl);
	  printf("h = %6.5f, t = %6.5f\n",
		 hh,tt);
	  printf("where the number of potentially bad slopes is %4d\n",badslopes);
	  printf("and the potential maximal intersection number is %4d\n",maxint);
	  printf("and the area is   ");
	  printDouble(Area);
	  printf("\n\n");
	  exit(1);
		  }
		}
      }	
	  }	
          printf("\nDone with the loops.\n");
		  printf("underflowtestexcept = %x \n", fetestexcept(FE_UNDERFLOW));
		  printf("overflowtestexcept = %x \n", fetestexcept(FE_OVERFLOW));

// The presence of underflow at any time in the program results in 
// underflowtestexcept = 10 being printed.
// The presence of overflow results in overflowtestexcept = 8 being printed.
// If neither underflow nor overflow occurs then the lines are printed out with "=0" in
// both cases.

		  
}