diff --git a/libsharp/sharp_core_inc.c b/libsharp/sharp_core_inc.c index 2c86463..fedf506 100644 --- a/libsharp/sharp_core_inc.c +++ b/libsharp/sharp_core_inc.c @@ -109,7 +109,7 @@ static inline int Y(rescale) (Tb * restrict lam1, Tb * restrict lam2, int did_scale=0; for (int i=0;iv[i]),vone); + Tv mask = vgt(vabs(lam2->v[i]),vload(sharp_ftol)); if (vanyTrue(mask)) { did_scale=1; @@ -126,20 +126,20 @@ static inline void Y(normalize) (Tb * restrict val, Tb * restrict scale) const Tv vfsmall=vload(sharp_fsmall), vfbig=vload(sharp_fbig); for (int i=0;iv[i]),vone); + Tv mask = vgt(vabs(val->v[i]),vload(sharp_ftol)); while (vanyTrue(mask)) { vmuleq(val->v[i],vblend(mask,vfsmall,vone)); vaddeq(scale->v[i],vblend(mask,vone,vzero)); - mask = vgt(vabs(val->v[i]),vone); + mask = vgt(vabs(val->v[i]),vload(sharp_ftol)); } - mask = vlt(vabs(val->v[i]),vfsmall); + mask = vlt(vabs(val->v[i]),vload(sharp_fsmall*sharp_ftol)); mask = vand(mask,vne(val->v[i],vzero)); while (vanyTrue(mask)) { vmuleq(val->v[i],vblend(mask,vfbig,vone)); vsubeq(scale->v[i],vblend(mask,vone,vzero)); - mask = vlt(vabs(val->v[i]),vfsmall); + mask = vlt(vabs(val->v[i]),vload(sharp_fsmall*sharp_ftol)); mask = vand(mask,vne(val->v[i],vzero)); } } diff --git a/libsharp/sharp_ylmgen_c.h b/libsharp/sharp_ylmgen_c.h index 144bc96..e935ce2 100644 --- a/libsharp/sharp_ylmgen_c.h +++ b/libsharp/sharp_ylmgen_c.h @@ -36,8 +36,9 @@ extern "C" { #endif -enum { sharp_minscale=-10, sharp_limscale=0, sharp_maxscale=2 }; -static const double sharp_fbig=0x1p+90,sharp_fsmall=0x1p-90; +enum { sharp_minscale=-1, sharp_limscale=1, sharp_maxscale=1 }; +static const double sharp_fbig=0x1p+450,sharp_fsmall=0x1p-450; +static const double sharp_ftol=0x1p-60; typedef struct { double f[2]; } sharp_ylmgen_dbl2; typedef struct { double f[3]; } sharp_ylmgen_dbl3;