192static 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);
535int 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);
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);
703 for (ild = 0; ild < ild_max; ild++)
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));
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;