Check cluster chto

cfastpt/utils_complex.c

@@ -65,59 +65,71 @@ Calculate g_l = exp( zln2 + lngamma( (l+nu)/2 + I*eta/2 ) - lngamma( (3+l-nu)/2
 
 void g_m_vals(double mu, double q_real, double *q_imag, double complex *gm, long N){
 	long i;
-	for(i=0; i<N; i++) {
-		gm[i] = cexp(lngamma_lanczos((mu+1.+q_real+I*q_imag[i])/2.) - lngamma_lanczos((mu+1.-q_real-I*q_imag[i])/2.) );		
-		// if(isnan(gl[i])) {printf("nan at l,nu,eta, = %lf %lg %lg %lg %lg\n", l,nu, eta[i], gamma_lanczos((l+z)/2.),gamma_lanczos((3.+l-z)/2.));exit(0);}
-	}
-}
-
-void gamma_ratios(double l, double nu, double *eta, double complex *gl, long N) {
-/* z = nu + I*eta
-Calculate g_l = exp( zln2 + lngamma( (l+nu)/2 + I*eta/2 ) - lngamma( (3+l-nu)/2 - I*eta/2 ) ) */
-	long i;
-	double complex z;
-	for(i=0; i<N; i++) {
-		z = nu+I*eta[i];
-		// gl[i] = cexp(z*log(2.) + clog(gamma_lanczos((l+z)/2.) ) - clog(gamma_lanczos((3.+l-z)/2.)));
-		gl[i] = cexp(lngamma_lanczos((l+z)/2.) - lngamma_lanczos((3.+l-z)/2.) );		
-		// if(isnan(gl[i])) {printf("nan at l,nu,eta, = %lf %lg %lg %lg %lg\n", l,nu, eta[i], gamma_lanczos((l+z)/2.),gamma_lanczos((3.+l-z)/2.));exit(0);}
-	}
-}
+	double complex plus;
+	double complex minus;
 
-void g_l(double l, double nu, double *eta, double complex *gl, long N) {
-/* z = nu + I*eta
-Calculate g_l = exp( zln2 + lngamma( (l+nu)/2 + I*eta/2 ) - lngamma( (3+l-nu)/2 - I*eta/2 ) ) */
-	long i;
-	double complex z;
 	for(i=0; i<N; i++) {
-		z = nu+I*eta[i];
-		// gl[i] = cexp(z*log(2.) + clog(gamma_lanczos((l+z)/2.) ) - clog(gamma_lanczos((3.+l-z)/2.)));
-		gl[i] = cexp(z*log(2.) + lngamma_lanczos((l+z)/2.) - lngamma_lanczos((3.+l-z)/2.) );		
+		plus = (mu+1.+q_real+I*q_imag[i])/2.;
+		minus= (mu+1.-q_real-I*q_imag[i])/2.;
+		if(fabs(q_imag[i])<200){
+			gm[i] = cexp(lngamma_lanczos(plus) - lngamma_lanczos(minus));
+		}else{
+			gm[i] = cexp((plus-0.5)*clog(plus) - (minus-0.5)*clog(minus) - q_real - I*q_imag[i] \
+					+1./12 *(1./plus - 1./minus) \
+					+1./360.*(1./cpow(minus,3) - 1./cpow(plus,3)) \
+					+1./1260*(1./cpow(plus,5) - 1./cpow(minus,5)));
+		}
 		// if(isnan(gl[i])) {printf("nan at l,nu,eta, = %lf %lg %lg %lg %lg\n", l,nu, eta[i], gamma_lanczos((l+z)/2.),gamma_lanczos((3.+l-z)/2.));exit(0);}
 	}
 }
 
-void g_l_1(double l, double nu, double *eta, double complex *gl1, long N) {
-/* z = nu + I*eta
-Calculate g_l_1 = exp(zln2 + lngamma( (l+nu-1)/2 + I*eta/2 ) - lngamma( (4+l-nu)/2 - I*eta/2 ) ) */
-	long i;
-	double complex z;
-	for(i=0; i<N; i++) {
-		z = nu+I*eta[i];
-		gl1[i] = -(z-1.)* cexp((z-1.)*log(2.) + lngamma_lanczos((l+z-1.)/2.) - lngamma_lanczos((4.+l-z)/2.));
-	}
-}
-
-void g_l_2(double l, double nu, double *eta, double complex *gl2, long N) {
-/* z = nu + I*eta
-Calculate g_l_1 = exp(zln2 + lngamma( (l+nu-2)/2 + I*eta/2 ) - lngamma( (5+l-nu)/2 - I*eta/2 ) ) */
-	long i;
-	double complex z;
-	for(i=0; i<N; i++) {
-		z = nu+I*eta[i];
-		gl2[i] = (z-1.)* (z-2.)* cexp((z-2.)*log(2.) + lngamma_lanczos((l+z-2.)/2.) - lngamma_lanczos((5.+l-z)/2.));
-	}
-}
+// void gamma_ratios(double l, double nu, double *eta, double complex *gl, long N) {
+// /* z = nu + I*eta
+// Calculate g_l = exp( zln2 + lngamma( (l+nu)/2 + I*eta/2 ) - lngamma( (3+l-nu)/2 - I*eta/2 ) ) */
+// 	long i;
+// 	double complex z;
+// 	for(i=0; i<N; i++) {
+// 		z = nu+I*eta[i];
+// 		// gl[i] = cexp(z*log(2.) + clog(gamma_lanczos((l+z)/2.) ) - clog(gamma_lanczos((3.+l-z)/2.)));
+// 		gl[i] = cexp(lngamma_lanczos((l+z)/2.) - lngamma_lanczos((3.+l-z)/2.) );		
+// 		// if(isnan(gl[i])) {printf("nan at l,nu,eta, = %lf %lg %lg %lg %lg\n", l,nu, eta[i], gamma_lanczos((l+z)/2.),gamma_lanczos((3.+l-z)/2.));exit(0);}
+// 	}
+// }
+
+// void g_l(double l, double nu, double *eta, double complex *gl, long N) {
+// /* z = nu + I*eta
+// Calculate g_l = exp( zln2 + lngamma( (l+nu)/2 + I*eta/2 ) - lngamma( (3+l-nu)/2 - I*eta/2 ) ) */
+// 	long i;
+// 	double complex z;
+// 	for(i=0; i<N; i++) {
+// 		z = nu+I*eta[i];
+// 		// gl[i] = cexp(z*log(2.) + clog(gamma_lanczos((l+z)/2.) ) - clog(gamma_lanczos((3.+l-z)/2.)));
+// 		gl[i] = cexp(z*log(2.) + lngamma_lanczos((l+z)/2.) - lngamma_lanczos((3.+l-z)/2.) );		
+// 		// if(isnan(gl[i])) {printf("nan at l,nu,eta, = %lf %lg %lg %lg %lg\n", l,nu, eta[i], gamma_lanczos((l+z)/2.),gamma_lanczos((3.+l-z)/2.));exit(0);}
+// 	}
+// }
+
+// void g_l_1(double l, double nu, double *eta, double complex *gl1, long N) {
+// /* z = nu + I*eta
+// Calculate g_l_1 = exp(zln2 + lngamma( (l+nu-1)/2 + I*eta/2 ) - lngamma( (4+l-nu)/2 - I*eta/2 ) ) */
+// 	long i;
+// 	double complex z;
+// 	for(i=0; i<N; i++) {
+// 		z = nu+I*eta[i];
+// 		gl1[i] = -(z-1.)* cexp((z-1.)*log(2.) + lngamma_lanczos((l+z-1.)/2.) - lngamma_lanczos((4.+l-z)/2.));
+// 	}
+// }
+
+// void g_l_2(double l, double nu, double *eta, double complex *gl2, long N) {
+// /* z = nu + I*eta
+// Calculate g_l_1 = exp(zln2 + lngamma( (l+nu-2)/2 + I*eta/2 ) - lngamma( (5+l-nu)/2 - I*eta/2 ) ) */
+// 	long i;
+// 	double complex z;
+// 	for(i=0; i<N; i++) {
+// 		z = nu+I*eta[i];
+// 		gl2[i] = (z-1.)* (z-2.)* cexp((z-2.)*log(2.) + lngamma_lanczos((l+z-2.)/2.) - lngamma_lanczos((5.+l-z)/2.));
+// 	}
+// }
 
 void c_window(double complex *out, double c_window_width, long halfN) {
 	// 'out' is (halfN+1) complex array

cfastpt/utils_complex.h

@@ -1,9 +1,6 @@
 #include <complex.h>
 #include <fftw3.h>
 void f_z(double z_real, double *z_imag, double complex *fz, long N);
-void g_l(double l, double nu, double *eta, double complex *gl, long N);
-void g_l_1(double l, double nu, double *eta, double complex *gl1, long N);
-void g_l_2(double l, double nu, double *eta, double complex *gl2, long N);
 
 void c_window(double complex *out, double c_window_width, long halfN);
 
@@ -12,9 +9,6 @@ void c_window(double complex *out, double c_window_width, long halfN);
 double complex gamma_lanczos(double complex z);
 double complex lngamma_lanczos(double complex z);
 
-void gamma_ratios(double l, double nu, double *eta, double complex *gl, long N);
-
-
 void g_m_vals(double mu, double q_real, double *q_imag, double complex *gm, long N);
 
 

cfftlog/cfftlog.c

@@ -309,3 +309,94 @@ void cfftlog_ells_increment(double *x, double *fx, long N, config *config, int*
 	free(out_ifft);
 	free(fb);
 }
+
+
+// same ell but multiple integrands with same x-range
+void cfftlog_multiple(double *x, double **fx, long N, long Nf, config *config, int ell, double *y, double **Fy) {
+
+	long N_original = N;
+	long N_pad = config->N_pad;
+	N += 2*N_pad;
+
+	if(N % 2) {printf("Please use even number of x !\n"); exit(0);}
+	long halfN = N/2;
+
+	double x0, y0;
+	x0 = x[0];
+
+	double dlnx;
+	dlnx = log(x[1]/x0);
+
+	// Only calculate the m>=0 part
+	double eta_m[halfN+1];
+	long i, j;
+	for(i=0; i<=halfN; i++) {eta_m[i] = 2*M_PI / dlnx / N * i;}
+
+	double complex gl[halfN+1];
+
+	switch(config->derivative) {
+		case 0: g_l_cfft((double)ell, config->nu, eta_m, gl, halfN+1); break;
+		case 1: g_l_1_cfft((double)ell, config->nu, eta_m, gl, halfN+1); break;
+		case 2: g_l_2_cfft((double)ell, config->nu, eta_m, gl, halfN+1); break;
+		default: printf("Integral Not Supported! Please choose config->derivative from [0,1,2].\n");
+	}
+	// printf("g2[0]: %.15e+I*(%.15e)\n", creal(g2[0]),cimag(g2[0]));
+
+	// calculate y arrays
+	for(i=0; i<N_original; i++) {y[i] = (ell+1.) / x[N_original-1-i];}
+	y0 = y[0];
+
+	// biased input func
+	double **fb;
+	fb = malloc(Nf* sizeof(double*));
+	for(j=0; j<Nf; j++) {
+		fb[j] = malloc(N* sizeof(double));
+
+		for(i=0; i<N_pad; i++) {
+			fb[j][i] = 0.;
+			fb[j][N-1-i] = 0.;
+		}
+		for(i=N_pad; i<N_pad+N_original; i++) {
+			fb[j][i] = fx[j][i-N_pad] / pow(x[i-N_pad], config->nu) ;
+		}
+	}
+
+	fftw_complex *out;
+	fftw_plan plan_forward, plan_backward;
+	out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (halfN+1) );
+	double *in_fb;
+	in_fb = malloc(sizeof(double) * N );
+	plan_forward = fftw_plan_dft_r2c_1d(N, in_fb, out, FFTW_ESTIMATE);
+
+	double *out_ifft;
+	out_ifft = malloc(sizeof(double) * N );
+	plan_backward = fftw_plan_dft_c2r_1d(N, out, out_ifft, FFTW_ESTIMATE);
+
+	for(j=0; j<Nf; j++){
+		for(i=0; i<N; i++) { in_fb[i] = fb[j][i];}
+		fftw_execute(plan_forward);
+
+		c_window_cfft(out, config->c_window_width, halfN);
+		// printf("out[1]:%.15e+i*(%.15e)\n", creal(out[1]), cimag(out[1]));
+
+		for(i=0; i<=halfN; i++) {
+			out[i] *= cpow(x0*y0/exp(2*N_pad*dlnx), -I*eta_m[i]) * gl[i] ;
+			out[i] = conj(out[i]);
+		}
+		// printf("out[1]:%.15e+i*(%.15e)\n", creal(out[1]), cimag(out[1]));
+
+		fftw_execute(plan_backward);
+
+		for(i=0; i<N_original; i++) {
+			Fy[j][i] = out_ifft[i-N_pad] * sqrt(M_PI) / (4.*N * pow(y[i], config->nu));
+		}
+	}
+
+	fftw_destroy_plan(plan_forward);
+	fftw_destroy_plan(plan_backward);
+	fftw_free(out);
+	free(out_ifft);
+	free(in_fb);
+	for(j=0; j<Nf; j++) { free(fb[j]); }
+	free(fb);
+}

cfftlog/cfftlog.h

@@ -11,4 +11,6 @@ void cfftlog(double *x, double *fx, long N, config *config, int ell, double *y,
 
 void cfftlog_ells(double *x, double *fx, long N, config *config, int* ell, long Nell, double **y, double **Fy);
 
-void cfftlog_ells_increment(double *x, double *fx, long N, config *config, int* ell, long Nell, double **y, double **Fy);
+void cfftlog_ells_increment(double *x, double *fx, long N, config *config, int* ell, long Nell, double **y, double **Fy);
+
+void cfftlog_multiple(double *x, double **fx, long N, long Nf, config *config, int ell, double *y, double **Fy);

cfftlog/utils_complex.c

@@ -48,15 +48,31 @@ double complex lngamma_lanczos_cfft(double complex z) {
 	return log(2*M_PI) /2.  + (z+0.5)*clog(t) -t + clog(x);
 }
 
+double complex ln_g_m_vals_cfft(double mu, double complex q) {
+/* similar routine as python version.
+use asymptotic expansion for large |mu+q| */
+	double complex asym_plus = (mu+1+ q)/2.;
+	double complex asym_minus= (mu+1- q)/2.;
+
+	return (asym_plus-0.5)*clog(asym_plus) - (asym_minus-0.5)*clog(asym_minus) - q \
+		+1./12 *(1./asym_plus - 1./asym_minus) \
+		+1./360.*(1./cpow(asym_minus,3) - 1./cpow(asym_plus,3)) \
+		+1./1260*(1./cpow(asym_plus,5) - 1./cpow(asym_minus,5));
+}
+
 void g_l_cfft(double l, double nu, double *eta, double complex *gl, long N) {
 /* z = nu + I*eta
 Calculate g_l = exp( zln2 + lngamma( (l+nu)/2 + I*eta/2 ) - lngamma( (3+l-nu)/2 - I*eta/2 ) ) */
 	long i;
 	double complex z;
 	for(i=0; i<N; i++) {
 		z = nu+I*eta[i];
-		// gl[i] = cexp(z*log(2.) + clog(gamma_lanczos((l+z)/2.) ) - clog(gamma_lanczos((3.+l-z)/2.)));
-		gl[i] = cexp(z*log(2.) + lngamma_lanczos_cfft((l+z)/2.) - lngamma_lanczos_cfft((3.+l-z)/2.) );		
+		if(l+fabs(eta[i])<200){
+			// gl[i] = cexp(z*log(2.) + clog(gamma_lanczos((l+z)/2.) ) - clog(gamma_lanczos((3.+l-z)/2.)));
+			gl[i] = cexp(z*log(2.) + lngamma_lanczos_cfft((l+z)/2.) - lngamma_lanczos_cfft((3.+l-z)/2.) );	
+		}else{
+			gl[i] = cexp(z*log(2.) + ln_g_m_vals_cfft(l+0.5, z-1.5));
+		}
 		// if(isnan(gl[i])) {printf("nan at l,nu,eta, = %lf %lg %lg %lg %lg\n", l,nu, eta[i], gamma_lanczos((l+z)/2.),gamma_lanczos((3.+l-z)/2.));exit(0);}
 	}
 }
@@ -68,7 +84,11 @@ Calculate g_l_1 = exp(zln2 + lngamma( (l+nu-1)/2 + I*eta/2 ) - lngamma( (4+l-nu)
 	double complex z;
 	for(i=0; i<N; i++) {
 		z = nu+I*eta[i];
-		gl1[i] = -(z-1.)* cexp((z-1.)*log(2.) + lngamma_lanczos_cfft((l+z-1.)/2.) - lngamma_lanczos_cfft((4.+l-z)/2.));
+		if(l+fabs(eta[i])<200){
+			gl1[i] = -(z-1.)* cexp((z-1.)*log(2.) + lngamma_lanczos_cfft((l+z-1.)/2.) - lngamma_lanczos_cfft((4.+l-z)/2.));
+		}else{
+			gl1[i] = -(z-1.)* cexp((z-1.)*log(2.) + ln_g_m_vals_cfft(l+0.5, z-2.5));
+		}
 	}
 }
 
@@ -79,7 +99,11 @@ Calculate g_l_1 = exp(zln2 + lngamma( (l+nu-2)/2 + I*eta/2 ) - lngamma( (5+l-nu)
 	double complex z;
 	for(i=0; i<N; i++) {
 		z = nu+I*eta[i];
-		gl2[i] = (z-1.)* (z-2.)* cexp((z-2.)*log(2.) + lngamma_lanczos_cfft((l+z-2.)/2.) - lngamma_lanczos_cfft((5.+l-z)/2.));
+		if(l+fabs(eta[i])<200){
+			gl2[i] = (z-1.)* (z-2.)* cexp((z-2.)*log(2.) + lngamma_lanczos_cfft((l+z-2.)/2.) - lngamma_lanczos_cfft((5.+l-z)/2.));
+		}else{
+			gl2[i] = (z-1.)* (z-2.)* cexp((z-2.)*log(2.) + ln_g_m_vals_cfft(l+0.5, z-3.5));
+		}
 	}
 }
 

cfftlog/utils_complex.h

@@ -10,4 +10,5 @@ void c_window_cfft(double complex *out, double c_window_width, long halfN);
 // void resample_fourier_gauss(double *k, double *fk, config *config);
 
 double complex gamma_lanczos_cfft(double complex z);
-double complex lngamma_lanczos_cfft(double complex z);
+double complex lngamma_lanczos_cfft(double complex z);
+double complex ln_g_m_vals_cfft(double mu, double complex q);