Commit 6a4661dc authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

(JEC) 20/9/16 use new implementation of the spherical bessel functions. Gain a...

(JEC) 20/9/16 use new implementation of the spherical bessel functions. Gain a factor 10 in speed. Could be less precise than the double precision default Boost implementation.
parent 0f185319
......@@ -3,7 +3,20 @@
#include "laguerreTransform.h"
namespace LagSHT {
//singleton init
BesselJImp* BesselJImp::instance = 0;
BesselJImp* BesselJImp::getInstance() {
if(!instance) {
instance = new BesselJImp();
return instance;
}
else {
return instance;
}
}//getInstance
void Bessel::BesselRoots() {
using namespace std;
......
......@@ -42,6 +42,220 @@ namespace LagSHT {
//// --------- Class Spherical Bessel utils -- //
////////////////////////////////////////////////////////////////
/* JEC 20/9/16
POLICY = 2 => Numerical Recipies-like function
= 1 => BOOST + policy reduite
= 0 => BOOST call avec r_4 arguments
*/
#define _GAMMA1_ 2.6789385347077476336556
#define _GAMMA2_ 1.3541179394264004169452
#define POLICY 2
// j_l(x) hub
class BesselJImp {
public:
//the unique ptr
static BesselJImp* getInstance();
//the code
inline virtual r_8 operator()(int el, r_8 x) const{
r_8 val=0;
#if POLICY == 2
bessel_j(el,x,&val);
#elif POLICY == 1
using namespace boost::math;
using namespace boost::math::policies;
typedef policy<digits10<5>, promote_float<false>, promote_double<false>> pol;
val = sph_bessel(el, (r_4)(x), pol()); //JEc 19mai16 flaot arg/return
#else
using namespace boost::math;
val = sph_bessel(el, (r_4)(x));
#endif
return val;
}
private:
//Forbid Ctor
BesselJImp() {}
//the ptr
static BesselJImp* instance;
//The local implementation
inline int bessel_j(
int l,
double x,
double * jl
) const {
double nu,nu2,beta,beta2;
double x2,sx,sx2,cx;
double cotb,cot3b,cot6b,secb,sec2b;
double trigarg,expterm,fl;
double l3,cosb;
fl = (double)l;
x2 = x*x;
/************* Use closed form for l<7 **********/
if (l < 7) {
sx=sin(x);
cx=cos(x);
if(l == 0) {
if (x > 0.1) *jl=sx/x;
else *jl=1.-x2/6.*(1.-x2/20.);
return 0;
}
if(l == 1) {
if (x > 0.2) *jl=(sx/x -cx)/x;
else *jl=x/3.*(1.-x2/10.*(1.-x2/28.));
return 0;
}
if (l == 2) {
if (x > 0.3) *jl=(-3.*cx/x-sx*(1.-3./x2))/x;
else *jl=x2/15.*(1.-x2/14.*(1.-x2/36.));
return 0;
}
if (l == 3) {
if (x > 0.4) *jl=(cx*(1.-15./x2)-sx*(6.-15./x2)/x)/x;
else *jl=x*x2/105.*(1.-x2/18.*(1.-x2/44.));
return 0;
}
if (l == 4) {
if (x > 0.6) *jl=(sx*(1.-45./x2+105./x2/x2) +cx*(10.-105./x2)/x)/x;
else *jl=x2*x2/945.*(1.-x2/22.*(1.-x2/52.));
return 0;
}
if (l == 5) {
if (x > 1) *jl=(sx*(15.-420./x2+945./x2/x2)/x -cx*(1.0-105./x2+945./x2/x2))/x;
else *jl=x2*x2*x/10395.*(1.-x2/26.*(1.-x2/60.));
return 0;
}
if (l == 6) {
if (x > 1) *jl=(sx*(-1.+(210.-(4725.-10395./x2)/x2)/x2)+
cx*(-21.+(1260.-10395./x2)/x2)/x)/x;
else *jl=x2*x2*x2/135135.*(1.-x2/30.*(1.-x2/68.));
return 0;
}
}
else {
if (x <= 1.e-40) {
*jl=0.0;
return 0;
}
nu= fl + 0.5;
nu2=nu*nu;
if ((x2/fl) < 0.5) {
*jl=exp(fl*log(x/nu/2.)+nu*(1-log(2.))-(1.-(1.-3.5/nu2)/nu2/30.)/12./nu)
/nu*(1.-x2/(4.*nu+4.)*(1.-x2/(8.*nu+16.)*(1.-x2/(12.*nu+36.))));
return 0;
}
if ((fl*fl/x) < 0.5) {
beta = x - M_PI/2.*(fl+1.);
*jl = (cos(beta)*(1.-(nu2-0.25)*(nu2-2.25)/8./x2*(1.-(nu2-6.25)*(nu2-12.25)/48./x2))
-sin(beta)*(nu2-0.25)/2./x* (1.-(nu2-2.25)*(nu2-6.25)/24./x2*(1.-(nu2-12.25)*(nu2-20.25)/80./x2)) )/x;
return 0;
}
l3 = pow(nu,0.325);
if (x < nu-1.31*l3) {
cosb=nu/x;
sx=sqrt(nu2-x2);
cotb=nu/sx;
secb=x/nu;
beta=log(cosb+sx/x);
cot3b=cotb*cotb*cotb;
cot6b=cot3b*cot3b;
sec2b=secb*secb;
expterm=((2.+3.*sec2b)*cot3b/24.
- ((4.+sec2b)*sec2b*cot6b/16.
+ ((16.-(1512.+(3654.+375.*sec2b)*sec2b)*sec2b)*cot3b/5760.
+ (32.+(288.+(232.+13.*sec2b)*sec2b)*sec2b)*sec2b*cot6b/128./nu)*cot6b/nu)/nu)/nu;
*jl=sqrt(cotb*cosb)/(2.*nu)*exp(-nu*beta+nu/cotb-expterm);
return 0;
}
if (x > nu+1.48*l3) {
cosb=nu/x;
sx=sqrt(x2-nu2);
cotb=nu/sx;
secb=x/nu;
beta=acos(cosb);
cot3b=cotb*cotb*cotb;
cot6b=cot3b*cot3b;
sec2b=secb*secb;
trigarg=nu/cotb-nu*beta-M_PI/4.
-((2.+3.*sec2b)*cot3b/24.
+(16.-(1512.+(3654.+375.*sec2b)*sec2b)*sec2b)*cot3b*cot6b/5760./nu2)/nu;
expterm=((4.+sec2b)*sec2b*cot6b/16.
-(32.+(288.+(232.+13.*sec2b)*sec2b)*sec2b)*sec2b*cot6b*cot6b/128./nu2)/nu2;
*jl=sqrt(cotb*cosb)/nu*exp(-expterm)*cos(trigarg);
return 0;
}
/* last possible case */
beta=x-nu;
beta2=beta*beta;
sx=6./x;
sx2=sx*sx;
secb=pow(sx,1./3.);
sec2b=secb*secb;
*jl=(_GAMMA1_*secb + beta*_GAMMA2_*sec2b
-(beta2/18.-1./45.)*beta*sx*secb*_GAMMA1_
-((beta2-1.)*beta2/36.+1./420.)*sx*sec2b*_GAMMA2_
+(((beta2/1620.-7./3240.)*beta2+1./648.)*beta2-1./8100.)*sx2*secb*_GAMMA1_
+(((beta2/4536.-1./810.)*beta2+19./11340.)*beta2-13./28350.)*beta*sx2*sec2b*_GAMMA2_
-((((beta2/349920.-1./29160.)*beta2+71./583200.)*beta2-121./874800.)*
beta2+7939./224532000.)*beta*sx2*sx*secb*_GAMMA1_)*sqrt(sx)/12./sqrt(M_PI);
return 0;
}
std::cout << "l or x out of bounds: " << l << ", " << x << std::endl;
return -1;
}//bessel_j
};//class
class Bessel {
......@@ -71,7 +285,9 @@ class Bessel {
// Sperical Bessel functions
static inline r_8 jn(int n, r_8 x) {
return boost::math::sph_bessel(n, x);
// return boost::math::sph_bessel(n, x);
BesselJImp* jn_imp = BesselJImp::getInstance();
return (*jn_imp)(n,x); //JEC 20/9/16 delegate the implementation
}
//Helper
......@@ -107,7 +323,15 @@ private:
Explicit prohibit Assignment Operator
*/
Bessel& operator=(const Bessel&);
}; //class
#undef _GAMMA1_
#undef _GAMMA2_
}//end namespace
#endif //LAGSHTBESSEL_SEEN
......@@ -68,8 +68,7 @@ Laguerre2Bessel::Laguerre2Bessel(BaseGeometry* sphere,
cout << "READ Jlnp......." << endl;
//read the Rmax-independant values of the jlnp and rescale
// stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
std::stringstream ss; ss << "limber0-jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
std::stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
std::string fname(jlnpDir+"/"+ss.str());
std::ifstream ifs(fname.c_str(), std::ifstream::in);
......@@ -161,7 +160,9 @@ void Laguerre2Bessel::ComputeJInteg(std::string clenshawDir, std::string jlnpDir
cout << "(JEC) Plan Inverse Init START" << endl;
//Unique Product plan
int nOrdProd= 2*nOrdGen+1; //order max of product is 2*nOrdGen but to get a power of 2 for the FFT I add 1.
// int nOrdProd= 2*nOrdGen+1; //order max of product is 2*nOrdGen but to get a power of 2 for the FFT I add 1.
int nOrdProd= 2*nOrdGen-1; //JEC 20/9/16 this is ok with -1
vector<r_8> VecDCT1Inv(nOrdProd+1,0.);
FFTPlanning planInv(nOrdProd,VecDCT1Inv);
......@@ -188,7 +189,7 @@ void Laguerre2Bessel::ComputeJInteg(std::string clenshawDir, std::string jlnpDir
tstack_pop("Common CC weights init...");
stringstream ss; ss << "v230316-jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
stringstream ss; ss << "jlnp-L"<<Lmax_<<"-N"<<Nmax_<<"-P"<<Pmax_<<".txt";
std::ofstream ofs;
string fname(jlnpDir+"/"+ss.str());
ofs.open (fname.c_str(), std::ofstream::out);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment