/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*

 Copyright (C) 2006 François du Vignaud

 This file is part of QuantLib, a free-software/open-source library

 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it

 under the terms of the QuantLib license.  You should have received a

 copy of the license along with this program; if not, please email

 <quantlib-dev@lists.sf.net>. The license is also available online at

 <https://www.quantlib.org/license.shtml>.

 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.

*/

#include <ql/math/integrals/kronrodintegral.hpp>

#include <ql/types.hpp>

namespace QuantLib {

    static Real rescaleError(Real err,

                             const Real resultAbs,

                             const Real resultAsc) {

        err = std::fabs(err) ;

        if (resultAsc != 0 && err != 0){

            Real scale = std::pow((200 * err / resultAsc), 1.5) ;

            if (scale < 1)

                err = resultAsc * scale ;

            else

                err = resultAsc ;

            }

        if (resultAbs > QL_MIN_POSITIVE_REAL / (50 * QL_EPSILON )){

            Real min_err = 50 * QL_EPSILON  * resultAbs ;

            if (min_err > err)

                err = min_err ;

            }

        return err ;

    }

    /* Gauss-Kronrod-Patterson quadrature coefficients for use in

    quadpack routine qng. These coefficients were calculated with

    101 decimal digit arithmetic by L. W. Fullerton, Bell Labs, Nov

    1981. */

    /* x1, abscissae common to the 10-, 21-, 43- and 87-point rule */

    static const Real x1[5] = {

        0.973906528517171720077964012084452,

        0.865063366688984510732096688423493,

        0.679409568299024406234327365114874,

        0.433395394129247190799265943165784,

        0.148874338981631210884826001129720

        } ;

    /* w10, weights of the 10-point formula */

    static const Real w10[5] = {

        0.066671344308688137593568809893332,

        0.149451349150580593145776339657697,

        0.219086362515982043995534934228163,

        0.269266719309996355091226921569469,

        0.295524224714752870173892994651338

        } ;

    /* x2, abscissae common to the 21-, 43- and 87-point rule */

    static const Real x2[5] = {

        0.995657163025808080735527280689003,

        0.930157491355708226001207180059508,

        0.780817726586416897063717578345042,

        0.562757134668604683339000099272694,

        0.294392862701460198131126603103866

        } ;

    /* w21a, weights of the 21-point formula for abscissae x1 */

    static const Real w21a[5] = {

        0.032558162307964727478818972459390,

        0.075039674810919952767043140916190,

        0.109387158802297641899210590325805,

        0.134709217311473325928054001771707,

        0.147739104901338491374841515972068

        } ;

    /* w21b, weights of the 21-point formula for abscissae x2 */

    static const Real w21b[6] = {

        0.011694638867371874278064396062192,

        0.054755896574351996031381300244580,

        0.093125454583697605535065465083366,

        0.123491976262065851077958109831074,

        0.142775938577060080797094273138717,

        0.149445554002916905664936468389821

        } ;

    /* x3, abscissae common to the 43- and 87-point rule */

    static const Real x3[11] = {

        0.999333360901932081394099323919911,

        0.987433402908088869795961478381209,

        0.954807934814266299257919200290473,

        0.900148695748328293625099494069092,

        0.825198314983114150847066732588520,

        0.732148388989304982612354848755461,

        0.622847970537725238641159120344323,

        0.499479574071056499952214885499755,

        0.364901661346580768043989548502644,

        0.222254919776601296498260928066212,

        0.074650617461383322043914435796506

        } ;

    /* w43a, weights of the 43-point formula for abscissae x1, x3 */

    static const Real w43a[10] = {

        0.016296734289666564924281974617663,

        0.037522876120869501461613795898115,

        0.054694902058255442147212685465005,

        0.067355414609478086075553166302174,

        0.073870199632393953432140695251367,

        0.005768556059769796184184327908655,

        0.027371890593248842081276069289151,

        0.046560826910428830743339154433824,

        0.061744995201442564496240336030883,

        0.071387267268693397768559114425516

        } ;

