laguerreBuilder.h 7.78 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 30
#include <iostream>

31 32 33 34
#include "lagsht_numbers.h"

using namespace std;

35

36
#define DEBUG 0
37

38 39 40
namespace LagSHT {


41
  /////////////////////////////////////////////////////////////////////////////////////////////
42 43 44 45 46 47 48 49 50
//// - -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
51 52 53 54 55 56 57 58 59 60 61 62 63 64
  LaguerreFuncQuad(int N, int alpha=2): N_(N), alpha_(alpha),
    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);
  }
65

66
  //! Destructor
67 68
  virtual ~LaguerreFuncQuad() {
  }
69 70
   
  //! Getters
71 72
  int GetOrder() const { return N_; } 
  int GetAlpha() const { return alpha_; }
73 74

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

77
  //! Compute  L_n^alpha(x) * exp(-x/2) * scale with  n = N_ 
78
  r_8 Value(r_8 x, r_8 scale=1.0) const;
79
  
80 81 82 83 84
  //! 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:

85 86 87 88 89 90 91 92 93
  /*! returns the amount by which x must be shifted to go closer to a root of
    L_N^alpha(x)
  */
  double newtoncorr (double x, int N) const {
    using namespace std;
    if(N<1) return 0.;
    
    double pnm1=0., pn=1.;
    double c1=alpha_-1., c2=alpha_-1.;
94 95
    for(int i=1;i<=N;i++)
      {
96 97 98 99 100 101 102 103 104
        double pnm2 = pnm1;
        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; }
105
      }
106
    double ppn=(N*pn-(N+alpha_)*pnm1)/x;
107
    return -pn/ppn;
108
   }//newtoncorr
109

110 111 112
  /*! returns L_N^alpha(x)*exp(-x/2)
   */
  double value(double x, int N) const {
113 114 115
    using namespace std;
    if (N<0) return 0;

116 117 118
    double log2expfac = -0.5*x*inv_ln2;
    int scale=int(log2expfac/large_exponent2);
    log2expfac-=scale*large_exponent2;
119

120 121
    double pnm1=0., pn=exp(ln2*log2expfac);
    double c1=alpha_-1., c2=alpha_-1.;
122 123
    for(int i=1;i<=N;i++)
      {
124 125 126 127 128 129 130
        double pnm2 = pnm1;
        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; }
131
      }
132 133
    return (scale<minscale) ? 0 : pn*cf[scale-minscale];
  }//value
134

135 136 137
  /*! returns L_n^alpha(x)*exp(-x/2) [n=0..N] in res.
   */
  void values (double x, int N, std::vector<double> &res) {
138 139 140 141 142
    using namespace std;
    if (N<0) { res.resize(0); return; }
    res.resize(N+1);
    res[0]=exp(-0.5*x);

143 144 145
    double log2expfac = -0.5*x*inv_ln2;
    int scale=int(log2expfac/large_exponent2);
    log2expfac-=scale*large_exponent2;
146

147 148
    double pnm1=0., pn=exp(ln2*log2expfac);
    double c1=alpha_-1., c2=alpha_-1.;
149 150
    for(int i=1;i<=N;i++)
      {
151 152 153 154 155 156 157 158
        double pnm2 = pnm1;
        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];
159
      }
160 161 162 163 164 165 166 167
  }//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)
  */
    void nodes_and_weights(int N, std::vector<double> &nodes,
			   std::vector<double> &weights) const {
168 169 170 171
    using namespace std;
    nodes.resize(N);
    weights.resize(N);

172 173
    double tmp1 = N+alpha_;
    for (int i=1;i<alpha_-1;i++) tmp1 *= (tmp1-1.);
174 175 176
    tmp1 /= (N+1.);

    double x=0.;
177
    // Approximation by Stroud & Secrest, Gaussian quadrature formulas, 1966
178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
    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;
	  double c0 = (1.+2.55*km1)/(1.9*km1)+1.26*km1*alpha_/(1.+3.5*km1);
	  double c1 = 1.+0.3*alpha_;
	  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)
	{
	  double xold=x;
	  x+=newtoncorr (x, N);
	  
	  if (done) break;
	  if (abs(x-xold)<max(1e-12,1e-12*abs(x))) done=true;
	}
      
      nodes[k] = x;
      
      double tmp0 = value(x, N+1);
      weights[k] = x*tmp1/(tmp0*tmp0);
    }//end for

    //verification 
    r_16 sumR = 0.;
    r_16 sumW = 0.;
    r_16 wip;
    for (int i=0;i<N;i++){
      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);
    }
    int alphaFact = 1;
    for(int i=1;i<=alpha_; i++) alphaFact *= i;
    
222 223

#if DEBUG >= 1  
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244
    r_16 sumRTh = (r_16)(N*(N+alpha_));
    r_16 sumWTh = (r_16)(alphaFact);

    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;
    }
245
#endif
246
    
247
  }//nodes_and_weights
248 249 250 251 252

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

253 254 255 256 257 258 259
private:
  enum { large_exponent2=90, minscale=-14, maxscale=14 };
  double ln2, inv_ln2, fsmall, fbig;
  std::vector<double> xi; 
  std::vector<double> cf; 


260 261 262 263 264 265
}; //LaguerreFuncQuad



} // Fin du namespace

266
#undef DEBUG
267 268

#endif //LAGUERREBUILDER_SEEN