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