/*
 Copyright (C) 2008 Richard Gomes

 This source code is release under the BSD License.

 This file is part of JQuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://jquantlib.org/

 JQuantLib is free software: you can redistribute it and/or modify it
 under the terms of the JQuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 . The license is also available online at
 .

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.

 JQuantLib is based on QuantLib. http://quantlib.org/
 When applicable, the original copyright notice follows this notice.
 */

package org.jquantlib.math.distributions;

import org.jquantlib.QL;
import org.jquantlib.Settings;
import org.jquantlib.lang.exceptions.LibraryException;
import org.jquantlib.math.Closeness;
import org.jquantlib.math.Constants;
import org.jquantlib.math.randomnumbers.InverseCumulative;


/**
 * The inverse normal cumulative distribution is a non-linear function for which
 * no closed-form solution exists. The function is continuous, monotonically increasing,
 * infinitely differentiable, and maps the open interval (0,1) to the whole real line.
 *
 * @see 
 * 		An algorithm for computing the inverse normal cumulative distribution function
 *
 * @author Dominik Holenstein
 */

public class InverseCumulativeNormal implements InverseCumulative {

    static final private String SIGMA_MUST_BE_POSITIVE = "sigma must be greater than 0.0";

    //
    // static final private fields
    //

    static final private double a1 = -3.969683028665376e+01;
    static final private double a2 =  2.209460984245205e+02;
    static final private double a3 = -2.759285104469687e+02;
    static final private double a4 =  1.383577518672690e+02;
    static final private double a5 = -3.066479806614716e+01;
    static final private double a6 =  2.506628277459239e+00;

    static final private double b1 = -5.447609879822406e+01;
    static final private double b2 =  1.615858368580409e+02;
    static final private double b3 = -1.556989798598866e+02;
    static final private double b4 =  6.680131188771972e+01;
    static final private double b5 = -1.328068155288572e+01;

    static final private double c1 = -7.784894002430293e-03;
    static final private double c2 = -3.223964580411365e-01;
    static final private double c3 = -2.400758277161838e+00;
    static final private double c4 = -2.549732539343734e+00;
    static final private double c5 =  4.374664141464968e+00;
    static final private double c6 =  2.938163982698783e+00;

    static final private double d1 =  7.784695709041462e-03;
    static final private double d2 =  3.224671290700398e-01;
    static final private double d3 =  2.445134137142996e+00;
    static final private double d4 =  3.754408661907416e+00;

    //
    // Limits of the approximation regions (break-points)
    //
    static final private double xlow = 0.02425;
    static final private double xhigh = 1.0 - xlow;


    //
    // protected fields
    //

    protected double average;
    protected double sigma;


    //
    // private fields
    //

    // refinement for higher precision
    private final boolean highPrecision;


    //
    // public constructors
    //

    public InverseCumulativeNormal() {
        this(0.0);
    }

    public InverseCumulativeNormal(final double average) {
        this(average, 1.0);
    }

    public InverseCumulativeNormal(final double average, final double sigma) {
        QL.require(sigma > 0.0 , SIGMA_MUST_BE_POSITIVE); // TODO: message
        this.average = average;
        this.sigma = sigma;
        this.highPrecision = new Settings().isRefineHighPrecisionUsingHalleysMethod();
    }


    //
    // implements UnaryFunction
    //

    /**
     * Computes the inverse cumulative normal distribution.
     * @param x
     * @returns average + z * sigma
     */
    @Override
    public double op(double x)/* @ReadOnly */{
        QL.require(sigma > 0.0 , SIGMA_MUST_BE_POSITIVE); // TODO: message

        double z;
        double r;

        if (x < 0.0 || x > 1.0) {
            // try to recover if due to numerical error
            if (Closeness.isCloseEnough(x, 1.0)) {
                x = 1.0;
            } else if (Math.abs(x) < Constants.QL_EPSILON) {
                x = 0.0;
            } else {
                throw new LibraryException(SIGMA_MUST_BE_POSITIVE); // QA:[RG]::verified
            }
        }

        if (x < xlow) {
            // Rational approximation for the lower region 0