    /* w43b, weights of the 43-point formula for abscissae x3 */

    static const Real w43b[12] = {

        0.001844477640212414100389106552965,

        0.010798689585891651740465406741293,

        0.021895363867795428102523123075149,

        0.032597463975345689443882222526137,

        0.042163137935191811847627924327955,

        0.050741939600184577780189020092084,

        0.058379395542619248375475369330206,

        0.064746404951445885544689259517511,

        0.069566197912356484528633315038405,

        0.072824441471833208150939535192842,

        0.074507751014175118273571813842889,

        0.074722147517403005594425168280423

        } ;

    /* x4, abscissae of the 87-point rule */

    static const Real x4[22] = {

        0.999902977262729234490529830591582,

        0.997989895986678745427496322365960,

        0.992175497860687222808523352251425,

        0.981358163572712773571916941623894,

        0.965057623858384619128284110607926,

        0.943167613133670596816416634507426,

        0.915806414685507209591826430720050,

        0.883221657771316501372117548744163,

        0.845710748462415666605902011504855,

        0.803557658035230982788739474980964,

        0.757005730685495558328942793432020,

        0.706273209787321819824094274740840,

        0.651589466501177922534422205016736,

        0.593223374057961088875273770349144,

        0.531493605970831932285268948562671,

        0.466763623042022844871966781659270,

        0.399424847859218804732101665817923,

        0.329874877106188288265053371824597,

        0.258503559202161551802280975429025,

        0.185695396568346652015917141167606,

        0.111842213179907468172398359241362,

        0.037352123394619870814998165437704

        } ;

    /* w87a, weights of the 87-point formula for abscissae x1, x2, x3 */

    static const Real w87a[21] = {

        0.008148377384149172900002878448190,

        0.018761438201562822243935059003794,

        0.027347451050052286161582829741283,

        0.033677707311637930046581056957588,

        0.036935099820427907614589586742499,

        0.002884872430211530501334156248695,

        0.013685946022712701888950035273128,

        0.023280413502888311123409291030404,

        0.030872497611713358675466394126442,

        0.035693633639418770719351355457044,

        0.000915283345202241360843392549948,

        0.005399280219300471367738743391053,

        0.010947679601118931134327826856808,

        0.016298731696787335262665703223280,

        0.021081568889203835112433060188190,

        0.025370969769253827243467999831710,

        0.029189697756475752501446154084920,

        0.032373202467202789685788194889595,

        0.034783098950365142750781997949596,

        0.036412220731351787562801163687577,

        0.037253875503047708539592001191226

        } ;

    /* w87b, weights of the 87-point formula for abscissae x4    */

    static const Real w87b[23] = {

        0.000274145563762072350016527092881,

        0.001807124155057942948341311753254,

        0.004096869282759164864458070683480,

        0.006758290051847378699816577897424,

        0.009549957672201646536053581325377,

        0.012329447652244853694626639963780,

        0.015010447346388952376697286041943,

        0.017548967986243191099665352925900,

        0.019938037786440888202278192730714,

        0.022194935961012286796332102959499,

        0.024339147126000805470360647041454,

        0.026374505414839207241503786552615,

        0.028286910788771200659968002987960,

        0.030052581128092695322521110347341,

        0.031646751371439929404586051078883,

        0.033050413419978503290785944862689,

        0.034255099704226061787082821046821,

        0.035262412660156681033782717998428,

        0.036076989622888701185500318003895,

        0.036698604498456094498018047441094,

        0.037120549269832576114119958413599,

        0.037334228751935040321235449094698,

        0.037361073762679023410321241766599

        } ;

     void GaussKronrodNonAdaptive::setRelativeAccuracy(Real relativeAccuracy) { 

        relativeAccuracy_ = relativeAccuracy;

    }

