package jsci.maths;

import jsci.maths.analysis.RealFunction;
import jsci.maths.analysis.RealFunction2D;

/**
* The numerical math library.
* This class cannot be subclassed or instantiated because all methods are static.
* @version 0.9
* @author Mark Hale
*/
public final class NumericalMath extends AbstractMath {
        private NumericalMath() {}

        /**
        * Calculates the roots of the quadratic equation
        * ax2+bx+c=0.
        * @return an array containing the two roots.
        */
        public static double[] solveQuadratic(final double a,final double b,final double c) {
                final double roots[]=new double[2];
                double q;
                if(b < 0.0)
                        q = (-b+Math.sqrt(b*b-4.0*a*c))/2.0;
                else
                        q = (-b-Math.sqrt(b*b-4.0*a*c))/2.0;
                roots[0] = q/a;
                roots[1] = c/q;
                return roots;
        }
        /**
        * Finds a root using the bisection method.
        * @param a lower bound.
        * @param b upper bound.
        */
        public static double bisection(Mapping func, double a, double b, final int maxIter, final double tol) throws MaximumIterationsExceededException {
                final int signa = ExtraMath.sign(func.map(a));
                final int signb = ExtraMath.sign(func.map(b));
                if(signa == signb)
                        throw new IllegalArgumentException("Bounds do not bracket a root.");
                double x;
                int n = 0;
                do {
                        x = (a + b)/2.0;
                        int signx = ExtraMath.sign(func.map(x));
                        if(signx == signa) {
                                a = x;
                        } else if(signx == signb) {
                                b = x;
                        } else {
                                a = x;
                                b = x;
                        }
                        if(++n > maxIter)
                                throw new MaximumIterationsExceededException("No convergence after "+maxIter+" iterations.");
                } while(Math.abs(a-b) > tol);
                return x;
        }
        /**
        * Finds a root using the false position method.
        * @param a lower bound.
        * @param b upper bound.
        */
        public static double falsePosition(Mapping func, double a, double b, final int maxIter, final double tol) throws MaximumIterationsExceededException {
                double fa = func.map(a);
                double fb = func.map(b);
                final int signa = ExtraMath.sign(fa);
                final int signb = ExtraMath.sign(fb);
                if(signa == signb)
                        throw new IllegalArgumentException("Bounds do not bracket a root.");
                double x;
                double delta;
                int n = 0;
                do {
                        x = a - (b-a)*fa/(fb-fa);
                        double fx = func.map(x);
                        int signx = ExtraMath.sign(fx);
                        if(signx == signa) {
                                delta = x-a;
                                a = x;
                                fa = fx;
                        } else if(signx == signb) {
                                delta = b-x;
                                b = x;
                                fb = fx;
                        } else {
                                delta = 0.0;
                        }
                        if(++n > maxIter)
                                throw new MaximumIterationsExceededException("No convergence after "+maxIter+" iterations.");
                } while(Math.abs(delta) > tol);
                return x;
        }
        /**
        * Finds a root using the Newton-Raphson method.
        * @param x initial guess.
        */
        public static double newtonRaphson(RealFunction func, double x, final int maxIter, final double tol) throws MaximumIterationsExceededException {
                RealFunction deriv = func.differentiate();
                double delta;
                int n = 0;
                do {
                        delta = -func.map(x)/deriv.map(x);
                        x += delta;
                        if(++n > maxIter)
                                throw new MaximumIterationsExceededException("No convergence after "+maxIter+" iterations.");
                } while(Math.abs(delta) > tol);
                return x;
        }
        /**
        * Uses the Euler method to solve an ODE.
        * @param y an array to be filled with y values, set y[0] to initial condition.
        * @param func dy/dt as a function of y.
        * @param dt step size.
        * @return y.
        */
        public static double[] euler(final double y[],final Mapping func,final double dt) {
                for(int i=0;iij=dfi/dxj.
        */
        public static double[][] differentiate(final MappingND func,final double x[],final double dx[]) {
                final double xplus[]=new double[x.length];
                final double xminus[]=new double[x.length];
                System.arraycopy(x,0,xplus,0,x.length);
                System.arraycopy(x,0,xminus,0,x.length);
                xplus[0]+=dx[0];
                xminus[0]-=dx[0];
                double funcdiff[]=ArrayMath.scalarMultiply(0.5/dx[0],ArrayMath.subtract(func.map(xplus),func.map(xminus)));
                final double diff[][]=new double[funcdiff.length][x.length];
                for(int i=0;i