From 8b9c5a794c3410fb383bd65579797775573beb86 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 19 Aug 2019 08:01:59 +0300 Subject: [PATCH] translations legacy code removal Former-commit-id: 98e883f9c7586208623acca05585f747dc57f911 --- qpms/cytranslations.pyx | 39 +-- qpms/qpms_cdefs.pxd | 4 - qpms/translations.c | 660 +--------------------------------------- qpms/translations.h | 41 --- 4 files changed, 4 insertions(+), 740 deletions(-) diff --git a/qpms/cytranslations.pyx b/qpms/cytranslations.pyx index ab41a64..801346a 100644 --- a/qpms/cytranslations.pyx +++ b/qpms/cytranslations.pyx @@ -43,7 +43,7 @@ In simple words, it works like this: #------------------------------------- - +# This one is probably not used anymore an can perhaps be removed: cdef void loop_D_iiiidddii_As_D_lllldddbl(char **args, np.npy_intp *dims, np.npy_intp *steps, void *data) nogil: cdef np.npy_intp i, n = dims[0] cdef void *func = (data)#[0] @@ -99,11 +99,6 @@ np.import_ufunc() # BTW, aren't there anonymous arrays in cython? -cdef np.PyUFuncGenericFunction trans_X_taylor_loop_func[1] -cdef void *trans_A_taylor_elementwise_funcs[1] -cdef void *trans_B_taylor_elementwise_funcs[1] - -trans_X_taylor_loop_func[0] = loop_D_iiiidddii_As_D_lllldddbl # types to be used for all of the single-type translation # coefficient retrieval ufuncs called like @@ -139,38 +134,6 @@ ufunc__get_both_coeff_types[9] = np.NPY_CDOUBLE ufunc__get_both_coeff_types[10] = np.NPY_CDOUBLE -trans_A_taylor_elementwise_funcs[0] = qpms_trans_single_A_Taylor_ext -trans_B_taylor_elementwise_funcs[0] = qpms_trans_single_B_Taylor_ext - -trans_A_Taylor = np.PyUFunc_FromFuncAndData( - trans_X_taylor_loop_func, # func - trans_A_taylor_elementwise_funcs, #data - ufunc__get_either_trans_coeff_types, # types - 1, # ntypes: number of supported input types - 9, # nin: number of input args - 1, # nout: number of output args - 0, # identity element, unused - "trans_A_Taylor", # name - """ - TODO computes the E-E or M-M translation coefficient in Taylor's normalisation - """, # doc - 0 # unused, for backward compatibility of numpy c api - ) - -trans_B_Taylor = np.PyUFunc_FromFuncAndData( - trans_X_taylor_loop_func, - trans_B_taylor_elementwise_funcs, - ufunc__get_either_trans_coeff_types, - 1, # number of supported input types - 9, # number of input args - 1, # number of output args - 0, # identity element, unused - "trans_B_Taylor", - """ - TODO computes the E-E or M-M translation coefficient in Taylor's normalisation - """, - 0 # unused - ) # --------------------------------------------- # Wrapper for the qpms_trans_calculator "class" diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 2474e20..b331110 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -265,10 +265,6 @@ cdef extern from "symmetries.h": cdef extern from "translations.h": - cdouble qpms_trans_single_A_Taylor_ext(int m, int n, int mu, int nu, - double r, double th, double ph, int r_ge_d, int J) nogil - cdouble qpms_trans_single_B_Taylor_ext(int m, int n, int mu, int nu, - double r, double th, double ph, int r_ge_d, int J) nogil struct qpms_trans_calculator: int lMax size_t nelem diff --git a/qpms/translations.c b/qpms/translations.c index 323309f..c8733e4 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -114,41 +114,6 @@ 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)); } -#if 0 //Apparently, this is duplicit to that in bessel.c (which also support complex x) -int qpms_sph_bessel_fill(qpms_bessel_t typ, int lmax, double x, complex double *result_array) { - int retval; - double tmparr[lmax+1]; - switch(typ) { - case QPMS_BESSEL_REGULAR: - retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr); - for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l]; - return retval; - break; - case QPMS_BESSEL_SINGULAR: //FIXME: is this precise enough? Would it be better to do it one-by-one? - retval = gsl_sf_bessel_yl_array(lmax,x,tmparr); - for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l]; - return retval; - break; - case QPMS_HANKEL_PLUS: - case QPMS_HANKEL_MINUS: - retval = gsl_sf_bessel_jl_steed_array(lmax, x, tmparr); - for (int l = 0; l <= lmax; ++l) result_array[l] = tmparr[l]; - if(retval) return retval; - retval = gsl_sf_bessel_yl_array(lmax, x, tmparr); - if (typ==QPMS_HANKEL_PLUS) - for (int l = 0; l <= lmax; ++l) result_array[l] += I * tmparr[l]; - else - for (int l = 0; l <= lmax; ++l) result_array[l] +=-I * tmparr[l]; - return retval; - break; - default: - abort(); - //return GSL_EDOM; - } - assert(0); -} -#endif - static inline double qpms_trans_normlogfac(qpms_normalisation_t norm, int m, int n, int mu, int nu) { return -0.5*(lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1)); @@ -237,109 +202,6 @@ complex double qpms_trans_single_A(qpms_normalisation_t norm, return presum * sum; } -complex double qpms_trans_single_A_Taylor(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); - - int qmax = gaunt_q_max(-m,n,mu,nu); // nemá tu být +m? - // N.B. -m !!!!!! - double a1q[qmax+1]; - int err; - gaunt_xu(-m,n,mu,nu,qmax,a1q,&err); - double a1q0 = a1q[0]; - if (err) abort(); - - double leg[gsl_sf_legendre_array_n(n+nu)]; - if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,leg)) abort(); - complex double bes[n+nu+1]; - if (qpms_sph_bessel_fill(J, n+nu, kdlj.r, bes)) abort(); - complex double sum = 0; - for(int q = 0; q <= qmax; ++q) { - int p = n+nu-2*q; - int Pp_order = mu-m; - //if(p < abs(Pp_order)) continue; // FIXME raději nastav lépe meze - assert(p >= abs(Pp_order)); - double a1q_n = a1q[q] / a1q0; - double Pp = leg[gsl_sf_legendre_array_index(p, abs(Pp_order))]; - if (Pp_order < 0) Pp *= min1pow(mu-m) * exp(lgamma(1+p+Pp_order)-lgamma(1+p-Pp_order)); - complex double zp = bes[p]; - complex double summandq = (n*(n+1) + nu*(nu+1) - p*(p+1)) * min1pow(q) * a1q_n * zp * Pp; - sum += summandq; // TODO KAHAN - } - - double exponent=(lgamma(2*n+1)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2) - +lgamma(n+nu+m-mu+1)-lgamma(n-m+1)-lgamma(nu+mu+1) - +lgamma(n+nu+1) - lgamma(2*(n+nu)+1)); - complex double presum = exp(exponent); - presum *= cexp(I*(mu-m)*kdlj.phi) * min1pow(m) * ipow(nu+n) / (4*n); - - // N.B. ipow(nu-n) is different from the general formula! - complex double prenormratio = ipow(nu-n) * sqrt(((2.*nu+1)/(2.*n+1))* exp( - lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1))); - return (presum / prenormratio) * sum; -} - -// [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) { - assert(0); // FIXME probably gives wrong values, do not use. - 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) { - for (int q = 0; q <= q2max; ++q) - a2q[q] = 0; - a2q0 = 1; - } - else { - 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(); - a3q0 = a3q[0]; - - double leg[gsl_sf_legendre_array_n(n+nu+1)]; - if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,-1,leg)) abort(); - complex double bes[n+nu+2]; - if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bes)) abort(); - - complex double sum = 0; - 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; - complex double zp_ = bes[p+1]; - int Pp_order_ = mu-m; - //if(p+1 < abs(Pp_order_)) continue; // FIXME raději nastav lépe meze - assert(p+1 >= abs(Pp_order_)); - double Pp_ = leg[gsl_sf_legendre_array_index(p+1, abs(Pp_order_))]; - if (Pp_order_ < 0) Pp_ *= min1pow(mu-m) * exp(lgamma(1+1+p+Pp_order_)-lgamma(1+1+p-Pp_order_)); - complex 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) * zp_ * Pp_); - sum += summandq; // TODO KAHAN - } - - double exponent=(lgamma(2*n+3)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2) - +lgamma(n+nu+m-mu+2)-lgamma(n-m+1)-lgamma(nu+mu+1) - +lgamma(n+nu+2) - lgamma(2*(n+nu)+3)); - complex double presum = exp(exponent); - presum *= cexp(I*(mu-m)*kdlj.phi) * min1pow(m) * ipow(nu+n+1) / ( - (4*n)*(n+1)*(n+m+1)); - - // Taylor normalisation v2, proven to be equivalent - complex double prenormratio = ipow(nu-n); - - return (presum / prenormratio) * sum; -} complex double qpms_trans_single_B(qpms_normalisation_t norm, int m, int n, int mu, int nu, sph_t kdlj, @@ -428,93 +290,6 @@ complex double qpms_trans_single_B(qpms_normalisation_t norm, return presum * sum; } -complex double qpms_trans_single_B_Taylor(int m, int n, int mu, int nu, sph_t kdlj, - bool r_ge_d, qpms_bessel_t J) { - assert(0); // FIXME probably gives wrong values, do not use. - if(r_ge_d) J = QPMS_BESSEL_REGULAR; - double costheta = cos(kdlj.theta); - - 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) { - for (int q = 0; q <= q2max; ++q) - a2q[q] = 0; - a2q0 = 1; - } - else { - 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(); - a3q0 = a3q[0]; - - double leg[gsl_sf_legendre_array_n(n+nu+1)]; - if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,-1,leg)) abort(); - complex double bes[n+nu+2]; - if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bes)) abort(); - - complex double sum = 0; - 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; - complex double zp_ = bes[p+1]; - int Pp_order_ = mu-m; - //if(p+1 < abs(Pp_order_)) continue; // FIXME raději nastav lépe meze - assert(p+1 >= abs(Pp_order_)); - double Pp_ = leg[gsl_sf_legendre_array_index(p+1, abs(Pp_order_))]; - if (Pp_order_ < 0) Pp_ *= min1pow(mu-m) * exp(lgamma(1+1+p+Pp_order_)-lgamma(1+1+p-Pp_order_)); - complex 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) * zp_ * Pp_); - sum += summandq; //TODO KAHAN - } - - double exponent=(lgamma(2*n+3)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2) - +lgamma(n+nu+m-mu+2)-lgamma(n-m+1)-lgamma(nu+mu+1) - +lgamma(n+nu+2) - lgamma(2*(n+nu)+3)); - complex double presum = exp(exponent); - presum *= cexp(I*(mu-m)*kdlj.phi) * min1pow(m) * ipow(nu+n+1) / ( - (4*n)*(n+1)*(n+m+1)); - - // Taylor normalisation v2, proven to be equivalent - // ipow(nu-n) is different from the new general formula!!! - complex double prenormratio = ipow(nu-n) * sqrt(((2.*nu+1)/(2.*n+1))* exp( - lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1))); - - return (presum / prenormratio) * sum; -} - -complex double qpms_trans_single_A_Taylor_ext(int m, int n, int mu, int nu, - double kdlj_r, double kdlj_theta, double kdlj_phi, int r_ge_d, int J) { - sph_t kdlj = {kdlj_r, kdlj_theta, kdlj_phi}; - return qpms_trans_single_A_Taylor(m,n,mu,nu,kdlj,r_ge_d,J); -} - -complex double qpms_trans_single_B_Taylor_ext(int m, int n, int mu, int nu, - double kdlj_r, double kdlj_theta, double kdlj_phi, int r_ge_d, int J) { - sph_t kdlj = {kdlj_r, kdlj_theta, kdlj_phi}; - return qpms_trans_single_B_Taylor(m,n,mu,nu,kdlj,r_ge_d,J); -} - -void qpms_trans_calculator_free(qpms_trans_calculator *c) { - free(c->A_multipliers[0]); - free(c->A_multipliers); - free(c->B_multipliers[0]); - free(c->B_multipliers); -#ifdef LATTICESUMS - qpms_ewald3_constants_free(e3c); -#endif -#ifdef LATTICESUMS_OLD - free(c->hct); -#endif - free(c->legendre0); - free(c); -} - static inline size_t qpms_trans_calculator_index_mnmunu(const qpms_trans_calculator *c, int m, int n, int mu, int nu){ return c->nelem * qpms_mn2y(m,n) + qpms_mn2y(mu,nu); @@ -530,7 +305,7 @@ static inline size_t qpms_trans_calculator_index_yyu(const qpms_trans_calculator static inline double fsq(double x) {return x * x; } -static void qpms_trans_calculator_multipliers_A_general( +static void qpms_trans_calculator_multipliers_A( qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) { assert(qmax == gaunt_q_max(-m,n,mu,nu)); @@ -589,7 +364,7 @@ static void qpms_trans_calculator_multipliers_A_general( // as in [Xu](61) -double cruzan_bfactor(int M, int n, int mu, int nu, int p) { +static double cruzan_bfactor(int M, int n, int mu, int nu, int p) { double logprefac = lgamma(n+M+1) - lgamma(n-M+1) + lgamma(nu+mu+1) - lgamma(nu-mu+1) + lgamma(p-M-mu+2) - lgamma(p+M+mu+2); logprefac *= 0.5; @@ -599,7 +374,7 @@ double cruzan_bfactor(int M, int n, int mu, int nu, int p) { } -void qpms_trans_calculator_multipliers_B_general( +void qpms_trans_calculator_multipliers_B( qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int Qmax){ // This is according to the Cruzan-type formula [Xu](59) @@ -658,192 +433,6 @@ void qpms_trans_calculator_multipliers_B_general( } } -/*static*/ void qpms_trans_calculator_multipliers_B_general_oldXu( - qpms_normalisation_t norm, - complex double *dest, int m, int n, int mu, int nu, int Qmax) { - assert(0); // FIXME probably gives wrong values, do not use. - 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 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) { - for (int q = 0; q <= q2max; ++q) - a2q[q] = 0; - a2q0 = 1; - } - else { - 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,q2max,a3q,&err); if (err) abort(); // FIXME this should probably go away - a3q0 = a3q[0]; - - - int csphase = qpms_normalisation_t_csphase(norm); //TODO FIXME use this - double normlogfac= qpms_trans_normlogfac(norm,m,n,mu,nu); - double normfac = qpms_trans_normfac(norm,m,n,mu,nu); - /// N<-M type coefficients w.r.t. Kristensson's convention. Csphase has been already taken into acct ^^^. - normfac *= qpms_normalisation_factor_M_noCS(norm, nu, mu) - / qpms_normalisation_factor_N_noCS(norm, n, m); - // TODO use csphase to modify normfac here!!!! - // normfac = xxx ? -normfac : normfac; - normfac *= min1pow(m);//different from old taylor - - - - double exponent=(lgamma(2*n+3)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2) - +lgamma(n+nu+m-mu+2)-lgamma(n-m+1)-lgamma(nu+mu+1) - +lgamma(n+nu+2) - lgamma(2*(n+nu)+3)) - +normlogfac; - complex double presum = exp(exponent); - presum *= I * ipow(nu+n) /*different from old Taylor */ * normfac / ( - (4*n)*(n+1)*(n+m+1)); - - 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; - int Pp_order_ = mu-m; - //if(p+1 < abs(Pp_order_)) continue; // FIXME raději nastav lépe meze - assert(p+1 >= abs(Pp_order_)); - double Ppfac = (Pp_order_ >= 0) ? 1 : - - min1pow(mu-m) * exp(lgamma(1+1+p+Pp_order_)-lgamma(1+1+p-Pp_order_)); - 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-BQ_OFFSET] = Ppfac * summandq * presum; - } -} - -//#if 0 -static void qpms_trans_calculator_multipliers_A_Taylor( - complex double *dest, int m, int n, int mu, int nu, int qmax) { - assert(qmax == gaunt_q_max(-m,n,mu,nu)); - double a1q[qmax+1]; - int err; - gaunt_xu(-m,n,mu,nu,qmax,a1q,&err); - if (err) abort(); - double a1q0 = a1q[0]; - - double exponent=(lgamma(2*n+1)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2) - +lgamma(n+nu+m-mu+1)-lgamma(n-m+1)-lgamma(nu+mu+1) - +lgamma(n+nu+1) - lgamma(2*(n+nu)+1)) - 0.5*( // ex-prenormratio - lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1)); - double presum = exp(exponent); - presum *= min1pow(m+n) * sqrt((2.*n+1)/(2.*nu+1)) / (4*n); - - for(int q = 0; q <= qmax; q++) { - int p = n+nu-2*q; - int Pp_order = mu - m; - assert(p >= abs(Pp_order)); - double a1q_n = a1q[q] / a1q0; - // Assuming non_normalized legendre polynomials! - double Ppfac = (Pp_order >= 0) ? 1 : - min1pow(mu-m) * exp(lgamma(1+p+Pp_order)-lgamma(1+p-Pp_order)); - double summandfac = (n*(n+1) + nu*(nu+1) - p*(p+1)) * min1pow(q) * a1q_n; - dest[q] = presum * summandfac * Ppfac; - // FIXME I might not need complex here - } -} -//#endif -#if 0 -static void qpms_trans_calculator_multipliers_A_Taylor( - complex double *dest, int m, int n, int mu, int nu, int qmax) { - assert(qmax == gaunt_q_max(-m,n,mu,nu)); - double a1q[qmax+1]; - int err; - gaunt_xu(-m,n,mu,nu,qmax,a1q,&err); - if (err) abort(); - double a1q0 = a1q[0]; - for(int q = 0; q <= qmax; ++q) { - int p = n+nu-2*q; - int Pp_order = mu-m; - //if(p < abs(Pp_order)) continue; // FIXME raději nastav lépe meze - assert(p >= abs(Pp_order)); - double a1q_n = a1q[q] / a1q0; - //double Pp = leg[gsl_sf_legendre_array_index(p, abs(Pp_order))]; - //complex double zp = bes[p]; - dest[q] = (n*(n+1) + nu*(nu+1) - p*(p+1)) * min1pow(q) * a1q_n /* * zp * Pp*/; - if (Pp_order < 0) dest[q] *= min1pow(mu-m) * exp(lgamma(1+p+Pp_order)-lgamma(1+p-Pp_order)); - //sum += summandq; - } - - double exponent=(lgamma(2*n+1)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2) - +lgamma(n+nu+m-mu+1)-lgamma(n-m+1)-lgamma(nu+mu+1) - +lgamma(n+nu+1) - lgamma(2*(n+nu)+1)); - complex double presum = exp(exponent); - presum *=/* cexp(I*(mu-m)*kdlj.phi) * */ min1pow(m) * ipow(nu+n) / (4*n); - - complex double prenormratio = ipow(nu-n) * sqrt(((2.*nu+1)/(2.*n+1))* exp( - lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1))); - //return (presum / prenormratio) * sum; - for(int q=0;q<=qmax;++q) dest[q] *= presum / prenormratio; -} -#endif - - - -static void qpms_trans_calculator_multipliers_B_Taylor( - complex double *dest, int m, int n, int mu, int nu, int Qmax) { - assert(0); // FIXME probably gives wrong values, do not use. - 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 - double a2q[q2max+1], a3q[q2max+1], a2q0, a3q0; - int err; - if (mu == nu) { - for (int q = 0; q <= q2max; ++q) - a2q[q] = 0; - a2q0 = 1; - } - else { - 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,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) - +lgamma(n+nu+m-mu+2)-lgamma(n-m+1)-lgamma(nu+mu+1) - +lgamma(n+nu+2) - lgamma(2*(n+nu)+3)) - 0.5 * ( - lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1) - -lgamma(nu+mu+1)); - complex double presum = exp(exponent); - presum *= I * (min1pow(m+n) *sqrt((2.*n+1)/(2.*nu+1)) / ( - (4*n)*(n+1)*(n+m+1))); - - 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; - int Pp_order_ = mu-m; - //if(p+1 < abs(Pp_order_)) continue; // FIXME raději nastav lépe meze - assert(p+1 >= abs(Pp_order_)); - double Ppfac = (Pp_order_ >= 0) ? 1 : - - min1pow(mu-m) * exp(lgamma(1+1+p+Pp_order_)-lgamma(1+1+p-Pp_order_)); - 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-BQ_OFFSET] = Ppfac * summandq * presum; - } -} - -int qpms_trans_calculator_multipliers_A(qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int qmax) { - qpms_trans_calculator_multipliers_A_general(norm, dest, m, n, mu, nu, qmax); - return 0; -} - -int qpms_trans_calculator_multipliers_B(qpms_normalisation_t norm, complex double *dest, int m, int n, int mu, int nu, int Qmax) { - qpms_trans_calculator_multipliers_B_general(norm, dest, m, n, mu, nu, Qmax); - return 0; -} - qpms_trans_calculator *qpms_trans_calculator_init (const int lMax, const qpms_normalisation_t normalisation) { TROPS_ONLY_EIMF_IMPLEMENTED(normalisation); @@ -1556,249 +1145,6 @@ int qpms_trans_calculator_e32_short_points_and_shift(const qpms_trans_calculator #endif // 0 #endif // LATTICESUMS32 -#ifdef LATTICESUMS_OLD - -int qpms_trans_calculator_get_shortrange_AB_arrays_buf(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, - size_t deststride, size_t srcstride, - sph_t kdlj, qpms_bessel_t J, - qpms_l_t lrcutoff, unsigned kappa, double cc, // regularisation params - complex double *bessel_buf, double *legendre_buf - ) { - assert(J == QPMS_HANKEL_PLUS); // support only J == 3 for now - if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) { - for (size_t i = 0; i < c->nelem; ++i) - for (size_t j = 0; j < c->nelem; ++j) { - *(Adest + i*srcstride + j*deststride) = NAN+I*NAN; - *(Bdest + i*srcstride + j*deststride) = NAN+I*NAN; - } - // TODO warn? different return value? - return 0; - } - { - double costheta = cos(kdlj.theta); - if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*c->lMax+1, - costheta,-1,legendre_buf)) abort(); - // if (qpms_sph_bessel_fill(J, 2*c->lMax+1, kdlj.r, bessel_buf)) abort(); // original - hankelparts_fill(NULL, bessel_buf, 2*c->lMax+1, lrcutoff, c->hct, kappa, cc, kdlj.r); - size_t desti = 0, srci = 0; - for (int n = 1; n <= c->lMax; ++n) for (int m = -n; m <= n; ++m) { - for (int nu = 1; nu <= c->lMax; ++nu) for (int mu = -nu; mu <= nu; ++mu) { - size_t assertindex = qpms_trans_calculator_index_mnmunu(c,m,n,mu,nu); - assert(assertindex == desti*c->nelem + srci); - *(Adest + deststride * desti + srcstride * srci) = - qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu, - kdlj,false,J,bessel_buf,legendre_buf); - *(Bdest + deststride * desti + srcstride * srci) = - qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu, - kdlj,false,J,bessel_buf,legendre_buf); - ++srci; - } - ++desti; - srci = 0; - } - return 0; - } -} - -int qpms_trans_calculator_get_shortrange_AB_buf_p(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, - int m, int n, int mu, int nu, sph_t kdlj, - qpms_bessel_t J, - qpms_l_t lrcutoff, unsigned kappa, double cc, // regularisation params - complex double *bessel_buf, double *legendre_buf) { - assert(J == QPMS_HANKEL_PLUS); // support only J == 3 for now - if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) { - *Adest = NAN+I*NAN; - *Bdest = NAN+I*NAN; - // TODO warn? different return value? - return 0; - } - { - double costheta = cos(kdlj.theta); - if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1, - costheta,-1,legendre_buf)) abort(); - //if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort(); // original - hankelparts_fill(NULL, bessel_buf, 2*c->lMax+1, lrcutoff, c->hct, kappa, cc, kdlj.r); - - *Adest = qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu, - kdlj,false,J,bessel_buf,legendre_buf); - *Bdest = qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu, - kdlj,false,J,bessel_buf,legendre_buf); - return 0; - } -} - -// Short-range parts of the translation coefficients -int qpms_trans_calculator_get_shortrange_AB_p(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, - qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, - qpms_bessel_t J /* Only J=3 valid for now */, - qpms_l_t lrcutoff, unsigned kappa, double cc) { - double leg[gsl_sf_legendre_array_n(2*c->lMax+1)]; - complex double bes[2*c->lMax+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth. - return qpms_trans_calculator_get_shortrange_AB_buf_p(c,Adest, Bdest,m,n,mu,nu,kdlj,J, - lrcutoff, kappa, cc, - bes, leg); -} - -int qpms_trans_calculator_get_shortrange_AB_arrays(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, - size_t deststride, size_t srcstride, - sph_t kdlj, qpms_bessel_t J /* Only J=3 valid for now */, - qpms_l_t lrcutoff, unsigned kappa, double cc) { - double leg[gsl_sf_legendre_array_n(c->lMax+c->lMax+1)]; - complex double bes[2*c->lMax+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth. - return qpms_trans_calculator_get_shortrange_AB_arrays_buf(c, - Adest, Bdest, deststride, srcstride, - kdlj, J, - lrcutoff, kappa, cc, - bes, leg); -} - -// Long-range parts -static inline complex double qpms_trans_calculator_get_2DFT_longrange_A_precalcbuf(const qpms_trans_calculator *c, - int m, int n, int mu, int nu, sph_t k_sph /* theta must be M_PI_2 */, - qpms_bessel_t J /* must be 3 for now */, - const complex double *lrhankel_recparts_buf) { - assert(J == QPMS_HANKEL_PLUS); - //assert(k_sph.theta == M_PI_2); CHECK IN ADVANCE INSTEAD - //assert(k_sph.r > 0); - size_t i = qpms_trans_calculator_index_mnmunu(c, m, n, mu, nu); - size_t qmax = c->A_multipliers[i+1] - c->A_multipliers[i] - 1; - assert(qmax == gaunt_q_max(-m,n,mu,nu)); - complex double sum, kahanc; - ckahaninit(&sum, &kahanc); - for(size_t q = 0; q <= qmax; ++q) { - int p = n+nu-2*q; - double Pp = c->legendre0[gsl_sf_legendre_array_index(p, abs(mu-m))]; - complex double zp = trindex_cd(lrhankel_recparts_buf, p)[abs(mu-m)]; // orig: bessel_buf[p]; - if (mu - m < 0) zp *= min1pow(mu-m); // DLMF 10.4.1 - complex double multiplier = c->A_multipliers[i][q]; - ckahanadd(&sum, &kahanc, Pp * zp * multiplier); - } - complex double eimf = cexp(I*(mu-m)*k_sph.phi); - return sum * eimf * ipow(mu-m); -} - -static inline complex double qpms_trans_calculator_get_2DFT_longrange_B_precalcbuf(const qpms_trans_calculator *c, - int m, int n, int mu, int nu, sph_t k_sph /* theta must be M_PI_2 */, - qpms_bessel_t J /* must be 3 for now */, - const complex double *lrhankel_recparts_buf) { - assert(J == QPMS_HANKEL_PLUS); - 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 - BQ_OFFSET); - assert(qmax == gauntB_Q_max(-m,n,mu,nu)); - complex double sum, kahanc; - ckahaninit(&sum, &kahanc); - for(int q = BQ_OFFSET; q <= qmax; ++q) { - int p = n+nu-2*q; - double Pp_ = c->legendre0[gsl_sf_legendre_array_index(p+1, abs(mu-m))]; - complex double zp_ = trindex_cd(lrhankel_recparts_buf, p+1)[abs(mu-m)]; // orig: bessel_buf[p+1]; - if (mu - m < 0) zp_ *= min1pow(mu-m); // DLMF 10.4.1 - complex double multiplier = c->B_multipliers[i][q-BQ_OFFSET]; - ckahanadd(&sum, &kahanc, Pp_ * zp_ * multiplier); - } - complex double eimf = cexp(I*(mu-m)*k_sph.phi); - return sum * eimf * ipow(mu-m); -} - -int qpms_trans_calculator_get_2DFT_longrange_AB_buf_p(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, - int m, int n, int mu, int nu, sph_t k_sph, - qpms_bessel_t J, - qpms_l_t lrk_cutoff, unsigned kappa, double cv, double k0, - complex double *lrhankel_recparts_buf) { - - assert (J == QPMS_HANKEL_PLUS); - assert(k_sph.theta == M_PI_2); - - { - //double costheta = cos(kdlj.theta); - //if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1, - // costheta,-1,legendre_buf)) abort(); - //if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort(); - lrhankel_recpart_fill(lrhankel_recparts_buf, 2*c->lMax+1 /* TODO n+nu+1 might be enough */, - lrk_cutoff, c->hct, kappa, cv, k0, k_sph.r); - *Adest = qpms_trans_calculator_get_2DFT_longrange_A_precalcbuf(c,m,n,mu,nu, - k_sph,J,lrhankel_recparts_buf); - *Bdest = qpms_trans_calculator_get_2DFT_longrange_B_precalcbuf(c,m,n,mu,nu, - k_sph,J,lrhankel_recparts_buf); - return 0; - } -} - -// Fourier transforms of the long-range parts of the translation coefficients -int qpms_trans_calculator_get_2DFT_longrange_AB_p(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, - qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t k_sph, - qpms_bessel_t J /* Only J=3 valid for now */, - qpms_l_t lrcutoff, unsigned kappa, double cv, double k0) { - int maxp = 2*c->lMax+1; // TODO this may not be needed here, n+nu+1 could be enough instead - complex double lrhankel_recpart[maxp * (maxp+1) / 2]; - return qpms_trans_calculator_get_2DFT_longrange_AB_buf_p(c, Adest, Bdest,m,n,mu,nu,k_sph, - J, lrcutoff, kappa, cv, k0, lrhankel_recpart); -} - -int qpms_trans_calculator_get_2DFT_longrange_AB_arrays_buf(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, - size_t deststride, size_t srcstride, - sph_t k_sph, qpms_bessel_t J /* must be 3 for now */, - qpms_l_t lrk_cutoff, unsigned kappa, double cv, double k0, - complex double *lrhankel_recparts_buf) { - assert(J == QPMS_HANKEL_PLUS); - assert(k_sph.theta == M_PI_2); -#if 0 - if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) { - for (size_t i = 0; i < c->nelem; ++i) - for (size_t j = 0; j < c->nelem; ++j) { - *(Adest + i*srcstride + j*deststride) = NAN+I*NAN; - *(Bdest + i*srcstride + j*deststride) = NAN+I*NAN; - } - // TODO warn? different return value? - return 0; - } -#endif - { - lrhankel_recpart_fill(lrhankel_recparts_buf, 2*c->lMax+1, - lrk_cutoff, c->hct, kappa, cv, k0, k_sph.r); - // if (qpms_sph_bessel_fill(J, 2*c->lMax+1, kdlj.r, bessel_buf)) abort(); - size_t desti = 0, srci = 0; - for (int n = 1; n <= c->lMax; ++n) for (int m = -n; m <= n; ++m) { - for (int nu = 1; nu <= c->lMax; ++nu) for (int mu = -nu; mu <= nu; ++mu) { - size_t assertindex = qpms_trans_calculator_index_mnmunu(c,m,n,mu,nu); - assert(assertindex == desti*c->nelem + srci); - *(Adest + deststride * desti + srcstride * srci) = - qpms_trans_calculator_get_2DFT_longrange_A_precalcbuf(c,m,n,mu,nu, - k_sph,J,lrhankel_recparts_buf); - *(Bdest + deststride * desti + srcstride * srci) = - qpms_trans_calculator_get_2DFT_longrange_B_precalcbuf(c,m,n,mu,nu, - k_sph,J,lrhankel_recparts_buf); - ++srci; - } - ++desti; - srci = 0; - } - return 0; - } -} - -// FIXME i get stack smashing error inside the following function if compiled with optimizations -int qpms_trans_calculator_get_2DFT_longrange_AB_arrays(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, - size_t deststride, size_t srcstride, - sph_t k_sph, qpms_bessel_t J /* Only J=3 valid for now */, - qpms_l_t lrcutoff, unsigned kappa, double cv, double k0) { - int maxp = 2*c->lMax+1; - complex double lrhankel_recpart[maxp * (maxp+1) / 2]; - return qpms_trans_calculator_get_2DFT_longrange_AB_arrays_buf(c, - Adest, Bdest, deststride, srcstride, k_sph, J, - lrcutoff, kappa, cv, k0, - lrhankel_recpart); -} - -#endif // LATTICESUMS_OLD - complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c, int m, int n, int mu, int nu, diff --git a/qpms/translations.h b/qpms/translations.h index f5974ca..893864f 100644 --- a/qpms/translations.h +++ b/qpms/translations.h @@ -40,19 +40,6 @@ #include "ewald.h" #endif -// TODO replace the xplicit "Taylor" functions with general, -// taking qpms_normalisation_t argument. -complex double qpms_trans_single_A_Taylor(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, - bool r_ge_d, qpms_bessel_t J); - -complex double qpms_trans_single_B_Taylor(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, - bool r_ge_d, qpms_bessel_t J); - -complex double qpms_trans_single_A_Taylor_ext(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, double kdlj_r, - double kdlj_th, double kdlj_phi, int r_ge_d, int J); - -complex double qpms_trans_single_B_Taylor_ext(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, double kdlj_r, - double kdlj_th, double kdlj_phi, int r_ge_d, int J); complex double qpms_trans_single_A(qpms_normalisation_t norm, qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, bool r_ge_d, qpms_bessel_t J); @@ -168,34 +155,6 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p( /// Workspace has to be at least 2 * c->neleme**2 long ); -#ifdef LATTICESUMS_OLD -// Short-range parts of the translation coefficients -int qpms_trans_calculator_get_shortrange_AB_p(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, - qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, - qpms_bessel_t J /* Only J=3 valid for now */, - qpms_l_t longrange_order_cutoff, unsigned kappa, double cc); -int qpms_trans_calculator_get_shortrange_AB_arrays(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, - size_t deststride, size_t srcstride, - sph_t kdlj, qpms_bessel_t J /* Only J=3 valid for now */, - qpms_l_t longrange_order_cutoff, unsigned kappa, double cc); - -// Fourier transforms of the long-range parts of the translation coefficients -int qpms_trans_calculator_get_2DFT_longrange_AB_p(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, - qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t k_sph, - qpms_bessel_t J /* Only J=3 valid for now */, - qpms_l_t longrange_order_cutoff, unsigned kappa, double cv, double k0); - -int qpms_trans_calculator_get_2DFT_longrange_AB_arrays(const qpms_trans_calculator *c, - complex double *Adest, complex double *Bdest, - size_t deststride, size_t srcstride, - sph_t k_sph, qpms_bessel_t J /* Only J=3 valid for now */, - qpms_l_t longrange_order_cutoff, unsigned kappa, double cv, double k0); -#endif // LATTICESUMS_OLD - - #ifdef LATTICESUMS32 // for the time being there are only those relatively low-level quick&dirty functions // according to what has been implemented from ewald.h;