NFFT Logo 3.2.2
lambda.c
00001 /*
00002  * Copyright (c) 2002, 2012 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 /* $Id: util.c 3483 2010-04-23 19:02:34Z keiner $ */
00020 
00021 #include "infft.h"
00022 
00023 /* Coefficients for Lanzcos's approximation to the Gamma function. Can be
00024  * regenerated with Mathematica from file lambda.nb. */
00025 
00026 #if defined(NFFT_LDOUBLE)
00027   #if LDBL_MANT_DIG > 64
00028     /* long double 128 bit wide */
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 /*    static const R denom[24] =
00058     {
00059       K(0.0),
00060       K(1124000727777607680000.0),
00061       K(4148476779335454720000.0),
00062       K(6756146673770930688000.0),
00063       K(6548684852703068697600.0),
00064       K(4280722865357147142912.0),
00065       K(2021687376910682741568.0),
00066       K(720308216440924653696.0),
00067       K(199321978221066137360.0),
00068       K(43714229649594412832.0),
00069       K(7707401101297361068.0),
00070       K(1103230881185949736.0),
00071       K(129006659818331295.0),
00072       K(12363045847086207.0),
00073       K(971250460939913.0),
00074       K(62382416421941.0),
00075       K(3256091103430.0),
00076       K(136717357942.0),
00077       K(4546047198.0),
00078       K(116896626.0),
00079       K(2240315.0),
00080       K(30107.0),
00081       K(253.0),
00082       K(1.0L)
00083     };*/
00084     static const R g = K(20.32098218798637390136718750000000000000);
00085   #elif LDBL_MANT_DIG == 64
00086     /* long double 96 bit wide */
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 /*    static const R denom[17] =
00109     {
00110       K(0.0),
00111       K(1307674368000.0),
00112       K(4339163001600.0),
00113       K(6165817614720.0),
00114       K(5056995703824.0),
00115       K(2706813345600.0),
00116       K(1009672107080.0),
00117       K(272803210680.0),
00118       K(54631129553.0),
00119       K(8207628000.0),
00120       K(928095740.0),
00121       K(78558480.0),
00122       K(4899622.0),
00123       K(218400.0),
00124       K(6580.0),
00125       K(120.0),
00126       K(1.0L)
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   /* float */
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 /*  static const R denom[6] =
00145   {
00146     K(0.0),
00147     K(24.0),
00148     K(50.0),
00149     K(35.0),
00150     K(10.0),
00151     K(1.0)
00152   };*/
00153   static const R g = K(1.428456135094165802001953125000000000000);
00154 #else
00155   /* double */
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 /*  static const R denom[13] =
00174   {
00175     K(0.0),
00176     K(39916800.0),
00177     K(120543840.0),
00178     K(150917976.0),
00179     K(105258076.0),
00180     K(45995730.0),
00181     K(13339535.0),
00182     K(2637558.0),
00183     K(357423.0),
00184     K(32670.0),
00185     K(1925.0),
00186     K(66.0),
00187     K(1.0)
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 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409