1 /*
2  *  This file is part of the MR utility library.
3  *
4  *  This code is free software; you can redistribute it and/or modify
5  *  it under the terms of the GNU General Public License as published by
6  *  the Free Software Foundation; either version 2 of the License, or
7  *  (at your option) any later version.
8  *
9  *  This code is distributed in the hope that it will be useful,
10  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
query_cb(struct Address * result,void * data)11  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  *  GNU General Public License for more details.
13  *
14  *  You should have received a copy of the GNU General Public License
15  *  along with this code; if not, write to the Free Software
16  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
17  */
18 
19 /** \file ducc0/math/gl_integrator.h
20  *  Functionality for Gauss-Legendre quadrature
21  *
22  *  \copyright Copyright (C) 2019-2021 Max-Planck-Society, Ignace Bogaert
23  *  \author Martin Reinecke
24  */
25 
26 #ifndef DUCC0_GL_INTEGRATOR_H
27 #define DUCC0_GL_INTEGRATOR_H
28 
29 #include <cmath>
30 #include <array>
31 #include <cstddef>
32 #include <utility>
33 #include <vector>
34 #include "ducc0/math/constants.h"
35 #include "ducc0/infra/error_handling.h"
36 
37 namespace ducc0 {
38 
39 namespace detail_gl_integrator {
40 
41 using namespace std;
42 
43 template<typename T> inline T one_minus_x2 (T x)
44   { if (x<0) x=-x; return (x>T(0.1)) ? (T(1)+x)*(T(1)-x) : T(1)-x*x; }
45 
46 pair<double, double> calc_gl_iterative(size_t n, size_t i)
47   {
48   using Tfloat = long double;
49   constexpr Tfloat eps=Tfloat(3e-14L);
50   const Tfloat dn = Tfloat(n);
51   const Tfloat t0 = Tfloat(1) - (1-Tfloat(1)/dn) / (Tfloat(8)*dn*dn);
52   const Tfloat t1 = Tfloat(1)/(Tfloat(4)*dn+Tfloat(2));
53   Tfloat x0 = cos(double(pi * ((i<<2)-1) * t1)) * t0;
54 
55   bool dobreak=false;
56   size_t j=0;
57   Tfloat dpdx;
58   while(1)
59     {
60     Tfloat P_1 = 1;
61     Tfloat P0 = x0;
62     Tfloat dx, x1;
63 
64     for (size_t k=2; k<=n; k++)
65       {
66       Tfloat P_2 = P_1;
67       P_1 = P0;
68       P0 = x0*P_1 + (Tfloat(k)-Tfloat(1))/Tfloat(k) * (x0*P_1-P_2);
69       }
70 
71     dpdx = (P_1 - x0*P0) * n / one_minus_x2(x0);
72 
73     /* Newton step */
74     x1 = x0 - P0/dpdx;
75     dx = x0-x1;
76     x0 = x1;
77     if (dobreak) break;
78 
79     if ((dx<0?-dx:dx)<=eps) dobreak=1;
80     MR_assert(++j<100, "convergence problem");
81     }
82 
83   return make_pair(double(x0), double(2./(one_minus_x2(x0)*dpdx*dpdx)));
84   }
85 
86 
87 // The next three functions are modified versions of the FastGL code
88 // by Ignace Bogaert. The code is available at
89 // https://sourceforge.net/projects/fastgausslegendrequadrature/
90 // A paper describing the algorithms is available at
91 // https://epubs.siam.org/doi/pdf/10.1137/140954969
92 
93 // This function computes the kth zero of the BesselJ(0,x)
94 double besseljzero(int k)
95   {
96   constexpr static array<double,20> JZ
97     {2.40482555769577276862163187933,  5.52007811028631064959660411281,
98      8.65372791291101221695419871266,  11.7915344390142816137430449119,
99      14.9309177084877859477625939974, 18.0710639679109225431478829756,
100      21.2116366298792589590783933505, 24.3524715307493027370579447632,
101      27.4934791320402547958772882346, 30.6346064684319751175495789269,
102      33.7758202135735686842385463467, 36.9170983536640439797694930633,
103      40.0584257646282392947993073740, 43.1997917131767303575240727287,
104      46.3411883716618140186857888791, 49.4826098973978171736027615332,
105      52.6240518411149960292512853804, 55.7655107550199793116834927735,
106      58.9069839260809421328344066346, 62.0484691902271698828525002646};
107 
108   if (k<=20) return JZ[k-1];
109 
110   double z = pi*(k-0.25);
111   double r = 1.0/z;
112   double r2 = r*r;
113   return z
114        + r*(0.125
115        +r2*(-0.807291666666666666666666666667e-1
116        +r2*(0.246028645833333333333333333333
117        +r2*(-1.82443876720610119047619047619
118        +r2*(25.3364147973439050099206349206
119        +r2*(-567.644412135183381139802038240
120        +r2*(18690.4765282320653831636345064
121        +r2*(-8.49353580299148769921876983660e5
122        +r2*5.09225462402226769498681286758e7))))))));
123   }
124 
125 // This function computes the square of BesselJ(1, BesselZero(0,k))
126 double besselj1squared(int k)
127   {
128   constexpr static array<double,21> J1
129     {0.269514123941916926139021992911 , 0.115780138582203695807812836182,
130      0.0736863511364082151406476811985, 0.0540375731981162820417749182758,
131      0.0426614290172430912655106063495, 0.0352421034909961013587473033648,
132      0.0300210701030546726750888157688, 0.0261473914953080885904584675399,
133      0.0231591218246913922652676382178, 0.0207838291222678576039808057297,
134      0.0188504506693176678161056800214, 0.0172461575696650082995240053542,
135      0.0158935181059235978027065594287, 0.0147376260964721895895742982592,
136      0.0137384651453871179182880484134, 0.0128661817376151328791406637228,
137      0.0120980515486267975471075438497, 0.0114164712244916085168627222986,
138      0.0108075927911802040115547286830, 0.0102603729262807628110423992790,
139      0.00976589713979105054059846736696};
140 
141   if (k<=21) return J1[k-1];
142 
143   double x = 1.0/(k-0.25);
144   double x2 = x*x;
145   return x * (0.202642367284675542887758926420
146        + x2*x2*(-0.303380429711290253026202643516e-3
147        + x2*(0.198924364245969295201137972743e-3
148        + x2*(-0.228969902772111653038747229723e-3
149        + x2*(0.433710719130746277915572905025e-3
150        + x2*(-0.123632349727175414724737657367e-2
151        + x2*(0.496101423268883102872271417616e-2
152        + x2*(-0.266837393702323757700998557826e-1
153        + x2*.185395398206345628711318848386))))))));
154   }
155 
156 // Compute a node-weight pair, with k limited to half the range
157 pair<double,double> calc_gl_bogaert(size_t n, size_t k0)
158   {
159   size_t k = ((2*k0-1)<=n) ? k0 : n-k0+1;
160   // First get the Bessel zero
161   double w = 1.0/(n+0.5);
162   double nu = besseljzero(k);
163   double theta = w*nu;
164   double x = theta*theta;
165 
166   // Get the asymptotic BesselJ(1,nu) squared
167   double B = besselj1squared(k);
168 
169   // Get the Chebyshev interpolants for the nodes...
170   double SF1T = (((((-1.29052996274280508473467968379e-12*x
171                      +2.40724685864330121825976175184e-10)*x
172                      -3.13148654635992041468855740012e-8)*x
173                      +0.275573168962061235623801563453e-5)*x
174                      -0.148809523713909147898955880165e-3)*x
175                      +0.416666666665193394525296923981e-2)*x
176                      -0.416666666666662959639712457549e-1;
177   double SF2T = (((((+2.20639421781871003734786884322e-9*x
178                      -7.53036771373769326811030753538e-8)*x
179                      +0.161969259453836261731700382098e-5)*x
180                      -0.253300326008232025914059965302e-4)*x
181                      +0.282116886057560434805998583817e-3)*x
182                      -0.209022248387852902722635654229e-2)*x
183                      +0.815972221772932265640401128517e-2;
184   double SF3T = (((((-2.97058225375526229899781956673e-8*x
185                      +5.55845330223796209655886325712e-7)*x
186                      -0.567797841356833081642185432056e-5)*x
187                      +0.418498100329504574443885193835e-4)*x
188                      -0.251395293283965914823026348764e-3)*x
189                      +0.128654198542845137196151147483e-2)*x
190                      -0.416012165620204364833694266818e-2;
191 
192   // ...and for the weights
193   double WSF1T = ((((((((-2.20902861044616638398573427475e-14*x
194                          +2.30365726860377376873232578871e-12)*x
195                          -1.75257700735423807659851042318e-10)*x
196                          +1.03756066927916795821098009353e-8)*x
197                          -4.63968647553221331251529631098e-7)*x
198                          +0.149644593625028648361395938176e-4)*x
199                          -0.326278659594412170300449074873e-3)*x
200                          +0.436507936507598105249726413120e-2)*x
201                          -0.305555555555553028279487898503e-1)*x
202                          +0.833333333333333302184063103900e-1;
203   double WSF2T = (((((((+3.63117412152654783455929483029e-12*x
204                         +7.67643545069893130779501844323e-11)*x
205                         -7.12912857233642220650643150625e-9)*x
206                         +2.11483880685947151466370130277e-7)*x
207                         -0.381817918680045468483009307090e-5)*x
208                         +0.465969530694968391417927388162e-4)*x
209                         -0.407297185611335764191683161117e-3)*x
210                         +0.268959435694729660779984493795e-2)*x
211                         -0.111111111111214923138249347172e-1;
212   double WSF3T = (((((((+2.01826791256703301806643264922e-9*x
213                         -4.38647122520206649251063212545e-8)*x
214                         +5.08898347288671653137451093208e-7)*x
215                         -0.397933316519135275712977531366e-5)*x
216                         +0.200559326396458326778521795392e-4)*x
217                         -0.422888059282921161626339411388e-4)*x
218                         -0.105646050254076140548678457002e-3)*x
219                         -0.947969308958577323145923317955e-4)*x
220                         +0.656966489926484797412985260842e-2;
221 
222   // Then refine with the paper expansions
223   double NuoSin = nu/sin(theta);
224   double BNuoSin = B*NuoSin;
225   double WInvSinc = w*w*NuoSin;
226   double WIS2 = WInvSinc*WInvSinc;
227 
228   // Finally compute the node and the weight
229   theta = w*(nu + theta * WInvSinc * (SF1T + WIS2*(SF2T + WIS2*SF3T)));
230   double Deno = BNuoSin + BNuoSin * WIS2*(WSF1T + WIS2*(WSF2T + WIS2*WSF3T));
231   double weight = (2.0*w)/Deno;
232   return make_pair((k==k0) ? cos(theta) : -cos(theta), weight);
233   }
234 
235 pair<double, double> calc_gl(size_t n, size_t k)
236   {
237   MR_assert(n>=k, "k must not be greater than n");
238   MR_assert(k>0, "k must be positive");
239   return (n<=100) ? calc_gl_iterative(n,k) : calc_gl_bogaert(n,k);
240   }
241 
242 /// Class for computing Gauss-Lgendre abscissas, weights and intgrals
243 class GL_Integrator
244   {
245   private:
246     size_t n_;
247     vector<double> x, w;
248 
249   public:
250     /// Creates an integrator for \a n abscissas
251     /** \note The \a nthreads parameter is obsolescent and ignored. */
252     GL_Integrator(size_t n, size_t /*nthreads*/=1)
253       : n_(n)
254       {
255       MR_assert(n>=1, "number of points must be at least 1");
256       size_t m = (n+1)>>1;
257       x.resize(m);
258       w.resize(m);
259       for (size_t i=0; i<m; ++i)
260         {
261         auto tmp = calc_gl(n, m-i);
262         x[i] = tmp.first;
263         w[i] = tmp.second;
264         }
265       }
266 
267     /// Returns the approximated integral of \a func in [-1;1]
268     template<typename Func> auto integrate(Func f) -> decltype(f(0.))
269       {
270       using T = decltype(f(0.));
271       T res=0;
272       size_t istart=0;
273       if (n_&1)
274         {
275         res = f(x[0])*w[0];
276         istart=1;
277         }
278       for (size_t i=istart; i<x.size(); ++i)
279         res += (f(x[i])+f(-x[i]))*w[i];
280       return res;
281       }
282 
283     /// Returns the approximated integral of \a func in [-1;1], where \a func
284     /// is symmetric with respect to x=0.
285     template<typename Func> auto integrateSymmetric(Func f) -> decltype(f(0.))
286       {
287       using T = decltype(f(0.));
288       T res=f(x[0])*w[0];
289       if (n_&1) res *= 0.5;
290       for (size_t i=1; i<x.size(); ++i)
291         res += f(x[i])*w[i];
292       return res*2;
293       }
294 
295     /// Returns the Gauss-Legendre abscissas.
296     vector<double> coords() const
297       {
298       vector<double> res(n_);
299       for (size_t i=0; i<x.size(); ++i)
300         {
301         res[i]=-x[x.size()-1-i];
302         res[n_-1-i] = x[x.size()-1-i];
303         }
304       return res;
305       }
306     /// Returns the non-negative Gauss-Legendre abscissas.
307     const vector<double> &coordsSymmetric() const
308       { return x; }
309 
310     /// Returns the Gauss-Legendre weights.
311     vector<double> weights() const
312       {
313       vector<double> res(n_);
314       for (size_t i=0; i<w.size(); ++i)
315         res[i]=res[n_-1-i]=w[w.size()-1-i];
316       return res;
317       }
318     /// Returns the Gauss-Legendre weights for the non-negative abscissas,
319     /// with an additional factor of 2 for positive abscissas.
320     vector<double> weightsSymmetric() const
321       {
322       auto res = w;
323       if (n_&1) res[0]*=0.5;
324       for (auto &v:res) v*=2;
325       return res;
326       }
327   };
328 
329 }
330 
331 using detail_gl_integrator::GL_Integrator;
332 
333 }
334 
335 #endif
336