*
 * Programmed by Jack W. Crenshaw, Crenshaw Associates
 * CIS #72325,1327
 *
 * This file defines certain fundamental math functions.
 * These functions replace functions max, min, abs, and
 * sqrt as defined in math.h.
 *
 * They are more general, and are defined via templates
 * so can be used without name changes for any numeric
 * type.
 *
 */


#include <iostream.h>
#include <math.h>

 #ifndef JMATH
 #define JMATH

 #define FALSE 0
 #define TRUE (!FALSE)


 // Return the maximum of two numbers

 template <class T>
 T max(T x, T y)
 {
		 return (x > y)? x: y;
 };


// Return the minimum of two numbers

 template <class T>
 T min(T x, T y)
 {
		 return (x < y)? x: y;
 };


// Return the absolute value of a number

 template <class T>
 T abs(T x)
 {
		 return (x > (T)0)? x: -x;
 };


// Return a number with the magnitude of the first
// argument, and the sign of the second

 template <class T>
 T signum(T x, T y)
 {
		 return (y > (T)0)? abs(x): -abs(x);
 };


// Return a number modulo the second argument

template <class T>
T mod(T x, T y){
		 if(y == (T)0) return x;
		 int i = (int)(x/y);
		 if(x*y < (T)0)--i;
		 x = x-((T) i)*y;
		 if(x==y)x -= y;
		 return x;
}

// return the square root of a number

 double square_root(double x);

 long double square_root(long double x);

#endif
/* jmath.cpp
 *
 * Programmed by Jack W. Crenshaw, Crenshaw Associates
 * CIS #72325,1327
 *
 * This file, with jmath.h, defines certain fundamental
 * math functions. These functions replace functions max,
 * fmax, min, fmin, abs, fabs, and sqrt as defined in math.h.
 *
 * They are more general, and are defined via templates
 * so can be used without name changes for any numeric
 * type.
 *
 */
#include <math.h>

// Return the square root of a number


float square_root(float x){
		 if(x < 0.0) return 0.0;
		 else return sqrt(x);
}


double square_root(double x){
		 if(x < 0.0) return 0.0;
		 else return sqrt(x);
}


long double square_root(long double x){
		 if(x < 0.0L) return 0.0L;
		 else return sqrtl(x);
}
/* jutils.h
 *
 * Copyright (c) 2004 Crenshaw Technologies
 *
 * This file contains some simple functions that I find useful
 * in writing and testing software.
 */

void pause(void);		 	 // use with Borland compiler to hold screen in view
void nl(void);		 		 // newline (function name follows Small C convention)
void error(char *s);    		 // show error and continue
void errorStop(char *s);		 // show error and terminate

/* constant.cpp
 *
 * Copyright (c) 2004 Crenshaw Technologies
 *
 * This file initiates a number of useful constants. It calculates the
 * values of the constants by use of executable assignments at run time.
 */

#include <cmath>

using namespace std;

// Pi and relatives
double Pi = 4*atan(1.0);				 // 180 degrees
double TwoPi = 2*Pi;		 		  	 // 360 degrees
double PiOver2 = Pi/2;		 		 //  90 degrees
double PiOver4 = Pi/4;		 		 //  45 degrees
double PiOver6 = Pi/6;		 		 //  30 degrees

// conversion factors
double RPD = Pi/180;
double DPR = 180/Pi;

// Square roots
double Root2 = sqrt(2.0);
double Root3 = sqrt(3.0);
double Root5 = sqrt(5.0);
double Golden = (Root5 + 1)/2;

// Trig values
double Sin30 = 0.5;
double Cos30 = Root3/2;
double Sin45 = Root2/2;

// Synonyms
#define Cos45 (Sin45)
#define Sin60 (Cos30)
#define Cos60 (Sin30)
#define pi    (Pi)
#define HalfPi (PiOver2)
#define halfpi (PiOver2)
#define Radians_Per_Degree (RPD)
#define Degrees_Per_Radian (DPR)


/* constant.cpp
 *
 * Copyright (c) 2004 Crenshaw Technologies
 *
 * This file initiates a number of useful constants. It calculates the
 * values of the constants by use of executable assignments at run time.
 */

// Pi and relatives
extern double Pi;		 		 	// 180 degrees
extern double TwoPi;		 		 	// 360 degrees
extern double PiOver2;		 		//  90 degrees
extern double PiOver4;		 		//  45 degrees
extern double PiOver6;		 		//  30 degrees

// conversion factors
extern double RPD;
extern double DPR;

// Square roots
extern double Root2;
extern double Root3;
extern double Root5;
extern double Golden;

// Trig values
extern double Sin30;
extern double Cos30;
extern double Sin45;

// Synonyms
#define Cos45 (Sin45)
#define Sin60 (Cos30)
#define Cos60 (Sin30)
#define pi    (Pi)
#define HalfPi (PiOver2)
#define halfpi (PiOver2)
#define Radians_Per_Degree (RPD)
#define Degrees_Per_Radian (DPR)




 /* jutils.cpp
 *
 * Copyright (c) 2004 Crenshaw Technologies
 *
 * This file contains some simple functions that I find useful
 * in writing and testing software.
 */

