laguerreBuilder.h 8.68 KB
Newer Older
1 2
#ifndef LAGUERREBUILDER_SEEN
#define LAGUERREBUILDER_SEEN
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
/*
 *  This file is part of LagSHT.
 *
 *  LagSHT is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  LagSHT is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with libsharp; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

/*
 *  LagSHT is being developed at the Linear Accelerateur Laboratory (LAL)
 *  91898 ORSAY CEDEX - FRANCE
 *  main author: J.E Campagne
 */
26

27
#include <cmath> //very important for abs definition
28
#include <vector>
29
#include <iostream>
30
#include <iomanip>      // std::setprecision
31

32 33 34 35
//JEC 12/1/16: alpha becomes double instead of int
//JEC (N+alpha)! -> Gamma[N+alpha+1]. Use Boost::tgamma function
#include "boost/math/special_functions/gamma.hpp"

36 37 38 39
#include "lagsht_numbers.h"

using namespace std;

40

41
#define DEBUG 0
42

43 44 45
namespace LagSHT {


46
/////////////////////////////////////////////////////////////////////////////////////////////
47 48 49 50 51 52 53 54 55
//// - -Class generalized Laguerre Function and Gauss-Laguerre Quadrature --------------   //
////   The Laguerre Function is scaled by exp(-r/2) compared to the Laguerre Polynomials   //
////   the exp(r_i) scales the weights compared to the Polynomial Gauss-Laguerre Quadrature//
//// (nb. an public heritance from LaguerrePolyQuad would made confusion)                  //
/////////////////////////////////////////////////////////////////////////////////////////////

class LaguerreFuncQuad   {
public:
  //! Creator
56
  LaguerreFuncQuad(int N, r_8 alpha=2): N_(N), alpha_(alpha),
57 58 59 60 61 62 63 64 65 66 67 68 69
    ln2(0.693147180559945309417232121458176),
    inv_ln2(1.4426950408889634073599246810018921),
    fsmall(std::ldexp(1.,-large_exponent2)),
    fbig(std::ldexp(1.,large_exponent2)),
    xi(N+2), cf(maxscale+1-minscale) {

    using namespace std;
    xi[0]=0.;
    for (size_t i=1; i<xi.size(); ++i)
      xi[i]=1./i;
    for (size_t i=0; i<cf.size(); ++i)
      cf[i] = ldexp(1.,(int(i)+minscale)*large_exponent2);
  }
70

71
  //! Destructor
72 73
  virtual ~LaguerreFuncQuad() {
  }
74 75
   
  //! Getters
76
  int GetOrder() const { return N_; } 
77
  r_8 GetAlpha() const { return alpha_; }
78 79

