npstat is hosted by Hepforge, IPPP Durham
NPStat  5.10.0
GaussHermiteQuadrature.hh
Go to the documentation of this file.
1 #ifndef NPSTAT_GAUSSHERMITEQUADRATURE_HH_
2 #define NPSTAT_GAUSSHERMITEQUADRATURE_HH_
3 
4 /*!
5 // \file GaussHermiteQuadrature.hh
6 //
7 // \brief Gauss-Hermite quadratures in long double precision
8 //
9 // Author: I. Volobouev
10 //
11 // April 2010
12 */
13 
14 #include <vector>
15 #include <utility>
16 #include <cassert>
17 
19 
20 namespace npstat {
21  /**
22  // Gauss-Hermite quadrature. Internally, operations are performed
23  // in long double precision.
24  */
26  {
27  public:
28  enum {interval_quadrature = 0};
29 
30  /**
31  // At the moment, the following numbers of points are supported:
32  // 16, 32, 64, 100, 128, 256, 512.
33  //
34  // If an unsupported number of points is given in the
35  // constructor, std::invalid_argument exception will be thrown.
36  //
37  // Note that, for the quadrature with 512 points, the
38  // arguments become as large as 31.4 and the weights become
39  // as small as 4.9e-430. Use of such weights with doubles is
40  // rather meaningless. Long double integrands are needed.
41  */
42  explicit GaussHermiteQuadrature(unsigned npoints);
43 
44  /** Return the number of quadrature points */
45  inline unsigned npoints() const {return npoints_;}
46 
47  /**
48  // The abscissae are returned for positive points only,
49  // so the buffer length should be at least npoints/2.
50  */
51  void getAbscissae(long double* abscissae, unsigned len) const;
52 
53  /** Return abscissae for all points */
54  void getAllAbscissae(long double* abscissae, unsigned len) const;
55 
56  /**
57  // The weights are returned for positive points only,
58  // so the buffer length should be at least npoints/2.
59  */
60  void getWeights(long double* weights, unsigned len) const;
61 
62  /** Return weights for all points */
63  void getAllWeights(long double* weights, unsigned len) const;
64 
65  /** Perform the quadrature */
66  template <typename FcnResult, typename FcnArg>
67  long double integrate(const Functor1<FcnResult,FcnArg>& fcn) const;
68 
69  /** Perform the quadrature with Gaussian density weight */
70  template <typename FcnResult, typename FcnArg>
71  long double integrateProb(long double mean, long double sigma,
72  const Functor1<FcnResult,FcnArg>& fcn) const;
73 
74  /** Check if the rule with the given number of points is supported */
75  static bool isAllowed(unsigned npoints);
76 
77  /** The complete set of allowed rules, in the increasing order */
78  static std::vector<unsigned> allowedNPonts();
79 
80  /**
81  // Minimum number of points, among the supported rules, which
82  // integrates a polynomial with the given degree exactly (with
83  // the appropriate weight). Returns 0 if the degree is out of range.
84  */
85  static unsigned minimalExactRule(unsigned polyDegree);
86 
87  private:
89 
90  const long double* a_;
91  const long double* w_;
92  mutable std::vector<std::pair<long double, long double> > buf_;
93  unsigned npoints_;
94 
95 #ifdef SWIG
96  public:
97  inline std::vector<double> abscissae2() const
98  {
99  return std::vector<double>(a_, a_+npoints_/2U);
100  }
101 
102  inline std::vector<double> weights2() const
103  {
104  return std::vector<double>(w_, w_+npoints_/2U);
105  }
106 
107  inline std::vector<double> allabscissae2() const
108  {
109  assert(a_);
110  std::vector<double> abscissae(npoints_);
111  const unsigned halfpoints = npoints_/2;
112  for (unsigned i=0; i<halfpoints; ++i)
113  abscissae[i] = -a_[halfpoints - 1U - i];
114  for (unsigned i=halfpoints; i<npoints_; ++i)
115  abscissae[i] = a_[i - halfpoints];
116  return abscissae;
117  }
118 
119  inline std::vector<double> allweights2() const
120  {
121  assert(w_);
122  std::vector<double> weights(npoints_);
123  const unsigned halfpoints = npoints_/2;
124  for (unsigned i=0; i<halfpoints; ++i)
125  weights[i] = w_[halfpoints - 1U - i];
126  for (unsigned i=halfpoints; i<npoints_; ++i)
127  weights[i] = w_[i - halfpoints];
128  return weights;
129  }
130 
131  template <typename FcnResult, typename FcnArg>
132  inline double integrate2(const Functor1<FcnResult,FcnArg>& fcn) const
133  {
134  return integrate(fcn);
135  }
136 
137  template <typename FcnResult, typename FcnArg>
138  inline double integrateProb2(const double mean, const double sigma,
139  const Functor1<FcnResult,FcnArg>& fcn) const
140  {
141  return integrateProb(mean, sigma, fcn);
142  }
143 #endif // SWIG
144  };
145 }
146 
147 #include "npstat/nm/GaussHermiteQuadrature.icc"
148 
149 #endif // NPSTAT_GAUSSHERMITEQUADRATURE_HH_
Interface definitions and concrete simple functors for a variety of functor-based calculations.
Definition: GaussHermiteQuadrature.hh:26
void getAllWeights(long double *weights, unsigned len) const
void getAbscissae(long double *abscissae, unsigned len) const
static unsigned minimalExactRule(unsigned polyDegree)
long double integrate(const Functor1< FcnResult, FcnArg > &fcn) const
void getWeights(long double *weights, unsigned len) const
long double integrateProb(long double mean, long double sigma, const Functor1< FcnResult, FcnArg > &fcn) const
GaussHermiteQuadrature(unsigned npoints)
void getAllAbscissae(long double *abscissae, unsigned len) const
static bool isAllowed(unsigned npoints)
unsigned npoints() const
Definition: GaussHermiteQuadrature.hh:45
static std::vector< unsigned > allowedNPonts()
Definition: AbsArrayProjector.hh:14
Definition: SimpleFunctors.hh:58