44 #define SYMBOL_ABEL_POISSON(k,h) (pow(h,k))
47 #define SYMBOL_SINGULARITY(k,h) ((2.0/(2*k+1))*pow(h,k))
52 #define KT_ABEL_POISSON (0)
54 #define KT_SINGULARITY (1)
56 #define KT_LOC_SUPP (2)
58 #define KT_GAUSSIAN (3)
61 enum pvalue {NO = 0, YES = 1, BOTH = 2};
63 static inline int scaled_modified_bessel_i_series(
const R x,
const R alpha,
64 const int nb,
const int ize, R *b)
66 const R enmten = K(4.0)*nfft_float_property(NFFT_R_MIN);
67 R tempa = K(1.0), empal = K(1.0) + alpha, halfx = K(0.0), tempb = K(0.0);
74 tempa = POW(halfx, alpha)/TGAMMA(empal);
79 if (K(1.0) < x + K(1.0))
82 b[0] = tempa + tempa*tempb/empal;
84 if (x != K(0.0) && b[0] == K(0.0))
92 R tempc = halfx, tover = (enmten + enmten)/x;
97 for (n = 1; n < nb; n++)
103 if (tempa <= tover*empal)
106 b[n] = tempa + tempa*tempb/empal;
108 if (b[n] == K(0.0) && n < ncalc)
113 for (n = 1; n < nb; n++)
119 static inline void scaled_modified_bessel_i_normalize(
const R x,
120 const R alpha,
const int nb,
const int ize, R *b,
const R sum_)
122 const R enmten = K(4.0)*nfft_float_property(NFFT_R_MIN);
128 sum = sum * TGAMMA(K(1.0) + alpha) * POW(x/K(2.0), -alpha);
138 for (n = 1; n <= nb; n++)
194 static int smbi(
const R x,
const R alpha,
const int nb,
const int ize, R *b)
216 const int nsig = MANT_DIG + 2;
217 const R enten = nfft_float_property(NFFT_R_MAX);
218 const R ensig = POW(K(10.0),(R)nsig);
219 const R rtnsig = POW(K(10.0),-CEIL((R)nsig/K(4.0)));
220 const R xlarge = K(1E4);
221 const R exparg = FLOOR(LOG(POW(K(R_RADIX),K(DBL_MAX_EXP-1))));
224 int l, n, nend, magx, nbmx, ncalc, nstart;
225 R p, em, en, sum, pold, test, empal, tempa, tempb, tempc, psave, plast, tover,
228 magx = LRINT(FLOOR(x));
231 if ( nb <= 0 || x < K(0.0) || alpha < K(0.0) || K(1.0) <= alpha
232 || ((ize != 1 || exparg < x) && (ize != 2 || xlarge < x)))
233 return (MIN(nb,0) - 1);
237 return scaled_modified_bessel_i_series(x,alpha,nb,ize,b);
245 en = (R) (n+n) + (alpha+alpha);
250 test = ensig + ensig;
252 if ((5*nsig) < (magx << 1))
255 test /= POW(K(1.585),(R)magx);
264 for (n = nstart; n <= nend; n++)
269 p = en*plast/x + pold;
287 p = en*plast/x + pold;
288 }
while (p <= K(1.0));
293 test = pold*plast*(K(0.5) - K(0.5)/(tempb * tempb))/ensig;
299 for (ncalc = nstart; ncalc <= nend; ncalc++)
303 psave = en*psavel/x + pold;
304 if (test < psave * psavel)
314 en = (R) (n+n) + (alpha+alpha);
317 test = FMAX(test,SQRT(plast*ensig)*SQRT(p+p));
327 p = en*plast/x + pold;
338 emp2al = em - K(1.0) + (alpha+alpha);
339 sum = tempa*empal*emp2al/em;
347 for (l = 1; l <= nend; ++l)
355 for (l = 1; l <= nend; ++l)
361 tempa = en*tempb/x + tempc;
372 sum = (sum + tempa*empal)*emp2al/em;
381 sum = sum + sum + tempa;
382 scaled_modified_bessel_i_normalize(x,alpha,nb,ize,b,sum);
389 b[n-1] = en*tempa/x + tempb;
393 sum = sum + sum + b[0];
394 scaled_modified_bessel_i_normalize(x,alpha,nb,ize,b,sum);
405 sum = (sum + b[n-1]*empal)*emp2al/em;
413 for (l = 1; l <= nend; ++l)
417 b[n-1] = en*b[n]/x + b[n+1];
425 sum = (sum + b[n-1]*empal)*emp2al/em;
430 b[0] = K(2.0)*empal*b[1]/x + b[2];
431 sum = sum + sum + b[0];
433 scaled_modified_bessel_i_normalize(x,alpha,nb,ize,b,sum);
451 static inline double innerProduct(
const double phi1,
const double theta1,
452 const double phi2,
const double theta2)
454 double pi2theta1 = K2PI*theta1, pi2theta2 = K2PI*theta2;
455 return (cos(pi2theta1)*cos(pi2theta2)
456 + sin(pi2theta1)*sin(pi2theta2)*cos(K2PI*(phi1-phi2)));
472 return (1.0/(K4PI))*((1.0-h)*(1.0+h))/pow(sqrt(1.0-2.0*h*x+h*h),3.0);
488 return (1.0/(K2PI))/sqrt(1.0-2.0*h*x+h*h);
507 return (x<=h)?(0.0):(pow((x-h),lambda));
524 return exp(2.0*sigma*(x-1.0));
537 int main (
int argc,
char **argv)
586 fftw_complex *prec = NULL;
601 fscanf(stdin,
"testcases=%d\n",&tc_max);
602 fprintf(stdout,
"%d\n",tc_max);
605 for (tc = 0; tc < tc_max; tc++)
608 fscanf(stdin,
"nfsft=%d\n",&use_nfsft);
609 fprintf(stdout,
"%d\n",use_nfsft);
613 fscanf(stdin,
"nfft=%d\n",&use_nfft);
614 fprintf(stdout,
"%d\n",use_nfft);
618 fscanf(stdin,
"cutoff=%d\n",&cutoff);
619 fprintf(stdout,
"%d\n",cutoff);
628 fscanf(stdin,
"fpt=%d\n",&use_fpt);
629 fprintf(stdout,
"%d\n",use_fpt);
631 fscanf(stdin,
"threshold=%lf\n",&threshold);
632 fprintf(stdout,
"%lf\n",threshold);
639 threshold = 1000000000000.0;
655 fscanf(stdin,
"kernel=%d\n",&kt);
656 fprintf(stdout,
"%d\n",kt);
659 fscanf(stdin,
"parameter_sets=%d\n",&ip_max);
660 fprintf(stdout,
"%d\n",ip_max);
663 p = (
double**)
nfft_malloc(ip_max*
sizeof(
double*));
668 fscanf(stdin,
"parameters=%d\n",&ipp_max);
669 fprintf(stdout,
"%d\n",ipp_max);
671 for (ip = 0; ip < ip_max; ip++)
674 p[ip] = (
double*)
nfft_malloc(ipp_max*
sizeof(
double));
677 for (ipp = 0; ipp < ipp_max; ipp++)
680 fscanf(stdin,
"%lf\n",&p[ip][ipp]);
681 fprintf(stdout,
"%lf\n",p[ip][ipp]);
686 fscanf(stdin,
"bandwidths=%d\n",&im_max);
687 fprintf(stdout,
"%d\n",im_max);
691 for (im = 0; im < im_max; im++)
694 fscanf(stdin,
"%d\n",&m[im]);
695 fprintf(stdout,
"%d\n",m[im]);
696 m_max = MAX(m_max,m[im]);
700 fscanf(stdin,
"node_sets=%d\n",&ild_max);
701 fprintf(stdout,
"%d\n",ild_max);
705 for (ild = 0; ild < ild_max; ild++)
711 fscanf(stdin,
"L=%d ",&ld[ild][0]);
712 fprintf(stdout,
"%d\n",ld[ild][0]);
713 l_max = MAX(l_max,ld[ild][0]);
716 fscanf(stdin,
"D=%d ",&ld[ild][1]);
717 fprintf(stdout,
"%d\n",ld[ild][1]);
718 d_max = MAX(d_max,ld[ild][1]);
721 fscanf(stdin,
"compare=%d ",&ld[ild][2]);
722 fprintf(stdout,
"%d\n",ld[ild][2]);
725 if (ld[ild][2] == YES)
728 fscanf(stdin,
"precomputed=%d\n",&ld[ild][3]);
729 fprintf(stdout,
"%d\n",ld[ild][3]);
733 fscanf(stdin,
"repetitions=%d\n",&ld[ild][4]);
734 fprintf(stdout,
"%d\n",ld[ild][4]);
737 if (ld[ild][3] == YES)
740 ld_max_prec = MAX(ld_max_prec,ld[ild][0]*ld[ild][1]);
742 l_max_prec = MAX(l_max_prec,ld[ild][0]);
755 b = (fftw_complex*)
nfft_malloc(l_max*
sizeof(fftw_complex));
756 eta = (
double*)
nfft_malloc(2*l_max*
sizeof(
double));
757 f_hat = (fftw_complex*)
nfft_malloc(NFSFT_F_HAT_SIZE(m_max)*
sizeof(fftw_complex));
758 a = (fftw_complex*)
nfft_malloc((m_max+1)*
sizeof(fftw_complex));
759 xi = (
double*)
nfft_malloc(2*d_max*
sizeof(
double));
760 f_m = (fftw_complex*)
nfft_malloc(d_max*
sizeof(fftw_complex));
761 f = (fftw_complex*)
nfft_malloc(d_max*
sizeof(fftw_complex));
764 if (precompute == YES)
766 prec = (fftw_complex*)
nfft_malloc(ld_max_prec*
sizeof(fftw_complex));
770 for (l = 0; l < l_max; l++)
772 b[l] = (((double)rand())/RAND_MAX) - 0.5;
773 eta[2*l] = (((double)rand())/RAND_MAX) - 0.5;
774 eta[2*l+1] = acos(2.0*(((
double)rand())/RAND_MAX) - 1.0)/(K2PI);
778 for (d = 0; d < d_max; d++)
780 xi[2*d] = (((double)rand())/RAND_MAX) - 0.5;
781 xi[2*d+1] = acos(2.0*(((
double)rand())/RAND_MAX) - 1.0)/(K2PI);
785 nfsft_precompute(m_max,threshold,
786 ((use_nfsft==NO)?(NFSFT_NO_FAST_ALGORITHM):(0U)), 0U);
789 for (ip = 0; ip < ip_max; ip++)
796 for (k = 0; k <= m_max; k++)
797 a[k] = SYMBOL_ABEL_POISSON(k,p[ip][0]);
803 for (k = 0; k <= m_max; k++)
804 a[k] = SYMBOL_SINGULARITY(k,p[ip][0]);
812 a[1] = ((p[ip][1]+1+p[ip][0])/(p[ip][1]+2.0))*a[0];
813 for (k = 2; k <= m_max; k++)
814 a[k] = (1.0/(k+p[ip][1]+1))*((2*k-1)*p[ip][0]*a[k-1] -
815 (k-p[ip][1]-2)*a[k-2]);
820 steed = (
double*)
nfft_malloc((m_max+1)*
sizeof(double));
821 smbi(2.0*p[ip][0],0.5,m_max+1,2,steed);
822 for (k = 0; k <= m_max; k++)
823 a[k] = K2PI*(sqrt(KPI/p[ip][0]))*steed[k];
830 for (k = 0; k <= m_max; k++)
831 a[k] *= (2*k+1)/(K4PI);
834 for (ild = 0; ild < ild_max; ild++)
837 if (ld[ild][2] != NO)
841 if (ld[ild][3] != NO)
846 rinc = l_max_prec-ld[ild][0];
849 for (d = 0; d < ld[ild][1]; d++)
852 for (l = 0; l < ld[ild][0]; l++)
856 temp =
innerProduct(eta[2*l],eta[2*l+1],xi[2*d],xi[2*d+1]);
896 for (i = 0; i < ld[ild][4]; i++)
902 rinc = l_max_prec-ld[ild][0];
910 constant = ((p[ip][1]+1)/(K2PI*pow(1-p[ip][0],p[ip][1]+1)));
913 for (d = 0; d < ld[ild][1]; d++)
919 for (l = 0; l < ld[ild][0]; l++)
920 f[d] += b[l]*(*ptr++);
932 for (d = 0; d < ld[ild][1]; d++)
938 for (l = 0; l < ld[ild][0]; l++)
939 f[d] += b[l]*(*ptr++);
952 t_dp = t_dp/((double)ld[ild][4]);
967 for (i = 0; i < ld[ild][4]; i++)
975 for (d = 0; d < ld[ild][1]; d++)
981 for (l = 0; l < ld[ild][0]; l++)
985 temp =
innerProduct(eta[2*l],eta[2*l+1],xi[2*d],xi[2*d+1]);
996 for (d = 0; d < ld[ild][1]; d++)
1002 for (l = 0; l < ld[ild][0]; l++)
1006 temp =
innerProduct(eta[2*l],eta[2*l+1],xi[2*d],xi[2*d+1]);
1017 constant = ((p[ip][1]+1)/(K2PI*pow(1-p[ip][0],p[ip][1]+1)));
1020 for (d = 0; d < ld[ild][1]; d++)
1026 for (l = 0; l < ld[ild][0]; l++)
1030 temp =
innerProduct(eta[2*l],eta[2*l+1],xi[2*d],xi[2*d+1]);
1044 for (d = 0; d < ld[ild][1]; d++)
1050 for (l = 0; l < ld[ild][0]; l++)
1054 temp =
innerProduct(eta[2*l],eta[2*l+1],xi[2*d],xi[2*d+1]);
1068 t_d = t_d/((double)ld[ild][4]);
1085 for (im = 0; im < im_max; im++)
1088 nfsft_init_guru(&plan_adjoint, m[im],ld[ild][0],
1089 ((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
1090 ((use_fpt!=0)?(0U):(NFSFT_USE_DPT)),
1091 PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
1092 FFT_OUT_OF_PLACE, cutoff);
1093 nfsft_init_guru(&plan,m[im],ld[ild][1],
1094 ((use_nfft!=0)?(0U):(NFSFT_USE_NDFT)) |
1095 ((use_fpt!=0)?(0U):(NFSFT_USE_DPT)),
1096 PRE_PHI_HUT | PRE_PSI | FFTW_INIT |
1099 plan_adjoint.
f_hat = f_hat;
1100 plan_adjoint.
x = eta;
1105 nfsft_precompute_x(&plan_adjoint);
1106 nfsft_precompute_x(&plan);
1109 if (use_nfsft == BOTH)
1118 for (i = 0; i < ld[ild][4]; i++)
1122 nfsft_adjoint_direct(&plan_adjoint);
1125 for (k = 0; k <= m[im]; k++)
1126 for (n = -k; n <= k; n++)
1127 f_hat[NFSFT_INDEX(k,n,&plan_adjoint)] *= a[k];
1130 nfsft_trafo_direct(&plan);
1139 t_fd = t_fd/((double)ld[ild][4]);
1142 if (ld[ild][2] != NO)
1145 err_fd =
X(error_l_infty_1_complex)(f, f_m, ld[ild][1], b,
1151 if (use_nfsft != NO)
1167 for (i = 0; i < ld[ild][4]; i++)
1170 if (use_nfsft != NO)
1173 nfsft_adjoint(&plan_adjoint);
1178 nfsft_adjoint_direct(&plan_adjoint);
1182 for (k = 0; k <= m[im]; k++)
1183 for (n = -k; n <= k; n++)
1184 f_hat[NFSFT_INDEX(k,n,&plan_adjoint)] *= a[k];
1187 if (use_nfsft != NO)
1195 nfsft_trafo_direct(&plan);
1202 if (use_nfsft != NO)
1208 if (use_nfsft != NO)
1211 t_f = t_f/((double)ld[ild][4]);
1216 t_fd = t_fd/((double)ld[ild][4]);
1220 if (ld[ild][2] != NO)
1223 if (use_nfsft != NO)
1226 err_f =
X(error_l_infty_1_complex)(f, f_m, ld[ild][1], b,
1232 err_fd =
X(error_l_infty_1_complex)(f, f_m, ld[ild][1], b,
1238 fprintf(stdout,
"%e\n%e\n%e\n%e\n%e\n%e\n\n",t_d,t_dp,t_fd,t_f,err_fd,
1242 nfsft_finalize(&plan_adjoint);
1243 nfsft_finalize(&plan);
1254 if (precompute == YES)
1269 for (ild = 0; ild < ild_max; ild++)
1277 for (ip = 0; ip < ip_max; ip++)
1283 return EXIT_SUCCESS;
double * x
the nodes for ,
static double gaussianKernel(const double x, const double sigma)
Evaluates the spherical Gaussian kernel at a node .
#define KT_SINGULARITY
Singularity kernel.
static double poissonKernel(const double x, const double h)
Evaluates the Poisson kernel at a node .
fftw_complex * f_hat
Fourier coefficients.
R nfft_elapsed_seconds(ticks t1, ticks t0)
Return number of elapsed seconds between two time points.
static int smbi(const R x, const R alpha, const int nb, const int ize, R *b)
Calculates the modified bessel function , possibly scaled by , for real non-negative with ...
static double innerProduct(const double phi1, const double theta1, const double phi2, const double theta2)
Computes the standard inner product between two vectors on the unit sphere given in spherical coord...
data structure for an NFSFT (nonequispaced fast spherical Fourier transform) plan with double precisi...
static double locallySupportedKernel(const double x, const double h, const double lambda)
Evaluates the locally supported kernel at a node .
static double singularityKernel(const double x, const double h)
Evaluates the singularity kernel at a node .
#define X(name)
Include header for C99 complex datatype.
#define KT_ABEL_POISSON
Abel-Poisson kernel.
void * nfft_malloc(size_t n)
pvalue
Enumeration type for yes/no/both-type parameters.
#define KT_LOC_SUPP
Locally supported kernel.
#define KT_GAUSSIAN
Gaussian kernel.