import java.text.NumberFormat;

import Jama.Matrix;


public class MathMethods {

	// numerical methods
	
	interface RFunc1d { // real function of 1 variable
		double f(double x);
	}

	// numerical integration using the trapezoid method
	public static double trapezoidIntegrate(RFunc1d fn, double xMin, double xMax, double dx) {
		double area = 0;
		double x = xMin;
		double fx = fn.f(x);
		boolean go = true;
		while (go) {
			double x2 = x + dx;
			if (x2 >= xMax) {
				x2 = xMax;
				dx = x2 - x;
				go = false;
			}
			double fx2 = fn.f(x2);
			double fAvg = (fx + fx2) / 2;
			area += dx * fAvg;
			x = x2;
			fx = fx2;
		}
		return area;
	}
	
	// estimated derivative using the symmetric difference quotient formula.  See https://en.wikipedia.org/wiki/Numerical_differentiation.
	public static RFunc1d symDiffQuotient(RFunc1d fn, double h) {
		return x -> (fn.f(x + h) - fn.f(x - h)) / (2 * h); 
	}
	
	// root finding with Newton-Raphson's method See https://en.wikipedia.org/wiki/Newton%27s_method .
	public static double newtonRaphsonRoot(RFunc1d fn, RFunc1d dfndx, double x0, double eps) {
		double x = x0;
		double fx = fn.f(x);
		while (Math.abs(fx) > eps) {
			double slopeAtX = dfndx.f(x);
			x -= fx / slopeAtX;
			fx = fn.f(x);
		}
		return x;
	}
	
	// root finding with Newton-Raphson's method using estimated derivatives
	public static double newtonRaphsonRoot(RFunc1d fn, double dx, double x0, double eps) {
		return newtonRaphsonRoot(fn, symDiffQuotient(fn, dx), x0, eps);
	}
	
	// For linear algebra operations, I recommend the use of JAMA (http://math.nist.gov/javanumerics/jama/)
	// In Eclipse, right-click the project, select Build Path -> Configure Build Path...
	// Select the Libraries tab. Click the Add JARs button. Select the JAMA .jar file that you downloaded.  Click OK.
	// Then uncomment the code below and Control-Shift-o to update imports.
	public static void linearAlgebraDemo() {		
		// This code is from http://math.nist.gov/javanumerics/jama/ :
		double[][] array = {{1.,2.,3.},{4.,5.,6.},{7.,8.,10.}};
		Matrix A = new Matrix(array);
		Matrix b = Matrix.random(3,1);
		Matrix x = A.solve(b);
		Matrix Residual = A.times(x).minus(b);
		double rnorm = Residual.normInf();
		
		// Print Results:
		System.out.println("A:");
		NumberFormat nf = NumberFormat.getNumberInstance();
		A.print(nf, 3);
		nf.setMinimumFractionDigits(4);
		System.out.println("times x:");
		x.print(nf, 8);
		System.out.println("equals b:");
		b.print(nf, 8);
		System.out.println("Residual:");
		Residual.print(nf,  8);
		System.out.println("Residual norm: " + rnorm);		
	}
	
	public static void testAll() {
		linearAlgebraDemo();
		
		RFunc1d sin = x -> Math.sin(x);
		RFunc1d cos = x -> Math.cos(x);
		double xMin = 0;
		double xMax = Math.PI;
		double dx = .001;
		System.out.printf("The integral of sin from %f to %f with stepsize %f is approximately %f.\n",
				xMin, xMax, dx, trapezoidIntegrate(sin, xMin, xMax, dx));
		System.out.printf("The integral of cos from %f to %f with stepsize %f is approximately %f.\n",
				xMin, xMax, dx, trapezoidIntegrate(cos, xMin, xMax, dx));
		RFunc1d dSin = 	symDiffQuotient(sin, dx);
		System.out.printf("The estimated derivative of sin at 0 using the symmetric difference quotient with h=%f is %f.\n",
				dx, dSin.f(0));		
		System.out.printf("The integral of the symmetric difference quotient of sin from %f to %f with stepsize %f is approximately %f.\n",
				xMin, xMax, dx, trapezoidIntegrate(dSin, xMin, xMax, dx));
		RFunc1d f = x -> (x - 2) * (x + 1); // x^2 - x - 2
		RFunc1d fPrime = x -> 2 * x - 1; 
		double x0 = 1;
		double epsilon = 1e-4;
		System.out.printf("The Newton-Raphson root of f(x) = x^2 - x - 2 found\n starting at x=%f with tolerance %f is %f.\n",
				x0, epsilon, newtonRaphsonRoot(f, fPrime, x0, epsilon));
		x0 = 0;
		System.out.printf("The Newton-Raphson root of f(x) = x^2 - x - 2 found\n starting at x=%f with tolerance %f is %f.\n",
				x0, epsilon, newtonRaphsonRoot(f, fPrime, x0, epsilon));
		System.out.println("The same results using derivatives estimated by the symmetric difference quotient:");
		x0 = 1;
		System.out.printf("The Newton-Raphson root of f(x) = x^2 - x - 2 found\n starting at x=%f with tolerance %f is %f.\n",
				x0, epsilon, newtonRaphsonRoot(f, dx, x0, epsilon));
		x0 = 0;
		System.out.printf("The Newton-Raphson root of f(x) = x^2 - x - 2 found\n starting at x=%f with tolerance %f is %f.\n",
				x0, epsilon, newtonRaphsonRoot(f, dx, x0, epsilon));
		double xExtreme = newtonRaphsonRoot(fPrime, dx, x0, epsilon);
		System.out.printf("The Newton-Raphson extremum of f(x) = x^2 - x - 2 found\n starting at x=%f with tolerance %f is f(%f)=%f.\n",
				x0, epsilon, xExtreme, f.f(xExtreme));
		f = x -> (x + 1) * (x - 1) * (x - 3); // f(x) = x^3 - 3x^2 - x + 3
		fPrime = symDiffQuotient(f, dx); // approx f'(x) = 3x^2 - 6x - 1, (-b+-sqrt(b*b-4*a*c))/2*a = (6+-sqrt(36+12))/6 = (6+-4sqrt(3))/6 = 1 +- 2*sqrt(3)/3 = -.1547005384, 2.1547005384 
		for (int initX = -1; initX <= 4; initX++) {
			xExtreme = newtonRaphsonRoot(fPrime, dx, x0, epsilon);
			x0 = initX;
			System.out.printf("The Newton-Raphson extremum of f(x) = x^3 - 3x^2 - x + 3 found\n starting at x=%f with tolerance %f is f(%f)=%f.\n",
					x0, epsilon, xExtreme, f.f(xExtreme));
		}
	}

	public static void main(String[] args) {
		testAll();
	}

}