00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00026 #include "infft.h"
00027
00028 #include <stdio.h>
00029 #include <sys/time.h>
00030 #include <sys/resource.h>
00031 #ifdef HAVE_MALLOC_H
00032 #include <malloc.h>
00033 #endif
00034 #include "cstripack.h"
00035
00036 #include <complex.h>
00037 #include "nfft3.h"
00038 #include "nfft3util.h"
00039
00043 double nfft_second(void)
00044 {
00045 static struct rusage temp;
00046 double foo, foo1;
00047
00048 getrusage(RUSAGE_SELF,&temp);
00049 foo = temp.ru_utime.tv_sec;
00050 foo1 = temp.ru_utime.tv_usec;
00051 return foo + (foo1/1000000.0);
00052 }
00053
00054 #ifdef HAVE_MALLOC_H
00055 int nfft_total_used_memory(void)
00056 {
00057 struct mallinfo m;
00058 m=mallinfo();
00059 return m.hblkhd + m.uordblks;
00060 }
00061 #else
00062 int nfft_total_used_memory()
00063 {
00064 fprintf(stderr,
00065 "\nWarning in util/total_used_memory: mallinfo() not available\n");
00066 return 0;
00067 }
00068 #endif
00069
00070 int nfft_int_2_pow(int a)
00071 {
00072 return (1U<< a);
00073 }
00074
00075 int nfft_ld(int m)
00076 {
00077 int l=0;
00078 int mm=m;
00079
00080 while(mm>0)
00081 {
00082 mm=(mm>>1);
00083 l++;
00084 }
00085 return (l-1);
00086 }
00087
00090 int nfft_next_power_of_2(int N)
00091 {
00092 int n,i,logn;
00093 int N_is_not_power_of_2=0;
00094
00095 if (N == 0)
00096 {
00097 return 1;
00098 }
00099 else
00100 {
00101 n=N;
00102 logn=0;
00103 while (n != 1)
00104 {
00105 if (n%2 == 1)
00106 {
00107 N_is_not_power_of_2=1;
00108 }
00109 n = n/2;
00110 logn++;
00111 }
00112
00113 if (!N_is_not_power_of_2)
00114 {
00115 logn--;
00116 }
00117
00118 for (i = 0; i <= logn; i++)
00119 {
00120 n = n*2;
00121 }
00122
00123 return n;
00124 }
00125 }
00126
00129 void nfft_next_power_of_2_exp(int N, int *N2, int *t)
00130 {
00131 int n,i,logn;
00132 int N_is_not_power_of_2=0;
00133
00134 if (N == 0)
00135 {
00136 *N2 = 1;
00137 *t = 0;
00138 }
00139 else
00140 {
00141 n=N;
00142 logn=0;
00143 while (n != 1)
00144 {
00145 if (n%2 == 1)
00146 {
00147 N_is_not_power_of_2=1;
00148 }
00149 n = n/2;
00150 logn++;
00151 }
00152
00153 if (!N_is_not_power_of_2)
00154 {
00155 logn--;
00156 }
00157
00158 for (i = 0; i <= logn; i++)
00159 {
00160 n = n*2;
00161 }
00162
00163 *N2 = n;
00164 *t = logn+1;
00165 }
00166 }
00167
00170 int nfft_prod_int(int *vec, int d)
00171 {
00172 int t, prod;
00173
00174 prod=1;
00175 for(t=0; t<d; t++)
00176 prod *= vec[t];
00177
00178 return prod;
00179 }
00180
00183 int nfct_prod_int(int *vec, int d)
00184 {
00185 return nfft_prod_int(vec, d);
00186 }
00187
00190 int nfst_prod_minus_a_int(int *vec, int a, int d)
00191 {
00192 int t, prod;
00193
00194 prod=1;
00195 for(t=0; t<d; t++)
00196 prod *= vec[t]-a;
00197
00198 return prod;
00199 }
00200
00201
00204 int nfft_plain_loop(int *idx,int *N,int d)
00205 {
00206 int t,sum;
00207
00208 sum=idx[0];
00209 for(t=1; t<d; t++)
00210 sum=sum*N[t]+idx[t];
00211
00212 return sum;
00213 }
00214
00217 double nfft_prod_real(double *vec,int d)
00218 {
00219 int t;
00220 double prod;
00221
00222 prod=1.0;
00223 for(t=0; t<d; t++)
00224 prod*=vec[t];
00225
00226 return prod;
00227 }
00228
00229 double nfft_sinc(double x)
00230 {
00231 if (fabs(x) < 1e-20)
00232 return(1.0);
00233 else
00234 return((double)(sin(x)/x));
00235 }
00236
00237 static void bspline_help(int k, double x, double *scratch, int j, int ug, int og,
00238 int r)
00239 {
00240 int i;
00241 int index;
00242 double a;
00244
00245 for(i=og+r-k+1, index=og; index>=ug; i--, index--)
00246 {
00247 a = ((double)(x - i)) / ((double)(k - j));
00248 scratch[index] = (1 - a) * scratch[index-1] + a * scratch[index];
00249 }
00250 }
00251
00255 double nfft_bspline(int k, double x, double *scratch)
00256 {
00257 double result_value;
00258 int r;
00259 int g1,g2;
00260 int j,index,ug,og;
00261 double a;
00263 result_value=0.0;
00264 if(0<x && x<k)
00265 {
00266
00267 if((k-x)<x) x=k-x;
00268
00269 r=(int)(ceil(x)-1.0);
00270
00271 for(index=0; index<k; index++)
00272 scratch[index]=0.0;
00273
00274 scratch[k-r-1]=1.0;
00275
00276
00277 g1 = r;
00278 g2 = k - 1 - r;
00279 ug = g2;
00280
00281
00282
00283 for(j=1, og=g2+1; j<=g1; j++, og++)
00284 {
00285 a = (x - r + k - 1 - og)/(k - j);
00286 scratch[og] = (1 - a) * scratch[og-1];
00287 bspline_help(k,x,scratch,j,ug+1,og-1,r);
00288 a = (x - r + k - 1 - ug)/(k - j);
00289 scratch[ug] = a * scratch[ug];
00290 }
00291 for(og-- ; j<=g2; j++)
00292 {
00293 bspline_help(k,x,scratch,j,ug+1,og,r);
00294 a = (x - r + k - 1 - ug)/(k - j);
00295 scratch[ug] = a * scratch[ug];
00296 }
00297 for( ; j<k; j++)
00298 {
00299 ug++;
00300 bspline_help(k,x,scratch,j,ug,og,r);
00301 }
00302 result_value = scratch[k-1];
00303 }
00304 return(result_value);
00305 }
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371 #define HAVE_LONG_DOUBLE 1
00372
00373
00374 #define RETSIGTYPE void
00375
00376
00377 #define STDC_HEADERS 1
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387 #define SIZEOF_INT 4
00388
00389
00390 #define HAVE_STRING_H 1
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401 #define DOMAIN 1
00402 #define SING 2
00403 #define OVERFLOW 3
00404 #define UNDERFLOW 4
00405 #define TLOSS 5
00406 #define PLOSS 6
00407
00408 #define EDOM 33
00409 #define ERANGE 34
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432 #define UNK 1
00433
00434
00435 #ifdef FLOAT_WORDS_BIGENDIAN
00436 #define BIGENDIAN 1
00437 #else
00438 #define BIGENDIAN 0
00439 #endif
00440
00441
00442
00443
00444
00445
00446
00447 #define VOLATILE
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459 #define XPD 0,
00460
00461
00462 #define DENORMAL 1
00463
00464
00465 #define INFINITIES 1
00466
00467
00468
00469 #define NANS 1
00470
00471
00472 #define MINUSZERO 1
00473
00474
00475
00476 #define ANSIC 1
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538 static double chbevl(double x, double array[], int n)
00539 {
00540 double b0, b1, b2, *p;
00541 int i;
00542
00543 p = array;
00544 b0 = *p++;
00545 b1 = 0.0;
00546 i = n - 1;
00547
00548 do
00549 {
00550 b2 = b1;
00551 b1 = b0;
00552 b0 = x * b1 - b2 + *p++;
00553 }
00554 while( --i );
00555
00556 return( 0.5*(b0-b2) );
00557 }
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606 #ifdef UNK
00607 static double A[] =
00608 {
00609 -4.41534164647933937950E-18,
00610 3.33079451882223809783E-17,
00611 -2.43127984654795469359E-16,
00612 1.71539128555513303061E-15,
00613 -1.16853328779934516808E-14,
00614 7.67618549860493561688E-14,
00615 -4.85644678311192946090E-13,
00616 2.95505266312963983461E-12,
00617 -1.72682629144155570723E-11,
00618 9.67580903537323691224E-11,
00619 -5.18979560163526290666E-10,
00620 2.65982372468238665035E-9,
00621 -1.30002500998624804212E-8,
00622 6.04699502254191894932E-8,
00623 -2.67079385394061173391E-7,
00624 1.11738753912010371815E-6,
00625 -4.41673835845875056359E-6,
00626 1.64484480707288970893E-5,
00627 -5.75419501008210370398E-5,
00628 1.88502885095841655729E-4,
00629 -5.76375574538582365885E-4,
00630 1.63947561694133579842E-3,
00631 -4.32430999505057594430E-3,
00632 1.05464603945949983183E-2,
00633 -2.37374148058994688156E-2,
00634 4.93052842396707084878E-2,
00635 -9.49010970480476444210E-2,
00636 1.71620901522208775349E-1,
00637 -3.04682672343198398683E-1,
00638 6.76795274409476084995E-1
00639 };
00640 #endif
00641
00642 #ifdef DEC
00643 static unsigned short A[] = {
00644 0121642,0162671,0004646,0103567,
00645 0022431,0115424,0135755,0026104,
00646 0123214,0023533,0110365,0156635,
00647 0023767,0033304,0117662,0172716,
00648 0124522,0100426,0012277,0157531,
00649 0025254,0155062,0054461,0030465,
00650 0126010,0131143,0013560,0153604,
00651 0026517,0170577,0006336,0114437,
00652 0127227,0162253,0152243,0052734,
00653 0027724,0142766,0061641,0160200,
00654 0130416,0123760,0116564,0125262,
00655 0031066,0144035,0021246,0054641,
00656 0131537,0053664,0060131,0102530,
00657 0032201,0155664,0165153,0020652,
00658 0132617,0061434,0074423,0176145,
00659 0033225,0174444,0136147,0122542,
00660 0133624,0031576,0056453,0020470,
00661 0034211,0175305,0172321,0041314,
00662 0134561,0054462,0147040,0165315,
00663 0035105,0124333,0120203,0162532,
00664 0135427,0013750,0174257,0055221,
00665 0035726,0161654,0050220,0100162,
00666 0136215,0131361,0000325,0041110,
00667 0036454,0145417,0117357,0017352,
00668 0136702,0072367,0104415,0133574,
00669 0037111,0172126,0072505,0014544,
00670 0137302,0055601,0120550,0033523,
00671 0037457,0136543,0136544,0043002,
00672 0137633,0177536,0001276,0066150,
00673 0040055,0041164,0100655,0010521
00674 };
00675 #endif
00676
00677 #ifdef IBMPC
00678 static unsigned short A[] = {
00679 0xd0ef,0x2134,0x5cb7,0xbc54,
00680 0xa589,0x977d,0x3362,0x3c83,
00681 0xbbb4,0x721e,0x84eb,0xbcb1,
00682 0x5eba,0x93f6,0xe6d8,0x3cde,
00683 0xfbeb,0xc297,0x5022,0xbd0a,
00684 0x2627,0x4b26,0x9b46,0x3d35,
00685 0x1af0,0x62ee,0x164c,0xbd61,
00686 0xd324,0xe19b,0xfe2f,0x3d89,
00687 0x6abc,0x7a94,0xfc95,0xbdb2,
00688 0x3c10,0xcc74,0x98be,0x3dda,
00689 0x9556,0x13ae,0xd4fe,0xbe01,
00690 0xcb34,0xa454,0xd903,0x3e26,
00691 0x30ab,0x8c0b,0xeaf6,0xbe4b,
00692 0x6435,0x9d4d,0x3b76,0x3e70,
00693 0x7f8d,0x8f22,0xec63,0xbe91,
00694 0xf4ac,0x978c,0xbf24,0x3eb2,
00695 0x6427,0xcba5,0x866f,0xbed2,
00696 0x2859,0xbe9a,0x3f58,0x3ef1,
00697 0x1d5a,0x59c4,0x2b26,0xbf0e,
00698 0x7cab,0x7410,0xb51b,0x3f28,
00699 0xeb52,0x1f15,0xe2fd,0xbf42,
00700 0x100e,0x8a12,0xdc75,0x3f5a,
00701 0xa849,0x201a,0xb65e,0xbf71,
00702 0xe3dd,0xf3dd,0x9961,0x3f85,
00703 0xb6f0,0xf121,0x4e9e,0xbf98,
00704 0xa32d,0xcea8,0x3e8a,0x3fa9,
00705 0x06ea,0x342d,0x4b70,0xbfb8,
00706 0x88c0,0x77ac,0xf7ac,0x3fc5,
00707 0xcd8d,0xc057,0x7feb,0xbfd3,
00708 0xa22a,0x9035,0xa84e,0x3fe5,
00709 };
00710 #endif
00711
00712 #ifdef MIEEE
00713 static unsigned short A[] = {
00714 0xbc54,0x5cb7,0x2134,0xd0ef,
00715 0x3c83,0x3362,0x977d,0xa589,
00716 0xbcb1,0x84eb,0x721e,0xbbb4,
00717 0x3cde,0xe6d8,0x93f6,0x5eba,
00718 0xbd0a,0x5022,0xc297,0xfbeb,
00719 0x3d35,0x9b46,0x4b26,0x2627,
00720 0xbd61,0x164c,0x62ee,0x1af0,
00721 0x3d89,0xfe2f,0xe19b,0xd324,
00722 0xbdb2,0xfc95,0x7a94,0x6abc,
00723 0x3dda,0x98be,0xcc74,0x3c10,
00724 0xbe01,0xd4fe,0x13ae,0x9556,
00725 0x3e26,0xd903,0xa454,0xcb34,
00726 0xbe4b,0xeaf6,0x8c0b,0x30ab,
00727 0x3e70,0x3b76,0x9d4d,0x6435,
00728 0xbe91,0xec63,0x8f22,0x7f8d,
00729 0x3eb2,0xbf24,0x978c,0xf4ac,
00730 0xbed2,0x866f,0xcba5,0x6427,
00731 0x3ef1,0x3f58,0xbe9a,0x2859,
00732 0xbf0e,0x2b26,0x59c4,0x1d5a,
00733 0x3f28,0xb51b,0x7410,0x7cab,
00734 0xbf42,0xe2fd,0x1f15,0xeb52,
00735 0x3f5a,0xdc75,0x8a12,0x100e,
00736 0xbf71,0xb65e,0x201a,0xa849,
00737 0x3f85,0x9961,0xf3dd,0xe3dd,
00738 0xbf98,0x4e9e,0xf121,0xb6f0,
00739 0x3fa9,0x3e8a,0xcea8,0xa32d,
00740 0xbfb8,0x4b70,0x342d,0x06ea,
00741 0x3fc5,0xf7ac,0x77ac,0x88c0,
00742 0xbfd3,0x7feb,0xc057,0xcd8d,
00743 0x3fe5,0xa84e,0x9035,0xa22a
00744 };
00745 #endif
00746
00747
00748
00749
00750
00751
00752
00753
00754 #ifdef UNK
00755 static double B[] =
00756 {
00757 -7.23318048787475395456E-18,
00758 -4.83050448594418207126E-18,
00759 4.46562142029675999901E-17,
00760 3.46122286769746109310E-17,
00761 -2.82762398051658348494E-16,
00762 -3.42548561967721913462E-16,
00763 1.77256013305652638360E-15,
00764 3.81168066935262242075E-15,
00765 -9.55484669882830764870E-15,
00766 -4.15056934728722208663E-14,
00767 1.54008621752140982691E-14,
00768 3.85277838274214270114E-13,
00769 7.18012445138366623367E-13,
00770 -1.79417853150680611778E-12,
00771 -1.32158118404477131188E-11,
00772 -3.14991652796324136454E-11,
00773 1.18891471078464383424E-11,
00774 4.94060238822496958910E-10,
00775 3.39623202570838634515E-9,
00776 2.26666899049817806459E-8,
00777 2.04891858946906374183E-7,
00778 2.89137052083475648297E-6,
00779 6.88975834691682398426E-5,
00780 3.36911647825569408990E-3,
00781 8.04490411014108831608E-1
00782 };
00783 #endif
00784
00785 #ifdef DEC
00786 static unsigned short B[] = {
00787 0122005,0066672,0123124,0054311,
00788 0121662,0033323,0030214,0104602,
00789 0022515,0170300,0113314,0020413,
00790 0022437,0117350,0035402,0007146,
00791 0123243,0000135,0057220,0177435,
00792 0123305,0073476,0144106,0170702,
00793 0023777,0071755,0017527,0154373,
00794 0024211,0052214,0102247,0033270,
00795 0124454,0017763,0171453,0012322,
00796 0125072,0166316,0075505,0154616,
00797 0024612,0133770,0065376,0025045,
00798 0025730,0162143,0056036,0001632,
00799 0026112,0015077,0150464,0063542,
00800 0126374,0101030,0014274,0065457,
00801 0127150,0077271,0125763,0157617,
00802 0127412,0104350,0040713,0120445,
00803 0027121,0023765,0057500,0001165,
00804 0030407,0147146,0003643,0075644,
00805 0031151,0061445,0044422,0156065,
00806 0031702,0132224,0003266,0125551,
00807 0032534,0000076,0147153,0005555,
00808 0033502,0004536,0004016,0026055,
00809 0034620,0076433,0142314,0171215,
00810 0036134,0146145,0013454,0101104,
00811 0040115,0171425,0062500,0047133
00812 };
00813 #endif
00814
00815 #ifdef IBMPC
00816 static unsigned short B[] = {
00817 0x8b19,0x54ca,0xadb7,0xbc60,
00818 0x9130,0x6611,0x46da,0xbc56,
00819 0x8421,0x12d9,0xbe18,0x3c89,
00820 0x41cd,0x0760,0xf3dd,0x3c83,
00821 0x1fe4,0xabd2,0x600b,0xbcb4,
00822 0xde38,0xd908,0xaee7,0xbcb8,
00823 0xfb1f,0xa3ea,0xee7d,0x3cdf,
00824 0xe6d7,0x9094,0x2a91,0x3cf1,
00825 0x629a,0x7e65,0x83fe,0xbd05,
00826 0xbb32,0xcf68,0x5d99,0xbd27,
00827 0xc545,0x0d5f,0x56ff,0x3d11,
00828 0xc073,0x6b83,0x1c8c,0x3d5b,
00829 0x8cec,0xfa26,0x4347,0x3d69,
00830 0x8d66,0x0317,0x9043,0xbd7f,
00831 0x7bf2,0x357e,0x0fd7,0xbdad,
00832 0x7425,0x0839,0x511d,0xbdc1,
00833 0x004f,0xabe8,0x24fe,0x3daa,
00834 0x6f75,0xc0f4,0xf9cc,0x3e00,
00835 0x5b87,0xa922,0x2c64,0x3e2d,
00836 0xd56d,0x80d6,0x5692,0x3e58,
00837 0x616e,0xd9cd,0x8007,0x3e8b,
00838 0xc586,0xc101,0x412b,0x3ec8,
00839 0x9e52,0x7899,0x0fa3,0x3f12,
00840 0x9049,0xa2e5,0x998c,0x3f6b,
00841 0x09cb,0xaca8,0xbe62,0x3fe9
00842 };
00843 #endif
00844
00845 #ifdef MIEEE
00846 static unsigned short B[] = {
00847 0xbc60,0xadb7,0x54ca,0x8b19,
00848 0xbc56,0x46da,0x6611,0x9130,
00849 0x3c89,0xbe18,0x12d9,0x8421,
00850 0x3c83,0xf3dd,0x0760,0x41cd,
00851 0xbcb4,0x600b,0xabd2,0x1fe4,
00852 0xbcb8,0xaee7,0xd908,0xde38,
00853 0x3cdf,0xee7d,0xa3ea,0xfb1f,
00854 0x3cf1,0x2a91,0x9094,0xe6d7,
00855 0xbd05,0x83fe,0x7e65,0x629a,
00856 0xbd27,0x5d99,0xcf68,0xbb32,
00857 0x3d11,0x56ff,0x0d5f,0xc545,
00858 0x3d5b,0x1c8c,0x6b83,0xc073,
00859 0x3d69,0x4347,0xfa26,0x8cec,
00860 0xbd7f,0x9043,0x0317,0x8d66,
00861 0xbdad,0x0fd7,0x357e,0x7bf2,
00862 0xbdc1,0x511d,0x0839,0x7425,
00863 0x3daa,0x24fe,0xabe8,0x004f,
00864 0x3e00,0xf9cc,0xc0f4,0x6f75,
00865 0x3e2d,0x2c64,0xa922,0x5b87,
00866 0x3e58,0x5692,0x80d6,0xd56d,
00867 0x3e8b,0x8007,0xd9cd,0x616e,
00868 0x3ec8,0x412b,0xc101,0xc586,
00869 0x3f12,0x0fa3,0x7899,0x9e52,
00870 0x3f6b,0x998c,0xa2e5,0x9049,
00871 0x3fe9,0xbe62,0xaca8,0x09cb
00872 };
00873 #endif
00874
00879 double nfft_i0(double x)
00880 {
00881 double y;
00882
00883 if( x < 0 )
00884 x = -x;
00885 if( x <= 8.0 )
00886 {
00887 y = (x/2.0) - 2.0;
00888 return( exp(x) * chbevl( y, A, 30 ) );
00889 }
00890
00891 return( exp(x) * chbevl( 32.0/x - 2.0, B, 25 ) / sqrt(x) );
00892 }
00893
00896 double nfft_dot_complex(double _Complex *x, int n)
00897 {
00898 int k;
00899 double dot;
00900
00901 for(k=0,dot=0; k<n; k++)
00902 dot+=conj(x[k])*x[k];
00903
00904 return dot;
00905 }
00906
00909 double nfft_dot_double(double *x, int n)
00910 {
00911 int k;
00912 double dot;
00913
00914 for(k=0,dot=0; k<n; k++)
00915 dot+=x[k]*x[k];
00916
00917 return dot;
00918 }
00919
00920
00923 double nfft_dot_w_complex(double _Complex *x, double *w, int n)
00924 {
00925 int k;
00926 double dot;
00927
00928 for(k=0,dot=0.0; k<n; k++)
00929 dot+=w[k]*conj(x[k])*x[k];
00930
00931 return dot;
00932 }
00933
00936 double nfft_dot_w_double(double *x, double *w, int n)
00937 {
00938 int k;
00939 double dot;
00940
00941 for(k=0,dot=0.0; k<n; k++)
00942 dot+=w[k]*x[k]*x[k];
00943
00944 return dot;
00945 }
00946
00947
00951 double nfft_dot_w_w2_complex(double _Complex *x, double *w, double *w2, int n)
00952 {
00953 int k;
00954 double dot;
00955
00956 for(k=0,dot=0.0; k<n; k++)
00957 dot+=w[k]*w2[k]*w2[k]*conj(x[k])*x[k];
00958
00959 return dot;
00960 }
00961
00965 double nfft_dot_w2_complex(double _Complex *x, double *w2, int n)
00966 {
00967 int k;
00968 double dot;
00969
00970 for(k=0,dot=0.0; k<n; k++)
00971 dot+=w2[k]*w2[k]*conj(x[k])*x[k];
00972
00973 return dot;
00974 }
00975
00978 void nfft_cp_complex(double _Complex *x, double _Complex *y, int n)
00979 {
00980 int k;
00981
00982 for(k=0;k<n;k++)
00983 x[k]=y[k];
00984 }
00985
00988 void nfft_cp_double(double *x, double *y, int n)
00989 {
00990 int k;
00991
00992 for(k=0;k<n;k++)
00993 x[k]=y[k];
00994 }
00995
00998 void nfft_cp_a_complex(double _Complex *x, double a, double _Complex *y, int n)
00999 {
01000 int k;
01001
01002 for(k=0;k<n;k++)
01003 x[k]=a*y[k];
01004 }
01005
01008 static void nfft_cp_a_double(double *x, double a, double *y, int n)
01009 {
01010 int k;
01011
01012 for(k=0;k<n;k++)
01013 x[k]=a*y[k];
01014 }
01015
01016
01019 void nfft_cp_w_complex(double _Complex *x, double *w, double _Complex *y, int n)
01020 {
01021 int k;
01022
01023 for(k=0;k<n;k++)
01024 x[k]=w[k]*y[k];
01025 }
01026
01029 void nfft_cp_w_double(double *x, double *w, double *y, int n)
01030 {
01031 int k;
01032
01033 for(k=0;k<n;k++)
01034 x[k]=w[k]*y[k];
01035 }
01036
01037
01038
01041 void nfft_upd_axpy_complex(double _Complex *x, double a, double _Complex *y, int n)
01042 {
01043 int k;
01044
01045 for(k=0;k<n;k++)
01046 x[k]=a*x[k]+y[k];
01047 }
01048
01051 void nfft_upd_axpy_double(double *x, double a, double *y, int n)
01052 {
01053 int k;
01054
01055 for(k=0;k<n;k++)
01056 x[k]=a*x[k]+y[k];
01057 }
01058
01059
01062 void nfft_upd_xpay_complex(double _Complex *x, double a, double _Complex *y, int n)
01063 {
01064 int k;
01065
01066 for(k=0;k<n;k++)
01067 x[k]+=a*y[k];
01068 }
01069
01072 void nfft_upd_xpay_double(double *x, double a, double *y, int n)
01073 {
01074 int k;
01075
01076 for(k=0;k<n;k++)
01077 x[k]+=a*y[k];
01078 }
01079
01080
01081
01084 void nfft_upd_axpby_complex(double _Complex *x, double a, double _Complex *y, double b, int n)
01085 {
01086 int k;
01087
01088 for(k=0;k<n;k++)
01089 x[k]=a*x[k]+b*y[k];
01090 }
01091
01094 void nfft_upd_axpby_double(double *x, double a, double *y, double b, int n)
01095 {
01096 int k;
01097
01098 for(k=0;k<n;k++)
01099 x[k]=a*x[k]+b*y[k];
01100 }
01101
01102
01105 void nfft_upd_xpawy_complex(double _Complex *x, double a, double *w, double _Complex *y, int n)
01106 {
01107 int k;
01108
01109 for(k=0;k<n;k++)
01110 x[k]+=a*w[k]*y[k];
01111 }
01112
01115 void nfft_upd_xpawy_double(double *x, double a, double *w, double *y, int n)
01116 {
01117 int k;
01118
01119 for(k=0;k<n;k++)
01120 x[k]+=a*w[k]*y[k];
01121 }
01122
01123
01124
01127 void nfft_upd_axpwy_complex(double _Complex *x, double a, double *w, double _Complex *y, int n)
01128 {
01129 int k;
01130
01131 for(k=0;k<n;k++)
01132 x[k]=a*x[k]+w[k]*y[k];
01133 }
01134
01137 void nfft_upd_axpwy_double(double *x, double a, double *w, double *y, int n)
01138 {
01139 int k;
01140
01141 for(k=0;k<n;k++)
01142 x[k]=a*x[k]+w[k]*y[k];
01143 }
01144
01145
01146 void nfft_fftshift_complex(double _Complex *x, int d, int* N)
01147 {
01148 int d_pre, d_act, d_post;
01149 int N_pre, N_act, N_post;
01150 int k_pre, k_act, k_post;
01151 int k,k_swap;
01152
01153 double _Complex x_swap;
01154
01155 for(d_act=0;d_act<d;d_act++)
01156 {
01157 for(d_pre=0, N_pre=1;d_pre<d_act;d_pre++)
01158 N_pre*=N[d_pre];
01159
01160 N_act=N[d_act];
01161
01162 for(d_post=d_act+1, N_post=1;d_post<d;d_post++)
01163 N_post*=N[d_post];
01164
01165 for(k_pre=0;k_pre<N_pre;k_pre++)
01166 for(k_act=0;k_act<N_act/2;k_act++)
01167 for(k_post=0;k_post<N_post;k_post++)
01168 {
01169 k=(k_pre*N_act+k_act)*N_post+k_post;
01170 k_swap=(k_pre*N_act+k_act+N_act/2)*N_post+k_post;
01171
01172 x_swap=x[k];
01173 x[k]=x[k_swap];
01174 x[k_swap]=x_swap;
01175 }
01176 }
01177 }
01178
01179 static double l_1_complex(double _Complex *x, double _Complex *y, int n)
01180 {
01181 int k;
01182 double l1;
01183
01184 if(y==NULL)
01185 for(k=0,l1=0; k<n; k++)
01186 l1+=cabs(x[k]);
01187 else
01188 for(k=0,l1=0; k<n; k++)
01189 l1+=cabs(x[k]-y[k]);
01190
01191 return l1;
01192 }
01193
01194 static double l_1_double(double *x, double *y, int n)
01195 {
01196 int k;
01197 double l1;
01198
01199 if(y==NULL)
01200 for(k=0,l1=0; k<n; k++)
01201 l1+=fabs(x[k]);
01202 else
01203 for(k=0,l1=0; k<n; k++)
01204 l1+=fabs(x[k]-y[k]);
01205
01206 return l1;
01207 }
01208
01209
01210
01211
01212 static double l_2_complex(double _Complex *x, double _Complex *y, int n)
01213 {
01214 int k;
01215 double l22;
01216
01217 if(y==NULL)
01218 for(k=0,l22=0; k<n; k++)
01219 l22+=conj(x[k])*x[k];
01220 else
01221 for(k=0,l22=0; k<n; k++)
01222 l22+=conj(x[k]-y[k])*(x[k]-y[k]);
01223
01224 return sqrt(l22);
01225 }
01226
01227 static double l_2_double(double *x, double *y, int n)
01228 {
01229 int k;
01230 double l22;
01231
01232 if(y==NULL)
01233 for(k=0,l22=0; k<n; k++)
01234 l22+=x[k]*x[k];
01235 else
01236 for(k=0,l22=0; k<n; k++)
01237 l22+=(x[k]-y[k])*(x[k]-y[k]);
01238
01239 return sqrt(l22);
01240 }
01241
01242 static double l_infty_complex(double _Complex *x, double _Complex *y, int n)
01243 {
01244 int k;
01245 double linfty;
01246
01247 if(y==NULL)
01248 for(k=0,linfty=0; k<n; k++)
01249 linfty=((linfty<cabs(x[k]))?cabs(x[k]):linfty);
01250 else
01251 for(k=0,linfty=0; k<n; k++)
01252 linfty=((linfty<cabs(x[k]-y[k]))?cabs(x[k]-y[k]):linfty);
01253
01254 return linfty;
01255 }
01256
01257
01258 static double l_infty_double(double *x, double *y, int n)
01259 {
01260 int k;
01261 double linfty;
01262
01263 if(y==NULL)
01264 for(k=0,linfty=0; k<n; k++)
01265 linfty=((linfty<fabs(x[k]))?fabs(x[k]):linfty);
01266 else
01267 for(k=0,linfty=0; k<n; k++)
01268 linfty=((linfty<fabs(x[k]-y[k]))?fabs(x[k]-y[k]):linfty);
01269
01270 return linfty;
01271 }
01272
01273
01274
01275
01276
01279 double nfft_error_l_infty_complex(double _Complex *x, double _Complex *y, int n)
01280 {
01281 return (l_infty_complex(x, y, n)/l_infty_complex(x, NULL, n));
01282 }
01283
01286 double nfft_error_l_infty_double(double *x, double *y, int n)
01287 {
01288 return (l_infty_double(x, y, n)/l_infty_double(x, NULL, n));
01289 }
01290
01291
01292
01295 double nfft_error_l_infty_1_complex(double _Complex *x, double _Complex *y, int n,
01296 double _Complex *z, int m)
01297 {
01298 return (l_infty_complex(x, y, n)/l_1_complex(z, NULL, m));
01299 }
01300
01303 double nfft_error_l_infty_1_double(double *x, double *y, int n,
01304 double *z, int m)
01305 {
01306 return (l_infty_double(x, y, n)/l_1_double(z, NULL, m));
01307 }
01308
01309
01310
01313 double nfft_error_l_2_complex(double _Complex *x, double _Complex *y, int n)
01314 {
01315 return (l_2_complex(x, y, n)/l_2_complex(x, NULL, n));
01316 }
01317
01320 double nfft_error_l_2_double(double *x, double *y, int n)
01321 {
01322 return (l_2_double(x, y, n)/l_2_double(x, NULL, n));
01323 }
01324
01325
01326
01329 void nfft_vpr_int(int *x, int n, char *text)
01330 {
01331 int k;
01332
01333 if(text!=NULL)
01334 {
01335 printf ("\n %s, adr=%p\n", text, (void*)x);
01336 for (k=0; k<n; k++)
01337 {
01338 if (k%8==0)
01339 printf("%6d.\t", k);
01340 printf("%d,", x[k]);
01341 if (k%8==7)
01342 printf("\n");
01343 }
01344 if (n%8!=0)
01345 printf("\n");
01346 }
01347 else
01348 for (k=0; k<n; k++)
01349 printf("%d,\n", x[k]);
01350 fflush(stdout);
01351 }
01352
01354 void X(vpr_double)(R *x, const int n, const char *text)
01355 {
01356 int k;
01357
01358 if (x == NULL)
01359 {
01360 printf("null pointer\n");
01361 fflush(stdout);
01362 exit(-1);
01363 }
01364
01365 if (text != NULL)
01366 {
01367 printf ("\n %s, adr=%p\n", text, (void*)x);
01368
01369 for (k = 0; k < n; k++)
01370 {
01371 if (k%8 == 0)
01372 printf("%6d.\t", k);
01373
01374 printf("%+.1" FE ",", x[k]);
01375
01376 if (k%8 == 7)
01377 printf("\n");
01378 }
01379
01380 if (n%8 != 0)
01381 printf("\n");
01382 }
01383 else
01384 for (k = 0; k < n; k++)
01385 printf("%+" FE ",\n", x[k]);
01386
01387 fflush(stdout);
01388 }
01389
01391 void X(vpr_complex)(C *x, const int n, const char *text)
01392 {
01393 int k;
01394
01395 if(text != NULL)
01396 {
01397 printf("\n %s, adr=%p\n", text, (void*)x);
01398 for (k = 0; k < n; k++)
01399 {
01400 if (k%4 == 0)
01401 printf("%6d.\t", k);
01402
01403 printf("%+.1" FE "%+.1" FE "i,", CREAL(x[k]), CIMAG(x[k]));
01404
01405 if (k%4==3)
01406 printf("\n");
01407 }
01408 if (n%4!=0)
01409 printf("\n");
01410 }
01411 else
01412 for (k = 0; k < n; k++)
01413 printf("%+" FE "%+" FE "i,\n", CREAL(x[k]), CIMAG(x[k]));
01414
01415 fflush(stdout);
01416 }
01417
01418 void X(vrand_unit_complex)(C *x, const int n)
01419 {
01420 int k;
01421
01422 for (k = 0; k < n; k++)
01423 x[k] = RAND + II*RAND;
01424 }
01425
01426 void X(vrand_shifted_unit_double)(R *x, const int n)
01427 {
01428 int k;
01429
01430 for (k = 0; k < n; k++)
01431 x[k] = RAND - K(0.5);
01432 }
01433
01435 void X(voronoi_weights_1d)(R *w, R *x, const int M)
01436 {
01437 int j;
01438
01439 w[0] = (x[1]-x[0])/K(2.0);
01440
01441 for(j = 1; j < M-1; j++)
01442 w[j] = (x[j+1]-x[j-1])/K(2.0);
01443
01444 w[M-1] = (x[M-1]-x[M-2])/K(2.0);
01445 }
01446
01447 void nfft_voronoi_weights_S2(double *w, double *xi, int M)
01448 {
01449 double *x;
01450 double *y;
01451 double *z;
01452 int j;
01453 int k;
01454 int el;
01455 int Mlocal = M;
01456 int lnew;
01457 int ier;
01458 int *list;
01459 int *lptr;
01460 int *lend;
01461 int *near;
01462 int *next;
01463 double *dist;
01464 int *ltri;
01465 int *listc;
01466 int nb;
01467 double *xc;
01468 double *yc;
01469 double *zc;
01470 double *rc;
01471 double *vr;
01472 int lp;
01473 int lpl;
01474 int kv;
01475 double a;
01476
01477
01478 x = (double*)nfft_malloc(M * sizeof(double));
01479 y = (double*)nfft_malloc(M * sizeof(double));
01480 z = (double*)nfft_malloc(M * sizeof(double));
01481
01482 list = (int*)nfft_malloc((6*M-12+1)*sizeof(int));
01483 lptr = (int*)nfft_malloc((6*M-12+1)*sizeof(int));
01484 lend = (int*)nfft_malloc((M+1)*sizeof(int));
01485 near = (int*)nfft_malloc((M+1)*sizeof(int));
01486 next = (int*)nfft_malloc((M+1)*sizeof(int));
01487 dist = (double*)nfft_malloc((M+1)*sizeof(double));
01488 ltri = (int*)nfft_malloc((6*M+1)*sizeof(int));
01489 listc = (int*)nfft_malloc((6*M-12+1)*sizeof(int));
01490 xc = (double*)nfft_malloc((2*M-4+1)*sizeof(double));
01491 yc = (double*)nfft_malloc((2*M-4+1)*sizeof(double));
01492 zc = (double*)nfft_malloc((2*M-4+1)*sizeof(double));
01493 rc = (double*)nfft_malloc((2*M-4+1)*sizeof(double));
01494 vr = (double*)nfft_malloc(3*(2*M-4+1)*sizeof(double));
01495
01496
01497
01498 for (k = 0; k < M; k++)
01499 {
01500 x[k] = sin(2.0*PI*xi[2*k+1])*cos(2.0*PI*xi[2*k]);
01501 y[k] = sin(2.0*PI*xi[2*k+1])*sin(2.0*PI*xi[2*k]);
01502 z[k] = cos(2.0*PI*xi[2*k+1]);
01503 }
01504
01505
01506 trmesh_(&Mlocal, x, y, z, list, lptr, lend, &lnew, near, next, dist, &ier);
01507
01508
01509 if (ier == 0)
01510 {
01511
01512 crlist_(&Mlocal, &Mlocal, x, y, z, list, lend, lptr, &lnew, ltri, listc, &nb, xc,
01513 yc, zc, rc, &ier);
01514
01515 if (ier == 0)
01516 {
01517
01518 for (k = 0; k < M; k++)
01519 {
01520
01521 lpl = lend[k];
01522 lp = lpl;
01523
01524 j = 0;
01525 vr[3*j] = x[k];
01526 vr[3*j+1] = y[k];
01527 vr[3*j+2] = z[k];
01528
01529 do
01530 {
01531 j++;
01532
01533 lp = lptr[lp-1];
01534 kv = listc[lp-1];
01535 vr[3*j] = xc[kv-1];
01536 vr[3*j+1] = yc[kv-1];
01537 vr[3*j+2] = zc[kv-1];
01538
01539 } while (lp != lpl);
01540
01541 a = 0;
01542 for (el = 0; el < j; el++)
01543 {
01544 a += areas_(vr, &vr[3*(el+1)],&vr[3*(((el+1)%j)+1)]);
01545 }
01546
01547 w[k] = a;
01548 }
01549 }
01550 }
01551
01552
01553 nfft_free(x);
01554 nfft_free(y);
01555 nfft_free(z);
01556
01557 nfft_free(list);
01558 nfft_free(lptr);
01559 nfft_free(lend);
01560 nfft_free(near);
01561 nfft_free(next);
01562 nfft_free(dist);
01563 nfft_free(ltri);
01564 nfft_free(listc);
01565 nfft_free(xc);
01566 nfft_free(yc);
01567 nfft_free(zc);
01568 nfft_free(rc);
01569 nfft_free(vr);
01570 }
01571
01576 R X(modified_fejer)(const int N, const int kk)
01577 {
01578 return (K(2.0)/((R)(N*N))*(K(1.0)-FABS(K(2.0)*kk+K(1.0))/((R)N)));
01579 }
01580
01582 R X(modified_jackson2)(const int N, const int kk)
01583 {
01584 int kj;
01585 const R n=(N/K(2.0)+K(1.0))/K(2.0);
01586 R result, k;
01587
01588 for (result = K(0.0), kj = kk; kj <= kk+1; kj++)
01589 {
01590 k = ABS(kj);
01591
01592 if(k/n < K(1.0))
01593 result += K(1.0) - (K(3.0)*k + K(6.0)*n*POW(k,K(2.0))
01594 - K(3.0)*POW(k,K(3.0)))/(K(2.0)*n*(K(2.0)*POW(n,K(2.0))+K(1.0)));
01595 else
01596 result+= (K(2.0)*n-k)*(POW(2*n-k,K(2.0))-K(1.0))/(K(2.0)
01597 *n*(K(2.0)*POW(n,K(2.0))+K(1.0)));
01598 }
01599
01600 return result;
01601 }
01602
01604 R X(modified_jackson4)(const int N, const int kk)
01605 {
01606 int kj;
01607 const R n = (N/K(2.0)+K(3.0))/K(4.0), normalisation = (K(2416.0)*POW(n,K(7.0))
01608 + K(1120.0)*POW(n,K(5.0)) + K(784.0)*POW(n,K(3.0)) + K(720.0)*n);
01609 R result, k;
01610
01611 for (result = K(0.0), kj = kk; kj <= kk + 1; kj++)
01612 {
01613 k = ABS(kj);
01614
01615 if (k/n < K(1.0))
01616 result += K(1.0) - (K(1260.0)*k + (K(1680.0)*POW(n, K(5.0))
01617 + K(2240.0)*POW(n, K(3.0)) + K(2940.0)*n)*POW(k, K(2.0))
01618 - K(1715.0)*POW(k, K(3.0)) - (K(560.0)*POW(n, K(3.0))
01619 + K(1400.0)*n)*POW(k, K(4.0)) + K(490.0)*POW(k, K(5.0))
01620 + K(140.0)*n*POW(k, K(6.0)) - K(35.0)*POW(k,K(7.0)))/normalisation;
01621
01622 if ((K(1.0) <= k/n) && (k/n < K(2.0)))
01623 result += ((K(2472.0)*POW(n, K(7.0)) + K(336.0)*POW(n, K(5.0))
01624 + K(3528.0)*POW(n, K(3.0)) - K(1296.0)*n) - (K(392.0)*POW(n, K(6.0))
01625 - K(3920.0)*POW(n, K(4.0)) + K(8232.0)*POW(n, K(2.0)) - K(756.0))*k
01626 - (K(504.0)*POW(n, K(5.0)) + K(10080.0)*POW(n, K(3.0))
01627 - K(5292.0)*n)*POW(k, K(2.0)) - (K(1960.0)*POW(n, K(4.0))
01628 - K(7840.0)*POW(n, K(2.0)) + K(1029.0))*POW(k, K(3.0))
01629 + (K(2520.0)*POW(n, K(3.0)) - K(2520.0)*n) * POW(k, K(4.0))
01630 - (K(1176.0)*POW(n, K(2.0)) - K(294.0)) * POW(k, K(5.0))
01631 + K(252.0)*n*POW(k, K(6.0)) - K(21.0)*POW(k, K(7.0)))/normalisation;
01632
01633 if ((K(2.0) <= k/n) && (k/n < K(3.0)))
01634 result += (-(K(1112.0)*POW(n, K(7.0)) - K(12880.0)*POW(n, K(5.0))
01635 + K(7448.0)*POW(n, K(3.0)) - K(720.0)*n) + (K(12152.0)*POW(n, K(6.0))
01636 - K(27440.0)*POW(n, K(4.0)) + K(8232.0)*POW(n, K(2.0)) - K(252.0))*k
01637 - (K(19320.0)*POW(n, K(5.0)) - K(21280.0)*POW(n, K(3.0))
01638 + K(2940.0)*n)*POW(k, K(2.0)) + (K(13720.0)*POW(n, K(4.0))
01639 - K(7840.0)*POW(n, K(2.0)) + K(343.0))*POW(k, K(3.0))
01640 - (K(5320.0)*POW(n, K(3.0)) - K(1400.0)*n)*POW(k, K(4.0))
01641 + (K(1176.0)*POW(n, K(2.0)) - K(98.0))*POW(k, K(5.0))
01642 - K(140.0)*n*POW(k, K(6.0)) + K(7.0) * POW(k, K(7.0)))/normalisation;
01643
01644 if ((K(3.0) <= k/n) && (k/n < K(4.0)))
01645 result += ((4*n-k)*(POW(4*n-k, K(2.0)) - K(1.0))*(POW(4*n-k, K(2.0))
01646 - K(4.0))*(POW(4*n-k, K(2.0)) - K(9.0)))/normalisation;
01647 }
01648
01649 return result;
01650 }
01651
01653 R X(modified_sobolev)(const R mu, const int kk)
01654 {
01655 R result;
01656 int kj, k;
01657
01658 for (result = K(0.0), kj = kk; kj <= kk+1; kj++)
01659 {
01660 k = ABS(kj);
01661 if (k == 0)
01662 result += K(1.0);
01663 else
01664 result += POW((double)k,-K(2.0)*mu);
01665 }
01666
01667 return result;
01668 }
01669
01671 R X(modified_multiquadric)(const R mu, const R c, const int kk)
01672 {
01673 R result;
01674 int kj, k;
01675
01676 for (result = K(0.0), kj = kk; kj <= kk+1; kj++)
01677 {
01678 k = ABS(kj);
01679 result += POW((double)(k*k + c*c), -mu);
01680 }
01681
01682 return result;
01683 }
01684
01685 static inline int scaled_modified_bessel_i_series(const R x, const R alpha,
01686 const int nb, const int ize, R *b)
01687 {
01688 const R enmten = K(4.0)*K(R_MIN);
01689 R tempa = K(1.0), empal = K(1.0) + alpha, halfx = K(0.0), tempb = K(0.0);
01690 int n, ncalc = nb;
01691
01692 if (enmten < x)
01693 halfx = x/K(2.0);
01694
01695 if (alpha != K(0.0))
01696 tempa = POW(halfx, alpha)/TGAMMA(empal);
01697
01698 if (ize == 2)
01699 tempa *= EXP(-x);
01700
01701 if (K(1.0) < x + K(1.0))
01702 tempb = halfx*halfx;
01703
01704 b[0] = tempa + tempa*tempb/empal;
01705
01706 if (x != K(0.0) && b[0] == K(0.0))
01707 ncalc = 0;
01708
01709 if (nb == 1)
01710 return ncalc;
01711
01712 if (K(0.0) < x)
01713 {
01714 R tempc = halfx, tover = (enmten + enmten)/x;
01715
01716 if (tempb != K(0.0))
01717 tover = enmten/tempb;
01718
01719 for (n = 1; n < nb; n++)
01720 {
01721 tempa /= empal;
01722 empal += K(1.0);
01723 tempa *= tempc;
01724
01725 if (tempa <= tover*empal)
01726 tempa = K(0.0);
01727
01728 b[n] = tempa + tempa*tempb/empal;
01729
01730 if (b[n] == K(0.0) && n < ncalc)
01731 ncalc = n;
01732 }
01733 }
01734 else
01735 for (n = 1; n < nb; n++)
01736 b[n] = K(0.0);
01737
01738 return ncalc;
01739 }
01740
01741 static inline void scaled_modified_bessel_i_normalize(const R x,
01742 const R alpha, const int nb, const int ize, R *b, const R sum_)
01743 {
01744 const R enmten = K(4.0)*K(R_MIN);
01745 R sum = sum_, tempa;
01746 int n;
01747
01748
01749 if (alpha != K(0.0))
01750 sum = sum * TGAMMA(K(1.0) + alpha) * POW(x/K(2.0), -alpha);
01751
01752 if (ize == 1)
01753 sum *= EXP(-x);
01754
01755 tempa = enmten;
01756
01757 if (K(1.0) < sum)
01758 tempa *= sum;
01759
01760 for (n = 1; n <= nb; n++)
01761 {
01762 if (b[n-1] < tempa)
01763 b[n-1] = K(0.0);
01764
01765 b[n-1] /= sum;
01766 }
01767 }
01768
01816 int nfft_smbi(const R x, const R alpha, const int nb, const int ize, R *b)
01817 {
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838 const int nsig = R_DIG + 2;
01839 const R enten = POW(K(10.0),K(R_MAX_10_EXP));
01840 const R ensig = POW(K(10.0),(R)nsig);
01841 const R rtnsig = POW(K(10.0),-CEIL((R)nsig/K(4.0)));
01842 const R xlarge = K(1E4);
01843 const R exparg = FLOOR(LOG(POW(K(R_RADIX),K(DBL_MAX_EXP-1))));
01844
01845
01846 int l, n, nend, magx, nbmx, ncalc, nstart;
01847 R p, em, en, sum, pold, test, empal, tempa, tempb, tempc, psave, plast, tover,
01848 emp2al, psavel;
01849
01850 magx = LRINT(FLOOR(x));
01851
01852
01853 if ( nb <= 0 || x < K(0.0) || alpha < K(0.0) || K(1.0) <= alpha
01854 || ((ize != 1 || exparg < x) && (ize != 2 || xlarge < x)))
01855 return (MIN(nb,0) - 1);
01856
01857
01858 if (x < rtnsig)
01859 return scaled_modified_bessel_i_series(x,alpha,nb,ize,b);
01860
01861 ncalc = nb;
01862
01863
01864 nbmx = nb - magx;
01865 n = magx + 1;
01866
01867 en = (R) (n+n) + (alpha+alpha);
01868 plast = K(1.0);
01869 p = en/x;
01870
01871
01872 test = ensig + ensig;
01873
01874 if ((5*nsig) < (magx << 1))
01875 test = SQRT(test*p);
01876 else
01877 test /= POW(K(1.585),(R)magx);
01878
01879 if (3 <= nbmx)
01880 {
01881
01882 tover = enten/ensig;
01883 nstart = magx+2;
01884 nend = nb - 1;
01885
01886 for (n = nstart; n <= nend; n++)
01887 {
01888 en += K(2.0);
01889 pold = plast;
01890 plast = p;
01891 p = en*plast/x + pold;
01892 if (p > tover)
01893 {
01894
01895
01896 tover = enten;
01897 p /= tover;
01898 plast /= tover;
01899 psave = p;
01900 psavel = plast;
01901 nstart = n + 1;
01902
01903 do
01904 {
01905 n++;
01906 en += K(2.0);
01907 pold = plast;
01908 plast = p;
01909 p = en*plast/x + pold;
01910 } while (p <= K(1.0));
01911
01912 tempb = en/x;
01913
01914
01915 test = pold*plast*(K(0.5) - K(0.5)/(tempb * tempb))/ensig;
01916 p = plast*tover;
01917 n--;
01918 en -= K(2.0);
01919 nend = MIN(nb,n);
01920
01921 for (ncalc = nstart; ncalc <= nend; ncalc++)
01922 {
01923 pold = psavel;
01924 psavel = psave;
01925 psave = en*psavel/x + pold;
01926 if (test < psave * psavel)
01927 break;
01928 }
01929
01930 ncalc--;
01931 goto L80;
01932 }
01933 }
01934
01935 n = nend;
01936 en = (R) (n+n) + (alpha+alpha);
01937
01938
01939 test = FMAX(test,SQRT(plast*ensig)*SQRT(p+p));
01940 }
01941
01942
01943 do
01944 {
01945 n++;
01946 en += K(2.0);
01947 pold = plast;
01948 plast = p;
01949 p = en*plast/x + pold;
01950 } while (p < test);
01951
01952
01953 L80:
01954 n++;
01955 en += K(2.0);
01956 tempb = K(0.0);
01957 tempa = K(1.0)/p;
01958 em = (R)(n-1);
01959 empal = em + alpha;
01960 emp2al = em - K(1.0) + (alpha+alpha);
01961 sum = tempa*empal*emp2al/em;
01962 nend = n-nb;
01963
01964 if (nend < 0)
01965 {
01966
01967 b[n-1] = tempa;
01968 nend = -nend;
01969 for (l = 1; l <= nend; ++l)
01970 b[n-1 + l] = K(0.0);
01971 }
01972 else
01973 {
01974 if (nend != 0)
01975 {
01976
01977 for (l = 1; l <= nend; ++l)
01978 {
01979 n--;
01980 en -= K(2.0);
01981 tempc = tempb;
01982 tempb = tempa;
01983 tempa = en*tempb/x + tempc;
01984 em -= K(1.0);
01985 emp2al -= K(1.0);
01986
01987 if (n == 1)
01988 break;
01989
01990 if (n == 2)
01991 emp2al = K(1.0);
01992
01993 empal -= K(1.0);
01994 sum = (sum + tempa*empal)*emp2al/em;
01995 }
01996 }
01997
01998
01999 b[n-1] = tempa;
02000
02001 if (nb <= 1)
02002 {
02003 sum = sum + sum + tempa;
02004 scaled_modified_bessel_i_normalize(x,alpha,nb,ize,b,sum);
02005 return ncalc;
02006 }
02007
02008
02009 n--;
02010 en -= 2.0;
02011 b[n-1] = en*tempa/x + tempb;
02012
02013 if (n == 1)
02014 {
02015 sum = sum + sum + b[0];
02016 scaled_modified_bessel_i_normalize(x,alpha,nb,ize,b,sum);
02017 return ncalc;
02018 }
02019
02020 em -= K(1.0);
02021 emp2al -= K(1.0);
02022
02023 if (n == 2)
02024 emp2al = K(1.0);
02025
02026 empal -= K(1.0);
02027 sum = (sum + b[n-1]*empal)*emp2al/em;
02028 }
02029
02030 nend = n - 2;
02031
02032 if (nend != 0)
02033 {
02034
02035 for (l = 1; l <= nend; ++l)
02036 {
02037 n--;
02038 en -= K(2.0);
02039 b[n-1] = en*b[n]/x + b[n+1];
02040 em -= K(1.0);
02041 emp2al -= K(1.0);
02042
02043 if (n == 2)
02044 emp2al = K(1.0);
02045
02046 empal -= K(1.0);
02047 sum = (sum + b[n-1]*empal)*emp2al/em;
02048 }
02049 }
02050
02051
02052 b[0] = K(2.0)*empal*b[1]/x + b[2];
02053 sum = sum + sum + b[0];
02054
02055 scaled_modified_bessel_i_normalize(x,alpha,nb,ize,b,sum);
02056 return ncalc;
02057 }
02058
02059
02060
02061 #if defined(NFFT_LDOUBLE)
02062 #if SIZEOF_LONG_DOUBLE == 16
02063 static const R lanzcos[22] =
02064 {
02065 K(6.992159592766260467081279017133763989752796675052e-10),
02066 K(2.8123669719981441131615873659455441181381798148252),
02067 K(-20.710427744802709693741006921011198460802137523364),
02068 K(68.91604215552677659752923125507927632926035487156),
02069 K(-137.06264623052621370808270752038118534584695692919),
02070 K(181.60793022155827236320731899262871987253012010211),
02071 K(-169.18458743339031548467251800473390300200144720615),
02072 K(114.00119644872831703640105977147389190644587171368),
02073 K(-56.31458215179027632548005714718150178553786189071),
02074 K(20.455656046361896023848762777393331001159409886467),
02075 K(-5.4335605693344340005752804679314670882396415124971),
02076 K(1.0410515397995220235480840944608088445246607731035),
02077 K(-0.14064836221797762435665594700765635768382988137516),
02078 K(0.01295767486030631796076454893320245106478793545067),
02079 K(-0.00077606781904462836224565012148234163222031403439914),
02080 K(0.000028229941587194672222791862511243708412223569786751),
02081 K(-5.6505789508859496251539126911727595443972191305184e-7),
02082 K(5.3642107416477670941393734212506122918693783684933e-9),
02083 K(-1.9054477424755607958337633536133218546180235027896e-11),
02084 K(1.6767560615458281885574906761496245907479292929211e-14),
02085 K(-1.6069159098936011461868918827916292610798729856278e-18),
02086 K(2.092833837734494231412616008353349620313749127668e-24)};
02087 static const R g = K(22.);
02088 #elif SIZEOF_LONG_DOUBLE == 12
02089 static const R lanzcos[16] =
02090 {
02091 K(0.9999999999999999999998151990326234090755978504801206971875240538249),
02092 K(42298.205314163579685884559965825876032898966930749727725308255),
02093 K(-139068.033365719940583314330900395851082150878479801707616019365),
02094 K(181526.445910064425982455538776735270618298440139068097226843822),
02095 K(-119998.518895134887942904852888927609503242910142130425166120992),
02096 K(42498.279665085565336840634629253119065594761580523677606897758),
02097 K(-7861.992932117389960417642813146578778617649628454257358410924),
02098 K(687.805169061068980481146382924730399157260064751317893043547),
02099 K(-23.0082756959823209227391335933370375125546087097458511646070585787792),
02100 K(0.1864230791562731665361253542681001772953257132731654277546996312946),
02101 K(-0.0001094217900494200137669846537980382942701406476771516598942134305),
02102 K(3.8119733228437385220001515981932509556399175176128943412654167182529e-10),
02103 K(-4.7545060695949573669792060004330071744683150710413248650156210893946e-10),
02104 K(2.3162039672883094408744543321924181587375284941017905602597030297461e-10),
02105 K(-5.469699080479773550184561017808390160175857819276162113300096555236e-11),
02106 K(5.761353296381615770264101368671565767532131044262156102602431704e-12)
02107 };
02108 static const R g = K(10.900511);
02109 #else
02110 #error Unsupported size of long double
02111 #endif
02112 #elif defined(NFFT_SINGLE)
02113 static const R lanzcos[5] =
02114 {
02115 0.032648921146387503785008589380436838603631779746936L,
02116 1.1886889953894899130227639793483334640624981688106L,
02117 -1.0684103609197946791553357282032624745327225621591L,
02118 0.18871201360609388869017219868989358810366207478127L,
02119 -0.002745022750941733257487411560167743487986763037079L
02120 };
02121 static const R g = K(4.3408820000000000000000000000000000000000000000000);
02122 #else
02123 static const R lanzcos[13] = {
02124 K(0.006061842346248906525783753964555936883222463665497),
02125 K(1.4256283548148038891120508241039168186977000466332),
02126 K(-2.1475341954959013287000186097849561222515097674559),
02127 K(0.95726691187623991672746418765475493079905679642641),
02128 K(-0.12868865928482566615781458898635486929956555429179),
02129 K(0.0030886434783642720644480000315021785318956671151758),
02130 K(-9.7849610972063267859543785146578805638938457810257e-7),
02131 K(-1.7768848405624834159045007110755286743551855068815e-8),
02132 K(1.9851561978478612433621720534889693635174319844665e-8),
02133 K(-1.2477453878744212747015501336682760715977546596628e-8),
02134 K(5.6096232232009359480838882000975403703116221325714e-9),
02135 K(-1.6133242708801833407734306938668413119166269573061e-9),
02136 K(2.191097792670207503670965523347118763341507005123e-10)};
02137 static const R lanzcosp[12] = {
02138 K(3.0627152859269464002981250710810845422307962730128e8),
02139 K(4.775603155294306744694136459666282665368397175884e8),
02140 K(3.3795733209388037978795882754036968577436873787256e8),
02141 K(1.4313371466803237581908862806209047500008516090763e8),
02142 K(4.026600756736916571780374444713775549949468190092e7),
02143 K(7.8910786118305655779880429845306886138919766854132e6),
02144 K(1.0980346455035079893427073585191634273991944222518e6),
02145 K(108372.57068540043287489056427702756001287882054205),
02146 K(7427.8953058166394494408791055208785413418341223703),
02147 K(336.44978500306352175687290720400758052504709702257),
02148 K(9.0582421437926674355635744116038966456459610306128),
02149 K(0.10976007071323978811079010281977761670657790941656)};
02150 static const R g = K(6.0246800407767295837402343750000000000000000000000);
02151 #endif
02152
02153
02154 static inline R csum(const R z)
02155 {
02156 const short n = sizeof(lanzcos)/sizeof(lanzcos[0]);
02157 R r = K(0.0), y = z + (R)(n-1);
02158 short i;
02159 for (i = n-1; i >= 1; i--, y -= K(1.0))
02160 r += lanzcos[i]/y;
02161 r += lanzcos[0];
02162 return r;
02163 }
02164
02165 static inline R csump(const R z)
02166 {
02167 const short n = sizeof(lanzcos)/sizeof(lanzcos[0]);
02168 R num = K(0.0), denom = K(1.0);
02169 short i;
02170 for (i = n-1; i >= 1; i--)
02171 {
02172 num = z * num + lanzcosp[i-1];
02173 denom *= z + (R)(i);
02174 }
02175 return lanzcos[0] + (num/denom);
02176 }
02177
02178
02179 static DK(KE,2.7182818284590452353602874713526624977572470937000);
02180
02181 R nfft_lambda(const R z, const R eps)
02182 {
02183 R d = K(1.0) - eps, zpg = z + g, emh = eps - K(0.5);
02184 return EXP(-LOG1P(d/(zpg+emh))*(z+emh)) * POW(KE/(zpg+K(0.5)),d)
02185 * (csump(z-d)/csump(z));
02186 }
02187
02188
02189
02190 R nfft_lambda2(const R mu, const R nu)
02191 {
02192 if (mu == K(0.0) || nu == K(0.0))
02193 return K(1.0);
02194 else
02195 return
02196 SQRT(
02197 POW((mu+nu+g+K(0.5))/(K(1.0)*(mu+g+K(0.5))),mu)
02198 * POW((mu+nu+g+K(0.5))/(K(1.0)*(nu+g+K(0.5))),nu)
02199 * SQRT(KE*(mu+nu+g+K(0.5))/((mu+g+K(0.5))*(nu+g+K(0.5))))
02200 * (csum(mu+nu)/(csum(mu)*csum(nu)))
02201 );
02202 }