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

// With the successful completion of the programs slopes1, slopes2, slopes3 and slopes4, 
// we may restrict to a compact parameter space:
// 1 \leq e_2 \leq 1.4
// e_2 \leq e_3 \leq 1.5
// e_3 \leq e_4 \leq 1.6
// e_2 \leq m \leq 2.5
// 1.7 \leq h \leq 4.0
// 0 \leq t \leq 0.5
// In this routine, we divide the parameter space into small regions,
// and use Tools 1, 2 and 3 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.

// 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);
Double		circRadiusSquared		( Double h, Double m, Double t);
Double		baseTriangleOne			( Double h, Double m, Double t);
Double		baseTriangleTwo			( Double h, Double m, Double t);
Double		baseTriangleThree		( Double h, Double m, Double t);
Double		offsetTriangleOne		( Double h, Double m, Double t);
Double		offsetTriangleTwo		( Double h, Double m, Double t);
Double		offsetTriangleThree		( Double h, Double m, Double t);
Double		altitudeTriangleOne		( Double h, Double m, Double t);
Double		altitudeTriangleTwo		( Double h, Double m, Double t);
Double		altitudeTriangleThree	( Double h, Double m, Double t);
Double		distanceSquared			( Double base, Double height, Double offset, Double e3);
bool		passe3CircRadius		( Double h, Double m, Double t, Double e2, Double e3);
bool		passTriangleOne			( Double h, Double m, Double t, Double e2, Double e3);
bool		passTriangleTwo			( Double h, Double m, Double t, Double e2, Double e3);
bool		passTriangleThree		( Double h, Double m, Double t, Double e2, Double e3);
bool		passTool3		( Double h, Double m, Double t, Double e2, Double e3);
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);
Double AreaBound( Double e2, Double e3, Double e4);

// We now define functions used in Tools 2 and 3 

Double circRadiusSquared( Double h, Double m, Double t)
// This is the square of the circumradius of the triangle with vertices
// O = (0,0), A = (m,0), B = (mt,h). (See Figure 10)
{
	Double result;
	result = ((t*t*m*m)+(h*h))*(((1-t)*(1-t)*m*m)+(h*h))/(4.0*h*h);
	return result;
}

Double baseTriangleOne( Double h, Double m, Double t)
// This is the length of AB.
{
	Double result;
	result = sqrt( (h*h) + (m*m*(1-t)*(1-t)) );
	return result;
}

Double baseTriangleTwo( Double h, Double m, Double t)
// This is the length of OB.
{
	Double result;
	result = sqrt( (h*h) + (m*m*t*t) );
	return result;
}

Double baseTriangleThree( Double h, Double m, Double t)
// This is the length of OA.
{
	return m;
}

Double offsetTriangleOne( Double h, Double m, Double t)
// If O' denotes the orthogonal projection of O onto AB,
// this is the length of AO'.
{
	Double result;
	result = (m*m*(1-t)) / baseTriangleOne(h,m,t);
	return result;
}

Double offsetTriangleTwo( Double h, Double m, Double t)
// If A' denotes the orthogonal projection of A onto OB,
// this is the length of OA'.
{
	Double result;
	result = (m*m*t) / baseTriangleTwo(h,m,t);
	return result;
}

Double offsetTriangleThree( Double h, Double m, Double t)
// If B' denotes the orthogonal projection of B onto OA,
// this is the length of OB'.
{
	Double result;
	result = m*t;
	return result;
}

Double altitudeTriangleOne( Double h, Double m, Double t)
// This is the length of the altitude of OAB, with base AB.
{
	Double result;
	result = h*m/baseTriangleOne(h,m,t);
	return result;
}

Double altitudeTriangleTwo( Double h, Double m, Double t)
// This is the length of the altitude of OAB, with base OB.
{
	Double result;
	result = h*m/baseTriangleTwo(h,m,t);
	return result;
}

Double altitudeTriangleThree( Double h, Double m, Double t)
// This is the length of the altitude of OAB, with base OA.
{
	return h;
}

Double distanceSquared( Double base, Double height, Double offset, Double e3)
// Let XYZ be a triangle with given base XY, height and offset, where
// the offset denotes distance between X and the orthogonal projection of
// Z onto XY. Suppose that e3 is definitely more than twice the
// length of XY. (Otherwise, this function evaluates the square root
// of a potentially negative number, and the routine ends with an error message.)
// Let P be the point that has distance e3 from both X and Y, and that
// lies on the Z side of XY. This function computes the square of
// the distance between P and Z.
// Note that if it is applied to the triangle OAB, then this
// is the square of the distance between O and one of the distinguished
// points in Figure 10.
{
	Double result;
	result = ((base/2.0) - offset)*((base/2.0) - offset)
	+(height - sqrt((e3*e3)-(base*base/4.0)))*(height - sqrt((e3*e3)-(base*base/4.0)));
	return result;
}


