MathUtils.hh

Go to the documentation of this file.

#ifndef NPSTAT_MATHUTILS_HH_
#define NPSTAT_MATHUTILS_HH_
 
/*!
// \file MathUtils.hh
//
// \brief Various simple mathematical utilities which did not end up inside
//        dedicated headers
//
// Author: I. Volobouev
//
// March 2010
*/
 
#include <utility>
#include <vector>
#include <cassert>
 
namespace npstat {
    /**
    // Solve the quadratic equation x*x + b*x + c == 0
    // in a numerically sound manner. Return the number of roots.
    */
    unsigned solveQuadratic(double b, double c,
                            double *x1, double *x2);
 
    /**
    // Find the real roots of the cubic:
    //   x**3 + p*x**2 + q*x + r = 0
    // The number of real roots is returned, and the roots are placed
    // into the "v3" array. Original code by Don Herbison-Evans (see
    // his article "Solving Quartics and Cubics for Graphics" in the
    // book "Graphics Gems V", page 3), with minimal adaptation for
    // this package by igv.
    */
    unsigned solveCubic(double p, double q, double r, double v3[3]);
 
    /**
    // Find an extremum of a parabola passing through the three given points.
    // The returned value is "true" for minimum and "false" for maximum.
    // std::invalid_argument will be thrown in case some of the x values
    // coincide or if the coordinates describe a straight line.
    */
    bool parabolicExtremum(double x0, double y0, double x1, double y1,
                           double x2, double y2, double* extremumCoordinate,
                           double* extremumValue);
 
    /**
    // Volume of the n-dimensional sphere of unit radius. Should be
    // multuplied by R^n to get the volume of the sphere with arbitrary radius.
    */
    double ndUnitSphereVolume(unsigned n);
 
    /**
    // Area of the sphere of unit radius embedded in the n-dimensional space.
    // Should be multuplied by R^(n-1) to get the area of the sphere with
    // arbitrary radius.
    */
    double ndUnitSphereArea(unsigned n);
 
    /**
    // Sum of polynomial series. The length of the
    // array of coefficients should be at least degree+1.
    // The highest degree coefficient is assumed to be
    // the last one in the "coeffs" array (0th degree
    // coefficient comes first).
    */
    template<typename Numeric>
    long double polySeriesSum(const Numeric *coeffs,
                              unsigned degree, long double x);
 
    /**
    // Sum and derivative of polynomial series. The length of the
    // array of coefficients should be at least degree+1.
    */
    template<typename Numeric>
    void polyAndDeriv(const Numeric *coeffs, unsigned degree,
                      long double x, long double *value, long double *deriv);
 
    /**
    // Coefficients for the integral of polynomial series.
    // The integration constant is set to 0. The length of the
    // array of coefficients should be at least degree+1,
    // and the length of the buffer for the integral coefficients
    // should be at least degree+2.
    */
    template<typename Numeric>
    void polyIntegralCoeffs(const Numeric* coeffs, unsigned degree,
                            Numeric* integralCoeffs);
 
    /**
    // Series using Legendre polynomials. 0th degree coefficient
    // comes first. Although any value of x can be specified, the result
    // is not going to be terribly meaningful in case |x| > 1.
    */
    template<typename Numeric>
    long double legendreSeriesSum(const Numeric *coeffs,
                                  unsigned degree, long double x);
 
    /**
    // Series for Hermite polynomials orthogonal with weight exp(-x*x/2).
    // These are sometimes called the "probabilists' Hermite polynomials".
    */
    template<typename Numeric>
    long double hermiteSeriesSumProb(const Numeric *coeffs,
                                     unsigned degree, long double x);
 
    /**
    // Series for Hermite polynomials orthogonal with weight exp(-x*x). These
    // are sometimes called the "physicists' Hermite polynomials". The weight
    // exp(-x*x) is also used in the "GaussHermiteQuadrature" class.
    */
    template<typename Numeric>
    long double hermiteSeriesSumPhys(const Numeric *coeffs,
                                     unsigned degree, long double x);
 
    /** Series for Gegenbauer polynomials. "lambda" is the parameter. */
    template<typename Numeric>
    long double gegenbauerSeriesSum(const Numeric *coeffs, unsigned degree,
                                    long double lambda, long double x);
 
    /**
    // Series for the Chebyshev polynomials of the first kind. The array
    // of coefficients must be at least degree+1 long.
    */
    template<typename Numeric>
    long double chebyshevSeriesSum(const Numeric *coeffs, unsigned degree,
                                   long double x);
 
    /**
    // Convert series for the Chebyshev polynomials of the first kind
    // into the series for monomials. The arrays of coefficients must be
    // at least degree+1 long.
    */
    template<typename Numeric1, typename Numeric2>
    void chebyshevMonomialCoeffs(const Numeric1 *coeffs, unsigned degree,
                                 Numeric2* monoCoeffs);
 
    /**
    // Series for the Chebyshev polynomials of the first kind
    // on the given interval [xmin, xmax]
    */
    template<typename Numeric>
    long double chebyshevSeriesSum(const Numeric *coeffs, unsigned degree,
                                   long double xmin, long double xmax,
                                   long double x);
 
