Page MenuHomeHEPForge
Files F11221715

No OneTemporary
Actions

View Options

Index: trunk/npstat/nm/ClassicalOrthoPolys1D.cc
===================================================================
--- trunk/npstat/nm/ClassicalOrthoPolys1D.cc (revision 936)
+++ trunk/npstat/nm/ClassicalOrthoPolys1D.cc (revision 937)
@@ -1,411 +1,405 @@
#include <cmath>
#include <cassert>
#include <stdexcept>
#include "npstat/nm/SpecialFunctions.hh"
#include "npstat/nm/GaussLegendreQuadrature.hh"
#include "npstat/nm/ClassicalOrthoPolys1D.hh"
#include "npstat/rng/permutation.hh"
// Calculate (-1)^k
static inline int powm1(const unsigned k)
{
return k % 2U ? -1 : 1;
}
namespace {
class PolyProdWithSquaredPoly : public npstat::Functor1<long double, long double>
{
public:
inline PolyProdWithSquaredPoly(const npstat::AbsClassicalOrthoPoly1D& fcn,
const unsigned deg1, const unsigned deg2)
: fcn_(fcn), degree1_(deg1), degree2_(deg2) {}
inline virtual ~PolyProdWithSquaredPoly() {}
inline virtual long double operator()(const long double& a) const
{
std::pair<long double,long double> p = fcn_.twopoly(
degree1_, degree2_, a);
return p.first * p.second * p.second;
}
private:
const npstat::AbsClassicalOrthoPoly1D& fcn_;
unsigned degree1_;
unsigned degree2_;
};
}
namespace npstat {
void LegendreOrthoPoly1D::monomialCoeffs(
const double* coeffs, const unsigned degree,
long double* monoCoeffs) const
{
assert(coeffs);
assert(monoCoeffs);
std::fill(monoCoeffs, monoCoeffs+(degree+1U), 0.0L);
for (unsigned deg=0U; deg<=degree; ++deg)
if (coeffs[deg])
{
const long double norm = coeffs[deg]*sqrtl(2.0L*deg + 1)/powl(2.0L, deg);
const unsigned kmax = deg/2U;
for (unsigned k=0; k<=kmax; ++k)
monoCoeffs[deg-2U*k] += norm*powm1(k)*ldfactorial(2U*(deg-k))/ldfactorial(k)/ldfactorial(deg-k)/ldfactorial(deg-2U*k);
}
}
void LegendreOrthoPoly1D::integralCoeffs(
const double* coeffs, const unsigned degree,
double* integCoeffs) const
{
const double sqr3 = 1.7320508075688773;
assert(coeffs);
assert(integCoeffs);
integCoeffs[0] = coeffs[0];
integCoeffs[1] = coeffs[0]/sqr3;
std::fill(integCoeffs+2U, integCoeffs+(degree+2U), 0.0);
for (unsigned deg=1U; deg<=degree; ++deg)
if (coeffs[deg])
{
const double tmp = coeffs[deg]/sqrt(2.0*deg + 1);
integCoeffs[deg+1U] += tmp/sqrt(2.0*deg + 3);
integCoeffs[deg-1U] -= tmp/sqrt(2.0*deg - 1);
}
}
void LegendreOrthoPoly1D::derivativeCoeffs(
const double* coeffs, const unsigned degree,
double* derivCoeffs) const
{
if (degree)
{
assert(coeffs);
assert(derivCoeffs);
std::fill(derivCoeffs, derivCoeffs+degree, 0.0);
for (unsigned deg=1U; deg<=degree; ++deg)
if (coeffs[deg])
{
const double tmp = coeffs[deg]*sqrt(2.0*deg + 1);
for (unsigned k=0U; k<deg; k+=2U)
derivCoeffs[deg-1U-k] += tmp*sqrt(2.0*(deg-k) - 1);
}
}
else
{
if (derivCoeffs)
derivCoeffs[0] = 0.0;
}
}
void LegendreOrthoPoly1D::squaredPolyCoeffs(
const unsigned degree, long double* expansionCoeffs) const
{
assert(expansionCoeffs);
expansionCoeffs[0] = 2.0;
if (degree)
{
const unsigned nIntegPt = GaussLegendreQuadrature::minimalExactRule(4U*degree);
assert(nIntegPt);
const GaussLegendreQuadrature glq(nIntegPt);
for (unsigned deg=1U; deg<=2U*degree; ++deg)
{
if (deg % 2U)
expansionCoeffs[deg] = 0.0;
else
{
PolyProdWithSquaredPoly prod(*this, deg, degree);
expansionCoeffs[deg] = glq.integrate(prod, -1.0, 1.0);
}
}
}
}
void ShiftedLegendreOrthoPoly1D::integralCoeffs(
const double* coeffs, const unsigned degree,
double* integCoeffs) const
{
const double twosqr3 = 3.4641016151377546;
assert(coeffs);
assert(integCoeffs);
integCoeffs[0] = coeffs[0]/2.0;
integCoeffs[1] = coeffs[0]/twosqr3;
std::fill(integCoeffs+2U, integCoeffs+(degree+2U), 0.0);
for (unsigned deg=1U; deg<=degree; ++deg)
if (coeffs[deg])
{
const double tmp = coeffs[deg]/2.0/sqrt(2.0*deg + 1);
integCoeffs[deg+1U] += tmp/sqrt(2.0*deg + 3);
integCoeffs[deg-1U] -= tmp/sqrt(2.0*deg - 1);
}
}
void ShiftedLegendreOrthoPoly1D::derivativeCoeffs(
const double* coeffs, const unsigned degree,
double* derivCoeffs) const
{
if (degree)
{
assert(coeffs);
assert(derivCoeffs);
std::fill(derivCoeffs, derivCoeffs+degree, 0.0);
for (unsigned deg=1U; deg<=degree; ++deg)
if (coeffs[deg])
{
const double tmp = 2.0*coeffs[deg]*sqrt(2.0*deg + 1);
for (unsigned k=0U; k<deg; k+=2U)
derivCoeffs[deg-1U-k] += tmp*sqrt(2.0*(deg-k) - 1);
}
}
else
{
if (derivCoeffs)
derivCoeffs[0] = 0.0;
}
}
void ShiftedLegendreOrthoPoly1D::squaredPolyCoeffs(
const unsigned degree, long double* expansionCoeffs) const
{
assert(expansionCoeffs);
expansionCoeffs[0] = 1.0;
if (degree)
{
const unsigned nIntegPt = GaussLegendreQuadrature::minimalExactRule(4U*degree);
assert(nIntegPt);
const GaussLegendreQuadrature glq(nIntegPt);
for (unsigned deg=1U; deg<=2U*degree; ++deg)
{
if (deg % 2U)
expansionCoeffs[deg] = 0.0;
else
{
PolyProdWithSquaredPoly prod(*this, deg, degree);
expansionCoeffs[deg] = glq.integrate(prod, 0.0, 1.0);
}
}
}
}
JacobiOrthoPoly1D::JacobiOrthoPoly1D(const double a, const double b)
: alpha_(a),
beta_(b)
{
if (!(alpha_ > -1.0 && beta_ > -1.0)) throw std::invalid_argument(
"In npstat::JacobiOrthoPoly1D constructor: invalid arguments");
norm_ = powl(2.0L, alpha_+beta_+1)*Gamma(alpha_+1)*Gamma(beta_+1)/
Gamma(alpha_+beta_+2);
}
long double JacobiOrthoPoly1D::weight(const long double x) const
{
if (x < -1.0L || x > 1.0L)
return 0.0L;
else
{
if (x == 1.0L && alpha_ < 0.0L) throw std::invalid_argument(
"In npstat::JacobiOrthoPoly1D::weight: invalid x=1 argument");
if (x == -1.0L && beta_ < 0.0L) throw std::invalid_argument(
"In npstat::JacobiOrthoPoly1D::weight: invalid x=-1 argument");
// Don't call "powl" if we can avoid it because this is a slow function
long double p;
if (alpha_ == beta_)
{
const long double tmp = 1.0L - x*x;
// Fourth power is a special case, often used in density estimation
if (alpha_ == 4.0)
{
const long double tmpsq = tmp*tmp;
p = tmpsq*tmpsq;
}
else if (alpha_ == 0.0)
p = 1.0L;
else if (alpha_ == 1.0)
p = tmp;
else if (alpha_ == 2.0)
p = tmp*tmp;
else if (alpha_ == 3.0)
p = tmp*tmp*tmp;
else if (alpha_ == 5.0)
{
const long double tmpsq = tmp*tmp;
p = tmpsq*tmpsq*tmp;
}
else if (alpha_ == 6.0)
{
const long double tmpsq = tmp*tmp;
p = tmpsq*tmpsq*tmpsq;
}
else
p = powl(tmp, alpha_);
}
else
p = powl(1.0L - x, alpha_)*powl(1.0L + x, beta_);
return p/norm_;
}
}
std::pair<long double,long double>
JacobiOrthoPoly1D::monicRecurrenceCoeffs(const unsigned k) const
{
long double a, b = 1.0L;
if (k)
{
const long double tmp = 2*k + alpha_ + beta_;
a = (beta_ - alpha_)*(beta_ + alpha_)/tmp/(tmp + 2);
if (k == 1U)
b = 4*k*(k+alpha_)*(k+beta_)/tmp/tmp/(tmp + 1);
else
b = 4*k*(k+alpha_)*(k+beta_)*(k+alpha_+beta_)/
tmp/tmp/(tmp + 1)/(tmp - 1);
}
else
a = (beta_ - alpha_)/(beta_ + alpha_ + 2);
return std::pair<long double,long double>(a, b);
}
std::pair<long double,long double>
JacobiOrthoPoly1D::recurrenceCoeffs(const unsigned k) const
{
const std::pair<long double,long double>& p = this->monicRecurrenceCoeffs(k);
return std::pair<long double,long double>(p.first, sqrtl(p.second));
}
void JacobiOrthoPoly1D::derivativeCoeffs(
const double* coeffs, const unsigned maxdeg, double* derivCoeffs) const
{
if (maxdeg)
{
assert(coeffs);
assert(derivCoeffs);
assert(derivCoeffs != coeffs);
const double aplusbplusone = alpha_ + beta_+ 1.0;
const double commonf = 0.25/(alpha_+1.0)/(beta_+1.0)*(aplusbplusone+1.0)*(aplusbplusone+2.0);
for (unsigned ideg=maxdeg; ideg>0U; --ideg)
derivCoeffs[ideg-1U] = coeffs[ideg]*sqrt(ideg*(aplusbplusone + ideg)*commonf);
}
}
long double HermiteProbOrthoPoly1D::weight(const long double x) const
{
const long double norm = 0.39894228040143267793994606L; // 1/sqrt(2 Pi)
return norm*expl(-0.5*x*x);
}
std::pair<long double,long double>
HermiteProbOrthoPoly1D::recurrenceCoeffs(const unsigned k) const
{
return std::pair<long double,long double>(0.0L, sqrtl(k));
}
std::pair<long double,long double>
HermiteProbOrthoPoly1D::monicRecurrenceCoeffs(const unsigned k) const
{
return std::pair<long double,long double>(0.0L, k);
}
std::pair<long double,long double>
LaguerreOrthoPoly1D::recurrenceCoeffs(const unsigned k) const
{
const long double kld = k;
return std::pair<long double,long double>(2*kld + 1, k > 0U ? kld : 1.0L);
}
std::pair<long double,long double>
LaguerreOrthoPoly1D::monicRecurrenceCoeffs(const unsigned k) const
{
const long double kld = k;
return std::pair<long double,long double>(2*kld + 1, k > 0U ? kld*kld : 1.0L);
}
GramOrthoPoly1D::GramOrthoPoly1D(const unsigned npt)
: npoints_(npt)
{
if (!npoints_) throw std::invalid_argument(
"In npstat::GramOrthoPoly1D constructor: argument must be positive");
}
long double GramOrthoPoly1D::weight(const long double x) const
{
// What can we do here? The weight is discrete.
return (x >= -1.0L && x <= 1.0L) ? 0.5L : 0.0L;
}
std::pair<long double,long double>
GramOrthoPoly1D::recurrenceCoeffs(const unsigned k) const
{
const std::pair<long double,long double>& p = this->monicRecurrenceCoeffs(k);
return std::pair<long double,long double>(p.first, sqrtl(p.second));
}
std::pair<long double,long double>
GramOrthoPoly1D::monicRecurrenceCoeffs(const unsigned k) const
{
long double b = 0.25L;
if (k)
{
const long double m = npoints_;
const long double n = k;
const long double msq = m*m;
const long double nsq = n*n;
b = 0.25L*nsq/msq*(msq - nsq)/(nsq - 0.25L);
}
return std::pair<long double,long double>(0.0L, b);
}
ShiftedGramOrthoPoly1D::ShiftedGramOrthoPoly1D(const unsigned npt)
: npoints_(npt)
{
if (!npoints_) throw std::invalid_argument(
"In npstat::ShiftedGramOrthoPoly1D constructor: argument must be positive");
}
long double ShiftedGramOrthoPoly1D::weight(const long double x) const
{
// What can we do here? The weight is discrete.
return (x >= 0.0L && x <= 1.0L) ? 1.0L : 0.0L;
}
std::pair<long double,long double>
ShiftedGramOrthoPoly1D::recurrenceCoeffs(const unsigned k) const
{
const std::pair<long double,long double>& p = this->monicRecurrenceCoeffs(k);
return std::pair<long double,long double>(p.first, sqrtl(p.second));
}
std::pair<long double,long double>
ShiftedGramOrthoPoly1D::monicRecurrenceCoeffs(const unsigned k) const
{
- long double b = 0.25L;
- if (k)
- {
- const long double m = npoints_;
- const long double n = k;
- const long double msq = m*m;
- const long double nsq = n*n;
- b = 0.0625L*nsq/msq*(msq - nsq)/(nsq - 0.25L);
- }
+ const long double n = npoints_;
+ const long double ksq = static_cast<long double>(k)*k;
+ const long double b = 0.0625L*ksq*(1.0L - ksq/n/n)/(ksq - 0.25L);
return std::pair<long double,long double>(0.5L, b);
}
}
Index: trunk/npstat/nm/ClassicalOrthoPolys1D.hh
===================================================================
--- trunk/npstat/nm/ClassicalOrthoPolys1D.hh (revision 936)
+++ trunk/npstat/nm/ClassicalOrthoPolys1D.hh (revision 937)
@@ -1,436 +1,436 @@
#ifndef NPSTAT_CLASSICALORTHOPOLYS1D_HH_
#define NPSTAT_CLASSICALORTHOPOLYS1D_HH_
/*!
// \file ClassicalOrthoPolys1D.hh
//
// \brief Orthonormal versions of some classical orthogonal polynomials
//
// Author: I. Volobouev
//
// July 2018
*/
#include <cmath>
#include <vector>
#include <algorithm>
#include "npstat/nm/AbsClassicalOrthoPoly1D.hh"
namespace npstat {
/** Orthonormal Legendre polynomials on [-1, 1] */
class LegendreOrthoPoly1D : public AbsClassicalOrthoPoly1D
{
public:
inline virtual LegendreOrthoPoly1D* clone() const
{return new LegendreOrthoPoly1D(*this);}
inline virtual ~LegendreOrthoPoly1D() {}
inline long double weight(const long double x) const
{return (x >= -1.0L && x <= 1.0L) ? 0.5L : 0.0L;}
inline double xmin() const {return -1.0L;}
inline double xmax() const {return 1.0L;}
/**
// Derive the series coefficients for the integral of the
// argument series from -1 to x. 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.
*/
void integralCoeffs(const double* coeffs, unsigned degree,
double* integCoeffs) const;
/**
// Derive the series coefficients for the derivative of the
// argument series. The array of coefficients must be at least
// degree+1 long, and the buffer for the derivative coefficients
// must be at least degree long.
*/
void derivativeCoeffs(const double* coeffs, unsigned degree,
double* derivCoeffs) const;
/**
// Derive the series coefficients for the monomial
// expansion. The arrays of coefficients must be
// at least degree+1 long.
*/
void monomialCoeffs(const double* coeffs, unsigned degree,
long double* monoCoeffs) const;
/**
// Expand a squared polynomial into polynomial series.
// The number of coefficients will be 2*degree + 1.
*/
void squaredPolyCoeffs(unsigned degree,
long double* expansionCoeffs) const;
inline virtual std::pair<long double,long double>
monicRecurrenceCoeffs(const unsigned k) const
{
- long double beta = 1.0L;
+ long double beta = 0.0L;
if (k)
beta = 1.0L/(4.0L - 1.0L/k/k);
return std::pair<long double,long double>(0.0L, beta);
}
inline virtual std::pair<long double,long double>
recurrenceCoeffs(const unsigned k) const
{
- long double sqb = 1.0L;
+ long double sqb = 0.0L;
if (k)
sqb = 1.0L/sqrtl(4.0L - 1.0L/k/k);
return std::pair<long double,long double>(0.0L, sqb);
}
#ifdef SWIG
public:
inline std::vector<double>
integralCoeffs2(const double *c, const unsigned degreep1) const
{
std::vector<double> tmp(degreep1+1U);
integralCoeffs(c, degreep1-1U, &tmp[0]);
return tmp;
}
inline std::vector<double>
derivativeCoeffs2(const double *c, const unsigned degreep1) const
{
std::vector<double> tmp(std::max(degreep1-1U, 1U));
derivativeCoeffs(c, degreep1-1U, &tmp[0]);
return tmp;
}
inline std::vector<double>
monomialCoeffs2(const double *c, const unsigned degreep1) const
{
std::vector<long double> ldtmp(degreep1+1U);
monomialCoeffs(c, degreep1-1U, &ldtmp[0]);
return std::vector<double>(ldtmp.begin(), ldtmp.end());
}
inline std::vector<double>
squaredPolyCoeffs2(const unsigned degree) const
{
std::vector<long double> tmp(2U*degree + 1U);
squaredPolyCoeffs(degree, &tmp[0]);
return std::vector<double>(tmp.begin(), tmp.end());
}
#endif
};
/** Orthonormal Legendre polynomials on [0, 1] (that is, shifted) */
class ShiftedLegendreOrthoPoly1D : public AbsClassicalOrthoPoly1D
{
public:
inline virtual ShiftedLegendreOrthoPoly1D* clone() const
{return new ShiftedLegendreOrthoPoly1D(*this);}
inline virtual ~ShiftedLegendreOrthoPoly1D() {}
inline long double weight(const long double x) const
{return (x >= 0.0L && x <= 1.0L) ? 1.0L : 0.0L;}
inline double xmin() const {return 0.0L;}
inline double xmax() const {return 1.0L;}
/**
// Derive the series coefficients for the integral of the
// argument series from 0 to x. 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.
*/
void integralCoeffs(const double* coeffs, unsigned degree,
double* integCoeffs) const;
/**
// Derive the series coefficients for the derivative of the
// argument series. The array of coefficients must be at least
// degree+1 long, and the buffer for the derivative coefficients
// must be at least degree long.
*/
void derivativeCoeffs(const double* coeffs, unsigned degree,
double* derivCoeffs) const;
/**
// Expand a squared polynomial into polynomial series.
// The number of coefficients will be 2*degree + 1.
*/
void squaredPolyCoeffs(unsigned degree,
long double* expansionCoeffs) const;
inline virtual std::pair<long double,long double>
monicRecurrenceCoeffs(const unsigned k) const
{
- long double beta = 1.0L;
+ long double beta = 0.0L;
if (k)
beta = 0.25L/(4.0L - 1.0L/k/k);
return std::pair<long double,long double>(0.5L, beta);
}
inline virtual std::pair<long double,long double>
recurrenceCoeffs(const unsigned k) const
{
- long double sqb = 1.0L;
+ long double sqb = 0.0L;
if (k)
sqb = 0.5L/sqrtl(4.0L - 1.0L/k/k);
return std::pair<long double,long double>(0.5L, sqb);
}
#ifdef SWIG
public:
inline std::vector<double>
integralCoeffs2(const double *c, const unsigned degreep1) const
{
std::vector<double> tmp(degreep1+1U);
integralCoeffs(c, degreep1-1U, &tmp[0]);
return tmp;
}
inline std::vector<double>
derivativeCoeffs2(const double *c, const unsigned degreep1) const
{
std::vector<double> tmp(std::max(degreep1-1U, 1U));
derivativeCoeffs(c, degreep1-1U, &tmp[0]);
return tmp;
}
inline std::vector<double>
squaredPolyCoeffs2(const unsigned degree) const
{
std::vector<long double> tmp(2U*degree + 1U);
squaredPolyCoeffs(degree, &tmp[0]);
return std::vector<double>(tmp.begin(), tmp.end());
}
#endif
};
/**
// The weight function for Jacobi polynomials is proportional
// to (1 - x)^alpha (1 + x)^beta
*/
class JacobiOrthoPoly1D : public AbsClassicalOrthoPoly1D
{
public:
JacobiOrthoPoly1D(double alpha, double beta);
inline virtual JacobiOrthoPoly1D* clone() const
{return new JacobiOrthoPoly1D(*this);}
inline virtual ~JacobiOrthoPoly1D() {}
inline double alpha() const {return alpha_;}
inline double beta() const {return beta_;}
long double weight(long double x) const;
inline double xmin() const {return -1.0L;}
inline double xmax() const {return 1.0L;}
/**
// The following function generates coefficients for derivatives
// of the Jacobi polynomial series. The derivative series should be
// evaluated using a JacobiOrthoPoly1D(alpha_+1, beta_+1) object and
// with maxdeg argument decreased by 1. The array "coeffs" should
// have at least maxdeg+1 elements, while the "derivativeCoeffs"
// buffer should have at least maxdeg elements.
*/
void derivativeCoeffs(const double* coeffs,
unsigned maxdeg,
double* derivCoeffs) const;
virtual std::pair<long double,long double>
monicRecurrenceCoeffs(unsigned k) const;
virtual std::pair<long double,long double>
recurrenceCoeffs(unsigned k) const;
private:
double alpha_;
double beta_;
long double norm_;
#ifdef SWIG
public:
inline std::vector<double>
derivativeCoeffs2(const double *c, const unsigned degreep1) const
{
std::vector<double> tmp(degreep1-1U);
derivativeCoeffs(c, degreep1-1U,
tmp.empty() ? (double *)0 : &tmp[0]);
return tmp;
}
#endif
};
/** Orthonormal Chebyshev polynomials of the 1st kind */
class ChebyshevOrthoPoly1st : public AbsClassicalOrthoPoly1D
{
public:
inline virtual ChebyshevOrthoPoly1st* clone() const
{return new ChebyshevOrthoPoly1st(*this);}
inline virtual ~ChebyshevOrthoPoly1st() {}
inline long double weight(const long double x) const
{
const long double Pi = 3.14159265358979323846264338L;
return (x > -1.0L && x < 1.0L) ? 1.0/sqrtl(1.0L - x*x)/Pi : 0.0L;
}
inline double xmin() const {return -1.0L;}
inline double xmax() const {return 1.0L;}
inline virtual std::pair<long double,long double>
monicRecurrenceCoeffs(const unsigned k) const
{
long double beta = 0.25L;
if (k == 0U)
beta = 3.14159265358979323846264338328L;
else if (k == 1U)
beta = 0.5L;
return std::pair<long double,long double>(0.0L, beta);
}
inline virtual std::pair<long double,long double>
recurrenceCoeffs(const unsigned k) const
{
long double sqb = 0.5L;
if (k == 0U)
sqb = 1.772453850905516027298167483L; // sqrt(Pi)
else if (k == 1U)
sqb = 0.7071067811865475244008443621L; // 1/sqrt(2)
return std::pair<long double,long double>(0.0L, sqb);
}
};
/** Orthonormal Chebyshev polynomials of the 2nd kind */
class ChebyshevOrthoPoly2nd : public AbsClassicalOrthoPoly1D
{
public:
inline virtual ChebyshevOrthoPoly2nd* clone() const
{return new ChebyshevOrthoPoly2nd(*this);}
inline virtual ~ChebyshevOrthoPoly2nd() {}
inline long double weight(const long double x) const
{
const long double PiOver2 = 1.57079632679489661923132169L;
return (x > -1.0L && x < 1.0L) ? sqrtl(1.0L - x*x)/PiOver2 : 0.0L;
}
inline double xmin() const {return -1.0L;}
inline double xmax() const {return 1.0L;}
inline virtual std::pair<long double,long double>
monicRecurrenceCoeffs(const unsigned k) const
{
long double beta = 0.25L;
if (!k)
beta = 1.57079632679489661923132169164L; // Pi/2
return std::pair<long double,long double>(0.0L, beta);
}
inline virtual std::pair<long double,long double>
recurrenceCoeffs(const unsigned k) const
{
long double sqb = 0.5L;
if (!k)
sqb = 1.25331413731550025120788264L; // sqrt(Pi/2)
return std::pair<long double,long double>(0.0L, sqb);
}
};
/** Orthonormal(!) probabilist's Hermite polynomials */
class HermiteProbOrthoPoly1D : public AbsClassicalOrthoPoly1D
{
public:
inline virtual HermiteProbOrthoPoly1D* clone() const
{return new HermiteProbOrthoPoly1D(*this);}
inline virtual ~HermiteProbOrthoPoly1D() {}
long double weight(long double x) const;
inline double xmin() const {return -xmax();}
inline double xmax() const {return 150.0;}
virtual std::pair<long double,long double>
monicRecurrenceCoeffs(unsigned k) const;
virtual std::pair<long double,long double>
recurrenceCoeffs(unsigned k) const;
};
/**
// Orthonormal Laguerre polynomials. Note that the definition
// here differs from the "standard" definition by the sign of
// all odd degree polynomials (we want to start with the monic
// representation).
*/
class LaguerreOrthoPoly1D : public AbsClassicalOrthoPoly1D
{
public:
inline virtual LaguerreOrthoPoly1D* clone() const
{return new LaguerreOrthoPoly1D(*this);}
inline virtual ~LaguerreOrthoPoly1D() {}
inline long double weight(const long double x) const
{return x < 0.0L ? 0.0L : expl(-x);}
inline double xmin() const {return 0.0;}
inline double xmax() const {return 11355.1;}
virtual std::pair<long double,long double>
monicRecurrenceCoeffs(unsigned k) const;
virtual std::pair<long double,long double>
recurrenceCoeffs(unsigned k) const;
};
/** Orthonormal(!) Gram polynomials */
class GramOrthoPoly1D : public AbsClassicalOrthoPoly1D
{
public:
GramOrthoPoly1D(unsigned npoints);
inline virtual GramOrthoPoly1D* clone() const
{return new GramOrthoPoly1D(*this);}
inline virtual ~GramOrthoPoly1D() {}
long double weight(long double x) const;
inline double xmin() const {return -1.0;}
inline double xmax() const {return 1.0;}
inline virtual unsigned maxDegree() const
{return npoints_ - 1U;}
virtual std::pair<long double,long double>
monicRecurrenceCoeffs(unsigned k) const;
virtual std::pair<long double,long double>
recurrenceCoeffs(unsigned k) const;
private:
unsigned npoints_;
};
/** Orthonormal(!) shifter Gram polynomials */
class ShiftedGramOrthoPoly1D : public AbsClassicalOrthoPoly1D
{
public:
ShiftedGramOrthoPoly1D(unsigned npoints);
inline virtual ShiftedGramOrthoPoly1D* clone() const
{return new ShiftedGramOrthoPoly1D(*this);}
inline virtual ~ShiftedGramOrthoPoly1D() {}
long double weight(long double x) const;
inline double xmin() const {return -1.0;}
inline double xmax() const {return 1.0;}
inline virtual unsigned maxDegree() const
{return npoints_ - 1U;}
virtual std::pair<long double,long double>
monicRecurrenceCoeffs(unsigned k) const;
virtual std::pair<long double,long double>
recurrenceCoeffs(unsigned k) const;
private:
unsigned npoints_;
};
}
#endif // NPSTAT_CLASSICALORTHOPOLYS1D_HH_

File Metadata

Mime Type
text/x-diff
Expires
Wed, May 14, 10:49 AM (1 d, 2 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5111260
Default Alt Text
(29 KB)

Event Timeline