// The following is Test 1 of Tool 3.
// It returns the value FALSE if the circumradius of OAB is definitely 
// less than e3. Otherwise, it returns the value TRUE.

bool passe3CircRadius( Double h, Double m, Double t, Double e2, Double e3)
{
	bool result;
	Double e3squared;
	e3squared = (e3 * e3);
	result = !(circRadiusSquared(h,m,t) < e3squared);
	return result;
}


// The following three tests form part of Test 2 of Tool 3.
// They check whether the six distinguished points in Figure 10
// have distance at least e2 from the origin O.
// Each test checks one of the distinguished points. (Note that we
// only need to check three such points, because for each
// distinguished point, its reflection in the origin is also a
// distinguished point.) If this distinguished point has distance
// definitely less than e2 from the origin, the function returns
// the value FALSE. Otherwise, it returns the value TRUE.

bool passTriangleOne( Double h, Double m, Double t, Double e2, Double e3)
{
	bool result;
	result = !( 
		distanceSquared( baseTriangleOne(h,m,t), altitudeTriangleOne(h,m,t), offsetTriangleOne(h,m,t), e3)
		< (e2*e2)
	);
	return result;
}

bool passTriangleTwo( Double h, Double m, Double t, Double e2, Double e3)
{
	bool result;
	result = !( 
		distanceSquared( baseTriangleTwo(h,m,t), altitudeTriangleTwo(h,m,t), offsetTriangleTwo(h,m,t), e3)
		< (e2*e2)
	);
	return result;
}

bool passTriangleThree( Double h, Double m, Double t, Double e2, Double e3)
{
	bool result;
	result = !( 
		distanceSquared( baseTriangleThree(h,m,t), altitudeTriangleThree(h,m,t), offsetTriangleThree(h,m,t), e3)
		< (e2*e2)
	);
	return result;
}

// The following is Tool 3. If Test 1 of Tool 3 or any of the checks
// that form Test 2 of Tool 3 returns the value TRUE, then this function
// returns the value TRUE. Otherwise, it returns the value FALSE.
// Note that the function passe3CircRadius is evaluated first.
// If it returns the value TRUE, the routine does not evaluate
// the functions used in Test 2 of Tool 3. On the other hand, if it returns
// the value FALSE, then the circumradius of OAB is definitely
// less than e3, and hence the function distanceSquared that is
// called in Test 2 does not generate an error.

bool passTool3( Double h, Double m, Double t, Double e2, Double e3)
{
	bool result;
	result = (passe3CircRadius(h,m,t,e2,e3) 
	|| passTriangleOne(h,m,t,e2,e3)
	|| passTriangleTwo(h,m,t,e2,e3)
	|| passTriangleThree(h,m,t,e2,e3));
	return result;
}

// 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 lower bound on the area of the cusp torus that is
// provided by Theorem 5.4, and Theorem 5.5 if applicable.

