NFFT 3.5.3alpha
lambda.c
1/*
2 * Copyright (c) 2002, 2017 Jens Keiner, Stefan Kunis, Daniel Potts
3 *
4 * This program is free software; you can redistribute it and/or modify it under
5 * the terms of the GNU General Public License as published by the Free Software
6 * Foundation; either version 2 of the License, or (at your option) any later
7 * version.
8 *
9 * This program is distributed in the hope that it will be useful, but WITHOUT
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12 * details.
13 *
14 * You should have received a copy of the GNU General Public License along with
15 * this program; if not, write to the Free Software Foundation, Inc., 51
16 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17 */
18
19#include "infft.h"
20
21/* Coefficients for Lanzcos's approximation to the Gamma function. Can be
22 * regenerated with Mathematica from file lambda.nb. */
23
24#if defined(NFFT_LDOUBLE)
25 #if LDBL_MANT_DIG > 64
26 /* long double 128 bit wide */
27 #define N 24
28 static const R num[24] =
29 {
30 K(3.035162425359883494754028782232869726547E21),
31 K(3.4967568944064301036001605717507506346E21),
32 K(1.9266526566893208886540195401514595829E21),
33 K(6.755170664882727663160830237424406199E20),
34 K(1.691728531049187527800862627495648317E20),
35 K(3.21979351672256057856444116302160246E19),
36 K(4.8378495427140832493758744745481812E18),
37 K(5.8843103809049324230843820398664955E17),
38 K(5.893958514163405862064178891925630E16),
39 K(4.919561837722192829918665308020810E15),
40 K(3.449165802442404074427531228315120E14),
41 K(2.041330296068782505988459692384726E13),
42 K(1.022234822943784007524609706893119E12),
43 K(4.33137871919821354846952908076307E10),
44 K(1.54921950559667418528481770869280E9),
45 K(4.6544421199876191938054157935810E7),
46 K(1.16527806807504975090675074910053E6),
47 K(24024.759267256769471083727721827),
48 K(400.96500811342195582435806376976),
49 K(5.2829901565447826961703902917085),
50 K(0.05289990244125101024092566765994),
51 K(0.0003783467106547406854542665695934),
52 K(1.7219414217921113919596660801124E-6),
53 K(3.747999317071488557713812635427084359354E-9)
54 };
55/* static const R denom[24] =
56 {
57 K(0.0),
58 K(1124000727777607680000.0),
59 K(4148476779335454720000.0),
60 K(6756146673770930688000.0),
61 K(6548684852703068697600.0),
62 K(4280722865357147142912.0),
63 K(2021687376910682741568.0),
64 K(720308216440924653696.0),
65 K(199321978221066137360.0),
66 K(43714229649594412832.0),
67 K(7707401101297361068.0),
68 K(1103230881185949736.0),
69 K(129006659818331295.0),
70 K(12363045847086207.0),
71 K(971250460939913.0),
72 K(62382416421941.0),
73 K(3256091103430.0),
74 K(136717357942.0),
75 K(4546047198.0),
76 K(116896626.0),
77 K(2240315.0),
78 K(30107.0),
79 K(253.0),
80 K(1.0L)
81 };*/
82 static const R g = K(20.32098218798637390136718750000000000000);
83 #elif LDBL_MANT_DIG == 64
84 /* long double 96 bit wide */
85 #define N 17
86 static const R num[17] =
87 {
88 K(2.715894658327717377557655133124376674911E12),
89 K(3.59017952609791210503852552872112955043E12),
90 K(2.22396659973781496931212735323581871017E12),
91 K(8.5694083451895624818099258668254858834E11),
92 K(2.2988587166874907293359744645339939547E11),
93 K(4.552617168754610815813502794395753410E10),
94 K(6.884887713165178784550917647709216425E9),
95 K(8.11048596140753186476028245385237278E8),
96 K(7.52139159654082231449961362311950170E7),
97 K(5.50924541722426515169752795795495283E6),
98 K(317673.536843541912671493184218236957),
99 K(14268.2798984503552014701437332033752),
100 K(489.361872040326367021390908360178781),
101 K(12.3894133003845444929588321786545861),
102 K(0.218362738950461496394157450728168315),
103 K(0.00239374952205844918669062799606398310),
104 K(0.00001229541408909435212800785616808830746135)
105 };
106/* static const R denom[17] =
107 {
108 K(0.0),
109 K(1307674368000.0),
110 K(4339163001600.0),
111 K(6165817614720.0),
112 K(5056995703824.0),
113 K(2706813345600.0),
114 K(1009672107080.0),
115 K(272803210680.0),
116 K(54631129553.0),
117 K(8207628000.0),
118 K(928095740.0),
119 K(78558480.0),
120 K(4899622.0),
121 K(218400.0),
122 K(6580.0),
123 K(120.0),
124 K(1.0L)
125 };*/
126 static const R g = K(12.22522273659706115722656250000000000000);
127 #else
128 #error Unsupported size of long double
129 #endif
130#elif defined(NFFT_SINGLE)
131 /* float */
132 #define N 6
133 static const R num[6] =
134 {
135 K(14.02614328749964766195705772850038393570),
136 K(43.74732405540314316089531289293124360129),
137 K(50.59547402616588964511581430025589038612),
138 K(26.90456680562548195593733429204228910299),
139 K(6.595765571169314946316366571954421695196),
140 K(0.6007854010515290065101128585795542383721)
141 };
142/* static const R denom[6] =
143 {
144 K(0.0),
145 K(24.0),
146 K(50.0),
147 K(35.0),
148 K(10.0),
149 K(1.0)
150 };*/
151 static const R g = K(1.428456135094165802001953125000000000000);
152#else
153 /* double */
154 #define N 13
155 static const R num[13] =
156 {
157 K(5.690652191347156388090791033559122686859E7),
158 K(1.037940431163445451906271053616070238554E8),
159 K(8.63631312881385914554692728897786842234E7),
160 K(4.33388893246761383477372374059053331609E7),
161 K(1.46055780876850680841416998279135921857E7),
162 K(3.48171215498064590882071018964774556468E6),
163 K(601859.61716810987866702265336993523025),
164 K(75999.293040145426498753034435989091371),
165 K(6955.9996025153761403563101155151989875),
166 K(449.944556906316811944685860765098840962),
167 K(19.5199278824761748284786096623565213621),
168 K(0.509841665565667618812517864480469450999),
169 K(0.006061842346248906525783753964555936883222)
170 };
171/* static const R denom[13] =
172 {
173 K(0.0),
174 K(39916800.0),
175 K(120543840.0),
176 K(150917976.0),
177 K(105258076.0),
178 K(45995730.0),
179 K(13339535.0),
180 K(2637558.0),
181 K(357423.0),
182 K(32670.0),
183 K(1925.0),
184 K(66.0),
185 K(1.0)
186 };*/
187 static const R g = K(6.024680040776729583740234375000000000000);
188#endif
189
190static inline R evaluate_rational(const R z_)
191{
192 R z = z_, s1, s2;
193 INT i;
194
195 if (z <= K(1.0))
196 {
197 s1 = num[N - 1];
198 s2 = K(1.0);
199 for (i = N - 2; i >= 0; --i)
200 {
201 s1 *= z;
202 s2 *= z + (R)(i);
203 s1 += num[i];
204 }
205 }
206 else
207 {
208 z = K(1.0)/z;
209 s1 = num[0];
210 s2 = K(1.0);
211 for (i = 1; i < N; ++i)
212 {
213 s1 *= z;
214 s2 *= K(1.0) + (R)(i-1) * z;
215 s1 += num[i];
216 }
217 }
218 return s1 / s2;
219}
220
221R Y(lambda)(const R z, const R eps)
222{
223 const R d = K(1.0) - eps, zpg = z + g, emh = eps - K(0.5);
224 return EXP(-LOG1P(d / (zpg + emh)) * (z + emh)) *
225 POW(KE / (zpg + K(0.5)),d) *
226 (evaluate_rational(z + eps) / evaluate_rational(z + K(1.0)));
227}
228
229R Y(lambda2)(const R mu, const R nu)
230{
231 if (mu == K(0.0))
232 return K(1.0);
233 else if (nu == K(0.0))
234 return K(1.0);
235 else
236 return
237 SQRT(
238 POW((mu + nu + g + K(0.5)) / (K(1.0) * (mu + g + K(0.5))), mu) *
239 POW((mu + nu + g + K(0.5)) / (K(1.0) * (nu + g + K(0.5))), nu) *
240 SQRT(KE * (mu + nu + g + K(0.5)) /
241 ((mu + g + K(0.5)) * (nu + g + K(0.5)))) *
242 (evaluate_rational(mu + nu + K(1.0)) /
243 (evaluate_rational(mu + K(1.0)) * evaluate_rational(nu + K(1.0))))
244 );
245}
Internal header file for auxiliary definitions and functions.