source: src/main/java/agents/anac/y2015/agentBuyogV2/flanagan/roots/RealRoot.java

/*
*   Class RealRoot
*
*   Contains methods for finding a real root
*
*   The function whose root is to be determined is supplied
*   by means of an interface, RealRootFunction,
*   if no derivative required
*
*   The function whose root is to be determined is supplied
*   by means of an interface, RealRootDerivFunction,
*   as is the first derivative if a derivative is required
*
*   WRITTEN BY: Dr Michael Thomas Flanagan
*
*   DATE:   18 May 2003
*   UPDATE: May 2003 - March 2008,  23-24 September 2008,  30 January 2010
*
*   DOCUMENTATION:
*   See Michael Thomas Flanagan's Java library on-line web page:
*   http://www.ee.ucl.ac.uk/~mflanaga/java/RealRoot.html
*   http://www.ee.ucl.ac.uk/~mflanaga/java/
*
*   Copyright (c) 2003 - 2010    Michael Thomas Flanagan
*
* Permission to use, copy and modify this software and its documentation for NON-COMMERCIAL purposes is granted, without fee,
* provided that an acknowledgement to the author, Dr Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies
* and associated documentation or publications.
*
* Redistributions of the source code of this source code, or parts of the source codes, must retain the above copyright notice, this list of conditions
* and the following disclaimer and requires written permission from the Michael Thomas Flanagan:
*
* Redistribution in binary form of all or parts of this class must reproduce the above copyright notice, this list of conditions and
* the following disclaimer in the documentation and/or other materials provided with the distribution and requires written permission from the Michael Thomas Flanagan:
*
* Dr Michael Thomas Flanagan makes no representations about the suitability or fitness of the software for any or for a particular purpose.
* Dr Michael Thomas Flanagan shall not be liable for any damages suffered as a result of using, modifying or distributing this software
* or its derivatives.
*
***************************************************************************************/

package agents.anac.y2015.agentBuyogV2.flanagan.roots;

import java.util.*;

import agents.anac.y2015.agentBuyogV2.flanagan.complex.Complex;
import agents.anac.y2015.agentBuyogV2.flanagan.complex.ComplexPoly;
import agents.anac.y2015.agentBuyogV2.flanagan.math.Fmath;


// RealRoot class
public class RealRoot{

    // INSTANCE VARIABLES

    private double root = Double.NaN;                           // root to be found
    private double tol = 1e-9;                                  // tolerance in determining convergence upon a root
    private int iterMax = 3000;                                 // maximum number of iterations allowed in root search
    private int iterN = 0;                                      // number of iterations taken in root search
    private double upperBound = 0;                              // upper bound for bisection and false position methods
    private double lowerBound = 0;                              // lower bound for bisection and false position methods
    private double estimate = 0;                                // estimate for Newton-Raphson method
    private int maximumBoundsExtension = 100;                   // number of times that the bounds may be extended
                                                                // by the difference separating them if the root is
                                                                // found not to be bounded
    private boolean noBoundExtensions = false;                  // = true if number of no extension to the  bounds allowed
    private boolean noLowerBoundExtensions = false;             // = true if number of no extension to the lower bound allowed
    private boolean noUpperBoundExtensions = false;             // = true if number of no extension to the upper bound allowed
    private boolean supressLimitReachedMessage = false;         // if true, supresses printing of the iteration limit reached message
    private boolean returnNaN = false;                          // if true exceptions resulting from a bound being NaN do not halt the prorgam but return NaN
                                                                // required by PsRandom and Stat classes calling RealRoot
    private boolean supressNaNmessage = false;                  // if = true the bound is NaN root returned as NaN message supressed

    // STATC VARIABLE

    private static int staticIterMax = 3000;                    // maximum number of iterations allowed in root search (static methods)

    private static int maximumStaticBoundsExtension = 100;      // number of times that the bounds may be extended
                                                                // by the difference separating them if the root is
                                                                // found not to be bounded (static methods)
    private static boolean noStaticBoundExtensions = false;     // = true if number of no extension to the  bounds allowed (static methods)
    private static boolean noStaticLowerBoundExtensions = false;// = true if number of no extension to the lower bound allowed (static methods)
    private static boolean noStaticUpperBoundExtensions = false;// = true if number of no extension to the upper bound allowed (static methods)
    private static boolean staticReturnNaN = false;             // if true exceptions resulting from a bound being NaN do not halt the prorgam but return NaN
                                                                // required by PsRandom and Stat classes calling RealRoot (static methods)
    private static double realTol = 1e-14;                      // tolerance as imag/real in deciding whether a root is real

    // CONSTRUCTOR
    public RealRoot(){
    }

    // INSTANCE METHODS

    // Set lower bound
    public void setLowerBound(double lower){
        this.lowerBound = lower;
    }

    // Set lower bound
    public void setUpperBound(double upper){
        this.upperBound = upper;
    }

    // Reset exception handling for NaN bound flag to true
    // when flag returnNaN = true exceptions resulting from a bound being NaN do not halt the prorgam but return NaN
    // required by PsRandom and Stat classes calling RealRoot
    public void resetNaNexceptionToTrue(){
        this.returnNaN = true;
    }

    // Reset exception handling for NaN bound flag to false
    // when flag returnNaN = false exceptions resulting from a bound being NaN  halts the prorgam
    // required by PsRandom and Stat classes calling RealRoot
    public void resetNaNexceptionToFalse(){
        this.returnNaN = false;
    }

    // Supress NaN bound message
    // if supressNaNmessage = true the bound is NaN root returned as NaN message supressed
    public void supressNaNmessage(){
        this.supressNaNmessage = true;
    }

    // Allow  NaN bound message
    // if supressNaNmessage = false the bound is NaN root returned as NaN message is written
    public void allowNaNmessage(){
        this.supressNaNmessage = false;
    }

    // Set estimate
    public void setEstimate(double estimate){
        this.estimate = estimate;
    }

    // Reset the default tolerance
    public void setTolerance(double tolerance){
        this.tol=tolerance;
    }

    // Get the default tolerance
    public double getTolerance(){
        return this.tol;
    }

    // Reset the maximum iterations allowed
    public void setIterMax(int imax){
        this.iterMax=imax;
    }

    // Get the maximum iterations allowed
    public int getIterMax(){
        return this.iterMax;
    }

    // Get the number of iterations taken
    public int getIterN(){
        return this.iterN;
    }

    // Reset the maximum number of bounds extensions
    public void setmaximumStaticBoundsExtension(int maximumBoundsExtension){
        this.maximumBoundsExtension=maximumBoundsExtension;
    }

    // Prevent extensions to the supplied bounds
    public void noBoundsExtensions(){
        this.noBoundExtensions = true;
        this.noLowerBoundExtensions = true;
        this.noUpperBoundExtensions = true;
    }

    // Prevent extension to the lower bound
    public void noLowerBoundExtension(){
        this.noLowerBoundExtensions = true;
        if(this.noUpperBoundExtensions)this.noBoundExtensions = true;
    }

    // Prevent extension to the upper bound
    public void noUpperBoundExtension(){
        this.noUpperBoundExtensions = true;
        if(this.noLowerBoundExtensions)this.noBoundExtensions = true;
    }

    // Supresses printing of the iteration limit reached message
    // USE WITH CARE - added only to accomadate a specific application using this class!!!!!
    public void supressLimitReachedMessage(){
        this.supressLimitReachedMessage = true;
    }