    /**
    // Generate Chebyshev series coefficients for a given functor on the
    // interval [xmin, xmax]. The array of coefficients must be at least
    // degree+1 long. The functor will be given a long double argument.
    */
    template<typename Functor, typename Numeric>
    void chebyshevSeriesCoeffs(const Functor& f,
                               long double xmin, long double xmax,
                               unsigned degree, Numeric *coeffs);
 
    /**
    // Converts Chebyshev series coefficients into coefficients
    // for the function integral. The 0's degree coefficient
    // will be chosen in such a way that the integral is 0
    // at xmin. The array of coefficients must be at least
    // degree+1 long, and the buffer for the integral coefficients
    // must be at least degree+2 long.
    */
    template<typename Numeric>
    void chebyshevIntegralCoeffs(const Numeric* coeffs, unsigned degree,
                                 long double xmin, long double xmax,
                                 Numeric* integralCoeffs);
 
    /**
    // Converts Chebyshev series coefficients into coefficients
    // for the function derivative. The array of coefficients
    // must be at least degree+1 long, and the buffer for the
    // derivative coefficients must be at least degree long.
    */
    template<typename Numeric>
    void chebyshevDerivativeCoeffs(const Numeric* coeffs, unsigned degree,
                                   long double xmin, long double xmax,
                                   Numeric* derivativeCoeffs);
 
#ifdef SWIG
    inline double polySeriesSum2(const double *c, const unsigned degreep1,
                                 const double x)
    {
        assert(degreep1);
        return polySeriesSum(c, degreep1-1U, x);
    }
 
    inline std::pair<double,double> polyAndDeriv2(
        const double *c, const unsigned degreep1,
        const double x)
    {
        assert(degreep1);
        long double p, d;
        polyAndDeriv(c, degreep1-1U, x, &p, &d);
        return std::pair<double,double>(p, d);
    }
 
    inline std::vector<double> polyIntegralCoeffs2(const double *c, const unsigned degreep1)
    {
        assert(degreep1);
        std::vector<double> tmp(degreep1+1U);
        polyIntegralCoeffs(c, degreep1-1U, &tmp[0]);
        return tmp;
    }
 
    inline double legendreSeriesSum2(const double *c, const unsigned degreep1,
                                     const double x)
    {
        assert(degreep1);
        return legendreSeriesSum(c, degreep1-1U, x);
    }
 
    inline double hermiteSeriesSumProb2(const double *c, const unsigned degreep1,
                                        const double x)
    {
        assert(degreep1);
        return hermiteSeriesSumProb(c, degreep1-1U, x);
    }
 
    inline double hermiteSeriesSumPhys2(const double *c, const unsigned degreep1,
                                        const double x)
    {
        assert(degreep1);
        return hermiteSeriesSumPhys(c, degreep1-1U, x);
    }
 
    inline double gegenbauerSeriesSum2(const double *c, const unsigned degreep1,
                                       const double lambda, const double x)
    {
        assert(degreep1);
        return gegenbauerSeriesSum(c, degreep1-1U, lambda, x);
    }
 
    inline double chebyshevSeriesSum2(const double *c, const unsigned degreep1,
                                      const double xmin, const double xmax,
                                      const double x)
    {
        assert(degreep1);
        return chebyshevSeriesSum(c, degreep1-1U, xmin, xmax, x);
    }
 
    inline std::vector<double> chebyshevIntegralCoeffs2(const double *c, const unsigned degreep1,
                                                        const double xmin, const double xmax)
    {
        assert(degreep1);
        std::vector<double> tmp(degreep1+1U);
        chebyshevIntegralCoeffs(c, degreep1-1U, xmin, xmax, &tmp[0]);
        return tmp;
    }
 
    inline std::vector<double> chebyshevDerivativeCoeffs2(const double *c, const unsigned degreep1,
                                                          const double xmin, const double xmax)
    {
        assert(degreep1);
        std::vector<double> tmp(degreep1-1U);
        chebyshevDerivativeCoeffs(c, degreep1-1U, xmin, xmax,
                                  tmp.empty() ? (double *)0 : &tmp[0]);
        return tmp;
    }
 
    inline std::vector<double> chebyshevMonomialCoeffs2(const double *c, const unsigned degreep1)
    {
        assert(degreep1);
        std::vector<double> tmp(degreep1);
        chebyshevMonomialCoeffs(c, degreep1-1U, &tmp[0]);
        return tmp;
    }
 
    template<typename Functor>
    inline std::vector<double> chebyshevSeriesCoeffs2(const Functor& f,
        const double xmin, const double xmax, const unsigned degree)
    {
        std::vector<double> tmp(degree+1U);
        chebyshevSeriesCoeffs(f, xmin, xmax, degree, &tmp[0]);
        return tmp;
    }
 
    inline std::vector<double> solveQuadratic2(const double b, const double c)
    {
        double x[2];
        const unsigned nRoots = solveQuadratic(b, c, &x[0], &x[1]);
        return std::vector<double>(x, x+nRoots);
    }
 
    inline std::vector<double> solveCubic2(const double p, const double q,
                                           const double r)
    {
        double v3[3];
        const unsigned nRoots = solveCubic(p, q, r, v3);
        return std::vector<double>(v3, v3+nRoots);
    }
#endif // SWIG
}
 
#include "npstat/nm/MathUtils.icc"
 
#endif // NPSTAT_MATHUTILS_HH_