"fix" the q ranges for the B coefficients in translations.c

Former-commit-id: 17084dd0b57cb4271c67d5d95b7e033c3968276f

qpms/test_translations_xu_multipliers.c

@@ -4,6 +4,11 @@
 //#include <math.h>
 #include <complex.h>
 
+#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))/2));
-					int Qmax_alt = MIN(n,MIN(nu,(n+nu+1-abs(mu-m))));
 					qpms_trans_calculator_multipliers_B_general(QPMS_NORMALISATION_XU_CS,
-							dest, m, n, mu, nu, Qmax);
+							dest, m, n, mu, nu, Qmax_alt);
-					for(int q = 0; q <= Qmax; ++q) {
+					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);

qpms/translations.c

@@ -9,9 +9,19 @@
 #include "assert_cython_workaround.h"
 #include <stdlib.h> //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);
+	// assert(Qmax == q2max);
-	// FIXME remove the q2max variable altogether, as it is probably equal
+	// FIXME is it safe to replace q2max with Qmax in gaunt_xu??
-	// to Qmax
 	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;
+	size_t qmax = c->B_multipliers[i+1] - c->B_multipliers[i] - (1 - BQ_OFFSET);
-	assert(qmax == gaunt_q_max(-m,n+1,mu,nu));
+	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);