HepLib
Lib3_GaussKronrod.h
Go to the documentation of this file.
1 // modified version from https://github.com/drjerry/quadpackpp
2 // GaussKronrod Rules
3 #pragma once
4 #include "mpreal.h"
5 #include <functional>
6 
7 /* error return codes */
8 #define SUCCESS 0
9 #define FAILURE 1
10 
11 namespace {
12  typedef mpfr::mpreal Real;
13  typedef const mpfr::mpreal & Real_t;
14 
15 
16  class GaussKronrod {
17  protected:
18  size_t m_; // Gauss-Legendre degree
19  size_t n_; // size of Gauss-Kronrod arrays
20  Real* xgk_; // Gauss-Kronrod abscissae
21  Real* wg_; // Gauss-Legendre weights
22  Real* wgk_; // Gauss-Kronrod weights
23 
24  private:
25  Real *coefs; // Chebyshev coefficients
26  Real *zeros; // zeros of Legendre polynomial
27 
28  void legendre_zeros();
29  void chebyshev_coefs();
30  void gauss_kronrod_abscissae();
31  void gauss_kronrod_weights();
32 
33  Real legendre_err(int deg, Real_t x, Real& err);
34  Real legendre_deriv(int deg, Real_t x);
35  Real chebyshev_series(Real_t x, Real& err);
36  Real chebyshev_series_deriv(Real_t x);
37 
38  public:
39  GaussKronrod(size_t m = 10);
40  ~GaussKronrod();
41  size_t size() { return n_; };
42  Real xgk(int k) { return (0 <= k && k < n_) ? xgk_[k] : Real(0); }
43  Real wgk(int k) { return (0 <= k && k < n_) ? wgk_[k] : Real(0); }
44  Real wg(int k) { return (0 <= k && k < n_/2) ? wg_[k] : Real(0); }
45  };
46 
47  GaussKronrod::GaussKronrod(size_t m) {
48  m_ = m;
49  n_ = m_ + 1;
50  xgk_ = new Real[n_];
51  wg_ = new Real[n_ / 2];
52  wgk_ = new Real[n_];
53  coefs = new Real[n_ + 1];
54  zeros = new Real[m_ + 2];
55 
56  legendre_zeros();
57  chebyshev_coefs();
58  gauss_kronrod_abscissae();
59  gauss_kronrod_weights();
60  }
61 
62  GaussKronrod::~GaussKronrod() {
63  delete[] xgk_;
64  delete[] wgk_;
65  delete[] wg_;
66  delete[] coefs;
67  delete[] zeros;
68  }
69 
70  void GaussKronrod::legendre_zeros() {
71  Real* temp = new Real[m_+1];
72  zeros[0] = Real(-1);
73  zeros[1] = Real(1);
74  Real delta, epsilon;
75 
76  for (int k = 1; k <= m_; ++k) {
77  // loop to locate zeros of P_k interlacing z_0,...,z_k
78  for (int j = 0; j < k; ++j) {
79  // Newton's method for P_k :
80  // initialize solver at midpoint of (z_j, z_{j+1})
81  delta = 1;
82  Real x_j = (zeros[j] + zeros[j+1]) / 2;
83  Real P_k = legendre_err(k, x_j, epsilon);
84  while (abs(P_k) > epsilon && abs(delta) > mpfr::machine_epsilon()) {
85  delta = P_k / legendre_deriv(k, x_j);
86  x_j -= delta;
87  P_k = legendre_err(k, x_j, epsilon);
88  }
89  temp[j] = x_j;
90  }
91 
92  // copy roots tmp_0,...,tmp_{k-1} to z_1,...z_k:
93  zeros[k+1] = zeros[k];
94  for (int j = 0; j < k; ++j) zeros[j+1] = temp[j];
95  }
96  delete[] temp;
97  }
98 
99  void GaussKronrod::chebyshev_coefs() {
100  size_t ell = (m_ + 1)/2;
101  Real* alpha = new Real[ell+1];
102  Real* f = new Real[ell+1];
103 
104  /* Care must be exercised in initalizing the constants in the definitions.
105  * Compilers interpret expressions like "(2*k + 1.0)/(k + 2.0)" as floating
106  * point precision, before casting to Real.
107  */
108  f[1] = Real(m_+1)/Real(2*m_ + 3);
109  alpha[0] = Real(1); // coefficient of T_{m+1}
110  alpha[1] = -f[1];
111 
112  for (int k = 2; k <= ell; ++k) {
113  f[k] = f[k-1] * (2*k - 1) * (m_ + k) / (k * (2*m_ + 2*k + 1));
114  alpha[k] = -f[k];
115  for (int i = 1; i < k; ++i) alpha[k] -= f[i] * alpha[k-i];
116  }
117 
118  for (int k = 0; k <= ell; ++k) {
119  coefs[m_ + 1 - 2*k] = alpha[k];
120  if (m_ >= 2*k) coefs[m_ - 2*k] = Real(0);
121  }
122 
123  delete[] alpha;
124  delete[] f;
125  }
126 
127  void GaussKronrod::gauss_kronrod_weights() {
128  Real err;
129  // Gauss-Legendre weights:
130  for (int k = 0; k < n_ / 2; ++k) {
131  Real x = xgk_[2*k + 1];
132  wg_[k] = (Real(-2) / ((m_ + 1) * legendre_deriv(m_, x) * legendre_err(m_+1, x, err)));
133  }
134 
135  /* The ratio of leading coefficients of P_n and T_{n+1} is computed
136  * from the recursive formulae for the respective polynomials.
137  */
138  Real F_m = Real(2) / Real(2*m_ + 1);
139  for (int k = 1; k <= m_; ++k) F_m *= (Real(2*k) / Real(2*k - 1));
140 
141  // Gauss-Kronrod weights:
142  for (size_t k = 0; k < n_; ++k) {
143  Real x = xgk_[k];
144  if (k % 2 == 0) {
145  wgk_[k] = F_m / (legendre_err(m_, x, err) * chebyshev_series_deriv(x));
146  } else {
147  wgk_[k] = (wg_[k/2] + F_m / (legendre_deriv(m_, x) * chebyshev_series(x, err)));
148  }
149  }
150  }
151 
152  void GaussKronrod::gauss_kronrod_abscissae() {
153  Real delta, epsilon;
154 
155  for (int k = 0; 2*k < n_; ++k) {
156  delta = 1;
157  // Newton's method for E_{n+1} :
158  Real x_k = (zeros[m_-k] + zeros[m_+1-k])/Real(2);
159  Real E = chebyshev_series(x_k, epsilon);
160  while (abs(E) > epsilon && abs(delta) > mpfr::machine_epsilon()) {
161  delta = E / chebyshev_series_deriv(x_k);
162  x_k -= delta;
163  E = chebyshev_series(x_k, epsilon);
164  }
165  xgk_[2*k] = x_k;
166  // copy adjacent Legendre-zero into the array:
167  if (2*k+1 < n_) xgk_[2*k+1] = zeros[m_-k];
168  }
169  }
170 
171  Real GaussKronrod::legendre_err(int n, Real_t x, Real& err) {
172  if (n == 0) {
173  err = Real(0);
174  return Real(1);
175  }
176  else if (n == 1) {
177  err = Real(0);
178  return x;
179  }
180 
181  Real P0 = Real(1), P1 = x, P2;
182  Real E0 = mpfr::machine_epsilon();
183  Real E1 = abs(x) * mpfr::machine_epsilon();
184  for (int k = 1; k < n; ++k) {
185  P2 = ((2*k + 1) * x * P1 - k * P0) / (k + 1);
186  err = ((2*k + 1) * abs(x) * E1 + k * E0) / (2*(k + 1));
187  P0 = P1; P1 = P2;
188  E0 = E1; E1 = err;
189  }
190  return P2;
191  }
192 
193  Real GaussKronrod::legendre_deriv(int n, Real_t x) {
194  if (n == 0) return Real(0);
195  else if (n == 1) return Real(1);
196 
197  Real P0 = Real(1), P1 = x, P2;
198  Real dP0 = Real(0), dP1 = Real(1), dP2;
199  for (int k = 1; k < n; ++k) {
200  P2 = ((2*k + 1) * x * P1 - k * P0) / (k + Real(1));
201  dP2 = (2*k + 1) * P1 + dP0;
202  P0 = P1; P1 = P2;
203  dP0 = dP1; dP1 = dP2;
204  }
205  return dP2;
206  }
207 
208  Real GaussKronrod::chebyshev_series(Real_t x, Real& err) {
209  Real d1(0), d2(0);
210  Real absc = abs(coefs[0]); // final term for truncation error
211  Real y2 = 2 * x; // linear term for Clenshaw recursion
212 
213  for (int k = n_; k >= 1; --k) {
214  Real temp = d1;
215  d1 = y2 * d1 - d2 + coefs[k];
216  d2 = temp;
217  absc += abs(coefs[k]);
218  }
219 
220  err = absc * mpfr::machine_epsilon();
221  return x * d1 - d2 + coefs[0]/2;
222  }
223 
224  Real GaussKronrod::chebyshev_series_deriv(Real_t x) {
225  Real d1(0), d2(0);
226  Real y2 = 2 * x; // linear term for Clenshaw recursion
227 
228  for (int k = n_; k >= 2; --k) {
229  Real temp = d1;
230  d1 = y2 * d1 - d2 + k * coefs[k];
231  d2 = temp;
232  }
233 
234  return y2 * d1 - d2 + coefs[1];
235  }
236 
237 }
Generated by doxygen 1.9.1