00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "infft.h"
00022
00023
00024
00025
00026 #if defined(NFFT_LDOUBLE)
00027 #if LDBL_MANT_DIG > 64
00028
00029 #define N 24
00030 static const R num[24] =
00031 {
00032 K(3.035162425359883494754028782232869726547E21),
00033 K(3.4967568944064301036001605717507506346E21),
00034 K(1.9266526566893208886540195401514595829E21),
00035 K(6.755170664882727663160830237424406199E20),
00036 K(1.691728531049187527800862627495648317E20),
00037 K(3.21979351672256057856444116302160246E19),
00038 K(4.8378495427140832493758744745481812E18),
00039 K(5.8843103809049324230843820398664955E17),
00040 K(5.893958514163405862064178891925630E16),
00041 K(4.919561837722192829918665308020810E15),
00042 K(3.449165802442404074427531228315120E14),
00043 K(2.041330296068782505988459692384726E13),
00044 K(1.022234822943784007524609706893119E12),
00045 K(4.33137871919821354846952908076307E10),
00046 K(1.54921950559667418528481770869280E9),
00047 K(4.6544421199876191938054157935810E7),
00048 K(1.16527806807504975090675074910053E6),
00049 K(24024.759267256769471083727721827),
00050 K(400.96500811342195582435806376976),
00051 K(5.2829901565447826961703902917085),
00052 K(0.05289990244125101024092566765994),
00053 K(0.0003783467106547406854542665695934),
00054 K(1.7219414217921113919596660801124E-6),
00055 K(3.747999317071488557713812635427084359354E-9)
00056 };
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084 static const R g = K(20.32098218798637390136718750000000000000);
00085 #elif LDBL_MANT_DIG == 64
00086
00087 #define N 17
00088 static const R num[17] =
00089 {
00090 K(2.715894658327717377557655133124376674911E12),
00091 K(3.59017952609791210503852552872112955043E12),
00092 K(2.22396659973781496931212735323581871017E12),
00093 K(8.5694083451895624818099258668254858834E11),
00094 K(2.2988587166874907293359744645339939547E11),
00095 K(4.552617168754610815813502794395753410E10),
00096 K(6.884887713165178784550917647709216425E9),
00097 K(8.11048596140753186476028245385237278E8),
00098 K(7.52139159654082231449961362311950170E7),
00099 K(5.50924541722426515169752795795495283E6),
00100 K(317673.536843541912671493184218236957),
00101 K(14268.2798984503552014701437332033752),
00102 K(489.361872040326367021390908360178781),
00103 K(12.3894133003845444929588321786545861),
00104 K(0.218362738950461496394157450728168315),
00105 K(0.00239374952205844918669062799606398310),
00106 K(0.00001229541408909435212800785616808830746135)
00107 };
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128 static const R g = K(12.22522273659706115722656250000000000000);
00129 #else
00130 #error Unsupported size of long double
00131 #endif
00132 #elif defined(NFFT_SINGLE)
00133
00134 #define N 6
00135 static const R num[6] =
00136 {
00137 K(14.02614328749964766195705772850038393570),
00138 K(43.74732405540314316089531289293124360129),
00139 K(50.59547402616588964511581430025589038612),
00140 K(26.90456680562548195593733429204228910299),
00141 K(6.595765571169314946316366571954421695196),
00142 K(0.6007854010515290065101128585795542383721)
00143 };
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153 static const R g = K(1.428456135094165802001953125000000000000);
00154 #else
00155
00156 #define N 13
00157 static const R num[13] =
00158 {
00159 K(5.690652191347156388090791033559122686859E7),
00160 K(1.037940431163445451906271053616070238554E8),
00161 K(8.63631312881385914554692728897786842234E7),
00162 K(4.33388893246761383477372374059053331609E7),
00163 K(1.46055780876850680841416998279135921857E7),
00164 K(3.48171215498064590882071018964774556468E6),
00165 K(601859.61716810987866702265336993523025),
00166 K(75999.293040145426498753034435989091371),
00167 K(6955.9996025153761403563101155151989875),
00168 K(449.944556906316811944685860765098840962),
00169 K(19.5199278824761748284786096623565213621),
00170 K(0.509841665565667618812517864480469450999),
00171 K(0.006061842346248906525783753964555936883222)
00172 };
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189 static const R g = K(6.024680040776729583740234375000000000000);
00190 #endif
00191
00192 static inline R evaluate_rational(const R z_)
00193 {
00194 R z = z_, s1, s2;
00195 int i;
00196
00197 if (z <= K(1.0))
00198 {
00199 s1 = num[N - 1];
00200 s2 = K(1.0);
00201 for (i = N - 2; i >= 0; --i)
00202 {
00203 s1 *= z;
00204 s2 *= (z + i);
00205 s1 += num[i];
00206 }
00207 }
00208 else
00209 {
00210 z = K(1.0)/z;
00211 s1 = num[0];
00212 s2 = K(1.0);
00213 for (i = 1; i < N; ++i)
00214 {
00215 s1 *= z;
00216 s2 *= K(1.0) + (i-1)*z;
00217 s1 += num[i];
00218 }
00219 }
00220 return s1 / s2;
00221 }
00222
00223 R X(lambda)(const R z, const R eps)
00224 {
00225 const R d = K(1.0) - eps, zpg = z + g, emh = eps - K(0.5);
00226 return EXP(-LOG1P(d / (zpg + emh)) * (z + emh)) *
00227 POW(KE / (zpg + K(0.5)),d) *
00228 (evaluate_rational(z + eps) / evaluate_rational(z + K(1.0)));
00229 }
00230
00231 R X(lambda2)(const R mu, const R nu)
00232 {
00233 if (mu == K(0.0))
00234 return K(1.0);
00235 else if (nu == K(0.0))
00236 return K(1.0);
00237 else
00238 return
00239 SQRT(
00240 POW((mu + nu + g + K(0.5)) / (K(1.0) * (mu + g + K(0.5))), mu) *
00241 POW((mu + nu + g + K(0.5)) / (K(1.0) * (nu + g + K(0.5))), nu) *
00242 SQRT(KE * (mu + nu + g + K(0.5)) /
00243 ((mu + g + K(0.5)) * (nu + g + K(0.5)))) *
00244 (evaluate_rational(mu + nu + K(1.0)) /
00245 (evaluate_rational(mu + K(1.0)) * evaluate_rational(nu + K(1.0))))
00246 );
00247 }