    Real GaussKronrodNonAdaptive::relativeAccuracy() const {

        return relativeAccuracy_;

    }

    GaussKronrodNonAdaptive::GaussKronrodNonAdaptive(Real absoluteAccuracy,

                                                     Size maxEvaluations,

                                                     Real relativeAccuracy)

    : Integrator(absoluteAccuracy, maxEvaluations),

      relativeAccuracy_(relativeAccuracy) {}

    Real

    GaussKronrodNonAdaptive::integrate(const std::function<Real (Real)>& f,

                                       Real a,

                                       Real b) const {

        Real result;

        //Size neval;

        Real fv1[5], fv2[5], fv3[5], fv4[5];

        Real savfun[21];  /* array of function values which have been computed */

        Real res10, res21, res43, res87;    /* 10, 21, 43 and 87 point results */

        Real err;

        Real resAbs; /* approximation to the integral of abs(f) */

        Real resasc; /* approximation to the integral of abs(f-i/(b-a)) */

        int k ;

        QL_REQUIRE(a<b, "b must be greater than a)");

        const Real halfLength = 0.5 * (b - a);

        const Real center = 0.5 * (b + a);

        const Real fCenter = f(center);

        // Compute the integral using the 10- and 21-point formula.

        res10 = 0;

        res21 = w21b[5] * fCenter;

        resAbs = w21b[5] * std::fabs(fCenter);

        for (k = 0; k < 5; k++) {

            Real abscissa = halfLength * x1[k];

            Real fval1 = f(center + abscissa);

            Real fval2 = f(center - abscissa);

            Real fval = fval1 + fval2;

            res10 += w10[k] * fval;

            res21 += w21a[k] * fval;

            resAbs += w21a[k] * (std::fabs(fval1) + std::fabs(fval2));

            savfun[k] = fval;

            fv1[k] = fval1;

            fv2[k] = fval2;

        }

        for (k = 0; k < 5; k++) {

            Real abscissa = halfLength * x2[k];

            Real fval1 = f(center + abscissa);

            Real fval2 = f(center - abscissa);

            Real fval = fval1 + fval2;

            res21 += w21b[k] * fval;

            resAbs += w21b[k] * (std::fabs(fval1) + std::fabs(fval2));

            savfun[k + 5] = fval;

            fv3[k] = fval1;

            fv4[k] = fval2;

        }

        result = res21 * halfLength;

        resAbs *= halfLength ;

        Real mean = 0.5 * res21;

        resasc = w21b[5] * std::fabs(fCenter - mean);

        for (k = 0; k < 5; k++)

            resasc += (w21a[k] * (std::fabs(fv1[k] - mean)

                        + std::fabs(fv2[k] - mean))

                        + w21b[k] * (std::fabs(fv3[k] - mean)

                        + std::fabs(fv4[k] - mean)));

        err = rescaleError ((res21 - res10) * halfLength, resAbs, resasc) ;

        resasc *= halfLength ;

        // test for convergence.

        if (err < absoluteAccuracy() || err < relativeAccuracy() * std::fabs(result)){

            setAbsoluteError(err);

            setNumberOfEvaluations(21);

            return result;

        }

        /* compute the integral using the 43-point formula. */

        res43 = w43b[11] * fCenter;

        for (k = 0; k < 10; k++)

            res43 += savfun[k] * w43a[k];

        for (k = 0; k < 11; k++){

            Real abscissa = halfLength * x3[k];

            Real fval = (f(center + abscissa)

                + f(center - abscissa));

            res43 += fval * w43b[k];

            savfun[k + 10] = fval;

            }

        // test for convergence.

        result = res43 * halfLength;

        err = rescaleError ((res43 - res21) * halfLength, resAbs, resasc);

       if (err < absoluteAccuracy() || err < relativeAccuracy() * std::fabs(result)){

            setAbsoluteError(err);

            setNumberOfEvaluations(43);

            return result;

        }

        /* compute the integral using the 87-point formula. */

        res87 = w87b[22] * fCenter;

        for (k = 0; k < 21; k++)

            res87 += savfun[k] * w87a[k];

        for (k = 0; k < 22; k++){

            Real abscissa = halfLength * x4[k];

            res87 += w87b[k] * (f(center + abscissa)

                + f(center - abscissa));

        }

        // test for convergence.

        result = res87 * halfLength ;

        err = rescaleError ((res87 - res43) * halfLength, resAbs, resasc);

        setAbsoluteError(err);

        setNumberOfEvaluations(87);

        return result;

    }