Double AreaBound( Double e2, Double e3, Double e4)
{
	Double AreaBound1;
	Double AreaBound2;
	Double result;
	Double Pi;
	Pi.value = 3.14159265358979323846;
	Pi.error = UU;
 
	AreaBound1 = Pi*e4*e4/2.0 + 2.0 *Pi*(e4/e2 - e4/2.0)*(e4/e2 - e4/2.0) + 
				2.0*Pi*(e4/e3 - e4/2.0)*(e4/e3 - e4/2.0) - 
				overlapApprox(e4/2.0,e4/2.0,e2) - overlapApprox(e4/2.0,e4/2.0,e3) - 
				2.0 * overlapApprox(e4/2.0,e4/e2 - e4/2.0,1.0/e2) - 
				2.0 * overlapApprox(e4/2.0,e4/e3 - e4/2.0,1.0/e3);
		
	AreaBound2 = Pi*e4*e4/2.0 + 2.0 *Pi*(e4/e2 - e4/2.0)*(e4/e2 - e4/2.0) + 
				2.0*Pi*(e4/e3 - e4/2.0)*(e4/e3 - e4/2.0) -
				2.0 * overlapApprox(e4/2.0,e4/e2 - e4/2.0, e3/e2) - 
				2.0 * overlapApprox(e4/e2 - e4/2.0, e4/e3 - e4/2.0, 1.0/(e2*e3)) - 
				2.0 * overlapApprox(e4/e3 - e4/2.0, e4/2.0, e2/e3) - 
				overlapApprox(e4/2.0,e4/2.0,e2) - 
				2.0 * overlapApprox(e4/2.0,e4/e2 - e4/2.0,1.0/e2);
				
	if (smallDsize(AreaBound1) < smallDsize(AreaBound2)) result = AreaBound1; 
		else result = AreaBound2;

	// We now improve this using Theorem 5.5, if applicable.

	if ((e2 + 0.9 > e4*e4) && (1.0/(e4*e2) - (e4/e2) + e4/2.0 - 0.1 > 0)){
	    result = result + 2.0*Pi*(1.0/(e4*e2) - (e4/e2) + e4/2.0)*(1.0/(e4*e2) - (e4/e2) + e4/2.0);
		result = result - 2.0*overlapApprox(e4/2.0, 1.0/(e4*e2) - (e4/e2) + e4/2.0, 1/e4);
		}
	  
	return result;
}
	

// The following is the main routine.

