gaussianquadratures.hpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
 
/*
 Copyright (C) 2005 Klaus Spanderen
 Copyright (C) 2005 Gary Kennedy
 
 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
 <http://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.
*/
 
#ifndef quantlib_gaussian_quadratures_hpp
#define quantlib_gaussian_quadratures_hpp
 
#include <ql/math/array.hpp>
#include <ql/math/integrals/integral.hpp>
#include <ql/math/integrals/gaussianorthogonalpolynomial.hpp>
 
namespace QuantLib {
    class GaussianOrthogonalPolynomial;
 
 
    class GaussianQuadrature {
      public:
        GaussianQuadrature(Size n,
                           const GaussianOrthogonalPolynomial& p);
 
#if defined(__GNUC__) && (__GNUC__ >= 7)
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wnoexcept-type"
#endif
 
        template <class F>
        Real operator()(const F& f) const {
            Real sum = 0.0;
            for (Integer i = Integer(order())-1; i >= 0; --i) {
                sum += w_[i] * f(x_[i]);
            }
            return sum;
        }
 
#if defined(__GNUC__) && (__GNUC__ >= 7)
#pragma GCC diagnostic pop
#endif
 
        Size order() const { return x_.size(); }
        const Array& weights() { return w_; }
        const Array& x()       { return x_; }
        
      protected:
        Array x_, w_;
    };
 
 
 
    class GaussLaguerreIntegration : public GaussianQuadrature {
      public:
        explicit GaussLaguerreIntegration(Size n, Real s = 0.0)
        : GaussianQuadrature(n, GaussLaguerrePolynomial(s)) {}
    };
 
 
    class GaussHermiteIntegration : public GaussianQuadrature {
      public:
        explicit GaussHermiteIntegration(Size n, Real mu = 0.0)
        : GaussianQuadrature(n, GaussHermitePolynomial(mu)) {}
    };
 
 
    class GaussJacobiIntegration : public GaussianQuadrature {
      public:
        GaussJacobiIntegration(Size n, Real alpha, Real beta)
        : GaussianQuadrature(n, GaussJacobiPolynomial(alpha, beta)) {}
    };
 
 
    class GaussHyperbolicIntegration : public GaussianQuadrature {
      public:
        explicit GaussHyperbolicIntegration(Size n)
        : GaussianQuadrature(n, GaussHyperbolicPolynomial()) {}
    };
 
 
    class GaussLegendreIntegration : public GaussianQuadrature {
      public:
        explicit GaussLegendreIntegration(Size n)
        : GaussianQuadrature(n, GaussJacobiPolynomial(0.0, 0.0)) {}
    };
 
 
    class GaussChebyshevIntegration : public GaussianQuadrature {
      public:
        explicit GaussChebyshevIntegration(Size n)
        : GaussianQuadrature(n, GaussJacobiPolynomial(-0.5, -0.5)) {}
    };
 
 
    class GaussChebyshev2ndIntegration : public GaussianQuadrature {
      public:
        explicit GaussChebyshev2ndIntegration(Size n)
      : GaussianQuadrature(n, GaussJacobiPolynomial(0.5, 0.5)) {}
    };
 
 
    class GaussGegenbauerIntegration : public GaussianQuadrature {
      public:
        GaussGegenbauerIntegration(Size n, Real lambda)
        : GaussianQuadrature(n, GaussJacobiPolynomial(lambda-0.5, lambda-0.5))
        {}
    };
 
 
    namespace detail {
        template <class Integration>
        class GaussianQuadratureIntegrator: public Integrator {
          public:
            explicit GaussianQuadratureIntegrator(Size n);
 
            ext::shared_ptr<Integration> getIntegration() const { return integration_; }
 
          private:
            Real integrate(const ext::function<Real (Real)>& f,
                                           Real a,
                                           Real b) const override;
 
            const ext::shared_ptr<Integration> integration_;
        };
    }
 
    typedef detail::GaussianQuadratureIntegrator<GaussLegendreIntegration>
        GaussLegendreIntegrator;
 
    typedef detail::GaussianQuadratureIntegrator<GaussChebyshevIntegration>
        GaussChebyshevIntegrator;
 
    typedef detail::GaussianQuadratureIntegrator<GaussChebyshev2ndIntegration>
        GaussChebyshev2ndIntegrator;
 
    class TabulatedGaussLegendre {
      public:
        explicit TabulatedGaussLegendre(Size n = 20) { order(n); }
        template <class F>
        Real operator() (const F& f) const {
            QL_ASSERT(w_ != nullptr, "Null weights");
            QL_ASSERT(x_ != nullptr, "Null abscissas");
            Size startIdx;
            Real val;
 
            const Size isOrderOdd = order_ & 1;
 
            if (isOrderOdd) {
              QL_ASSERT((n_>0), "assume at least 1 point in quadrature");
              val = w_[0]*f(x_[0]);
              startIdx=1;
            } else {
              val = 0.0;
              startIdx=0;
            }
 
            for (Size i=startIdx; i<n_; ++i) {
                val += w_[i]*f( x_[i]);
                val += w_[i]*f(-x_[i]);
            }
            return val;
        }
 
        void order(Size);
        Size order() const { return order_; }
 
      private:
        Size order_;
 
        const Real* w_;
        const Real* x_;
        Size  n_;
 
        static const Real w6[3];
        static const Real x6[3];
        static const Size n6;
 
        static const Real w7[4];
        static const Real x7[4];
        static const Size n7;
 
        static const Real w12[6];
        static const Real x12[6];
        static const Size n12;
 
        static const Real w20[10];
        static const Real x20[10];
        static const Size n20;
    };
 
}
 
#endif