    Real

    GaussKronrodAdaptive::integrate(const std::function<Real (Real)>& f,

                                    Real a,

                                    Real b) const {

        return integrateRecursively(f, a, b, absoluteAccuracy());

    }

    // weights for 7-point Gauss-Legendre integration

    // (only 4 values out of 7 are given as they are symmetric)

    static const Real g7w[] = { 0.417959183673469,

                                0.381830050505119,

                                0.279705391489277,

                                0.129484966168870 };

    // weights for 15-point Gauss-Kronrod integration

    static const Real k15w[] = { 0.209482141084728,

                                 0.204432940075298,

                                 0.190350578064785,

                                 0.169004726639267,

                                 0.140653259715525,

                                 0.104790010322250,

                                 0.063092092629979,

                                 0.022935322010529 };

    // abscissae (evaluation points)

    // for 15-point Gauss-Kronrod integration

    static const Real k15t[] = { 0.000000000000000,

                                 0.207784955007898,

                                 0.405845151377397,

                                 0.586087235467691,

                                 0.741531185599394,

                                 0.864864423359769,

                                 0.949107912342758,

                                 0.991455371120813 };

    Real GaussKronrodAdaptive::integrateRecursively(

                                    const std::function<Real (Real)>& f,

                                    Real a,

                                    Real b,

                                    Real tolerance) const {

            Real halflength = (b - a) / 2;

            Real center = (a + b) / 2;

            Real g7; // will be result of G7 integral

            Real k15; // will be result of K15 integral

            Real t, fsum; // t (abscissa) and f(t)

            Real fc = f(center);

            g7 = fc * g7w[0];

            k15 = fc * k15w[0];

            // calculate g7 and half of k15

            Integer j, j2;

            for (j = 1, j2 = 2; j < 4; j++, j2 += 2) {

                t = halflength * k15t[j2];

                fsum = f(center - t) + f(center + t);

                g7  += fsum * g7w[j];

                k15 += fsum * k15w[j2];

            }

            // calculate other half of k15

            for (j2 = 1; j2 < 8; j2 += 2) {

                t = halflength * k15t[j2];

                fsum = f(center - t) + f(center + t);

                k15 += fsum * k15w[j2];

            }

            // multiply by (a - b) / 2

            g7 = halflength * g7;

            k15 = halflength * k15;

            // 15 more function evaluations have been used

            increaseNumberOfEvaluations(15);

            // error is <= k15 - g7

            // if error is larger than tolerance then split the interval

            // in two and integrate recursively

            if (std::fabs(k15 - g7) < tolerance) {

                return k15;

            } else {

                QL_REQUIRE(numberOfEvaluations()+30 <=

                           maxEvaluations(),

                           "maximum number of function evaluations "

                           "exceeded");

                return integrateRecursively(f, a, center, tolerance/2)

                    + integrateRecursively(f, center, b, tolerance/2);

            }

        }

    GaussKronrodAdaptive::GaussKronrodAdaptive(Real absoluteAccuracy,

                                               Size maxEvaluations)

    : Integrator(absoluteAccuracy, maxEvaluations) {

        QL_REQUIRE(maxEvaluations >= 15,

                   "required maxEvaluations (" << maxEvaluations <<

                   ") not allowed. It must be >= 15");

    }

}