  //! Compute Weights and Nodes of the Gauss-Laguerre quadrature
80
  void QuadWeightNodes(vector<r_8>& nodes, vector<r_8>& weights) const;
81

82
  //! Compute  L_n^alpha(x) * exp(-x/2) * scale with  n = N_ 
83
  r_8 Value(r_8 x, r_8 scale=1.0) const;
84
  
85 86 87 88 89
  //! Compute  L_n^alpha(x) * exp(-x/2) * scale with  n=0,...,N_ (included)
  void Values(r_8 x, vector<r_8>& vals, r_8 scale=1.0);  

protected:

90 91 92
  /*! returns the amount by which x must be shifted to go closer to a root of
    L_N^alpha(x)
  */
93
  r_8 newtoncorr (r_8 x, int N) const {
94 95 96
    using namespace std;
    if(N<1) return 0.;
    
97 98
    r_8 pnm1=0., pn=1.;
    r_8 c1=alpha_-1., c2=alpha_-1.;
99 100
    for(int i=1;i<=N;i++)
      {
101
        r_8 pnm2 = pnm1;
102 103 104 105 106 107 108 109
        pnm1 = pn;
        c1+=2.;
        c2+=1.;
        pn = ((c1-x)*pnm1 - c2*pnm2)*xi[i];
        // if necessary, scale down
        // the function result is independent of the absolute scale 
        while (abs(pn)>=fbig)
          { pn*=fsmall; pnm1*=fsmall; }
110
      }
111
    r_8 ppn=(N*pn-(N+alpha_)*pnm1)/x;
112
    return -pn/ppn;
113
   }//newtoncorr
114

115 116
  //JEC 12/1/16 add a normalization factor 
  /*! returns L_N^alpha(x)*exp(-x/2)*norm
117
   */
118
  r_8 value(r_8 x, int N, r_8 norm = 1.0) const {
119 120 121
    using namespace std;
    if (N<0) return 0;

122
    r_8 log2expfac = -0.5*x*inv_ln2;
123 124
    int scale=int(log2expfac/large_exponent2);
    log2expfac-=scale*large_exponent2;
125

126 127
    r_8 pnm1=0., pn=exp(ln2*log2expfac)*norm;
    r_8 c1=alpha_-1., c2=alpha_-1.;
128 129
    for(int i=1;i<=N;i++)
      {
130
        r_8 pnm2 = pnm1;
131 132 133 134 135 136
        pnm1 = pn;
        c1+=2.;
        c2+=1.;
        pn = ((c1-x)*pnm1 - c2*pnm2)*xi[i];
        while (abs(pn)>=fbig) // if necessary, scale down and count rescalings
          { pn*=fsmall; pnm1*=fsmall; ++scale; }
137
      }
138 139
    return (scale<minscale) ? 0 : pn*cf[scale-minscale];
  }//value
140

141 142
  //JEC 12/1/16 add a normalization factor 
  /*! returns L_n^alpha(x)*exp(-x/2)*norm [n=0..N] in res.
143
   */
144
  void values (r_8 x, int N, std::vector<r_8> &res, r_8 norm=1.0) {
145 146 147
    using namespace std;
    if (N<0) { res.resize(0); return; }
    res.resize(N+1);
148
    res[0]=exp(-0.5*x)*norm;
149

150
    r_8 log2expfac = -0.5*x*inv_ln2;
151 152
    int scale=int(log2expfac/large_exponent2);
    log2expfac-=scale*large_exponent2;
153

154 155
    r_8 pnm1=0., pn=exp(ln2*log2expfac)*norm;
    r_8 c1=alpha_-1., c2=alpha_-1.;
156 157
    for(int i=1;i<=N;i++)
      {
158
        r_8 pnm2 = pnm1;
159 160 161 162 163 164 165
        pnm1 = pn;
        c1+=2.;
        c2+=1.;
        pn = ((c1-x)*pnm1 - c2*pnm2)*xi[i];
        while (abs(pn)>=fbig) // if necessary, scale down and count rescalings
          { pn*=fsmall; pnm1*=fsmall; ++scale; }
        res[i] = (scale<minscale) ? 0 : pn*cf[scale-minscale];
166
      }
167 168 169 170 171
  }//values

