From 4b7797edafeec0a9f9540ba32cd8702913941c0a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 7 Mar 2018 21:19:21 +0200 Subject: [PATCH] "fix" the q ranges for the B coefficients in translations.c Former-commit-id: 17084dd0b57cb4271c67d5d95b7e033c3968276f --- qpms/test_translations_xu_multipliers.c | 17 +++-- qpms/translations.c | 84 ++++++++++++++++--------- 2 files changed, 66 insertions(+), 35 deletions(-) diff --git a/qpms/test_translations_xu_multipliers.c b/qpms/test_translations_xu_multipliers.c index 43e678b..2aaee54 100644 --- a/qpms/test_translations_xu_multipliers.c +++ b/qpms/test_translations_xu_multipliers.c @@ -4,6 +4,11 @@ //#include #include +#ifdef QPMS_PACKED_B_MULTIPLIERS +#define BQ_OFFSET 1 +#else +#define BQ_OFFSET 0 +#endif void qpms_trans_calculator_multipliers_B_general( qpms_normalisation_t norm, @@ -13,6 +18,7 @@ int lMax=13; #define MIN(x,y) ((x)<(y)?(x):(y)) + // Python test: Qmax(M, n, mu, nu) = floor(min(n,nu,(n+nu+1-abs(M+mu))/2)) // q in IntegerRange(1, Qmax(-m,n,mu,nu)) @@ -24,15 +30,14 @@ int main() { for(int nu = 1; nu <= lMax; ++nu) for(int m = -n; m <= n; ++m) for(int mu = -nu; mu <= nu; ++mu){ - int Qmax = gaunt_q_max(-m, n+1, mu, nu); - int Qmax_alt = MIN(n,MIN(nu,(n+nu+1-abs(mu-m)))); + int Qmax_alt = MIN(n,MIN(nu,(n+nu+1-abs(mu-m))/2)); qpms_trans_calculator_multipliers_B_general(QPMS_NORMALISATION_XU_CS, - dest, m, n, mu, nu, Qmax); - for(int q = 0; q <= Qmax; ++q) { + dest, m, n, mu, nu, Qmax_alt); + for(int q = BQ_OFFSET; q <= Qmax_alt; ++q) { // int p = n + nu - 2*q; - int tubig = cabs(dest[q]) > 1e-8; + int tubig = cabs(dest[q-BQ_OFFSET]) > 1e-8; printf("%.16g + %.16g*I, // %d, %d, %d, %d, %d,%s\n", - creal(dest[q]), cimag(dest[q]), m, n, mu, nu, q, + creal(dest[q-BQ_OFFSET]), cimag(dest[q-BQ_OFFSET]), m, n, mu, nu, q, q > Qmax_alt ? (tubig?" //tubig":" //tu") : ""); } fflush(stdout); diff --git a/qpms/translations.c b/qpms/translations.c index 442ce8e..bafdbb7 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -9,9 +9,19 @@ #include "assert_cython_workaround.h" #include //abort() +// if defined, the pointer B_multipliers[y] corresponds to the q = 1 element; +// otherwise, it corresponds to the q = 0 element, which should be identically zero +#ifdef QPMS_PACKED_B_MULTIPLIERS +#define BQ_OFFSET 1 +#else +#define BQ_OFFSET 0 +#endif + + /* * References: - * [1] Yu-Lin Xu, Journal of Computational Physics 127, 285–298 (1996) + * [Xu_old] Yu-Lin Xu, Journal of Computational Physics 127, 285–298 (1996) + * [Xu] Yu-Lin Xu, Journal of Computational Physics 139, 137–165 (1998) */ /* @@ -40,6 +50,19 @@ double qpms_legendreD0(int m, int n) { return -2 * qpms_legendre0(m, n); } + +static inline int imin(int x, int y) { + return x > y ? y : x; +} + +// The uppermost value of q index for the B coefficient terms from [Xu](60). +// N.B. this is different from [Xu_old](79) due to the n vs. n+1 difference. +// However, the trailing terms in [Xu_old] are analytically zero (although +// the numerical values will carry some non-zero rounding error). +static inline int gauntB_Q_max(int M, int n, int mu, int nu) { + return imin(n, imin(nu, (n+nu+1-abs(M+mu))/2)); +} + int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double *result_array) { int retval; double tmparr[lmax+1]; @@ -202,14 +225,16 @@ complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kd return (presum / prenormratio) * sum; } -// [1], eq. (83) +// [Xu_old], eq. (83) complex double qpms_trans_single_B_Xu(int m, int n, int mu, int nu, sph_t kdlj, bool r_ge_d, qpms_bessel_t J) { if(r_ge_d) J = QPMS_BESSEL_REGULAR; double costheta = cos(kdlj.theta); + // TODO Qmax cleanup: can I replace Qmax with realQmax??? int q2max = gaunt_q_max(-m-1,n+1,mu+1,nu); int Qmax = gaunt_q_max(-m,n+1,mu,nu); + int realQmax = gauntB_Q_max(-m, n, mu, nu); double a2q[q2max+1], a3q[Qmax+1], a2q0, a3q0; int err; if (mu == nu) { @@ -230,7 +255,7 @@ complex double qpms_trans_single_B_Xu(int m, int n, int mu, int nu, sph_t kdlj, if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bes)) abort(); complex double sum = 0; - for (int q = 0; q <= Qmax; ++q) { + for (int q = 0; q <= realQmax; ++q) { int p = n+nu-2*q; double a2q_n = a2q[q]/a2q0; double a3q_n = a3q[q]/a3q0; @@ -267,6 +292,7 @@ complex double qpms_trans_single_B(qpms_normalisation_t norm, int q2max = gaunt_q_max(-m-1,n+1,mu+1,nu); int Qmax = gaunt_q_max(-m,n+1,mu,nu); + int realQmax = gauntB_Q_max(-m,n,mu,nu); double a2q[q2max+1], a3q[Qmax+1], a2q0, a3q0; int err; if (mu == nu) { @@ -288,7 +314,7 @@ complex double qpms_trans_single_B(qpms_normalisation_t norm, if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bes)) abort(); complex double sum = 0; - for (int q = 0; q <= Qmax; ++q) { + for (int q = 0; q <= realQmax; ++q) { int p = n+nu-2*q; double a2q_n = a2q[q]/a2q0; double a3q_n = a3q[q]/a3q0; @@ -327,6 +353,7 @@ complex double qpms_trans_single_B_Taylor(int m, int n, int mu, int nu, sph_t kd int q2max = gaunt_q_max(-m-1,n+1,mu+1,nu); int Qmax = gaunt_q_max(-m,n+1,mu,nu); + int realQmax = gauntB_Q_max(-m,n,mu,nu); double a2q[q2max+1], a3q[Qmax+1], a2q0, a3q0; int err; if (mu == nu) { @@ -347,7 +374,7 @@ complex double qpms_trans_single_B_Taylor(int m, int n, int mu, int nu, sph_t kd if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bes)) abort(); complex double sum = 0; - for (int q = 0; q <= Qmax; ++q) { + for (int q = 0; q <= realQmax; ++q) { int p = n+nu-2*q; double a2q_n = a2q[q]/a2q0; double a3q_n = a3q[q]/a3q0; @@ -453,14 +480,13 @@ static void qpms_trans_calculator_multipliers_A_general( } -static void qpms_trans_calculator_multipliers_B_general( +/*static*/ void qpms_trans_calculator_multipliers_B_general( qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int Qmax) { - assert(Qmax == gaunt_q_max(-m,n+1,mu,nu)); + assert(Qmax == gauntB_Q_max(-m,n,mu,nu)); int q2max = gaunt_q_max(-m-1,n+1,mu+1,nu); - assert(Qmax == q2max); - // FIXME remove the q2max variable altogether, as it is probably equal - // to Qmax + // assert(Qmax == q2max); + // FIXME is it safe to replace q2max with Qmax in gaunt_xu?? double a2q[q2max+1], a3q[Qmax+1], a2q0, a3q0; int err; if (mu == nu) { @@ -472,7 +498,7 @@ static void qpms_trans_calculator_multipliers_B_general( gaunt_xu(-m-1,n+1,mu+1,nu,q2max,a2q,&err); if (err) abort(); a2q0 = a2q[0]; } - gaunt_xu(-m,n+1,mu,nu,Qmax,a3q,&err); if (err) abort(); + gaunt_xu(-m,n+1,mu,nu,q2max,a3q,&err); if (err) abort(); // FIXME this should probably go away a3q0 = a3q[0]; @@ -494,7 +520,7 @@ static void qpms_trans_calculator_multipliers_B_general( presum *= I * ipow(nu+n) /*different from old Taylor */ * normfac / ( (4*n)*(n+1)*(n+m+1)); - for (int q = 0; q <= Qmax; ++q) { + for (int q = BQ_OFFSET; q <= Qmax; ++q) { int p = n+nu-2*q; double a2q_n = a2q[q]/a2q0; double a3q_n = a3q[q]/a3q0; @@ -507,7 +533,7 @@ static void qpms_trans_calculator_multipliers_B_general( double summandq = ((2*(n+1)*(nu-mu)*a2q_n -(-nu*(nu+1) - n*(n+3) - 2*mu*(n+1)+p*(p+3))* a3q_n) *min1pow(q)); - dest[q] = Ppfac * summandq * presum; + dest[q-BQ_OFFSET] = Ppfac * summandq * presum; } } @@ -581,12 +607,12 @@ static void qpms_trans_calculator_multipliers_A_Taylor( static void qpms_trans_calculator_multipliers_B_Taylor( complex double *dest, int m, int n, int mu, int nu, int Qmax) { - assert(Qmax == gaunt_q_max(-m,n+1,mu,nu)); + assert(Qmax == gauntB_Q_max(-m,n,mu,nu)); int q2max = gaunt_q_max(-m-1,n+1,mu+1,nu); - assert(Qmax == q2max); + //assert(Qmax == q2max); // FIXME remove the q2max variable altogether, as it is probably equal // to Qmax - double a2q[q2max+1], a3q[Qmax+1], a2q0, a3q0; + double a2q[q2max+1], a3q[q2max+1], a2q0, a3q0; int err; if (mu == nu) { for (int q = 0; q <= q2max; ++q) @@ -597,7 +623,7 @@ static void qpms_trans_calculator_multipliers_B_Taylor( gaunt_xu(-m-1,n+1,mu+1,nu,q2max,a2q,&err); if (err) abort(); a2q0 = a2q[0]; } - gaunt_xu(-m,n+1,mu,nu,Qmax,a3q,&err); if (err) abort(); + gaunt_xu(-m,n+1,mu,nu,q2max,a3q,&err); if (err) abort(); a3q0 = a3q[0]; double exponent=(lgamma(2*n+3)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2) @@ -609,7 +635,7 @@ static void qpms_trans_calculator_multipliers_B_Taylor( presum *= I * (min1pow(m+n) *sqrt((2.*n+1)/(2.*nu+1)) / ( (4*n)*(n+1)*(n+m+1))); - for (int q = 0; q <= Qmax; ++q) { + for (int q = BQ_OFFSET; q <= Qmax; ++q) { int p = n+nu-2*q; double a2q_n = a2q[q]/a2q0; double a3q_n = a3q[q]/a3q0; @@ -622,7 +648,7 @@ static void qpms_trans_calculator_multipliers_B_Taylor( double summandq = ((2*(n+1)*(nu-mu)*a2q_n -(-nu*(nu+1) - n*(n+3) - 2*mu*(n+1)+p*(p+3))* a3q_n) *min1pow(q)); - dest[q] = Ppfac * summandq * presum; + dest[q-BQ_OFFSET] = Ppfac * summandq * presum; } } @@ -645,17 +671,17 @@ int qpms_trans_calculator_multipliers_A(qpms_normalisation_t norm, complex doubl assert(0); } -int qpms_trans_calculator_multipliers_B(qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) { +int qpms_trans_calculator_multipliers_B(qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int Qmax) { switch (qpms_normalisation_t_normonly(norm)) { case QPMS_NORMALISATION_TAYLOR: #ifdef USE_SEPARATE_TAYLOR - qpms_trans_calculator_multipliers_B_Taylor(dest,m,n,mu,nu,qmax); + qpms_trans_calculator_multipliers_B_Taylor(dest,m,n,mu,nu,Qmax); return 0; break; #endif case QPMS_NORMALISATION_XU: case QPMS_NORMALISATION_KRISTENSSON: - qpms_trans_calculator_multipliers_B_general(norm, dest, m, n, mu, nu, qmax); + qpms_trans_calculator_multipliers_B_general(norm, dest, m, n, mu, nu, Qmax); return 0; break; default: @@ -705,14 +731,14 @@ qpms_trans_calculator int m, n, mu, nu; qpms_y2mn_p(y,&m,&n); qpms_y2mn_p(yu,&mu,&nu); - qmaxsum += 1 + ( + qmaxsum += (1 - BQ_OFFSET) + ( qmaxes[qpms_trans_calculator_index_yyu(c,y,yu)] - = gaunt_q_max(-m,n+1,mu,nu)); + = gauntB_Q_max(-m,n,mu,nu)); } c->B_multipliers[0] = malloc(qmaxsum * sizeof(complex double)); // calculate multiplier beginnings for(size_t i = 0; i < SQ(c->nelem); ++i) - c->B_multipliers[i+1] = c->B_multipliers[i] + qmaxes[i] + 1; + c->B_multipliers[i+1] = c->B_multipliers[i] + qmaxes[i] + (1 - BQ_OFFSET); // calculate the multipliers for(size_t y = 0; y < c->nelem; ++y) for(size_t yu = 0; yu < c->nelem; ++yu) { @@ -782,14 +808,14 @@ static inline complex double qpms_trans_calculator_get_B_precalcbuf(const qpms_t bool r_ge_d, qpms_bessel_t J, const complex double *bessel_buf, const double *legendre_buf) { size_t i = qpms_trans_calculator_index_mnmunu(c, m, n, mu, nu); - size_t qmax = c->B_multipliers[i+1] - c->B_multipliers[i] - 1; - assert(qmax == gaunt_q_max(-m,n+1,mu,nu)); + size_t qmax = c->B_multipliers[i+1] - c->B_multipliers[i] - (1 - BQ_OFFSET); + assert(qmax == gauntB_Q_max(-m,n,mu,nu)); complex double sum = 0; - for(int q = 0; q <= qmax; ++q) { + for(int q = BQ_OFFSET; q <= qmax; ++q) { int p = n+nu-2*q; double Pp_ = legendre_buf[gsl_sf_legendre_array_index(p+1, abs(mu-m))]; complex double zp_ = bessel_buf[p+1]; - complex double multiplier = c->B_multipliers[i][q]; + complex double multiplier = c->B_multipliers[i][q-BQ_OFFSET]; sum += Pp_ * zp_ * multiplier; } complex double eimf = cexp(I*(mu-m)*kdlj.phi);