int main()
{
  double ee2, ee3, ee4, mmtl, tt, hh;
  Double e2, e3, e4, mtl, Area, Area0, t, h, critSlopeLength;
  double widthe2, widthe3, widthm, widtht, widthh;
  Double Pi;
  int badslopes, maxint;
   	
  printf("We analyze exceptional slopes over a 6-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.  */

  // We now repeat the same set of code three times, with the only difference
  // being the range of e_2 that is considered. In the first part, we let e_2
  // vary from 1.0 to slightly more than 1.2, and the width of each e_2 interval 
  // and e_3 interval is 1/1024.  In the second part, e_2 increases to slightly 
  // more than 1.28, but now the e_2 and e_3 interval sizes are much smaller:
  // 1/8192 and 1/4096 respectively. In the final part, e_2 increases up to 1.4,
  // but now the interval sizes have increased to 1/512. The reason for doing this
  // is of course to reduce running time.
 
  
  widthe2 = 1.0/2048.0;
  widthe3 = 1.0/2048.0;
  widthm = 1.0/512.0;
  widtht = 1.0/256.0;
  widthh = 1.0/256.0;
  Pi.value = 3.14159265358979323846;
  Pi.error = UU;
  
  for (ee2 = 1 + 0.0/32.0 + widthe2; ee2 < 1 + 206.0/1024.0; ee2 = ee2 + 2.0 * widthe2){
	
	e2.value = ee2;
	e2.error = widthe2;
	
	// The following variable critSlopeLength is taken from Theorem 4.1.
    // Note that e2 is definitely less than sqrt(2);
	
	critSlopeLength = (Pi * e2)/psArcSin(e2/2.0);
	
	printf("e2 = %6.5f  critSlopeLength = %6.5f  \n",ee2,critSlopeLength.value);
	
    for (ee3 = ee2; ee3 < 1 + 1.0/2.0 + widthe3; ee3  = ee3 + 2.0 * widthe3){
	  e3.value = ee3;
	  e3.error = widthe3;

	  // We now compute a lower bound on cusp area, which assumes that e4 = e3,
	  // and store it in the variable Area0. Later, we will compare this with lower bounds
	  // on cusp area, when e4 > e3.
		
	Area0 = AreaBound( e2, e3, e3);
	
		for (ee4 = ee3; ee4 < 1 + 19.0/32.0 + widthe3; ee4 = ee4 + 2.0 * widthe3){
		
		e4.value = ee4;
		e4.error = widthe3;
		
		Area = AreaBound( e2, e3, e4);
		
	// In the following, we compare the lowest possible values of Area0 and
	// and Area. If the latter (when e4 is at least e3) is bigger than the former (when
	// e4 = e3), we can go to the next value of e4.
	// Note that when e4 = e3, then we do not pass to the next value of e4, 
        // and so in this case, we continue to further tests.

	if (smallDsize(Area0) < smallDsize(Area)){
		continue;
		}
		
      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;
				
//			In the following, we implement Tools 2 and 3.

			if ( circRadiusSquared( h, mtl, t) < e2*e2)
				continue;	
			if ( !passTool3(h, mtl, t, e2, e3)) 
				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, e3 = %6.5f, e4 = %6.5f, mtl = %6.5f\n",
		 ee2,ee3,e4.value,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);
		  }
		}
      }	
	  }	
    }
  }
		  


  widthe2 = 1.0/16384.0;
  widthe3 = 1.0/8192.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 + 206.0/1024.0 + widthe2; ee2 < 1 + 9.0/32.0; ee2 = ee2 + 2.0 * widthe2){
	
	e2.value = ee2;
	e2.error = widthe2;

	// The following variable critSlopeLength is taken from Theorem 4.1.
    // Note that e2 is definitely less than sqrt(2);
	
	critSlopeLength = (Pi * e2)/psArcSin(e2/2.0);
	
	printf("e2 = %6.5f  critSlopeLength = %6.5f  \n",ee2,critSlopeLength.value);
	
    for (ee3 = ee2; ee3 < 1 + 1.0/2.0 + widthe3; ee3  = ee3 + 2.0 * widthe3){
	  e3.value = ee3;
	  e3.error = widthe3;

	  // We now compute a lower bound on cusp area, which assumes that e4 = e3,
	  // and store it in the variable Area0. Later, we will compare this with lower bounds
	  // on cusp area, when e4 > e3.
		
	Area0 = AreaBound( e2, e3, e3);
	
		for (ee4 = ee3; ee4 < 1 + 19.0/32.0 + widthe3; ee4 = ee4 + 2.0 * widthe3){
		
		e4.value = ee4;
		e4.error = widthe3;
		
		Area = AreaBound( e2, e3, e4);
		
	// In the following, we compare the lowest possible values of Area0 and
	// and Area. If the latter (when e4 is at least e3) is bigger than the former (when
	// e4 = e3), we can go to the next value of e4.
	// Note that when e4 = e3, then we do not pass to the next value of e4, 
        // and so in this case, we continue to further tests.

	if (smallDsize(Area0) < smallDsize(Area)){
		continue;
		}
		
      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;
				
//			In the following, we implement Tools 2 and 3.

			if ( circRadiusSquared( h, mtl, t) < e2*e2)
				continue;	
			if ( !passTool3(h, mtl, t, e2, e3)) 
				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, e3 = %6.5f, e4 = %6.5f, mtl = %6.5f\n",
		 ee2,ee3,e4.value,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);
		  }
		}
      }	
	  }	
    }
  }
		  


  widthe2 = 1.0/1024.0;
  widthe3 = 1.0/1024.0;
  widthm = 1.0/512.0;
  widtht = 1.0/256.0;
  widthh = 1.0/256.0;
  Pi.value = 3.14159265358979323846;
  Pi.error = UU;
  
  for (ee2 = 1 + 9.0/32.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.
    // Note that e2 is definitely less than sqrt(2);
	
	critSlopeLength = (Pi * e2)/psArcSin(e2/2.0);
	
	printf("e2 = %6.5f  critSlopeLength = %6.5f  \n",ee2,critSlopeLength.value);
	
    for (ee3 = ee2; ee3 < 1 + 1.0/2.0 + widthe3; ee3  = ee3 + 2.0 * widthe3){
	  e3.value = ee3;
	  e3.error = widthe3;

	  // We now compute a lower bound on cusp area, which assumes that e4 = e3,
	  // and store it in the variable Area0. Later, we will compare this with lower bounds
	  // on cusp area, when e4 > e3.
		
	Area0 = AreaBound( e2, e3, e3);
	
		for (ee4 = ee3; ee4 < 1 + 19.0/32.0 + widthe3; ee4 = ee4 + 2.0 * widthe3){
		
		e4.value = ee4;
		e4.error = widthe3;
		
		Area = AreaBound( e2, e3, e4);
		
	// In the following, we compare the lowest possible values of Area0 and
	// and Area. If the latter (when e4 is at least e3) is bigger than the former (when
	// e4 = e3), we can go to the next value of e4.
	// Note that when e4 = e3, then we do not pass to the next value of e4, 
        // and so in this case, we continue to further tests.

	if (smallDsize(Area0) < smallDsize(Area)){
		continue;
		}
		
      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;
				
//			In the following, we implement Tools 2 and 3.

			if ( circRadiusSquared( h, mtl, t) < e2*e2)
				continue;	
			if ( !passTool3(h, mtl, t, e2, e3)) 
				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, e3 = %6.5f, e4 = %6.5f, mtl = %6.5f\n",
		 ee2,ee3,e4.value,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.


  
}