Skip to content

GitLab

Commit ac3eff0f authored by Martin Reinecke's avatar Martin Reinecke
Browse files

warning patrol

parent f514ac14
Expand all Hide whitespace changes
Inline Side-by-side
libsharp2/sharp.cc
View file @ ac3eff0f
......@@ -65,18 +65,18 @@ static void get_chunk_info (int ndata, int nmult, int &nchunks, int &chunksize)
nchunks = (ndata+chunksize-1)/chunksize;
}
MRUTIL_NOINLINE int sharp_get_mlim (int lmax, int spin, double sth, double cth)
MRUTIL_NOINLINE size_t sharp_get_mlim (size_t lmax, size_t spin, double sth, double cth)
{
double ofs=lmax*0.01;
if (ofs<100.) ofs=100.;
double b = -2*spin*fabs(cth);
double b = -2*double(spin)*abs(cth);
double t1 = lmax*sth+ofs;
double c = double(spin)*spin-t1*t1;
double discr = b*b-4*c;
if (discr<=0) return lmax;
double res=(-b+sqrt(discr))/2.;
if (res>lmax) res=lmax;
return int(res+0.5);
return size_t(res+0.5);
}
struct ringhelper
......@@ -90,7 +90,7 @@ struct ringhelper
ringhelper() : length(0) {}
void update(int nph, int mmax, double phi0)
{
norot = (fabs(phi0)<1e-14);
norot = (abs(phi0)<1e-14);
if (!(norot))
if ((mmax!=s_shift-1) || (!approx(phi0,phi0_,1e-12)))
{
......@@ -617,7 +617,7 @@ MRUTIL_NOINLINE void sharp_job::execute()
{
size_t llim=chunk*chunksize, ulim=min(llim+chunksize,ginfo->pair.size());
vector<int> ispair(ulim-llim);
vector<int> mlim(ulim-llim);
vector<size_t> mlim(ulim-llim);
vector<double> cth(ulim-llim), sth(ulim-llim);
for (size_t i=0; i<ulim-llim; ++i)
{
......@@ -669,7 +669,7 @@ void sharp_job::build_common (sharp_jobtype type,
if (type==SHARP_Yt) type=SHARP_MAP2ALM;
if (type==SHARP_WY) { type=SHARP_ALM2MAP; flags|=SHARP_USE_WEIGHTS; }
MR_assert((spin>=0)&&(spin<=alm_info.lmax), "bad spin");
MR_assert(spin<=alm_info.lmax, "bad spin");
this->type = type;
this->spin = spin;
nmaps = (type==SHARP_ALM2MAP_DERIV1) ? 2 : ((spin>0) ? 2 : 1);
......
libsharp2/sharp.h
View file @ ac3eff0f
......@@ -178,8 +178,8 @@ void sharp_set_nchunks_max(int new_nchunks_max);
/*! \} */
int sharp_get_mlim (int lmax, int spin, double sth, double cth);
int sharp_veclen(void);
size_t sharp_get_mlim (size_t lmax, size_t spin, double sth, double cth);
size_t sharp_veclen(void);
const char *sharp_architecture(void);
#endif
libsharp2/sharp_core.cc
View file @ ac3eff0f
......@@ -32,10 +32,10 @@
#undef ARCH
using t_inner_loop = void (*) (sharp_job &job, const int *ispair,
const double *cth_, const double *sth_, int llim, int ulim,
sharp_Ylmgen &gen, int mi, const int *mlim);
using t_veclen = int (*) (void);
using t_max_nvec = int (*) (int spin);
const double *cth_, const double *sth_, size_t llim, size_t ulim,
sharp_Ylmgen &gen, size_t mi, const size_t *mlim);
using t_veclen = size_t (*) (void);
using t_max_nvec = size_t (*) (size_t spin);
using t_architecture = const char *(*) (void);
static t_inner_loop inner_loop_ = nullptr;
......@@ -63,10 +63,10 @@ static int XCONCATX2(have,arch)(void) \
} \
\
void XCONCATX2(inner_loop,arch) (sharp_job &job, const int *ispair, \
const double *cth_, const double *sth_, int llim, int ulim, \
sharp_Ylmgen &gen, int mi, const int *mlim); \
int XCONCATX2(sharp_veclen,arch) (void); \
int XCONCATX2(sharp_max_nvec,arch) (int spin); \
const double *cth_, const double *sth_, size_t llim, size_t ulim, \
sharp_Ylmgen &gen, size_t mi, const size_t *mlim); \
size_t XCONCATX2(sharp_veclen,arch) (void); \
size_t XCONCATX2(sharp_max_nvec,arch) (size_t spin); \
const char *XCONCATX2(sharp_architecture,arch) (void);
#if (!defined(__APPLE__))
......@@ -108,14 +108,14 @@ DECL2(avx)
#pragma GCC visibility push(hidden)
void inner_loop (sharp_job &job, const int *ispair,const double *cth,
const double *sth, int llim, int ulim, sharp_Ylmgen &gen, int mi,
const int *mlim)
const double *sth, size_t llim, size_t ulim, sharp_Ylmgen &gen, size_t mi,
const size_t *mlim)
{
if (!inner_loop_) assign_funcs();
inner_loop_(job, ispair, cth, sth, llim, ulim, gen, mi, mlim);
}
int sharp_max_nvec(int spin)
size_t sharp_max_nvec(size_t spin)
{
if (!max_nvec_) assign_funcs();
return max_nvec_(spin);
......@@ -123,7 +123,7 @@ int sharp_max_nvec(int spin)
#pragma GCC visibility pop
int sharp_veclen(void)
size_t sharp_veclen(void)
{
if (!veclen_) assign_funcs();
return veclen_();
......
libsharp2/sharp_core_inc.cc
View file @ ac3eff0f
This diff is collapsed.
libsharp2/sharp_geomhelpers.cc
View file @ ac3eff0f
......@@ -34,22 +34,21 @@
using namespace std;
unique_ptr<sharp_geom_info> sharp_make_subset_healpix_geom_info (int nside, int stride, int nrings,
const int *rings, const double *weight)
unique_ptr<sharp_geom_info> sharp_make_subset_healpix_geom_info (size_t nside, ptrdiff_t stride, size_t nrings,
const size_t *rings, const double *weight)
{
const double pi=3.141592653589793238462643383279502884197;
ptrdiff_t npix=(ptrdiff_t)nside*nside*12;
ptrdiff_t ncap=2*(ptrdiff_t)nside*(nside-1);
ptrdiff_t npix=ptrdiff_t(nside*nside*12);
ptrdiff_t ncap=2*ptrdiff_t(nside*(nside-1));
vector<double> theta(nrings), weight_(nrings), phi0(nrings);
vector<size_t> nph(nrings);
vector<ptrdiff_t> stride_(nrings);
vector<ptrdiff_t> ofs(nrings);
vector<ptrdiff_t> stride_(nrings), ofs(nrings);
ptrdiff_t curofs=0, checkofs; /* checkofs used for assertion introduced when adding rings arg */
for (int m=0; m<nrings; ++m)
for (size_t m=0; m<nrings; ++m)
{
int ring = (rings==nullptr)? (m+1) : rings[m];
ptrdiff_t northring = (ring>2*nside) ? 4*nside-ring : ring;
auto ring = (rings==nullptr)? (m+1) : rings[m];
size_t northring = (ring>2*nside) ? 4*nside-ring : ring;
stride_[m] = stride;
if (northring < nside)
{
......@@ -88,14 +87,14 @@ unique_ptr<sharp_geom_info> sharp_make_subset_healpix_geom_info (int nside, int
return make_unique<sharp_geom_info>(nrings, nph.data(), ofs.data(), stride_.data(), phi0.data(), theta.data(), weight_.data());
}
unique_ptr<sharp_geom_info> sharp_make_weighted_healpix_geom_info (int nside, int stride,
unique_ptr<sharp_geom_info> sharp_make_weighted_healpix_geom_info (size_t nside, ptrdiff_t stride,
const double *weight)
{
return sharp_make_subset_healpix_geom_info(nside, stride, 4 * nside - 1, nullptr, weight);
return sharp_make_subset_healpix_geom_info(nside, stride, 4*nside-1, nullptr, weight);
}
unique_ptr<sharp_geom_info> sharp_make_gauss_geom_info (int nrings, int nphi, double phi0,
int stride_lon, int stride_lat)
unique_ptr<sharp_geom_info> sharp_make_gauss_geom_info (size_t nrings, size_t nphi, double phi0,
ptrdiff_t stride_lon, ptrdiff_t stride_lat)
{
const double pi=3.141592653589793238462643383279502884197;
......@@ -107,12 +106,12 @@ unique_ptr<sharp_geom_info> sharp_make_gauss_geom_info (int nrings, int nphi, do
mr::GL_Integrator integ(nrings);
auto theta = integ.coords();
auto weight = integ.weights();
for (int m=0; m<nrings; ++m)
for (size_t m=0; m<nrings; ++m)
{
theta[m] = acos(-theta[m]);
nph[m]=nphi;
phi0_[m]=phi0;
ofs[m]=(ptrdiff_t)m*stride_lat;
ofs[m]=ptrdiff_t(m*stride_lat);
stride_[m]=stride_lon;
weight[m]*=2*pi/nphi;
}
......@@ -121,18 +120,17 @@ unique_ptr<sharp_geom_info> sharp_make_gauss_geom_info (int nrings, int nphi, do
}
/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
unique_ptr<sharp_geom_info> sharp_make_fejer1_geom_info (int nrings, int ppring, double phi0,
int stride_lon, int stride_lat)
unique_ptr<sharp_geom_info> sharp_make_fejer1_geom_info (size_t nrings, size_t ppring, double phi0,
ptrdiff_t stride_lon, ptrdiff_t stride_lat)
{
const double pi=3.141592653589793238462643383279502884197;
vector<double> theta(nrings), weight(nrings), phi0_(nrings);
vector<size_t> nph(nrings);
vector<ptrdiff_t> stride_(nrings);
vector<ptrdiff_t> ofs(nrings);
vector<ptrdiff_t> stride_(nrings), ofs(nrings);
weight[0]=2.;
for (int k=1; k<=(nrings-1)/2; ++k)
for (size_t k=1; k<=(nrings-1)/2; ++k)
{
weight[2*k-1]=2./(1.-4.*k*k)*cos((k*pi)/nrings);
weight[2*k ]=2./(1.-4.*k*k)*sin((k*pi)/nrings);
......@@ -140,14 +138,14 @@ unique_ptr<sharp_geom_info> sharp_make_fejer1_geom_info (int nrings, int ppring,
if ((nrings&1)==0) weight[nrings-1]=0.;
mr::r2r_fftpack({size_t(nrings)}, {sizeof(double)}, {sizeof(double)}, {0}, false, false, weight.data(), weight.data(), 1.);
for (int m=0; m<(nrings+1)/2; ++m)
for (size_t m=0; m<(nrings+1)/2; ++m)
{
theta[m]=pi*(m+0.5)/nrings;
theta[nrings-1-m]=pi-theta[m];
nph[m]=nph[nrings-1-m]=ppring;
phi0_[m]=phi0_[nrings-1-m]=phi0;
ofs[m]=(ptrdiff_t)m*stride_lat;
ofs[nrings-1-m]=(ptrdiff_t)((nrings-1-m)*stride_lat);
ofs[m]=ptrdiff_t(m*stride_lat);
ofs[nrings-1-m]=ptrdiff_t((nrings-1-m)*stride_lat);
stride_[m]=stride_[nrings-1-m]=stride_lon;
weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(nrings*nph[m]);
}
......@@ -156,34 +154,33 @@ unique_ptr<sharp_geom_info> sharp_make_fejer1_geom_info (int nrings, int ppring,
}
/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
unique_ptr<sharp_geom_info> sharp_make_cc_geom_info (int nrings, int ppring, double phi0,
int stride_lon, int stride_lat)
unique_ptr<sharp_geom_info> sharp_make_cc_geom_info (size_t nrings, size_t ppring, double phi0,
ptrdiff_t stride_lon, ptrdiff_t stride_lat)
{
const double pi=3.141592653589793238462643383279502884197;
vector<double> theta(nrings), weight(nrings,0.), phi0_(nrings);
vector<size_t> nph(nrings);
vector<ptrdiff_t> stride_(nrings);
vector<ptrdiff_t> ofs(nrings);
vector<ptrdiff_t> stride_(nrings), ofs(nrings);
int n=nrings-1;
size_t n=nrings-1;
double dw=-1./(n*n-1.+(n&1));
weight[0]=2.+dw;
for (int k=1; k<=(n/2-1); ++k)
for (size_t k=1; k<=(n/2-1); ++k)
weight[2*k-1]=2./(1.-4.*k*k) + dw;
weight[2*(n/2)-1]=(n-3.)/(2*(n/2)-1) -1. -dw*((2-(n&1))*n-1);
mr::r2r_fftpack({size_t(n)}, {sizeof(double)}, {sizeof(double)}, {0}, false, false, weight.data(), weight.data(), 1.);
weight[n]=weight[0];
for (int m=0; m<(nrings+1)/2; ++m)
for (size_t m=0; m<(nrings+1)/2; ++m)
{
theta[m]=pi*m/(nrings-1.);
if (theta[m]<1e-15) theta[m]=1e-15;
theta[nrings-1-m]=pi-theta[m];
nph[m]=nph[nrings-1-m]=ppring;
phi0_[m]=phi0_[nrings-1-m]=phi0;
ofs[m]=(ptrdiff_t)m*stride_lat;
ofs[nrings-1-m]=(ptrdiff_t)((nrings-1-m)*stride_lat);
ofs[m]=ptrdiff_t(m*stride_lat);
ofs[nrings-1-m]=ptrdiff_t((nrings-1-m)*stride_lat);
stride_[m]=stride_[nrings-1-m]=stride_lon;
weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(n*nph[m]);
}
......@@ -192,33 +189,32 @@ unique_ptr<sharp_geom_info> sharp_make_cc_geom_info (int nrings, int ppring, dou
}
/* Weights from Waldvogel 2006: BIT Numerical Mathematics 46, p. 195 */
unique_ptr<sharp_geom_info> sharp_make_fejer2_geom_info (int nrings, int ppring, double phi0,
int stride_lon, int stride_lat)
unique_ptr<sharp_geom_info> sharp_make_fejer2_geom_info (size_t nrings, size_t ppring, double phi0,
ptrdiff_t stride_lon, ptrdiff_t stride_lat)
{
const double pi=3.141592653589793238462643383279502884197;
vector<double> theta(nrings), weight(nrings+1, 0.), phi0_(nrings);
vector<size_t> nph(nrings);
vector<ptrdiff_t> stride_(nrings);
vector<ptrdiff_t> ofs(nrings);
vector<ptrdiff_t> stride_(nrings), ofs(nrings);
int n=nrings+1;
size_t n=nrings+1;
weight[0]=2.;
for (int k=1; k<=(n/2-1); ++k)
for (size_t k=1; k<=(n/2-1); ++k)
weight[2*k-1]=2./(1.-4.*k*k);
weight[2*(n/2)-1]=(n-3.)/(2*(n/2)-1) -1.;
mr::r2r_fftpack({size_t(n)}, {sizeof(double)}, {sizeof(double)}, {0}, false, false, weight.data(), weight.data(), 1.);
for (int m=0; m<nrings; ++m)
for (size_t m=0; m<nrings; ++m)
weight[m]=weight[m+1];
for (int m=0; m<(nrings+1)/2; ++m)
for (size_t m=0; m<(nrings+1)/2; ++m)
{
theta[m]=pi*(m+1)/(nrings+1.);
theta[nrings-1-m]=pi-theta[m];
nph[m]=nph[nrings-1-m]=ppring;
phi0_[m]=phi0_[nrings-1-m]=phi0;
ofs[m]=(ptrdiff_t)m*stride_lat;
ofs[nrings-1-m]=(ptrdiff_t)((nrings-1-m)*stride_lat);
ofs[m]=ptrdiff_t(m*stride_lat);
ofs[nrings-1-m]=ptrdiff_t((nrings-1-m)*stride_lat);
stride_[m]=stride_[nrings-1-m]=stride_lon;
weight[m]=weight[nrings-1-m]=weight[m]*2*pi/(n*nph[m]);
}
......@@ -226,23 +222,22 @@ unique_ptr<sharp_geom_info> sharp_make_fejer2_geom_info (int nrings, int ppring,
return make_unique<sharp_geom_info> (nrings, nph.data(), ofs.data(), stride_.data(), phi0_.data(), theta.data(), weight.data());
}
unique_ptr<sharp_geom_info> sharp_make_mw_geom_info (int nrings, int ppring, double phi0,
int stride_lon, int stride_lat)
unique_ptr<sharp_geom_info> sharp_make_mw_geom_info (size_t nrings, size_t ppring, double phi0,
ptrdiff_t stride_lon, ptrdiff_t stride_lat)
{
const double pi=3.141592653589793238462643383279502884197;
vector<double> theta(nrings), phi0_(nrings);
vector<size_t> nph(nrings);
vector<ptrdiff_t> stride_(nrings);
vector<ptrdiff_t> ofs(nrings);
vector<ptrdiff_t> stride_(nrings), ofs(nrings);
for (int m=0; m<nrings; ++m)
for (size_t m=0; m<nrings; ++m)
{
theta[m]=pi*(2.*m+1.)/(2.*nrings-1.);
if (theta[m]>pi-1e-15) theta[m]=pi-1e-15;
nph[m]=ppring;
phi0_[m]=phi0;
ofs[m]=(ptrdiff_t)m*stride_lat;
ofs[m]=ptrdiff_t(m*stride_lat);
stride_[m]=stride_lon;
}
......
libsharp2/sharp_geomhelpers.h
View file @ ac3eff0f
......@@ -40,21 +40,21 @@
\note if \a weight is a null pointer, all weights are assumed to be 1.
\note if \a rings is a null pointer, take all rings
\ingroup geominfogroup */
std::unique_ptr<sharp_geom_info> sharp_make_subset_healpix_geom_info (int nside, int stride, int nrings,
const int *rings, const double *weight);
std::unique_ptr<sharp_geom_info> sharp_make_subset_healpix_geom_info (size_t nside, ptrdiff_t stride, size_t nrings,
const size_t *rings, const double *weight);
/*! Creates a geometry information describing a HEALPix map with an
Nside parameter \a nside. \a weight contains the relative ring
weights and must have \a 2*nside entries.
\note if \a weight is a null pointer, all weights are assumed to be 1.
\ingroup geominfogroup */
std::unique_ptr<sharp_geom_info> sharp_make_weighted_healpix_geom_info (int nside, int stride,
std::unique_ptr<sharp_geom_info> sharp_make_weighted_healpix_geom_info (size_t nside, ptrdiff_t stride,
const double *weight);
/*! Creates a geometry information describing a HEALPix map with an
Nside parameter \a nside.
\ingroup geominfogroup */
static inline std::unique_ptr<sharp_geom_info> sharp_make_healpix_geom_info (int nside, int stride)
static inline std::unique_ptr<sharp_geom_info> sharp_make_healpix_geom_info (size_t nside, ptrdiff_t stride)
{ return sharp_make_weighted_healpix_geom_info (nside, stride, nullptr); }
/*! Creates a geometry information describing a Gaussian map with \a nrings
......@@ -64,8 +64,8 @@ static inline std::unique_ptr<sharp_geom_info> sharp_make_healpix_geom_info (int
difference between the two start pixels in consecutive iso-latitude rings
is \a stride_lat.
\ingroup geominfogroup */
std::unique_ptr<sharp_geom_info> sharp_make_gauss_geom_info (int nrings, int nphi, double phi0,
int stride_lon, int stride_lat);
std::unique_ptr<sharp_geom_info> sharp_make_gauss_geom_info (size_t nrings, size_t nphi, double phi0,
ptrdiff_t stride_lon, ptrdiff_t stride_lat);
/*! Creates a geometry information describing an ECP map with \a nrings
iso-latitude rings and \a nphi pixels per ring. The azimuth of the first
......@@ -80,13 +80,13 @@ std::unique_ptr<sharp_geom_info> sharp_make_gauss_geom_info (int nrings, int nph
\a pi-0.5*(pi/nrings). There are no pixel centers at the poles.
\note This grid corresponds to Fejer's first rule.
\ingroup geominfogroup */
std::unique_ptr<sharp_geom_info> sharp_make_fejer1_geom_info (int nrings, int nphi, double phi0,
int stride_lon, int stride_lat);
std::unique_ptr<sharp_geom_info> sharp_make_fejer1_geom_info (size_t nrings, size_t nphi, double phi0,
ptrdiff_t stride_lon, ptrdiff_t stride_lat);
/*! Old name for sharp_make_fejer1_geom_info()
\ingroup geominfogroup */
static inline std::unique_ptr<sharp_geom_info> sharp_make_ecp_geom_info (int nrings, int nphi, double phi0,
int stride_lon, int stride_lat)
static inline std::unique_ptr<sharp_geom_info> sharp_make_ecp_geom_info (size_t nrings, size_t nphi, double phi0,
ptrdiff_t stride_lon, ptrdiff_t stride_lat)
{
return sharp_make_fejer1_geom_info (nrings, nphi, phi0, stride_lon, stride_lat);
}
......@@ -103,8 +103,8 @@ static inline std::unique_ptr<sharp_geom_info> sharp_make_ecp_geom_info (int nri
is \a 0 and that of the last ring is \a pi.
\note This grid corresponds to Clenshaw-Curtis integration.
\ingroup geominfogroup */
std::unique_ptr<sharp_geom_info> sharp_make_cc_geom_info (int nrings, int ppring, double phi0,
int stride_lon, int stride_lat);
std::unique_ptr<sharp_geom_info> sharp_make_cc_geom_info (size_t nrings, size_t ppring, double phi0,
ptrdiff_t stride_lon, ptrdiff_t stride_lat);
/*! Creates a geometry information describing an ECP map with \a nrings
iso-latitude rings and \a nphi pixels per ring. The azimuth of the first
......@@ -118,8 +118,8 @@ std::unique_ptr<sharp_geom_info> sharp_make_cc_geom_info (int nrings, int ppring
is \a pi/(nrings+1) and that of the last ring is \a pi-pi/(nrings+1).
\note This grid corresponds to Fejer's second rule.
\ingroup geominfogroup */
std::unique_ptr<sharp_geom_info> sharp_make_fejer2_geom_info (int nrings, int ppring, double phi0,
int stride_lon, int stride_lat);
std::unique_ptr<sharp_geom_info> sharp_make_fejer2_geom_info (size_t nrings, size_t ppring, double phi0,
ptrdiff_t stride_lon, ptrdiff_t stride_lat);
/*! Creates a geometry information describing a map with \a nrings
iso-latitude rings and \a nphi pixels per ring. The azimuth of the first
......@@ -134,7 +134,7 @@ std::unique_ptr<sharp_geom_info> sharp_make_fejer2_geom_info (int nrings, int pp
\note This is the grid introduced by McEwen & Wiaux 2011.
\note This function does \e not define any quadrature weights.
\ingroup geominfogroup */
std::unique_ptr<sharp_geom_info> sharp_make_mw_geom_info (int nrings, int ppring, double phi0,
int stride_lon, int stride_lat);
std::unique_ptr<sharp_geom_info> sharp_make_mw_geom_info (size_t nrings, size_t ppring, double phi0,
ptrdiff_t stride_lon, ptrdiff_t stride_lat);
#endif
libsharp2/sharp_internal.h
View file @ ac3eff0f
......@@ -67,9 +67,9 @@ struct sharp_job
};
void inner_loop (sharp_job &job, const int *ispair,const double *cth,
const double *sth, int llim, int ulim, sharp_Ylmgen &gen, int mi,
const int *mlim);
const double *sth, size_t llim, size_t ulim, sharp_Ylmgen &gen, size_t mi,
const size_t *mlim);
int sharp_max_nvec(int spin);
size_t sharp_max_nvec(size_t spin);
#endif
libsharp2/sharp_ylmgen.cc
View file @ ac3eff0f
......@@ -42,13 +42,12 @@ static inline void normalize (double &val, int &scale, double xfmax)
while (abs(val)<xfmax*sharp_fsmall) { val*=sharp_fbig; --scale; }
}
sharp_Ylmgen::sharp_Ylmgen (int l_max, int m_max, int spin)
sharp_Ylmgen::sharp_Ylmgen (size_t l_max, size_t m_max, size_t spin)
{
constexpr double inv_sqrt4pi = 0.2820947917738781434740397257803862929220;
lmax = l_max;
mmax = m_max;
MR_assert(spin>=0,"incorrect spin: must be nonnegative");
MR_assert(l_max>=spin,"incorrect l_max: must be >= spin");
MR_assert(l_max>=m_max,"incorrect l_max: must be >= m_max");
s = spin;
......@@ -56,30 +55,30 @@ sharp_Ylmgen::sharp_Ylmgen (int l_max, int m_max, int spin)
"bad value for min/maxscale");
cf.resize(sharp_maxscale-sharp_minscale+1);
cf[-sharp_minscale]=1.;
for (int m=-sharp_minscale-1; m>=0; --m)
cf[m]=cf[m+1]*sharp_fsmall;
for (int m=-sharp_minscale+1; m<(sharp_maxscale-sharp_minscale+1); ++m)
cf[m]=cf[m-1]*sharp_fbig;
for (int sc=-sharp_minscale-1; sc>=0; --sc)
cf[sc]=cf[sc+1]*sharp_fsmall;
for (int sc=-sharp_minscale+1; sc<(sharp_maxscale-sharp_minscale+1); ++sc)
cf[sc]=cf[sc-1]*sharp_fbig;
powlimit.resize(m_max+spin+1);
powlimit[0]=0.;
constexpr double ln2 = 0.6931471805599453094172321214581766;
constexpr double expo=-400*ln2;
for (int m=1; m<=m_max+spin; ++m)
powlimit[m]=exp(expo/m);
for (size_t i=1; i<=m_max+spin; ++i)
powlimit[i]=exp(expo/i);
m = -1;
m = ~size_t(0);
if (spin==0)
{
mfac.resize(mmax+1);
mfac[0] = inv_sqrt4pi;
for (int m=1; m<=mmax; ++m)
mfac[m] = mfac[m-1]*sqrt((2*m+1.)/(2*m));
for (size_t i=1; i<=mmax; ++i)
mfac[i] = mfac[i-1]*sqrt((2*i+1.)/(2*i));
root.resize(2*lmax+8);
iroot.resize(2*lmax+8);
for (int m=0; m<2*lmax+8; ++m)
for (size_t i=0; i<2*lmax+8; ++i)
{
root[m] = sqrt(m);
iroot[m] = (m==0) ? 0. : 1./root[m];
root[i] = sqrt(i);
iroot[i] = (i==0) ? 0. : 1./root[i];
}
eps.resize(lmax+4);
alpha.resize(lmax/2+2);
......@@ -87,35 +86,35 @@ sharp_Ylmgen::sharp_Ylmgen (int l_max, int m_max, int spin)
}
else
{
m=mlo=mhi=-1234567890;
m=mlo=mhi=~size_t(0);
coef.resize(lmax+3);
for (int m=0; m<lmax+3; ++m)
coef[m].a=coef[m].b=0.;
for (size_t i=0; i<lmax+3; ++i)
coef[i].a=coef[i].b=0.;
alpha.resize(lmax+3);
inv.resize(lmax+2);
inv[0]=0;
for (int m=1; m<lmax+2; ++m) inv[m]=1./m;
for (size_t i=1; i<lmax+2; ++i) inv[i]=1./i;
flm1.resize(2*lmax+3);
flm2.resize(2*lmax+3);
for (int m=0; m<2*lmax+3; ++m)
for (size_t i=0; i<2*lmax+3; ++i)
{
flm1[m] = sqrt(1./(m+1.));
flm2[m] = sqrt(m/(m+1.));
flm1[i] = sqrt(1./(i+1.));
flm2[i] = sqrt(i/(i+1.));
}
prefac.resize(mmax+1);
fscale.resize(mmax+1);
vector<double> fac(2*lmax+1);
vector<int> facscale(2*lmax+1);
fac[0]=1; facscale[0]=0;
for (int m=1; m<2*lmax+1; ++m)
for (size_t i=1; i<2*lmax+1; ++i)
{
fac[m]=fac[m-1]*sqrt(m);
facscale[m]=facscale[m-1];
normalize(fac[m],facscale[m],sharp_fbighalf);
fac[i]=fac[i-1]*sqrt(i);
facscale[i]=facscale[i-1];
normalize(fac[i],facscale[i],sharp_fbighalf);
}
for (int m=0; m<=mmax; ++m)
for (size_t i=0; i<=mmax; ++i)
{
int mlo_=s, mhi_=m;
int mlo_=s, mhi_=i;
if (mhi_<mlo_) swap(mhi_,mlo_);
double tfac=fac[2*mhi_]/fac[mhi_+mlo_];
int tscale=facscale[2*mhi_]-facscale[mhi_+mlo_];
......@@ -123,28 +122,27 @@ sharp_Ylmgen::sharp_Ylmgen (int l_max, int m_max, int spin)
tfac/=fac[mhi_-mlo_];
tscale-=facscale[mhi_-mlo_];
normalize(tfac,tscale,sharp_fbighalf);
prefac[m]=tfac;
fscale[m]=tscale;
prefac[i]=tfac;
fscale[i]=tscale;
}
}
}
void sharp_Ylmgen::prepare (int m)
void sharp_Ylmgen::prepare (size_t m_)
{
if (m==this->m) return;
MR_assert(m>=0,"incorrect m");
this->m = m;
if (m_==m) return;
m = m_;
if (s==0)
{
eps[m] = 0.;
for (int l=m+1; l<lmax+4; ++l)
for (size_t l=m+1; l<lmax+4; ++l)
eps[l] = root[l+m]*root[l-m]*iroot[2*l+1]*iroot[2*l-1];
alpha[0] = 1./eps[m+1];
alpha[1] = eps[m+1]/(eps[m+2]*eps[m+3]);
for (int il=1, l=m+2; l<lmax+1; ++il, l+=2)
for (size_t il=1, l=m+2; l<lmax+1; ++il, l+=2)
alpha[il+1]= ((il&1) ? -1 : 1) / (eps[l+2]*eps[l+3]*alpha[il]);