  /*! computes nodes and weights for Gauss-Laguerre integration of order N
    NOTE: the weights are correct for modified Laguerre polynomials, i.e.
    L_N^alpha(x)*exp(-x/2)
172
    The nodes/weights are indexed i:0,...,N-1
173
  */
174 175
    void nodes_and_weights(int N, std::vector<r_8> &nodes,
			   std::vector<r_8> &weights) const {
176 177 178 179
    using namespace std;
    nodes.resize(N);
    weights.resize(N);

180 181 182 183 184 185 186 187 188 189
    /* JEC 12/1/16 alpha is r_8. tmp1 = (N+alpha)!/N!/(N+1)^2 = Gamma(N+alpha+1)/Gamma(N+1)/(N+1)^2
      r_8 tmp1 = N+alpha_;
      for (int i=1;i<alpha_-1;i++) tmp1 *= (tmp1-1.);
      tmp1 /= (N+1.);

      Use the ratio of Gamma functions as in Boost library: Todo see if one should tune the Policy
    */

    r_8 tmp1 = boost::math::tgamma_ratio((r_8)N+alpha_+1.0, r_8(N)+1.0);
    tmp1 /= (N+1.)*(N+1.);
190

191
    r_8 x=0.;
192
    // Approximation by Stroud & Secrest, Gaussian quadrature formulas, 1966
193 194 195 196 197 198 199 200
    for(int k=0; k<N; k++) { 
      if (k==0)
	x = (1.0+alpha_)*(3.+0.92*alpha_)/(1.+2.4*N+1.8*alpha_);
      else if (k==1)
	x += (15.+6.25*alpha_)/(1.+0.9*alpha_+2.5*N);
      else
	{
	  int km1 = k-1;
201 202
	  r_8 c0 = (1.+2.55*km1)/(1.9*km1)+1.26*km1*alpha_/(1.+3.5*km1);
	  r_8 c1 = 1.+0.3*alpha_;
203 204 205 206 207 208 209 210
	  x += c0*(x-nodes[k-2])/c1;
	}
      
      // Refinement using Newton's method
      // Once we notice that we are close, do one more step and then exit
      bool done=false;
      while (true)
	{
211
	  r_8 xold=x;
212 213 214 215 216 217 218 219
	  x+=newtoncorr (x, N);
	  
	  if (done) break;
	  if (abs(x-xold)<max(1e-12,1e-12*abs(x))) done=true;
	}
      
      nodes[k] = x;
      
220
      r_8 tmp0 = value(x, N+1);
221 222 223
      weights[k] = x*tmp1/(tmp0*tmp0);
    }//end for

224
#if DEBUG >= 1  
225 226 227 228 229
    //verification 
    r_16 sumR = 0.;
    r_16 sumW = 0.;
    r_16 wip;
    for (int i=0;i<N;i++){
230 231
      
      cout << "root["<<i<<"]: " << setprecision(16) << nodes[i] << " LagN(ri) = " << Value(nodes[i]) << endl; 
232

233 234 235 236 237
      sumR += nodes[i];
      wip = exp(log(abs(weights[i]))-nodes[i]); if(weights[i]<0) wip *= -1.0;
      sumW += wip;
      //    sumW += weights[i]*exp(-nodes[i]/tau);
    }
238

239 240 241 242 243 244 245
    /* JEC 12/1/16 alpha ! = alphaFact = Gamma(alpha+1)
      int alphaFact = 1;
      for(int i=1;i<=alpha_; i++) alphaFact *= i;
      r_16 sumRTh = (r_16)(N*(N+alpha_));
      r_16 sumWTh = (r_16)(alphaFact);
    */
    
246
    r_16 sumRTh = (r_16)(N*(N+alpha_));
247 248
    r_16 sumWTh = boost::math::tgamma((r_16)(alpha_+1.0));

249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267

    cout << "Sum roots = " << sumR << " -> diff with theory: " << sumR - sumRTh << endl;
    cout << "Sum weights = " << sumW << " -> diff with theory: " << sumW - sumWTh << endl;
    
    if (fabs(sumR-sumRTh)>(0.000001*sumRTh)) {
      cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n" 
	   << "WARNING!!!! verify LaguerrePolyQuad sum of nodes not accurate: "
	   << " current value " << sumR << " theoretical value : " << sumRTh 
	   << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n" 
	   << endl;
    }
    
    if (fabs(sumW-sumWTh)>(0.000001*sumWTh)) {
      cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n" 
	   << "WARNING!!!! verify LaguerrePolyQuad sum of weights not accurate: "
	   << " current value " << sumW << " theoretical value : " << sumWTh 
	   << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n" 
	   << endl;
    }
268
#endif
269
    
270
  }//nodes_and_weights
271 272 273

protected:
  int N_; //!< Order of the polynomial
274
  r_8 alpha_; //!< second parameter of the generalized Laguerre polynomial 
275

276 277
private:
  enum { large_exponent2=90, minscale=-14, maxscale=14 };
278 279 280
  r_8 ln2, inv_ln2, fsmall, fbig;
  std::vector<r_8> xi; 
  std::vector<r_8> cf; 
281 282


283 284 285 286 287 288
}; //LaguerreFuncQuad



} // Fin du namespace

289
#undef DEBUG
290 291

#endif //LAGUERREBUILDER_SEEN