    // Combined bisection and Inverse Quadratic Interpolation method
    // bounds already entered
        public double brent(RealRootFunction g){
        return this.brent(g, this.lowerBound, this.upperBound);
    }

    // Combined bisection and Inverse Quadratic Interpolation method
    // bounds supplied as arguments
        public double brent(RealRootFunction g, double lower, double upper){
            this.lowerBound = lower;
            this.upperBound = upper;

            // check upper>lower
            if(upper==lower)throw new IllegalArgumentException("upper cannot equal lower");

        boolean testConv = true;    // convergence test: becomes false on convergence
        this.iterN = 0;
        double temp = 0.0D;

        if(upper<lower){
                temp = upper;
                upper = lower;
                lower = temp;
            }

            // calculate the function value at the estimate of the higher bound to x
            double fu = g.function(upper);
            // calculate the function value at the estimate of the lower bound of x
            double fl = g.function(lower);
            if(Double.isNaN(fl)){
                if(this.returnNaN){
                    if(!this.supressNaNmessage)System.out.println("Realroot: brent: lower bound returned NaN as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("lower bound returned NaN as the function value");
                }
            }
        if(Double.isNaN(fu)){
                if(this.returnNaN){
                    if(!this.supressNaNmessage)System.out.println("Realroot: brent: upper bound returned NaN as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("upper bound returned NaN as the function value");
                }
            }

        // check that the root has been bounded and extend bounds if not and extension allowed
        boolean testBounds = true;
        int numberOfBoundsExtension = 0;
        double initialBoundsDifference = (upper - lower)/2.0D;
        while(testBounds){
            if(fu*fl<=0.0D){
                testBounds=false;
            }
            else{
                if(this.noBoundExtensions){
                    String message = "RealRoot.brent: root not bounded and no extension to bounds allowed\n";
                    message += "NaN returned";
                    if(!this.supressNaNmessage)System.out.println(message);
                    return Double.NaN;
                }
                else{
                    numberOfBoundsExtension++;
                    if(numberOfBoundsExtension>this.maximumBoundsExtension){
                        String message = "RealRoot.brent: root not bounded and maximum number of extension to bounds allowed, " + this.maximumBoundsExtension + ", exceeded\n";
                        message += "NaN returned";
                        if(!this.supressNaNmessage)System.out.println(message);
                        return Double.NaN;
                    }
                    if(!this.noLowerBoundExtensions){
                        lower -= initialBoundsDifference;
                        fl = g.function(lower);
                    }
                    if(!this.noUpperBoundExtensions){
                        upper += initialBoundsDifference;
                        fu = g.function(upper);
                    }
                }
            }
        }

            // check initial values for true root value
            if(fl==0.0D){
                this.root=lower;
                testConv = false;
            }
            if(fu==0.0D){
                this.root=upper;
                testConv = false;
            }

            // Function at mid-point of initial estimates
        double mid=(lower+upper)/2.0D;  // mid point (bisect) or new x estimate (Inverse Quadratic Interpolation)
        double lastMidB = mid;          // last succesful mid point
        double fm = g.function(mid);
        double diff = mid-lower;        // difference between successive estimates of the root
        double fmB = fm;                // last succesful mid value function value
        double lastMid=mid;
        boolean lastMethod = true;      // true; last method = Inverse Quadratic Interpolation, false; last method = bisection method
        boolean nextMethod = true;      // true; next method = Inverse Quadratic Interpolation, false; next method = bisection method

            // search
            double rr=0.0D, ss=0.0D, tt=0.0D, pp=0.0D, qq=0.0D; // interpolation variables
            while(testConv){
                // test for convergence
                if(fm==0.0D || Math.abs(diff)<this.tol){
                    testConv=false;
                    if(fm==0.0D){
                        this.root=lastMid;
                    }
                    else{
                        if(Math.abs(diff)<this.tol)this.root=mid;
                    }
                }
                else{
                    lastMethod=nextMethod;
                    // test for succesfull inverse quadratic interpolation
                    if(lastMethod){
                        if(mid<lower || mid>upper){
                            // inverse quadratic interpolation failed
                            nextMethod=false;
                        }
                        else{
                            fmB=fm;
                            lastMidB=mid;
                        }
                    }
                    else{
                        nextMethod=true;
                    }
                        if(nextMethod){
                            // inverse quadratic interpolation
                            fl=g.function(lower);
                        fm=g.function(mid);
                        fu=g.function(upper);
                        rr=fm/fu;
                        ss=fm/fl;
                        tt=fl/fu;
                        pp=ss*(tt*(rr-tt)*(upper-mid)-(1.0D-rr)*(mid-lower));
                        qq=(tt-1.0D)*(rr-1.0D)*(ss-1.0D);
                        lastMid=mid;
                        diff=pp/qq;
                        mid=mid+diff;
                    }
                    else{
                        // Bisection procedure
                        fm=fmB;
                        mid=lastMidB;
                        if(fm*fl>0.0D){
                            lower=mid;
                            fl=fm;
                        }
                        else{
                            upper=mid;
                            fu=fm;
                        }
                        lastMid=mid;
                        mid=(lower+upper)/2.0D;
                        fm=g.function(mid);
                        diff=mid-lastMid;
                        fmB=fm;
                        lastMidB=mid;
                    }
                }
            this.iterN++;
            if(this.iterN>this.iterMax){
                if(!this.supressLimitReachedMessage){
                    if(!this.supressNaNmessage)System.out.println("Class: RealRoot; method: brent; maximum number of iterations exceeded - root at this point, " + Fmath.truncate(mid, 4) + ", returned");
                    if(!this.supressNaNmessage)System.out.println("Last mid-point difference = " + Fmath.truncate(diff, 4) +  ", tolerance = " + this.tol);
                }
                this.root = mid;
                testConv = false;
            }
        }
        return this.root;
    }

    // bisection method
    // bounds already entered
        public double bisect(RealRootFunction g){
            return this.bisect(g, this.lowerBound, this.upperBound);
        }

    // bisection method
        public double bisect(RealRootFunction g, double lower, double upper){
            this.lowerBound = lower;
            this.upperBound = upper;

            // check upper>lower
            if(upper==lower)throw new IllegalArgumentException("upper cannot equal lower");
            if(upper<lower){
            double temp = upper;
                upper = lower;
                lower = temp;
            }

        boolean testConv = true;    // convergence test: becomes false on convergence
        this.iterN = 0;             // number of iterations
        double diff = 1e300;        // abs(difference between the last two successive mid-pint x values)

            // calculate the function value at the estimate of the higher bound to x
            double fu = g.function(upper);
            // calculate the function value at the estimate of the lower bound of x
            double fl = g.function(lower);
            if(Double.isNaN(fl)){
                if(this.returnNaN){
                    if(!this.supressNaNmessage)System.out.println("RealRoot: bisect: lower bound returned NaN as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("lower bound returned NaN as the function value");
                }
            }
        if(Double.isNaN(fu)){
                if(this.returnNaN){
                    if(!this.supressNaNmessage)System.out.println("RealRoot: bisect: upper bound returned NaN as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("upper bound returned NaN as the function value");
                }
            }
        // check that the root has been bounded and extend bounds if not and extension allowed
        boolean testBounds = true;
        int numberOfBoundsExtension = 0;
        double initialBoundsDifference = (upper - lower)/2.0D;
        while(testBounds){
            if(fu*fl<=0.0D){
                testBounds=false;
            }
            else{
                if(this.noBoundExtensions){
                    String message = "RealRoot.bisect: root not bounded and no extension to bounds allowed\n";
                    message += "NaN returned";
                    if(!this.supressNaNmessage)System.out.println(message);
                    return Double.NaN;

                }
                else{
                    numberOfBoundsExtension++;
                    if(numberOfBoundsExtension>this.maximumBoundsExtension){
                        String message = "RealRoot.bisect: root not bounded and maximum number of extension to bounds allowed, " + this.maximumBoundsExtension + ", exceeded\n";
                        message += "NaN returned";
                        if(!this.supressNaNmessage)System.out.println(message);
                        return Double.NaN;
                    }
                    if(!this.noLowerBoundExtensions){
                        lower -= initialBoundsDifference;
                        fl = g.function(lower);
                    }
                    if(!this.noUpperBoundExtensions){
                        upper += initialBoundsDifference;
                        fu = g.function(upper);
                    }
                }
            }
        }

            // check initial values for true root value
            if(fl==0.0D){
                this.root=lower;
                testConv = false;
            }
            if(fu==0.0D){
                this.root=upper;
                testConv = false;
            }

            // start search
        double mid = (lower+upper)/2.0D;    // mid-point
        double lastMid = 1e300;             // previous mid-point
        double fm = g.function(mid);
        while(testConv){
            if(fm==0.0D || diff<this.tol){
                testConv=false;
                this.root=mid;
            }
            if(fm*fl>0.0D){
                lower = mid;
                fl=fm;
            }
            else{
                upper = mid;
                fu=fm;
            }
            lastMid = mid;
            mid = (lower+upper)/2.0D;
            fm = g.function(mid);
            diff = Math.abs(mid-lastMid);
            this.iterN++;
            if(this.iterN>this.iterMax){
                if(!this.supressLimitReachedMessage){
                    if(!this.supressNaNmessage)System.out.println("Class: RealRoot; method: bisect; maximum number of iterations exceeded - root at this point, " + Fmath.truncate(mid, 4) +  ", returned");
                    if(!this.supressNaNmessage)System.out.println("Last mid-point difference = " + Fmath.truncate(diff, 4) +  ", tolerance = " + this.tol);
                }
                this.root = mid;
                testConv = false;
            }
        }
        return this.root;
    }

    // false position  method
    // bounds already entered
        public double falsePosition(RealRootFunction g){
            return this.falsePosition(g, this.lowerBound, this.upperBound);
    }

    // false position  method
        public double falsePosition(RealRootFunction g, double lower, double upper){
            this.lowerBound = lower;
            this.upperBound = upper;

            // check upper>lower
            if(upper==lower)throw new IllegalArgumentException("upper cannot equal lower");
            if(upper<lower){
                double temp = upper;
                upper = lower;
                lower = temp;
            }

        boolean testConv = true;    // convergence test: becomes false on convergence
        this.iterN = 0;             // number of iterations
        double diff = 1e300;        // abs(difference between the last two successive mid-pint x values)

            // calculate the function value at the estimate of the higher bound to x
            double fu = g.function(upper);
            // calculate the function value at the estimate of the lower bound of x
            double fl = g.function(lower);
            if(Double.isNaN(fl)){
                if(this.returnNaN){
                    if(!this.supressNaNmessage)System.out.println("RealRoot: fals: ePositionlower bound returned NaN as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("lower bound returned NaN as the function value");
                }
            }
        if(Double.isNaN(fu)){
                if(this.returnNaN){
                    if(!this.supressNaNmessage)System.out.println("RealRoot: falsePosition: upper bound returned NaN as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("upper bound returned NaN as the function value");
                }
            }

        // check that the root has been bounded and extend bounds if not and extension allowed
        boolean testBounds = true;
        int numberOfBoundsExtension = 0;
        double initialBoundsDifference = (upper - lower)/2.0D;
        while(testBounds){
            if(fu*fl<=0.0D){
                testBounds=false;
            }
            else{
                if(this.noBoundExtensions){
                    String message = "RealRoot.falsePosition: root not bounded and no extension to bounds allowed\n";
                    message += "NaN returned";
                    if(!this.supressNaNmessage)System.out.println(message);
                    return Double.NaN;
                }
                else{
                    numberOfBoundsExtension++;
                    if(numberOfBoundsExtension>this.maximumBoundsExtension){
                        String message = "RealRoot.falsePosition: root not bounded and maximum number of extension to bounds allowed, " + this.maximumBoundsExtension + ", exceeded\n";
                        message += "NaN returned";
                        if(!this.supressNaNmessage)System.out.println(message);
                        return Double.NaN;
                    }
                    if(!this.noLowerBoundExtensions){
                        lower -= initialBoundsDifference;
                        fl = g.function(lower);
                    }
                    if(!this.noUpperBoundExtensions){
                        upper += initialBoundsDifference;
                        fu = g.function(upper);
                    }
                }
            }
        }

            // check initial values for true root value
            if(fl==0.0D){
                this.root=lower;
                testConv = false;
            }
            if(fu==0.0D){
                this.root=upper;
                testConv = false;
            }

            // start search
        double mid = lower+(upper-lower)*Math.abs(fl)/(Math.abs(fl)+Math.abs(fu));    // mid-point
        double lastMid = 1e300;             // previous mid-point
        double fm = g.function(mid);
        while(testConv){
            if(fm==0.0D || diff<this.tol){
                testConv=false;
                this.root=mid;
            }
            if(fm*fl>0.0D){
                lower = mid;
                fl=fm;
            }
            else{
                upper = mid;
                fu=fm;
            }
            lastMid = mid;
            mid = lower+(upper-lower)*Math.abs(fl)/(Math.abs(fl)+Math.abs(fu));    // mid-point
            fm = g.function(mid);
            diff = Math.abs(mid-lastMid);
            this.iterN++;
            if(this.iterN>this.iterMax){
                if(!this.supressLimitReachedMessage){
                    if(!this.supressNaNmessage)System.out.println("Class: RealRoot; method: falsePostion; maximum number of iterations exceeded - root at this point, " + Fmath.truncate(mid, 4) +  ", returned");
                    if(!this.supressNaNmessage)System.out.println("Last mid-point difference = " + Fmath.truncate(diff, 4) +  ", tolerance = " + this.tol);
                }
                this.root = mid;
                testConv = false;
            }
        }
        return this.root;
    }

    // Combined bisection and Newton Raphson method
    // bounds already entered
        public double bisectNewtonRaphson(RealRootDerivFunction g){
            return this.bisectNewtonRaphson(g, this.lowerBound, this.upperBound);
        }

    // Combined bisection and Newton Raphson method
    // default accuracy used
        public double bisectNewtonRaphson(RealRootDerivFunction g, double lower, double upper){
            this.lowerBound = lower;
            this.upperBound = upper;

            // check upper>lower
            if(upper==lower)throw new IllegalArgumentException("upper cannot equal lower");

        boolean testConv = true;    // convergence test: becomes false on convergence
        this.iterN = 0;             // number of iterations
        double temp = 0.0D;

        if(upper<lower){
                temp = upper;
                upper = lower;
                lower = temp;
            }

            // calculate the function value at the estimate of the higher bound to x
            double[] f = g.function(upper);
            double fu=f[0];
            // calculate the function value at the estimate of the lower bound of x
            f = g.function(lower);
            double fl=f[0];
            if(Double.isNaN(fl)){
                if(this.returnNaN){
                    if(!this.supressNaNmessage)System.out.println("RealRoot: bisectNewtonRaphson: lower bound returned NaN as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("lower bound returned NaN as the function value");
                }
            }
        if(Double.isNaN(fu)){
                if(this.returnNaN){
                    if(!this.supressNaNmessage)System.out.println("RealRoot: bisectNewtonRaphson: upper bound returned NaN as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("upper bound returned NaN as the function value");
                }
            }

        // check that the root has been bounded and extend bounds if not and extension allowed
        boolean testBounds = true;
        int numberOfBoundsExtension = 0;
        double initialBoundsDifference = (upper - lower)/2.0D;
        while(testBounds){
            if(fu*fl<=0.0D){
                testBounds=false;
            }
            else{
                if(this.noBoundExtensions){
                    String message = "RealRoot.bisectNewtonRaphson: root not bounded and no extension to bounds allowed\n";
                    message += "NaN returned";
                    if(!this.supressNaNmessage)System.out.println(message);
                    return Double.NaN;
                }
                else{
                    numberOfBoundsExtension++;
                    if(numberOfBoundsExtension>this.maximumBoundsExtension){
                        String message = "RealRoot.bisectNewtonRaphson: root not bounded and maximum number of extension to bounds allowed, " + this.maximumBoundsExtension + ", exceeded\n";
                        message += "NaN returned";
                        if(!this.supressNaNmessage)System.out.println(message);
                        return Double.NaN;
                    }
                    if(!this.noLowerBoundExtensions){
                        lower -= initialBoundsDifference;
                        f = g.function(lower);
                        fl = f[0];
                    }
                    if(!this.noUpperBoundExtensions){
                        upper += initialBoundsDifference;
                        f = g.function(upper);
                        fu = f[0];
                    }
                }
            }
        }

            // check initial values for true root value
            if(fl==0.0D){
                this.root=lower;
                testConv = false;
            }
            if(fu==0.0D){
                this.root=upper;
                testConv = false;
            }

            // Function at mid-point of initial estimates
        double mid=(lower+upper)/2.0D;      // mid point (bisect) or new x estimate (Newton-Raphson)
        double lastMidB = mid;              // last succesful mid point
        f = g.function(mid);
        double diff = f[0]/f[1];            // difference between successive estimates of the root
        double fm = f[0];
        double fmB = fm;                    // last succesful mid value function value
        double lastMid=mid;
        mid = mid-diff;
        boolean lastMethod = true;          // true; last method = Newton Raphson, false; last method = bisection method
        boolean nextMethod = true;          // true; next method = Newton Raphson, false; next method = bisection method

            // search
            while(testConv){
                // test for convergence
                if(fm==0.0D || Math.abs(diff)<this.tol){
                    testConv=false;
                    if(fm==0.0D){
                        this.root=lastMid;
                    }
                    else{
                        if(Math.abs(diff)<this.tol)this.root=mid;
                    }
                }
                else{
                    lastMethod=nextMethod;
                    // test for succesfull Newton-Raphson
                    if(lastMethod){
                        if(mid<lower || mid>upper){
                            // Newton Raphson failed
                            nextMethod=false;
                        }
                        else{
                            fmB=fm;
                            lastMidB=mid;
                        }
                    }
                    else{
                        nextMethod=true;
                    }
                        if(nextMethod){
                            // Newton-Raphson procedure
                        f=g.function(mid);
                        fm=f[0];
                        diff=f[0]/f[1];
                        lastMid=mid;
                        mid=mid-diff;
                    }
                    else{
                        // Bisection procedure
                        fm=fmB;
                        mid=lastMidB;
                        if(fm*fl>0.0D){
                            lower=mid;
                            fl=fm;
                        }
                        else{
                            upper=mid;
                            fu=fm;
                        }
                        lastMid=mid;
                        mid=(lower+upper)/2.0D;
                        f=g.function(mid);
                        fm=f[0];
                        diff=mid-lastMid;
                        fmB=fm;
                        lastMidB=mid;
                    }
                }
            this.iterN++;
            if(this.iterN>this.iterMax){
                if(!this.supressLimitReachedMessage){
                    if(!this.supressNaNmessage)System.out.println("Class: RealRoot; method: bisectNewtonRaphson; maximum number of iterations exceeded - root at this point, " + Fmath.truncate(mid, 4) +  ", returned");
                    if(!this.supressNaNmessage)System.out.println("Last mid-point difference = " + Fmath.truncate(diff, 4) +  ", tolerance = " + this.tol);
                }
                this.root = mid;
                testConv = false;
            }
        }
        return this.root;
    }

    // Newton Raphson method
    // estimate already entered
        public double newtonRaphson(RealRootDerivFunction g){
            return this.newtonRaphson(g, this.estimate);

        }

    // Newton Raphson method
        public double newtonRaphson(RealRootDerivFunction g, double x){
            this.estimate = x;
        boolean testConv = true;    // convergence test: becomes false on convergence
        this.iterN = 0;             // number of iterations
        double diff = 1e300;        // difference between the last two successive mid-pint x values

            // calculate the function and derivative value at the initial estimate  x
            double[] f = g.function(x);
            if(Double.isNaN(f[0])){
                if(this.returnNaN){
                    if(!this.supressNaNmessage)System.out.println("RealRoot: newtonRaphson: NaN returned as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("NaN returned as the function value");
                }
            }
        if(Double.isNaN(f[1])){
                if(this.returnNaN){
                    if(!this.supressNaNmessage)System.out.println("RealRoot: newtonRaphson: NaN returned as the derivative function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("NaN returned as the derivative function value");
                }
            }


            // search
        while(testConv){
            diff = f[0]/f[1];
            if(f[0]==0.0D || Math.abs(diff)<this.tol){
                this.root = x;
                testConv=false;
            }
            else{
                x -= diff;
                f = g.function(x);
                    if(Double.isNaN(f[0]))throw new ArithmeticException("NaN returned as the function value");
                    if(Double.isNaN(f[1]))throw new ArithmeticException("NaN returned as the derivative function value");
                    if(Double.isNaN(f[0])){
                        if(this.returnNaN){
                            if(!this.supressNaNmessage)System.out.println("RealRoot: bisect: NaN as the function value - NaN returned as root");
                            return Double.NaN;
                        }
                        else{
                            throw new ArithmeticException("NaN as the function value");
                        }
                    }
                if(Double.isNaN(f[1])){
                        if(this.returnNaN){
                            if(!this.supressNaNmessage)System.out.println("NaN as the function value - NaN returned as root");
                            return Double.NaN;
                        }
                        else{
                            throw new ArithmeticException("NaN as the function value");
                        }
                    }
            }
            this.iterN++;
            if(this.iterN>this.iterMax){
                if(!this.supressLimitReachedMessage){
                    if(!this.supressNaNmessage)System.out.println("Class: RealRoot; method: newtonRaphson; maximum number of iterations exceeded - root at this point, " + Fmath.truncate(x, 4) + ", returned");
                    if(!this.supressNaNmessage)System.out.println("Last mid-point difference = " + Fmath.truncate(diff, 4) +  ", tolerance = " + this.tol);
                }
                this.root = x;
                testConv = false;
            }
        }
        return this.root;
    }

    // STATIC METHODS

    // Reset the maximum iterations allowed  for static methods
    public void setStaticIterMax(int imax){
        RealRoot.staticIterMax = imax;
    }

    // Get the maximum iterations allowed  for static methods
    public int getStaticIterMax(){
        return RealRoot.staticIterMax;
    }

    // Reset the maximum number of bounds extensions for static methods
    public void setStaticMaximumStaticBoundsExtension(int maximumBoundsExtension){
        RealRoot.maximumStaticBoundsExtension = maximumBoundsExtension;
    }

    // Prevent extensions to the supplied bounds for static methods
    public void noStaticBoundsExtensions(){
        RealRoot.noStaticBoundExtensions = true;
        RealRoot.noStaticLowerBoundExtensions = true;
        RealRoot.noStaticUpperBoundExtensions = true;
    }

    // Prevent extension to the lower bound for static methods
    public void noStaticLowerBoundExtension(){
        RealRoot.noStaticLowerBoundExtensions = true;
        if(RealRoot.noStaticUpperBoundExtensions)RealRoot.noStaticBoundExtensions = true;
    }

    // Prevent extension to the upper bound for static methods
    public void noStaticUpperBoundExtension(){
        RealRoot.noStaticUpperBoundExtensions = true;
        if(RealRoot.noStaticLowerBoundExtensions)RealRoot.noStaticBoundExtensions = true;
    }

    // Reset exception handling for NaN bound flag to true for static methods
    // when flag returnNaN = true exceptions resulting from a bound being NaN do not halt the prorgam but return NaN
    // required by PsRandom and Stat classes calling RealRoot
    public void resetStaticNaNexceptionToTrue(){
        this.staticReturnNaN = true;
    }

    // Reset exception handling for NaN bound flag to false for static methods
    // when flag returnNaN = false exceptions resulting from a bound being NaN  halts the prorgam
    // required by PsRandom and Stat classes calling RealRoot
    public void resetStaticNaNexceptionToFalse(){
        this.staticReturnNaN= false;
    }




    // Combined bisection and Inverse Quadratic Interpolation method
    // bounds supplied as arguments
        public static double brent(RealRootFunction g, double lower, double upper, double tol){
            // check upper>lower
            if(upper==lower)throw new IllegalArgumentException("upper cannot equal lower");

        double root = Double.NaN;
        boolean testConv = true;    // convergence test: becomes false on convergence
        int iterN = 0;
        double temp = 0.0D;

        if(upper<lower){
                temp = upper;
                upper = lower;
                lower = temp;
            }

            // calculate the function value at the estimate of the higher bound to x
            double fu = g.function(upper);
            // calculate the function value at the estimate of the lower bound of x
            double fl = g.function(lower);
            if(Double.isNaN(fl)){
                if(RealRoot.staticReturnNaN){
                    System.out.println("Realroot: brent: lower bound returned NaN as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("lower bound returned NaN as the function value");
                }
            }
        if(Double.isNaN(fu)){
                if(RealRoot.staticReturnNaN){
                    System.out.println("Realroot: brent: upper bound returned NaN as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("upper bound returned NaN as the function value");
                }
            }

        // check that the root has been bounded and extend bounds if not and extension allowed
        boolean testBounds = true;
        int numberOfBoundsExtension = 0;
        double initialBoundsDifference = (upper - lower)/2.0D;
        while(testBounds){
            if(fu*fl<=0.0D){
                testBounds=false;
            }
            else{
                if(RealRoot.noStaticBoundExtensions){
                    String message = "RealRoot.brent: root not bounded and no extension to bounds allowed\n";
                    message += "NaN returned";
                    System.out.println(message);
                    return Double.NaN;
                }
                else{
                    numberOfBoundsExtension++;
                    if(numberOfBoundsExtension>RealRoot.maximumStaticBoundsExtension){
                        String message = "RealRoot.brent: root not bounded and maximum number of extension to bounds allowed, " + RealRoot.maximumStaticBoundsExtension + ", exceeded\n";
                        message += "NaN returned";
                        System.out.println(message);
                        return Double.NaN;
                    }
                    if(!RealRoot.noStaticLowerBoundExtensions){
                        lower -= initialBoundsDifference;
                        fl = g.function(lower);
                    }
                    if(!RealRoot.noStaticUpperBoundExtensions){
                        upper += initialBoundsDifference;
                        fu = g.function(upper);
                    }
                }
            }
        }

            // check initial values for true root value
            if(fl==0.0D){
                root=lower;
                testConv = false;
            }
            if(fu==0.0D){
                root=upper;
                testConv = false;
            }

            // Function at mid-point of initial estimates
        double mid=(lower+upper)/2.0D;      // mid point (bisect) or new x estimate (Inverse Quadratic Interpolation)
        double lastMidB = mid;              // last succesful mid point
        double fm = g.function(mid);
        double diff = mid-lower;            // difference between successive estimates of the root
        double fmB = fm;                    // last succesful mid value function value
        double lastMid=mid;
        boolean lastMethod = true;          // true; last method = Inverse Quadratic Interpolation, false; last method = bisection method
        boolean nextMethod = true;          // true; next method = Inverse Quadratic Interpolation, false; next method = bisection method

            // search
            double rr=0.0D, ss=0.0D, tt=0.0D, pp=0.0D, qq=0.0D; // interpolation variables
            while(testConv){
                // test for convergence
                if(fm==0.0D || Math.abs(diff)<tol){
                    testConv=false;
                    if(fm==0.0D){
                        root=lastMid;
                    }
                    else{
                        if(Math.abs(diff)<tol)root=mid;
                    }
                }
                else{
                    lastMethod=nextMethod;
                    // test for succesfull inverse quadratic interpolation
                    if(lastMethod){
                        if(mid<lower || mid>upper){
                            // inverse quadratic interpolation failed
                            nextMethod=false;
                        }
                        else{
                            fmB=fm;
                            lastMidB=mid;
                        }
                    }
                    else{
                        nextMethod=true;
                    }
                        if(nextMethod){
                            // inverse quadratic interpolation
                            fl=g.function(lower);
                        fm=g.function(mid);
                        fu=g.function(upper);
                        rr=fm/fu;
                        ss=fm/fl;
                        tt=fl/fu;
                        pp=ss*(tt*(rr-tt)*(upper-mid)-(1.0D-rr)*(mid-lower));
                        qq=(tt-1.0D)*(rr-1.0D)*(ss-1.0D);
                        lastMid=mid;
                        diff=pp/qq;
                        mid=mid+diff;
                    }
                    else{
                        // Bisection procedure
                        fm=fmB;
                        mid=lastMidB;
                        if(fm*fl>0.0D){
                            lower=mid;
                            fl=fm;
                        }
                        else{
                            upper=mid;
                            fu=fm;
                        }
                        lastMid=mid;
                        mid=(lower+upper)/2.0D;
                        fm=g.function(mid);
                        diff=mid-lastMid;
                        fmB=fm;
                        lastMidB=mid;
                    }
                }
            iterN++;
            if(iterN>RealRoot.staticIterMax){
                System.out.println("Class: RealRoot; method: brent; maximum number of iterations exceeded - root at this point, " + Fmath.truncate(mid, 4) + ", returned");
                System.out.println("Last mid-point difference = " + Fmath.truncate(diff, 4) +  ", tolerance = " + tol);
                root = mid;
                testConv = false;
            }
        }
        return root;
    }


    // bisection method
    // tolerance supplied
        public static double bisect(RealRootFunction g, double lower, double upper, double tol){

            // check upper>lower
            if(upper==lower)throw new IllegalArgumentException("upper cannot equal lower");
            if(upper<lower){
            double temp = upper;
                upper = lower;
                lower = temp;
            }

            double root = Double.NaN;   // variable to hold the returned root
        boolean testConv = true;    // convergence test: becomes false on convergence
        int iterN = 0;              // number of iterations
        double diff = 1e300;        // abs(difference between the last two successive mid-pint x values)

            // calculate the function value at the estimate of the higher bound to x
            double fu = g.function(upper);
            // calculate the function value at the estimate of the lower bound of x
            double fl = g.function(lower);
            if(Double.isNaN(fl)){
                if(RealRoot.staticReturnNaN){
                    System.out.println("RealRoot: bisect: lower bound returned NaN as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("lower bound returned NaN as the function value");
                }
            }
        if(Double.isNaN(fu)){
                if(RealRoot.staticReturnNaN){
                    System.out.println("RealRoot: bisect: upper bound returned NaN as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("upper bound returned NaN as the function value");
                }
            }
        // check that the root has been bounded and extend bounds if not and extension allowed
        boolean testBounds = true;
        int numberOfBoundsExtension = 0;
        double initialBoundsDifference = (upper - lower)/2.0D;
        while(testBounds){
            if(fu*fl<=0.0D){
                testBounds = false;
            }
            else{
                if(RealRoot.noStaticBoundExtensions){
                    String message = "RealRoot.bisect: root not bounded and no extension to bounds allowed\n";
                    message += "NaN returned";
                    System.out.println(message);
                    return Double.NaN;

                }
                else{
                    numberOfBoundsExtension++;
                    if(numberOfBoundsExtension>RealRoot.maximumStaticBoundsExtension){
                        String message = "RealRoot.bisect: root not bounded and maximum number of extension to bounds allowed, " + RealRoot.maximumStaticBoundsExtension + ", exceeded\n";
                        message += "NaN returned";
                        System.out.println(message);
                        return Double.NaN;
                    }
                    if(!RealRoot.noStaticLowerBoundExtensions){
                        lower -= initialBoundsDifference;
                        fl = g.function(lower);
                    }
                    if(!RealRoot.noStaticUpperBoundExtensions){
                        upper += initialBoundsDifference;
                        fu = g.function(upper);
                    }
                }
            }
        }

            // check initial values for true root value
            if(fl==0.0D){
                root=lower;
                testConv = false;
            }
            if(fu==0.0D){
                root=upper;
                testConv = false;
            }

            // start search
        double mid = (lower+upper)/2.0D;    // mid-point
        double lastMid = 1e300;             // previous mid-point
        double fm = g.function(mid);
        while(testConv){
            if(fm==0.0D || diff<tol){
                testConv=false;
                root=mid;
            }
            if(fm*fl>0.0D){
                lower = mid;
                fl=fm;
            }
            else{
                upper = mid;
                fu=fm;
            }
            lastMid = mid;
            mid = (lower+upper)/2.0D;
            fm = g.function(mid);
            diff = Math.abs(mid-lastMid);
            iterN++;
            if(iterN>RealRoot.staticIterMax){
                System.out.println("Class: RealRoot; method: bisect; maximum number of iterations exceeded - root at this point, " + Fmath.truncate(mid, 4) +  ", returned");
                System.out.println("Last mid-point difference = " + Fmath.truncate(diff, 4) +  ", tolerance = " + tol);
                root = mid;
                testConv = false;
            }
        }
        return root;
    }






    // false position  method
    // tolerance supplied
        public static double falsePosition(RealRootFunction g, double lower, double upper, double tol){

        // check upper>lower
            if(upper==lower)throw new IllegalArgumentException("upper cannot equal lower");
            if(upper<lower){
                double temp = upper;
                upper = lower;
                lower = temp;
            }

            double root = Double.NaN;   // variable to hold the returned root
        boolean testConv = true;    // convergence test: becomes false on convergence
        int iterN = 0;              // number of iterations
        double diff = 1e300;        // abs(difference between the last two successive mid-pint x values)

            // calculate the function value at the estimate of the higher bound to x
            double fu = g.function(upper);
            // calculate the function value at the estimate of the lower bound of x
            double fl = g.function(lower);
            if(Double.isNaN(fl)){
                if(RealRoot.staticReturnNaN){
                    System.out.println("RealRoot: fals: ePositionlower bound returned NaN as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("lower bound returned NaN as the function value");
                }
            }
        if(Double.isNaN(fu)){
                if(RealRoot.staticReturnNaN){
                    System.out.println("RealRoot: falsePosition: upper bound returned NaN as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("upper bound returned NaN as the function value");
                }
            }

        // check that the root has been bounded and extend bounds if not and extension allowed
        boolean testBounds = true;
        int numberOfBoundsExtension = 0;
        double initialBoundsDifference = (upper - lower)/2.0D;
        while(testBounds){
            if(fu*fl<=0.0D){
                testBounds=false;
            }
            else{
                if(RealRoot.noStaticBoundExtensions){
                    String message = "RealRoot.falsePosition: root not bounded and no extension to bounds allowed\n";
                    message += "NaN returned";
                    System.out.println(message);
                    return Double.NaN;
                }
                else{
                    numberOfBoundsExtension++;
                    if(numberOfBoundsExtension>RealRoot.maximumStaticBoundsExtension){
                        String message = "RealRoot.falsePosition: root not bounded and maximum number of extension to bounds allowed, " + RealRoot.maximumStaticBoundsExtension + ", exceeded\n";
                        message += "NaN returned";
                        System.out.println(message);
                        return Double.NaN;
                    }
                    if(!RealRoot.noStaticLowerBoundExtensions){
                        lower -= initialBoundsDifference;
                        fl = g.function(lower);
                    }
                    if(!RealRoot.noStaticUpperBoundExtensions){
                        upper += initialBoundsDifference;
                        fu = g.function(upper);
                    }
                }
            }
        }

            // check initial values for true root value
            if(fl==0.0D){
                root=lower;
                testConv = false;
            }
            if(fu==0.0D){
                root=upper;
                testConv = false;
            }

            // start search
        double mid = lower+(upper-lower)*Math.abs(fl)/(Math.abs(fl)+Math.abs(fu));    // mid-point
        double lastMid = 1e300;             // previous mid-point
        double fm = g.function(mid);
        while(testConv){
            if(fm==0.0D || diff<tol){
                testConv=false;
                root=mid;
            }
            if(fm*fl>0.0D){
                lower = mid;
                fl=fm;
            }
            else{
                upper = mid;
                fu=fm;
            }
            lastMid = mid;
            mid = lower+(upper-lower)*Math.abs(fl)/(Math.abs(fl)+Math.abs(fu));    // mid-point
            fm = g.function(mid);
            diff = Math.abs(mid-lastMid);
            iterN++;
            if(iterN>RealRoot.staticIterMax){
                System.out.println("Class: RealRoot; method: falsePostion; maximum number of iterations exceeded - root at this point, " + Fmath.truncate(mid, 4) +  ", returned");
                System.out.println("Last mid-point difference = " + Fmath.truncate(diff, 4) +  ", tolerance = " + tol);
                root = mid;
                testConv = false;
            }
        }
        return root;
    }






    // Combined bisection and Newton Raphson method
    // tolerance supplied
        public static double bisectNewtonRaphson(RealRootDerivFunction g, double lower, double upper, double tol){

            // check upper>lower
            if(upper==lower)throw new IllegalArgumentException("upper cannot equal lower");

        double root = Double.NaN;
        boolean testConv = true;    // convergence test: becomes false on convergence
        int iterN = 0;              // number of iterations
        double temp = 0.0D;

        if(upper<lower){
                temp = upper;
                upper = lower;
                lower = temp;
            }

            // calculate the function value at the estimate of the higher bound to x
            double[] f = g.function(upper);
            double fu=f[0];
            // calculate the function value at the estimate of the lower bound of x
            f = g.function(lower);
            double fl=f[0];
            if(Double.isNaN(fl)){
                if(RealRoot.staticReturnNaN){
                    System.out.println("RealRoot: bisectNewtonRaphson: lower bound returned NaN as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("lower bound returned NaN as the function value");
                }
            }
        if(Double.isNaN(fu)){
                if(RealRoot.staticReturnNaN){
                    System.out.println("RealRoot: bisectNewtonRaphson: upper bound returned NaN as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("upper bound returned NaN as the function value");
                }
            }

        // check that the root has been bounded and extend bounds if not and extension allowed
        boolean testBounds = true;
        int numberOfBoundsExtension = 0;
        double initialBoundsDifference = (upper - lower)/2.0D;
        while(testBounds){
            if(fu*fl<=0.0D){
                testBounds=false;
            }
            else{
                if(RealRoot.noStaticBoundExtensions){
                    String message = "RealRoot.bisectNewtonRaphson: root not bounded and no extension to bounds allowed\n";
                    message += "NaN returned";
                    System.out.println(message);
                    return Double.NaN;
                }
                else{
                    numberOfBoundsExtension++;
                    if(numberOfBoundsExtension>RealRoot.maximumStaticBoundsExtension){
                        String message = "RealRoot.bisectNewtonRaphson: root not bounded and maximum number of extension to bounds allowed, " + RealRoot.maximumStaticBoundsExtension + ", exceeded\n";
                        message += "NaN returned";
                        System.out.println(message);
                        return Double.NaN;
                    }
                    if(!RealRoot.noStaticLowerBoundExtensions){
                        lower -= initialBoundsDifference;
                        f = g.function(lower);
                        fl = f[0];
                    }
                    if(!RealRoot.noStaticUpperBoundExtensions){
                        upper += initialBoundsDifference;
                        f = g.function(upper);
                        fu = f[0];
                    }
                }
            }
        }

            // check initial values for true root value
            if(fl==0.0D){
                root=lower;
                testConv = false;
            }
            if(fu==0.0D){
                root=upper;
                testConv = false;
            }

            // Function at mid-point of initial estimates
        double mid=(lower+upper)/2.0D;      // mid point (bisect) or new x estimate (Newton-Raphson)
        double lastMidB = mid;              // last succesful mid point
        f = g.function(mid);
        double diff = f[0]/f[1];            // difference between successive estimates of the root
        double fm = f[0];
        double fmB = fm;                    // last succesful mid value function value
        double lastMid=mid;
        mid = mid-diff;
        boolean lastMethod = true;          // true; last method = Newton Raphson, false; last method = bisection method
        boolean nextMethod = true;          // true; next method = Newton Raphson, false; next method = bisection method

            // search
            while(testConv){
                // test for convergence
                if(fm==0.0D || Math.abs(diff)<tol){
                    testConv=false;
                    if(fm==0.0D){
                        root=lastMid;
                    }
                    else{
                        if(Math.abs(diff)<tol)root=mid;
                    }
                }
                else{
                    lastMethod=nextMethod;
                    // test for succesfull Newton-Raphson
                    if(lastMethod){
                        if(mid<lower || mid>upper){
                            // Newton Raphson failed
                            nextMethod=false;
                        }
                        else{
                            fmB=fm;
                            lastMidB=mid;
                        }
                    }
                    else{
                        nextMethod=true;
                    }
                        if(nextMethod){
                            // Newton-Raphson procedure
                        f=g.function(mid);
                        fm=f[0];
                        diff=f[0]/f[1];
                        lastMid=mid;
                        mid=mid-diff;
                    }
                    else{
                        // Bisection procedure
                        fm=fmB;
                        mid=lastMidB;
                        if(fm*fl>0.0D){
                            lower=mid;
                            fl=fm;
                        }
                        else{
                            upper=mid;
                            fu=fm;
                        }
                        lastMid=mid;
                        mid=(lower+upper)/2.0D;
                        f=g.function(mid);
                        fm=f[0];
                        diff=mid-lastMid;
                        fmB=fm;
                        lastMidB=mid;
                    }
                }
            iterN++;
            if(iterN>RealRoot.staticIterMax){
                System.out.println("Class: RealRoot; method: bisectNewtonRaphson; maximum number of iterations exceeded - root at this point, " + Fmath.truncate(mid, 4) +  ", returned");
                System.out.println("Last mid-point difference = " + Fmath.truncate(diff, 4) +  ", tolerance = " + tol);
                root = mid;
                testConv = false;
            }
        }
        return root;
    }





    // Newton Raphson method
        public static double newtonRaphson(RealRootDerivFunction g, double x, double tol){
            double root = Double.NaN;
        boolean testConv = true;    // convergence test: becomes false on convergence
        int iterN = 0;              // number of iterations
        double diff = 1e300;        // difference between the last two successive mid-pint x values

            // calculate the function and derivative value at the initial estimate  x
            double[] f = g.function(x);
            if(Double.isNaN(f[0])){
                if(RealRoot.staticReturnNaN){
                    System.out.println("RealRoot: newtonRaphson: NaN returned as the function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("NaN returned as the function value");
                }
            }
        if(Double.isNaN(f[1])){
                if(RealRoot.staticReturnNaN){
                    System.out.println("RealRoot: newtonRaphson: NaN returned as the derivative function value - NaN returned as root");
                    return Double.NaN;
                }
                else{
                    throw new ArithmeticException("NaN returned as the derivative function value");
                }
            }


            // search
        while(testConv){
            diff = f[0]/f[1];
            if(f[0]==0.0D || Math.abs(diff)<tol){
                root = x;
                testConv=false;
            }
            else{
                x -= diff;
                f = g.function(x);
                    if(Double.isNaN(f[0]))throw new ArithmeticException("NaN returned as the function value");
                    if(Double.isNaN(f[1]))throw new ArithmeticException("NaN returned as the derivative function value");
                    if(Double.isNaN(f[0])){
                        if(RealRoot.staticReturnNaN){
                            System.out.println("RealRoot: NewtonRaphson: NaN as the function value - NaN returned as root");
                            return Double.NaN;
                        }
                        else{
                            throw new ArithmeticException("NaN as the function value");
                        }
                    }
                if(Double.isNaN(f[1])){
                        if(RealRoot.staticReturnNaN){
                            System.out.println("NaN as the function value - NaN returned as root");
                            return Double.NaN;
                        }
                        else{
                            throw new ArithmeticException("NaN as the function value");
                        }
                    }
            }
            iterN++;
            if(iterN>RealRoot.staticIterMax){
                System.out.println("Class: RealRoot; method: newtonRaphson; maximum number of iterations exceeded - root at this point, " + Fmath.truncate(x, 4) + ", returned");
                System.out.println("Last mid-point difference = " + Fmath.truncate(diff, 4) +  ", tolerance = " + tol);
                root = x;
                testConv = false;
            }
        }
        return root;
    }

    // ROOTS OF A QUADRATIC EQUATION
    // c + bx + ax^2 = 0
    // roots returned in root[]
    // 4ac << b*b accomodated by these methods
    // roots returned as Double in an ArrayList if roots are real
    // roots returned as Complex in an ArrayList if any root is complex
    public static ArrayList<Object> quadratic(double c, double b, double a){

            ArrayList<Object> roots = new ArrayList<Object>(2);

            double bsquared = b*b;
            double fourac = 4.0*a*c;
            if(bsquared<fourac){
                Complex[] croots = ComplexPoly.quadratic(c, b, a);
                roots.add("complex");
                roots.add(croots);
            }
            else{
                double[] droots = new double[2];
                double  bsign = Fmath.sign(b);
                double qsqrt = Math.sqrt(bsquared - fourac);
                qsqrt = -0.5*(b + bsign*qsqrt);
                droots[0] = qsqrt/a;
                droots[1] = c/qsqrt;
                roots.add("real");
                roots.add(droots);
            }
            return roots;
    }


    // ROOTS OF A CUBIC EQUATION
    // a + bx + cx^2 + dx^3 = 0
    // roots returned as Double in an ArrayList if roots are real
    // roots returned as Complex in an ArrayList if any root is complex
    public static ArrayList<Object> cubic(double a, double b, double c, double d){

            ArrayList<Object> roots = new ArrayList<Object>(2);

            double aa = c/d;
            double bb = b/d;
            double cc = a/d;
            double bigQ = (aa*aa - 3.0*bb)/9.0;
            double bigQcubed = bigQ*bigQ*bigQ;
            double bigR = (2.0*aa*aa*aa - 9.0*aa*bb + 27.0*cc)/54.0;
            double bigRsquared = bigR*bigR;

            if(bigRsquared>=bigQcubed){
                Complex[] croots = ComplexPoly.cubic(a, b, c, d);
                roots.add("complex");
                roots.add(croots);
            }
            else{
                double[] droots = new double[3];
                double theta = Math.acos(bigR/Math.sqrt(bigQcubed));
                double aover3 = aa/3.0;
                double qterm = -2.0*Math.sqrt(bigQ);

                droots[0] = qterm*Math.cos(theta/3.0) - aover3;
                droots[1] = qterm*Math.cos((theta + 2.0*Math.PI)/3.0) - aover3;
                droots[2] = qterm*Math.cos((theta - 2.0*Math.PI)/3.0) - aover3;
                roots.add("real");
                roots.add(droots);
            }
            return roots;
    }

          // ROOTS OF A POLYNOMIAL
        // For general details of root searching and a discussion of the rounding errors
        // see Numerical Recipes, The Art of Scientific Computing
        // by W H Press, S A Teukolsky, W T Vetterling & B P Flannery
        // Cambridge University Press,   http://www.nr.com/

        // Calculate the roots  of a real polynomial
        // initial root estimate is zero [for deg>3]
        // roots are not olished [for deg>3]
        public static ArrayList<Object> polynomial(double[] coeff){
                boolean polish=true;
                double estx = 0.0;
                return RealRoot.polynomial(coeff, polish, estx);
        }

        // Calculate the roots  of a real polynomial
        // initial root estimate is zero [for deg>3]
        // roots are polished [for deg>3]
        public static ArrayList<Object> polynomial(double[] coeff, boolean polish){
                double estx = 0.0;
                return RealRoot.polynomial (coeff, polish, estx);
        }

        // Calculate the roots  of a real polynomial
        // initial root estimate is estx [for deg>3]
        // roots are not polished [for deg>3]
        public static ArrayList<Object> polynomial(double[] coeff, double estx){
                boolean polish=true;
                return RealRoot.polynomial(coeff, polish, estx);
        }

        // Calculate the roots  of a real polynomial
        // initial root estimate is estx [for deg>3]
        // roots are polished [for deg>3]
        public static ArrayList<Object> polynomial (double[] coeff, boolean polish, double estx){

                int nCoeff = coeff.length;
                if(nCoeff<2)throw new IllegalArgumentException("a minimum of two coefficients is required");
                ArrayList<Object> roots = new ArrayList<Object>(nCoeff);
                boolean realRoots = true;

                // check for zero roots
                int nZeros=0;
                int ii=0;
                boolean testZero=true;
                while(testZero){
                    if(coeff[ii]==0.0){
                        nZeros++;
                        ii++;
                    }
                    else{
                        testZero=false;
                    }
                }

                // Repack coefficients
                int nCoeffWz = nCoeff - nZeros;
                double[] coeffWz = new double[nCoeffWz];
                if(nZeros>0){
                    for(int i=0; i<nCoeffWz; i++)coeffWz[i] = coeff[i+nZeros];
                }
                else{
                    for(int i=0; i<nCoeffWz; i++)coeffWz[i] = coeff[i];
                }

                // Calculate non-zero roots
                ArrayList<Object> temp = new ArrayList<Object>(2);
                double[] cdreal = null;
                switch(nCoeffWz){
                        case 0:
                        case 1: break;
                        case 2: temp.add("real");
                                double[] dtemp = {-coeffWz[0]/coeffWz[1]};
                                temp.add(dtemp);
                                break;
                        case 3: temp = RealRoot.quadratic(coeffWz[0],coeffWz[1],coeffWz[2]);
                                if(((String)temp.get(0)).equals("complex"))realRoots = false;
                                break;
                        case 4: temp = RealRoot.cubic(coeffWz[0],coeffWz[1],coeffWz[2], coeffWz[3]);
                                if(((String)temp.get(0)).equals("complex"))realRoots = false;
                                break;
                        default: ComplexPoly cp = new ComplexPoly(coeffWz);
                                Complex[] croots = cp.roots(polish, new Complex(estx, 0.0));
                                cdreal = new double[nCoeffWz-1];
                                int counter = 0;
                                for(int i=0; i<(nCoeffWz-1); i++){
                                    if(croots[i].getImag()/croots[i].getReal()<RealRoot.realTol){
                                        cdreal[i] = croots[i].getReal();
                                        counter++;
                                    }
                                }
                                if(counter==(nCoeffWz-1)){
                                    temp.add("real");
                                    temp.add(cdreal);
                                }
                                else{
                                    temp.add("complex");
                                    temp.add(croots);
                                    realRoots = false;
                                }
                }

                // Pack roots into returned ArrayList
                if(nZeros==0){
                    roots = temp;
                }
                else{
                    if(realRoots){
                        double[] dtemp1 = new double[nCoeff-1];
                        double[] dtemp2 = (double[])temp.get(1);
                        for(int i=0; i<nCoeffWz-1; i++)dtemp1[i] = dtemp2[i];
                        for(int i=0; i<nZeros; i++)dtemp1[i+nCoeffWz-1] = 0.0;
                        roots.add("real");
                        roots.add(dtemp1);
                    }
                    else{
                        Complex[] dtemp1 = Complex.oneDarray(nCoeff-1);
                        Complex[] dtemp2 = (Complex[])temp.get(1);
                        for(int i=0; i<nCoeffWz-1; i++)dtemp1[i] = dtemp2[i];
                        for(int i=0; i<nZeros; i++)dtemp1[i+nCoeffWz-1] = new Complex(0.0, 0.0);
                        roots.add("complex");
                        roots.add(dtemp1);
                    }
                }

                return roots;
        }

        // Reset the criterion for deciding a that a root, calculated as Complex, is real
        // Default option; imag/real <1e-14
        // this method allows thew value of 1e-14 to be reset
        public void resetRealTest(double ratio){
            RealRoot.realTol = ratio;
        }

}