42 #define SYMBOL_ABEL_POISSON(k,h) (pow(h,k))
45 #define SYMBOL_SINGULARITY(k,h) ((2.0/(2*k+1))*pow(h,k))
50 #define KT_ABEL_POISSON (0)
52 #define KT_SINGULARITY (1)
54 #define KT_LOC_SUPP (2)
56 #define KT_GAUSSIAN (3)
59 enum pvalue {NO = 0, YES = 1, BOTH = 2};
61 static inline int scaled_modified_bessel_i_series(
const R x,
const R alpha,
62 const int nb,
const int ize, R *b)
64 const R enmten = K(4.0)*nfft_float_property(NFFT_R_MIN);
65 R tempa = K(1.0), empal = K(1.0) + alpha, halfx = K(0.0), tempb = K(0.0);
72 tempa = POW(halfx, alpha)/TGAMMA(empal);
77 if (K(1.0) < x + K(1.0))
80 b[0] = tempa + tempa*tempb/empal;
82 if (x != K(0.0) && b[0] == K(0.0))
90 R tempc = halfx, tover = (enmten + enmten)/x;
95 for (n = 1; n < nb; n++)
101 if (tempa <= tover*empal)
104 b[n] = tempa + tempa*tempb/empal;
106 if (b[n] == K(0.0) && n < ncalc)
111 for (n = 1; n < nb; n++)
117 static inline void scaled_modified_bessel_i_normalize(
const R x,
118 const R alpha,
const int nb,
const int ize, R *b,
const R sum_)
120 const R enmten = K(4.0)*nfft_float_property(NFFT_R_MIN);
126 sum = sum * TGAMMA(K(1.0) + alpha) * POW(x/K(2.0), -alpha);
136 for (n = 1; n <= nb; n++)
192 static int smbi(
const R x,
const R alpha,
const int nb,
const int ize, R *b)
214 const int nsig = MANT_DIG + 2;
215 const R enten = nfft_float_property(NFFT_R_MAX);
216 const R ensig = POW(K(10.0),(R)nsig);
217 const R rtnsig = POW(K(10.0),-CEIL((R)nsig/K(4.0)));
218 const R xlarge = K(1E4);
219 const R exparg = FLOOR(LOG(POW(K(R_RADIX),K(DBL_MAX_EXP-1))));
222 int l, n, nend, magx, nbmx, ncalc, nstart;
223 R p, em, en, sum, pold, test, empal, tempa, tempb, tempc, psave, plast, tover,
226 magx = LRINT(FLOOR(x));
229 if ( nb <= 0 || x < K(0.0) || alpha < K(0.0) || K(1.0) <= alpha
230 || ((ize != 1 || exparg < x) && (ize != 2 || xlarge < x)))
231 return (MIN(nb,0) - 1);
235 return scaled_modified_bessel_i_series(x,alpha,nb,ize,b);
243 en = (R) (n+n) + (alpha+alpha);
248 test = ensig + ensig;
250 if ((5*nsig) < (magx << 1))
253 test /= POW(K(1.585),(R)magx);
262 for (n = nstart; n <= nend; n++)
267 p = en*plast/x + pold;
285 p = en*plast/x + pold;
286 }
while (p <= K(1.0));
291 test = pold*plast*(K(0.5) - K(0.5)/(tempb * tempb))/ensig;
297 for (ncalc = nstart; ncalc <= nend; ncalc++)
301 psave = en*psavel/x + pold;
302 if (test < psave * psavel)
312 en = (R) (n+n) + (alpha+alpha);
315 test = FMAX(test,SQRT(plast*ensig)*SQRT(p+p));
325 p = en*plast/x + pold;
336 emp2al = em - K(1.0) + (alpha+alpha);
337 sum = tempa*empal*emp2al/em;
345 for (l = 1; l <= nend; ++l)
353 for (l = 1; l <= nend; ++l)
359 tempa = en*tempb/x + tempc;
370 sum = (sum + tempa*empal)*emp2al/em;
379 sum = sum + sum + tempa;
380 scaled_modified_bessel_i_normalize(x,alpha,nb,ize,b,sum);
387 b[n-1] = en*tempa/x + tempb;
391 sum = sum + sum + b[0];
392 scaled_modified_bessel_i_normalize(x,alpha,nb,ize,b,sum);
403 sum = (sum + b[n-1]*empal)*emp2al/em;
411 for (l = 1; l <= nend; ++l)
415 b[n-1] = en*b[n]/x + b[n+1];
423 sum = (sum + b[n-1]*empal)*emp2al/em;
428 b[0] = K(2.0)*empal*b[1]/x + b[2];
429 sum = sum + sum + b[0];
431 scaled_modified_bessel_i_normalize(x,alpha,nb,ize,b,sum);
449 static inline double innerProduct(
const double phi1,
const double theta1,
450 const double phi2,
const double theta2)
452 double pi2theta1 = K2PI*theta1, pi2theta2 = K2PI*theta2;
453 return (cos(pi2theta1)*cos(pi2theta2)
454 + sin(pi2theta1)*sin(pi2theta2)*cos(K2PI*(phi1-phi2)));
470 return (1.0/(K4PI))*((1.0-h)*(1.0+h))/pow(sqrt(1.0-2.0*h*x+h*h),3.0);
486 return (1.0/(K2PI))/sqrt(1.0-2.0*h*x+h*h);
505 return (x<=h)?(0.0):(pow((x-h),lambda));
522 return exp(2.0*sigma*(x-1.0));
535 int main (
int argc,
char **argv)
584 fftw_complex *prec = NULL;
599 fscanf(stdin,
"testcases=%d\n",&tc_max);
600 fprintf(stdout,
"%d\n",tc_max);
603 for (tc = 0; tc < tc_max; tc++)
606 fscanf(stdin,
"nfsft=%d\n",&use_nfsft);
607 fprintf(stdout,
"%d\n",use_nfsft);
611 fscanf(stdin,
"nfft=%d\n",&use_nfft);
612 fprintf(stdout,
"%d\n",use_nfft);
616 fscanf(stdin,
"cutoff=%d\n",&cutoff);
617 fprintf(stdout,
"%d\n",cutoff);
626 fscanf(stdin,
"fpt=%d\n",&use_fpt);
627 fprintf(stdout,
"%d\n",use_fpt);
629 fscanf(stdin,
"threshold=%lf\n",&threshold);
630 fprintf(stdout,
"%lf\n",threshold);
637 threshold = 1000000000000.0;
653 fscanf(stdin,
"kernel=%d\n",&kt);
654 fprintf(stdout,
"%d\n",kt);
657 fscanf(stdin,
"parameter_sets=%d\n",&ip_max);
658 fprintf(stdout,
"%d\n",ip_max);
661 p = (
double**) nfft_malloc(ip_max*
sizeof(
double*));
666 fscanf(stdin,
"parameters=%d\n",&ipp_max);
667 fprintf(stdout,
"%d\n",ipp_max);
669 for (ip = 0; ip < ip_max; ip++)
672 p[ip] = (
double*) nfft_malloc(ipp_max*
sizeof(
double));
675 for (ipp = 0; ipp < ipp_max; ipp++)
678 fscanf(stdin,
"%lf\n",&p[ip][ipp]);
679 fprintf(stdout,
"%lf\n",p[ip][ipp]);
684 fscanf(stdin,
"bandwidths=%d\n",&im_max);
685 fprintf(stdout,
"%d\n",im_max);
686 m = (
int*) nfft_malloc(im_max*
sizeof(
int));
689 for (im = 0; im < im_max; im++)
692 fscanf(stdin,
"%d\n",&m[im]);
693 fprintf(stdout,
"%d\n",m[im]);
694 m_max = MAX(m_max,m[im]);
698 fscanf(stdin,
"node_sets=%d\n",&ild_max);
699 fprintf(stdout,
"%d\n",ild_max);
700 ld = (
int**) nfft_malloc(ild_max*
sizeof(
int*));
703 for (ild = 0; ild < ild_max; ild++)
706 ld[ild] = (
int*) nfft_malloc(5*
sizeof(
int));
709 fscanf(stdin,
"L=%d ",&ld[ild][0]);
710 fprintf(stdout,
"%d\n",ld[ild][0]);
711 l_max = MAX(l_max,ld[ild][0]);
714 fscanf(stdin,
"D=%d ",&ld[ild][1]);
715 fprintf(stdout,
"%d\n",ld[ild][1]);
716 d_max = MAX(d_max,ld[ild][1]);
719 fscanf(stdin,
"compare=%d ",&ld[ild][2]);
720 fprintf(stdout,
"%d\n",ld[ild][2]);
723 if (ld[ild][2] == YES)
726 fscanf(stdin,
"precomputed=%d\n",&ld[ild][3]);
727 fprintf(stdout,
"%d\n",ld[ild][3]);
731 fscanf(stdin,
"repetitions=%d\n",&ld[ild][4]);
732 fprintf(stdout,
"%d\n",ld[ild][4]);
735 if (ld[ild][3] == YES)
738 ld_max_prec = MAX(ld_max_prec,ld[ild][0]*ld[ild][1]);
740 l_max_prec = MAX(l_max_prec,ld[ild][0]);
753 b = (fftw_complex*) nfft_malloc(l_max*
sizeof(fftw_complex));
754 eta = (
double*) nfft_malloc(2*l_max*
sizeof(
double));
755 f_hat = (fftw_complex*) nfft_malloc(
NFSFT_F_HAT_SIZE(m_max)*
sizeof(fftw_complex));
756 a = (fftw_complex*) nfft_malloc((m_max+1)*
sizeof(fftw_complex));
757 xi = (
double*) nfft_malloc(2*d_max*
sizeof(
double));
758 f_m = (fftw_complex*) nfft_malloc(d_max*
sizeof(fftw_complex));
759 f = (fftw_complex*) nfft_malloc(d_max*
sizeof(fftw_complex));
762 if (precompute == YES)
764 prec = (fftw_complex*) nfft_malloc(ld_max_prec*
sizeof(fftw_complex));
768 for (l = 0; l < l_max; l++)
770 b[l] = (((double)rand())/RAND_MAX) - 0.5;
771 eta[2*l] = (((double)rand())/RAND_MAX) - 0.5;
772 eta[2*l+1] = acos(2.0*(((
double)rand())/RAND_MAX) - 1.0)/(K2PI);
776 for (d = 0; d < d_max; d++)
778 xi[2*d] = (((double)rand())/RAND_MAX) - 0.5;
779 xi[2*d+1] = acos(2.0*(((
double)rand())/RAND_MAX) - 1.0)/(K2PI);
787 for (ip = 0; ip < ip_max; ip++)
794 for (k = 0; k <= m_max; k++)
795 a[k] = SYMBOL_ABEL_POISSON(k,p[ip][0]);
801 for (k = 0; k <= m_max; k++)
802 a[k] = SYMBOL_SINGULARITY(k,p[ip][0]);
810 a[1] = ((p[ip][1]+1+p[ip][0])/(p[ip][1]+2.0))*a[0];
811 for (k = 2; k <= m_max; k++)
812 a[k] = (1.0/(k+p[ip][1]+1))*((2*k-1)*p[ip][0]*a[k-1] -
813 (k-p[ip][1]-2)*a[k-2]);
818 steed = (
double*) nfft_malloc((m_max+1)*
sizeof(double));
819 smbi(2.0*p[ip][0],0.5,m_max+1,2,steed);
820 for (k = 0; k <= m_max; k++)
821 a[k] = K2PI*(sqrt(KPI/p[ip][0]))*steed[k];
828 for (k = 0; k <= m_max; k++)
829 a[k] *= (2*k+1)/(K4PI);
832 for (ild = 0; ild < ild_max; ild++)
835 if (ld[ild][2] != NO)
839 if (ld[ild][3] != NO)
844 rinc = l_max_prec-ld[ild][0];
847 for (d = 0; d < ld[ild][1]; d++)
850 for (l = 0; l < ld[ild][0]; l++)
854 temp =
innerProduct(eta[2*l],eta[2*l+1],xi[2*d],xi[2*d+1]);
894 for (i = 0; i < ld[ild][4]; i++)
900 rinc = l_max_prec-ld[ild][0];
908 constant = ((p[ip][1]+1)/(K2PI*pow(1-p[ip][0],p[ip][1]+1)));
911 for (d = 0; d < ld[ild][1]; d++)
917 for (l = 0; l < ld[ild][0]; l++)
918 f[d] += b[l]*(*ptr++);
930 for (d = 0; d < ld[ild][1]; d++)
936 for (l = 0; l < ld[ild][0]; l++)
937 f[d] += b[l]*(*ptr++);
950 t_dp = t_dp/((double)ld[ild][4]);
965 for (i = 0; i < ld[ild][4]; i++)
973 for (d = 0; d < ld[ild][1]; d++)
979 for (l = 0; l < ld[ild][0]; l++)
983 temp =
innerProduct(eta[2*l],eta[2*l+1],xi[2*d],xi[2*d+1]);
994 for (d = 0; d < ld[ild][1]; d++)
1000 for (l = 0; l < ld[ild][0]; l++)
1004 temp =
innerProduct(eta[2*l],eta[2*l+1],xi[2*d],xi[2*d+1]);
1015 constant = ((p[ip][1]+1)/(K2PI*pow(1-p[ip][0],p[ip][1]+1)));
1018 for (d = 0; d < ld[ild][1]; d++)
1024 for (l = 0; l < ld[ild][0]; l++)
1028 temp =
innerProduct(eta[2*l],eta[2*l+1],xi[2*d],xi[2*d+1]);
1042 for (d = 0; d < ld[ild][1]; d++)
1048 for (l = 0; l < ld[ild][0]; l++)
1052 temp =
innerProduct(eta[2*l],eta[2*l+1],xi[2*d],xi[2*d+1]);
1066 t_d = t_d/((double)ld[ild][4]);
1083 for (im = 0; im < im_max; im++)
1086 nfsft_init_guru(&plan_adjoint, m[im],ld[ild][0],
1091 nfsft_init_guru(&plan,m[im],ld[ild][1],
1096 plan_adjoint.f_hat = f_hat;
1097 plan_adjoint.x = eta;
1102 nfsft_precompute_x(&plan_adjoint);
1103 nfsft_precompute_x(&plan);
1106 if (use_nfsft == BOTH)
1115 for (i = 0; i < ld[ild][4]; i++)
1119 nfsft_adjoint_direct(&plan_adjoint);
1122 for (k = 0; k <= m[im]; k++)
1123 for (n = -k; n <= k; n++)
1127 nfsft_trafo_direct(&plan);
1136 t_fd = t_fd/((double)ld[ild][4]);
1139 if (ld[ild][2] != NO)
1142 err_fd =
X(error_l_infty_1_complex)(f, f_m, ld[ild][1], b,
1148 if (use_nfsft != NO)
1164 for (i = 0; i < ld[ild][4]; i++)
1167 if (use_nfsft != NO)
1175 nfsft_adjoint_direct(&plan_adjoint);
1179 for (k = 0; k <= m[im]; k++)
1180 for (n = -k; n <= k; n++)
1184 if (use_nfsft != NO)
1192 nfsft_trafo_direct(&plan);
1199 if (use_nfsft != NO)
1205 if (use_nfsft != NO)
1208 t_f = t_f/((double)ld[ild][4]);
1213 t_fd = t_fd/((double)ld[ild][4]);
1217 if (ld[ild][2] != NO)
1220 if (use_nfsft != NO)
1223 err_f =
X(error_l_infty_1_complex)(f, f_m, ld[ild][1], b,
1229 err_fd =
X(error_l_infty_1_complex)(f, f_m, ld[ild][1], b,
1235 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,
1251 if (precompute == YES)
1266 for (ild = 0; ild < ild_max; ild++)
1274 for (ip = 0; ip < ip_max; ip++)
1280 return EXIT_SUCCESS;
#define KT_ABEL_POISSON
Abel-Poisson kernel.
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 locallySupportedKernel(const double x, const double h, const double lambda)
Evaluates the locally supported kernel at a node .
static double gaussianKernel(const double x, const double sigma)
Evaluates the spherical Gaussian kernel at a node .
#define KT_GAUSSIAN
Gaussian kernel.
int main(int argc, char **argv)
The main program.
static double singularityKernel(const double x, const double h)
Evaluates the singularity kernel at a node .
static double poissonKernel(const double x, const double h)
Evaluates the Poisson kernel at a node .
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...
#define KT_SINGULARITY
Singularity kernel.
pvalue
Enumeration type for yes/no/both-type parameters.
#define KT_LOC_SUPP
Locally supported kernel.
*We expand this macro for each supported precision * X
R nfft_elapsed_seconds(ticks t1, ticks t0)
Return number of elapsed seconds between two time points.
void nfsft_trafo(nfsft_plan *plan)
void nfsft_adjoint(nfsft_plan *plan)
#define NFSFT_INDEX(k, n, plan)
#define NFSFT_NO_FAST_ALGORITHM
void nfsft_finalize(nfsft_plan *plan)
void nfsft_precompute(int N, double kappa, unsigned int nfsft_flags, unsigned int fpt_flags)
#define NFSFT_F_HAT_SIZE(N)
Internal header file for auxiliary definitions and functions.
Header file for the nfft3 library.