#include <math.h>
#include <fstream.h>
#include <iomanip.h>

// Chapter 3.1, p.93
// Simple bisection
template <class Type>
void bisection1 (ofstream& out, Type f(Type),Type a, Type b, int n_max)
{
	int n;
	Type c, fa, fb, w;
  
	fa = f(a);
	fb = f(b);

	if ((fa*fb) >= 0)
		return;
	else
	{
		out << setw(5) << "n" << setw(15) << "c" << setw(15) << "f(c)" << endl;
		for (n = 0; n <= n_max; n++)
		{
			c = (a + b)/2;
			w = f(c);
			out << setw(5) << n << setw(15) << c << setw(15) << w << endl;
			if ((w*fa) == 0)
				break;
			else
			{
				if ((w*fa) < 0)
				{
					b = c;
					fb = w;
				}
				else
				{
					a = c;
					fa = w;
				}
			}
		}
	}
} 

template <class Type>
int neg(Type x)
{
	return x<0 ? 1:0;
}

// Chapter 3.1, p.94
// Bisection with tests for convergence
template <class Type>
void bisection2 (ofstream& out, Type f(Type),Type a, Type b, int n_max, Type epsilon)
{
	int n;
	Type c, fa, fb, fc, error;
	int sign(Type );

	fa = f(a);
	fb = f(b);
	n = 0;

	if (neg(fa) == neg(fb))
	{
		out << "a = " << a << "   b = " << b;
		out << "   f(a) = " << fa << "   f(b) = " << fb << endl;
		out << "Function has the same signs at a and b" << endl;
		return;
	}
	error = b - a;

	out << setw(5) << "n" << setw(15) << "c" << setw(15) << "f(c)";
	out << setw(15) << "error" << endl;
	for (n = 0; n <= n_max; n++)
	{
		error = error / 2;
		c = a + error;
		fc = f(c);
		out << setw(5) << n << setw(15) << c << setw(15) << fc;
		out << setw(15) << error << endl;
		if (fabs(error) < epsilon)
		{
			out << "Convergence" << endl;
			return;
		}
		if (neg(fa) != neg(fc))
		{
			b = c;
			fb = fc;
		}
		else
		{
			a = c;
			fa = fc;
		}
	}
}

// Chapter 3.1, p.94-95
// Recursive solution for Bisection
template <class Type>
void bisection3 (ofstream& out, Type f(Type),Type a,Type b,Type fa,Type fb,int n_max, Type epsilon, int n)
{
	Type c,fc,error;
	inline int sign(Type );

	error = (b - a) / 2;
	c = a + error;
	fc = f(c);

	out << setw(5) << n << setw(15) << c << setw(15) << fc;
	out << setw(15) << error << endl;
	if (fabs(error) < epsilon)
	{
		out << "Convergence" << endl;
		return;
	}
	n = n + 1;
	if (n <= n_max)
	{
		if (neg(fa) != neg(fc))
			bisection3(out,f,a,c,fa,fc,n_max,epsilon,n);
		else
			bisection3(out,f,c,b,fc,fb,n_max,epsilon,n);
	}
}

// Chapter 3.2, p.104-105
// Newton's method
template <class Type>
void newton(ofstream& out, Type f(Type), Type deriv_f(Type), Type x, int n_max, Type epsilon, Type delta)
{
	int n = 0;
	Type d,fx, fp;

	fx = f(x);
	out << setw(5) << "n" << setw(15) << "c" << setw(15) << "f(c)" << endl;
	out << setw(5) << n << setw(15) << x << setw(15) << fx << endl;

	for (n = 1; n <= n_max; n++)
	{
		fp = deriv_f(x);
		if (fabs(fp) < delta)
		{
			cout << "Near root" << endl;
			return;
		}
		d = fx / fp;
		x = x - d;
		fx = f(x);
		out << setw(5) << n << setw(15) << x << setw(15) << fx << endl;
		if (fabs(d) < epsilon)
		{
			out << "Convergence" << endl;
			return;
		}
	}
}

template <class Type>
void swap(Type& first, Type& second)
{
	Type temp;
	temp = first;
	first = second;
	second = temp;
}

// Chapter 3.3, p.124
// Secant Method
template <class Type>
void secant(ofstream& out, Type f(Type), Type a, Type b, int n_max, Type epsilon)
{
	int n;
	Type temp_a, temp_b, fa, fb, d;

	fa = f(a);
	fb = f(b);

	temp_a = a;
	temp_b = b;
	if ( fabs(fa) > fabs(fb))
	{
		swap (temp_a,temp_b);
		swap (fa,fb);
	}

	out << setw(5) << "n" << setw(15) << "x" << setw(15) << "f(x)" << endl;
	out << setw(5) << "0" << setw(15) << temp_a << setw(15) << fa << endl;
	out << setw(5) << "1" << setw(15) << temp_b << setw(15) << fb << endl;

	for (n = 2 ; n <= n_max; n++)
	{
		if ( fabs(fa) > fabs(fb))
		{
			swap(temp_a,temp_b);
			swap(fa,fb);
		}

		d = (temp_b - temp_a) / (fb - fa);
		temp_b = temp_a;
		fb = fa;
		d = d * fa;
		if (fabs(d) < epsilon)
		{
			out << "Convergence" << endl;
			return;
		}
		temp_a = temp_a - d;
		fa = f(temp_a);
		out << setw(5) << n << setw(15) << temp_a << setw(15) << fa << endl;
	}
}