![]() |
NFFT
3.3.1
|
00001 /* 00002 * Copyright (c) 2002, 2016 Jens Keiner, Stefan Kunis, Daniel Potts 00003 * 00004 * This program is free software; you can redistribute it and/or modify it under 00005 * the terms of the GNU General Public License as published by the Free Software 00006 * Foundation; either version 2 of the License, or (at your option) any later 00007 * version. 00008 * 00009 * This program is distributed in the hope that it will be useful, but WITHOUT 00010 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 00011 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 00012 * details. 00013 * 00014 * You should have received a copy of the GNU General Public License along with 00015 * this program; if not, write to the Free Software Foundation, Inc., 51 00016 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 00017 */ 00018 00019 #include "infft.h" 00020 00021 /* Coefficients for Lanzcos's approximation to the Gamma function. Can be 00022 * regenerated with Mathematica from file lambda.nb. */ 00023 00024 #if defined(NFFT_LDOUBLE) 00025 #if LDBL_MANT_DIG > 64 00026 /* long double 128 bit wide */ 00027 #define N 24 00028 static const R num[24] = 00029 { 00030 K(3.035162425359883494754028782232869726547E21), 00031 K(3.4967568944064301036001605717507506346E21), 00032 K(1.9266526566893208886540195401514595829E21), 00033 K(6.755170664882727663160830237424406199E20), 00034 K(1.691728531049187527800862627495648317E20), 00035 K(3.21979351672256057856444116302160246E19), 00036 K(4.8378495427140832493758744745481812E18), 00037 K(5.8843103809049324230843820398664955E17), 00038 K(5.893958514163405862064178891925630E16), 00039 K(4.919561837722192829918665308020810E15), 00040 K(3.449165802442404074427531228315120E14), 00041 K(2.041330296068782505988459692384726E13), 00042 K(1.022234822943784007524609706893119E12), 00043 K(4.33137871919821354846952908076307E10), 00044 K(1.54921950559667418528481770869280E9), 00045 K(4.6544421199876191938054157935810E7), 00046 K(1.16527806807504975090675074910053E6), 00047 K(24024.759267256769471083727721827), 00048 K(400.96500811342195582435806376976), 00049 K(5.2829901565447826961703902917085), 00050 K(0.05289990244125101024092566765994), 00051 K(0.0003783467106547406854542665695934), 00052 K(1.7219414217921113919596660801124E-6), 00053 K(3.747999317071488557713812635427084359354E-9) 00054 }; 00055 /* static const R denom[24] = 00056 { 00057 K(0.0), 00058 K(1124000727777607680000.0), 00059 K(4148476779335454720000.0), 00060 K(6756146673770930688000.0), 00061 K(6548684852703068697600.0), 00062 K(4280722865357147142912.0), 00063 K(2021687376910682741568.0), 00064 K(720308216440924653696.0), 00065 K(199321978221066137360.0), 00066 K(43714229649594412832.0), 00067 K(7707401101297361068.0), 00068 K(1103230881185949736.0), 00069 K(129006659818331295.0), 00070 K(12363045847086207.0), 00071 K(971250460939913.0), 00072 K(62382416421941.0), 00073 K(3256091103430.0), 00074 K(136717357942.0), 00075 K(4546047198.0), 00076 K(116896626.0), 00077 K(2240315.0), 00078 K(30107.0), 00079 K(253.0), 00080 K(1.0L) 00081 };*/ 00082 static const R g = K(20.32098218798637390136718750000000000000); 00083 #elif LDBL_MANT_DIG == 64 00084 /* long double 96 bit wide */ 00085 #define N 17 00086 static const R num[17] = 00087 { 00088 K(2.715894658327717377557655133124376674911E12), 00089 K(3.59017952609791210503852552872112955043E12), 00090 K(2.22396659973781496931212735323581871017E12), 00091 K(8.5694083451895624818099258668254858834E11), 00092 K(2.2988587166874907293359744645339939547E11), 00093 K(4.552617168754610815813502794395753410E10), 00094 K(6.884887713165178784550917647709216425E9), 00095 K(8.11048596140753186476028245385237278E8), 00096 K(7.52139159654082231449961362311950170E7), 00097 K(5.50924541722426515169752795795495283E6), 00098 K(317673.536843541912671493184218236957), 00099 K(14268.2798984503552014701437332033752), 00100 K(489.361872040326367021390908360178781), 00101 K(12.3894133003845444929588321786545861), 00102 K(0.218362738950461496394157450728168315), 00103 K(0.00239374952205844918669062799606398310), 00104 K(0.00001229541408909435212800785616808830746135) 00105 }; 00106 /* static const R denom[17] = 00107 { 00108 K(0.0), 00109 K(1307674368000.0), 00110 K(4339163001600.0), 00111 K(6165817614720.0), 00112 K(5056995703824.0), 00113 K(2706813345600.0), 00114 K(1009672107080.0), 00115 K(272803210680.0), 00116 K(54631129553.0), 00117 K(8207628000.0), 00118 K(928095740.0), 00119 K(78558480.0), 00120 K(4899622.0), 00121 K(218400.0), 00122 K(6580.0), 00123 K(120.0), 00124 K(1.0L) 00125 };*/ 00126 static const R g = K(12.22522273659706115722656250000000000000); 00127 #else 00128 #error Unsupported size of long double 00129 #endif 00130 #elif defined(NFFT_SINGLE) 00131 /* float */ 00132 #define N 6 00133 static const R num[6] = 00134 { 00135 K(14.02614328749964766195705772850038393570), 00136 K(43.74732405540314316089531289293124360129), 00137 K(50.59547402616588964511581430025589038612), 00138 K(26.90456680562548195593733429204228910299), 00139 K(6.595765571169314946316366571954421695196), 00140 K(0.6007854010515290065101128585795542383721) 00141 }; 00142 /* static const R denom[6] = 00143 { 00144 K(0.0), 00145 K(24.0), 00146 K(50.0), 00147 K(35.0), 00148 K(10.0), 00149 K(1.0) 00150 };*/ 00151 static const R g = K(1.428456135094165802001953125000000000000); 00152 #else 00153 /* double */ 00154 #define N 13 00155 static const R num[13] = 00156 { 00157 K(5.690652191347156388090791033559122686859E7), 00158 K(1.037940431163445451906271053616070238554E8), 00159 K(8.63631312881385914554692728897786842234E7), 00160 K(4.33388893246761383477372374059053331609E7), 00161 K(1.46055780876850680841416998279135921857E7), 00162 K(3.48171215498064590882071018964774556468E6), 00163 K(601859.61716810987866702265336993523025), 00164 K(75999.293040145426498753034435989091371), 00165 K(6955.9996025153761403563101155151989875), 00166 K(449.944556906316811944685860765098840962), 00167 K(19.5199278824761748284786096623565213621), 00168 K(0.509841665565667618812517864480469450999), 00169 K(0.006061842346248906525783753964555936883222) 00170 }; 00171 /* static const R denom[13] = 00172 { 00173 K(0.0), 00174 K(39916800.0), 00175 K(120543840.0), 00176 K(150917976.0), 00177 K(105258076.0), 00178 K(45995730.0), 00179 K(13339535.0), 00180 K(2637558.0), 00181 K(357423.0), 00182 K(32670.0), 00183 K(1925.0), 00184 K(66.0), 00185 K(1.0) 00186 };*/ 00187 static const R g = K(6.024680040776729583740234375000000000000); 00188 #endif 00189 00190 static inline R evaluate_rational(const R z_) 00191 { 00192 R z = z_, s1, s2; 00193 INT i; 00194 00195 if (z <= K(1.0)) 00196 { 00197 s1 = num[N - 1]; 00198 s2 = K(1.0); 00199 for (i = N - 2; i >= 0; --i) 00200 { 00201 s1 *= z; 00202 s2 *= z + (R)(i); 00203 s1 += num[i]; 00204 } 00205 } 00206 else 00207 { 00208 z = K(1.0)/z; 00209 s1 = num[0]; 00210 s2 = K(1.0); 00211 for (i = 1; i < N; ++i) 00212 { 00213 s1 *= z; 00214 s2 *= K(1.0) + (R)(i-1) * z; 00215 s1 += num[i]; 00216 } 00217 } 00218 return s1 / s2; 00219 } 00220 00221 R Y(lambda)(const R z, const R eps) 00222 { 00223 const R d = K(1.0) - eps, zpg = z + g, emh = eps - K(0.5); 00224 return EXP(-LOG1P(d / (zpg + emh)) * (z + emh)) * 00225 POW(KE / (zpg + K(0.5)),d) * 00226 (evaluate_rational(z + eps) / evaluate_rational(z + K(1.0))); 00227 } 00228 00229 R Y(lambda2)(const R mu, const R nu) 00230 { 00231 if (mu == K(0.0)) 00232 return K(1.0); 00233 else if (nu == K(0.0)) 00234 return K(1.0); 00235 else 00236 return 00237 SQRT( 00238 POW((mu + nu + g + K(0.5)) / (K(1.0) * (mu + g + K(0.5))), mu) * 00239 POW((mu + nu + g + K(0.5)) / (K(1.0) * (nu + g + K(0.5))), nu) * 00240 SQRT(KE * (mu + nu + g + K(0.5)) / 00241 ((mu + g + K(0.5)) * (nu + g + K(0.5)))) * 00242 (evaluate_rational(mu + nu + K(1.0)) / 00243 (evaluate_rational(mu + K(1.0)) * evaluate_rational(nu + K(1.0)))) 00244 ); 00245 }