#include <iostream>
#include <conio>

using namespace std;

typedef unsigned int uint;

// let the customer see the screen <!> (use with Borland compilers only)
void pause(void){
		 cout << "Press any key to continue ...";
		 while(!kbhit());
}

// make a new line (function name follows Small C)
void nl(void){
		 cout << endl;
}

// show error and continue
void error(char *s){
		 cout << s << endl;
		 pause();
}

// show error and terminate
void errorStop(char *s){
		 error(s);
		 exit(1);
}//---------------------------------------------------------------------------

#ifndef MullerTestH
#define MullerTestH
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------


#pragma hdrstop
#include <iostream>
#include <cmath>

#include "MullerTest.h"
#include "constant.h"
#include "jutils.h"
#include "jmath.h"


using namespace std;

// The TWBRF test function
double f(double x){
		 return(4*sin(x)+exp(-x/6));
//		 return x-1;
}

// This function evaluates a new point, sets the y range,
// and tests for convergence

double get_y(double x, double (*f)(double), double eps,
		 		 		 		 		 double &ymax, double &ymin, bool &converged){
		 double y = f(x);
		 ymax = max(ymax, y);
		 ymin = min(ymin, y);
		 converged = (abs(y) < eps*(ymax-ymin));
		 return y;
}



// The iterator
// Based upon the old IBM subroutine rtmi.for
// It uses alternating bisection and inverse
// parabolic interpolation
// restructured by Jack Crenshaw 03/2004
double iterate(double (*f)(double), double x1, double x3, double eps, int imax){

		 double x2, y2;		 		 		 		 		 		 		 // the middle point in bisections
		 double xm, ym;           		 		 		 // the estimated root in interpolations
		 double xeps = eps*abs(x3-x1);		 // convergence criterion in x
		 double ymin = 1.0e-300;		 		 		 		 // set up for criterion in y
		 double ymax = -ymin;
		 bool converged = false;

		 // evaluate function at limits
		 double y1 = get_y(x1, f, eps, ymax, ymin, converged);
		 if(y1 == 0) return x1;		 		 		 // only exact zero accepted here
		 double y3 = get_y(x3, f, eps, ymax, ymin, converged);
		 if(y3 == 0) return x3;		 		 		 // and here

		 // initial points must straddle root
		 if(y1 * y3 > 0)
		 		 errorStop("Iterate: Initial bounds do not straddle root");

		 // begin main bisection loop
		 for(int i=0; i<imax; ++i){
		 		 x2 = (x1+x3)/2;		 		 		 		   // bisect x range
		 		 if(abs(x3 - x1) < eps){
		 		 		 // If we're bisecting, we stop when the interval gets small
		 		 		 return x2;
		 		 }
		 		 y2 = get_y(x2, f, eps, ymax, ymin, converged);

		 		 // perhaps we nailed the root?
		 		 if(converged){
		 		 		 return x2;
		 		 }

		 		 // relabel points to keep root between x1 and x2
		 		 if(y1*y2>0){
		 		 		 swap(x1, x3);
		 		 		 swap(y1, y3);
		 		 }

		 		 // here's where we try parabolic interpolation.
		 		 // there are two criteria for accepting the point.
		 		 // if any one of them fails, just bisect again.

		 		 // note that y21 and y31 cannot be zero here
		 		 // but y32 might be.
		 		 double y21=y2-y1;
		 		 double y32=y3-y2;
		 		 if(y32 != 0){
		 		 		 double y31=y3-y1;
		 		 		 double b=(x2-x1)/y21;
		 		 		 double c=(y21-y32)/(y32*y31);

		 		 		 // Test for the RTMI condition
		 		 		 if(y3*y31 >= 2.0*y2*y21){
		 		 		 		 xm=x1-b*y1*(1.0-c*y2);
		 		 		 		 ym = get_y(xm, f, eps, ymax, ymin, converged);

		 		 		 		 // perhaps we nailed the root?
		 		 		 		 if(converged){
		 		 		 		 		  return xm;
		 		 		 		 }
		 		 		 		 // Relabel the points 9as needed
		 		 		 		 // to keep root between x1 and x2
		 		 		 		 if(ym*y1<0.0){
		 		 		 		 		 x2=xm;
		 		 		 		 		 y2=ym;
		 		 		 		 }
		 		 		 		 else{
		 		 		 		 		 x1=xm;
		 		 		 		 		 y1=ym;
		 		 		 		 }
		 		 		 }
		 		 }
		 		 // we didn't do parabolic interpolation.
		 		 // just relabel points and continue.
		 		 x3 = x2;
		 		 y3 = y2;
		 }
// end of bisection loop. If we get here, we did not converge.
		 error("Iterate: Failed to converge in imax tries.");
		 return x2;
}

//---------------------------------------------------------------------------

#pragma package(smart_init)

int main(int argc, char *argv[]){
		 cout << iterate(f, 0, 6, 1.0e-8, 40) << endl;
		 pause();
		 return 0;
}