Alex Bikfalvi
SimStream Documentation
Gamma.cpp
00001 #include "Headers.h" 00002 #include "Gamma.h" 00003 #include "Psi.h" 00004 #include "Chebyshev.h" 00005 #include "Exp.h" 00006 00007 static struct {int n; double f; long i; } fact_table[MATH_SF_FACT_NMAX + 1] = { 00008 { 0, 1.0, 1L }, 00009 { 1, 1.0, 1L }, 00010 { 2, 2.0, 2L }, 00011 { 3, 6.0, 6L }, 00012 { 4, 24.0, 24L }, 00013 { 5, 120.0, 120L }, 00014 { 6, 720.0, 720L }, 00015 { 7, 5040.0, 5040L }, 00016 { 8, 40320.0, 40320L }, 00017 00018 { 9, 362880.0, 362880L }, 00019 { 10, 3628800.0, 3628800L }, 00020 { 11, 39916800.0, 39916800L }, 00021 { 12, 479001600.0, 479001600L }, 00022 00023 { 13, 6227020800.0, 0 }, 00024 { 14, 87178291200.0, 0 }, 00025 { 15, 1307674368000.0, 0 }, 00026 { 16, 20922789888000.0, 0 }, 00027 { 17, 355687428096000.0, 0 }, 00028 { 18, 6402373705728000.0, 0 }, 00029 { 19, 121645100408832000.0, 0 }, 00030 { 20, 2432902008176640000.0, 0 }, 00031 { 21, 51090942171709440000.0, 0 }, 00032 { 22, 1124000727777607680000.0, 0 }, 00033 { 23, 25852016738884976640000.0, 0 }, 00034 { 24, 620448401733239439360000.0, 0 }, 00035 { 25, 15511210043330985984000000.0, 0 }, 00036 { 26, 403291461126605635584000000.0, 0 }, 00037 { 27, 10888869450418352160768000000.0, 0 }, 00038 { 28, 304888344611713860501504000000.0, 0 }, 00039 { 29, 8841761993739701954543616000000.0, 0 }, 00040 { 30, 265252859812191058636308480000000.0, 0 }, 00041 { 31, 8222838654177922817725562880000000.0, 0 }, 00042 { 32, 263130836933693530167218012160000000.0, 0 }, 00043 { 33, 8683317618811886495518194401280000000.0, 0 }, 00044 { 34, 2.95232799039604140847618609644e38, 0 }, 00045 { 35, 1.03331479663861449296666513375e40, 0 }, 00046 { 36, 3.71993326789901217467999448151e41, 0 }, 00047 { 37, 1.37637530912263450463159795816e43, 0 }, 00048 { 38, 5.23022617466601111760007224100e44, 0 }, 00049 { 39, 2.03978820811974433586402817399e46, 0 }, 00050 { 40, 8.15915283247897734345611269600e47, 0 }, 00051 { 41, 3.34525266131638071081700620534e49, 0 }, 00052 { 42, 1.40500611775287989854314260624e51, 0 }, 00053 { 43, 6.04152630633738356373551320685e52, 0 }, 00054 { 44, 2.65827157478844876804362581101e54, 0 }, 00055 { 45, 1.19622220865480194561963161496e56, 0 }, 00056 { 46, 5.50262215981208894985030542880e57, 0 }, 00057 { 47, 2.58623241511168180642964355154e59, 0 }, 00058 { 48, 1.24139155925360726708622890474e61, 0 }, 00059 { 49, 6.08281864034267560872252163321e62, 0 }, 00060 { 50, 3.04140932017133780436126081661e64, 0 }, 00061 { 51, 1.55111875328738228022424301647e66, 0 }, 00062 { 52, 8.06581751709438785716606368564e67, 0 }, 00063 { 53, 4.27488328406002556429801375339e69, 0 }, 00064 { 54, 2.30843697339241380472092742683e71, 0 }, 00065 { 55, 1.26964033536582759259651008476e73, 0 }, 00066 { 56, 7.10998587804863451854045647464e74, 0 }, 00067 { 57, 4.05269195048772167556806019054e76, 0 }, 00068 { 58, 2.35056133128287857182947491052e78, 0 }, 00069 { 59, 1.38683118545689835737939019720e80, 0 }, 00070 { 60, 8.32098711274139014427634118320e81, 0 }, 00071 { 61, 5.07580213877224798800856812177e83, 0 }, 00072 { 62, 3.14699732603879375256531223550e85, 0 }, 00073 { 63, 1.982608315404440064116146708360e87, 0 }, 00074 { 64, 1.268869321858841641034333893350e89, 0 }, 00075 { 65, 8.247650592082470666723170306800e90, 0 }, 00076 { 66, 5.443449390774430640037292402480e92, 0 }, 00077 { 67, 3.647111091818868528824985909660e94, 0 }, 00078 { 68, 2.480035542436830599600990418570e96, 0 }, 00079 { 69, 1.711224524281413113724683388810e98, 0 }, 00080 { 70, 1.197857166996989179607278372170e100, 0 }, 00081 { 71, 8.504785885678623175211676442400e101, 0 }, 00082 { 72, 6.123445837688608686152407038530e103, 0 }, 00083 { 73, 4.470115461512684340891257138130e105, 0 }, 00084 { 74, 3.307885441519386412259530282210e107, 0 }, 00085 { 75, 2.480914081139539809194647711660e109, 0 }, 00086 { 76, 1.885494701666050254987932260860e111, 0 }, 00087 { 77, 1.451830920282858696340707840860e113, 0 }, 00088 { 78, 1.132428117820629783145752115870e115, 0 }, 00089 { 79, 8.946182130782975286851441715400e116, 0 }, 00090 { 80, 7.156945704626380229481153372320e118, 0 }, 00091 { 81, 5.797126020747367985879734231580e120, 0 }, 00092 { 82, 4.753643337012841748421382069890e122, 0 }, 00093 { 83, 3.945523969720658651189747118010e124, 0 }, 00094 { 84, 3.314240134565353266999387579130e126, 0 }, 00095 { 85, 2.817104114380550276949479442260e128, 0 }, 00096 { 86, 2.422709538367273238176552320340e130, 0 }, 00097 { 87, 2.107757298379527717213600518700e132, 0 }, 00098 { 88, 1.854826422573984391147968456460e134, 0 }, 00099 { 89, 1.650795516090846108121691926250e136, 0 }, 00100 { 90, 1.485715964481761497309522733620e138, 0 }, 00101 { 91, 1.352001527678402962551665687590e140, 0 }, 00102 { 92, 1.243841405464130725547532432590e142, 0 }, 00103 { 93, 1.156772507081641574759205162310e144, 0 }, 00104 { 94, 1.087366156656743080273652852570e146, 0 }, 00105 { 95, 1.032997848823905926259970209940e148, 0 }, 00106 { 96, 9.916779348709496892095714015400e149, 0 }, 00107 { 97, 9.619275968248211985332842594960e151, 0 }, 00108 { 98, 9.426890448883247745626185743100e153, 0 }, 00109 { 99, 9.332621544394415268169923885600e155, 0 }, 00110 { 100, 9.33262154439441526816992388563e157, 0 }, 00111 { 101, 9.42594775983835942085162312450e159, 0 }, 00112 { 102, 9.61446671503512660926865558700e161, 0 }, 00113 { 103, 9.90290071648618040754671525458e163, 0 }, 00114 { 104, 1.02990167451456276238485838648e166, 0 }, 00115 { 105, 1.08139675824029090050410130580e168, 0 }, 00116 { 106, 1.146280563734708354534347384148e170, 0 }, 00117 { 107, 1.226520203196137939351751701040e172, 0 }, 00118 { 108, 1.324641819451828974499891837120e174, 0 }, 00119 { 109, 1.443859583202493582204882102460e176, 0 }, 00120 { 110, 1.588245541522742940425370312710e178, 0 }, 00121 { 111, 1.762952551090244663872161047110e180, 0 }, 00122 { 112, 1.974506857221074023536820372760e182, 0 }, 00123 { 113, 2.231192748659813646596607021220e184, 0 }, 00124 { 114, 2.543559733472187557120132004190e186, 0 }, 00125 { 115, 2.925093693493015690688151804820e188, 0 }, 00126 { 116, 3.393108684451898201198256093590e190, 0 }, 00127 { 117, 3.96993716080872089540195962950e192, 0 }, 00128 { 118, 4.68452584975429065657431236281e194, 0 }, 00129 { 119, 5.57458576120760588132343171174e196, 0 }, 00130 { 120, 6.68950291344912705758811805409e198, 0 }, 00131 { 121, 8.09429852527344373968162284545e200, 0 }, 00132 { 122, 9.87504420083360136241157987140e202, 0 }, 00133 { 123, 1.21463043670253296757662432419e205, 0 }, 00134 { 124, 1.50614174151114087979501416199e207, 0 }, 00135 { 125, 1.88267717688892609974376770249e209, 0 }, 00136 { 126, 2.37217324288004688567714730514e211, 0 }, 00137 { 127, 3.01266001845765954480997707753e213, 0 }, 00138 { 128, 3.85620482362580421735677065923e215, 0 }, 00139 { 129, 4.97450422247728744039023415041e217, 0 }, 00140 { 130, 6.46685548922047367250730439554e219, 0 }, 00141 { 131, 8.47158069087882051098456875820e221, 0 }, 00142 { 132, 1.11824865119600430744996307608e224, 0 }, 00143 { 133, 1.48727070609068572890845089118e226, 0 }, 00144 { 134, 1.99294274616151887673732419418e228, 0 }, 00145 { 135, 2.69047270731805048359538766215e230, 0 }, 00146 { 136, 3.65904288195254865768972722052e232, 0 }, 00147 { 137, 5.01288874827499166103492629211e234, 0 }, 00148 { 138, 6.91778647261948849222819828311e236, 0 }, 00149 { 139, 9.61572319694108900419719561353e238, 0 }, 00150 { 140, 1.34620124757175246058760738589e241, 0 }, 00151 { 141, 1.89814375907617096942852641411e243, 0 }, 00152 { 142, 2.69536413788816277658850750804e245, 0 }, 00153 { 143, 3.85437071718007277052156573649e247, 0 }, 00154 { 144, 5.55029383273930478955105466055e249, 0 }, 00155 { 145, 8.04792605747199194484902925780e251, 0 }, 00156 { 146, 1.17499720439091082394795827164e254, 0 }, 00157 { 147, 1.72724589045463891120349865931e256, 0 }, 00158 { 148, 2.55632391787286558858117801578e258, 0 }, 00159 { 149, 3.80892263763056972698595524351e260, 0 }, 00160 { 150, 5.71338395644585459047893286526e262, 0 }, 00161 { 151, 8.62720977423324043162318862650e264, 0 }, 00162 { 152, 1.31133588568345254560672467123e267, 0 }, 00163 { 153, 2.00634390509568239477828874699e269, 0 }, 00164 { 154, 3.08976961384735088795856467036e271, 0 }, 00165 { 155, 4.78914290146339387633577523906e273, 0 }, 00166 { 156, 7.47106292628289444708380937294e275, 0 }, 00167 { 157, 1.17295687942641442819215807155e278, 0 }, 00168 { 158, 1.85327186949373479654360975305e280, 0 }, 00169 { 159, 2.94670227249503832650433950735e282, 0 }, 00170 { 160, 4.71472363599206132240694321176e284, 0 }, 00171 { 161, 7.59070505394721872907517857094e286, 0 }, 00172 { 162, 1.22969421873944943411017892849e289, 0 }, 00173 { 163, 2.00440157654530257759959165344e291, 0 }, 00174 { 164, 3.28721858553429622726333031164e293, 0 }, 00175 { 165, 5.42391066613158877498449501421e295, 0 }, 00176 { 166, 9.00369170577843736647426172359e297, 0 }, 00177 { 167, 1.50361651486499904020120170784e300, 0 }, 00178 { 168, 2.52607574497319838753801886917e302, 0 }, 00179 { 169, 4.26906800900470527493925188890e304, 0 }, 00180 { 170, 7.25741561530799896739672821113e306, 0 }, 00181 00182 /* 00183 { 171, 1.24101807021766782342484052410e309, 0 }, 00184 { 172, 2.13455108077438865629072570146e311, 0 }, 00185 { 173, 3.69277336973969237538295546352e313, 0 }, 00186 { 174, 6.42542566334706473316634250653e315, 0 }, 00187 { 175, 1.12444949108573632830410993864e318, 0 }, 00188 { 176, 1.97903110431089593781523349201e320, 0 }, 00189 { 177, 3.50288505463028580993296328086e322, 0 }, 00190 { 178, 6.23513539724190874168067463993e324, 0 }, 00191 { 179, 1.11608923610630166476084076055e327, 0 }, 00192 { 180, 2.00896062499134299656951336898e329, 0 }, 00193 { 181, 3.63621873123433082379081919786e331, 0 }, 00194 { 182, 6.61791809084648209929929094011e333, 0 }, 00195 { 183, 1.21107901062490622417177024204e336, 0 }, 00196 { 184, 2.22838537954982745247605724535e338, 0 }, 00197 { 185, 4.12251295216718078708070590390e340, 0 }, 00198 { 186, 7.66787409103095626397011298130e342, 0 }, 00199 { 187, 1.43389245502278882136241112750e345, 0 }, 00200 { 188, 2.69571781544284298416133291969e347, 0 }, 00201 { 189, 5.09490667118697324006491921822e349, 0 }, 00202 { 190, 9.68032267525524915612334651460e351, 0 }, 00203 { 191, 1.84894163097375258881955918429e354, 0 }, 00204 { 192, 3.54996793146960497053355363384e356, 0 }, 00205 { 193, 6.85143810773633759312975851330e358, 0 }, 00206 { 194, 1.32917899290084949306717315158e361, 0 }, 00207 { 195, 2.59189903615665651148098764559e363, 0 }, 00208 { 196, 5.08012211086704676250273578535e365, 0 }, 00209 { 197, 1.00078405584080821221303894971e368, 0 }, 00210 { 198, 1.98155243056480026018181712043e370, 0 }, 00211 { 199, 3.94328933682395251776181606966e372, 0 }, 00212 { 200, 7.88657867364790503552363213932e374, 0 } 00213 */ 00214 }; 00215 00216 /* Chebyshev coefficients for Gamma*(3/4(t+1)+1/2), -1<t<1 00217 */ 00218 static double gstar_a_data[30] = { 00219 2.16786447866463034423060819465, 00220 -0.05533249018745584258035832802, 00221 0.01800392431460719960888319748, 00222 -0.00580919269468937714480019814, 00223 0.00186523689488400339978881560, 00224 -0.00059746524113955531852595159, 00225 0.00019125169907783353925426722, 00226 -0.00006124996546944685735909697, 00227 0.00001963889633130842586440945, 00228 -6.3067741254637180272515795142e-06, 00229 2.0288698405861392526872789863e-06, 00230 -6.5384896660838465981983750582e-07, 00231 2.1108698058908865476480734911e-07, 00232 -6.8260714912274941677892994580e-08, 00233 2.2108560875880560555583978510e-08, 00234 -7.1710331930255456643627187187e-09, 00235 2.3290892983985406754602564745e-09, 00236 -7.5740371598505586754890405359e-10, 00237 2.4658267222594334398525312084e-10, 00238 -8.0362243171659883803428749516e-11, 00239 2.6215616826341594653521346229e-11, 00240 -8.5596155025948750540420068109e-12, 00241 2.7970831499487963614315315444e-12, 00242 -9.1471771211886202805502562414e-13, 00243 2.9934720198063397094916415927e-13, 00244 -9.8026575909753445931073620469e-14, 00245 3.2116773667767153777571410671e-14, 00246 -1.0518035333878147029650507254e-14, 00247 3.4144405720185253938994854173e-15, 00248 -1.0115153943081187052322643819e-15 00249 }; 00250 00251 00252 /* Chebyshev coefficients for 00253 * x^2(Gamma*(x) - 1 - 1/(12x)), x = 4(t+1)+2, -1 < t < 1 00254 */ 00255 static double gstar_b_data[] = { 00256 0.0057502277273114339831606096782, 00257 0.0004496689534965685038254147807, 00258 -0.0001672763153188717308905047405, 00259 0.0000615137014913154794776670946, 00260 -0.0000223726551711525016380862195, 00261 8.0507405356647954540694800545e-06, 00262 -2.8671077107583395569766746448e-06, 00263 1.0106727053742747568362254106e-06, 00264 -3.5265558477595061262310873482e-07, 00265 1.2179216046419401193247254591e-07, 00266 -4.1619640180795366971160162267e-08, 00267 1.4066283500795206892487241294e-08, 00268 -4.6982570380537099016106141654e-09, 00269 1.5491248664620612686423108936e-09, 00270 -5.0340936319394885789686867772e-10, 00271 1.6084448673736032249959475006e-10, 00272 -5.0349733196835456497619787559e-11, 00273 1.5357154939762136997591808461e-11, 00274 -4.5233809655775649997667176224e-12, 00275 1.2664429179254447281068538964e-12, 00276 -3.2648287937449326771785041692e-13, 00277 7.1528272726086133795579071407e-14, 00278 -9.4831735252566034505739531258e-15, 00279 -2.3124001991413207293120906691e-15, 00280 2.8406613277170391482590129474e-15, 00281 -1.7245370321618816421281770927e-15, 00282 8.6507923128671112154695006592e-16, 00283 -3.9506563665427555895391869919e-16, 00284 1.6779342132074761078792361165e-16, 00285 -6.0483153034414765129837716260e-17 00286 }; 00287 00288 00289 /* coefficients for gamma=7, kmax=8 Lanczos method */ 00290 static double lanczos_7_c[9] = { 00291 0.99999999999980993227684700473478, 00292 676.520368121885098567009190444019, 00293 -1259.13921672240287047156078755283, 00294 771.3234287776530788486528258894, 00295 -176.61502916214059906584551354, 00296 12.507343278686904814458936853, 00297 -0.13857109526572011689554707, 00298 9.984369578019570859563e-6, 00299 1.50563273514931155834e-7 00300 }; 00301 00302 00303 /* Lanczos method for real x > 0; 00304 * gamma=7, truncated at 1/(z+8) 00305 * [J. SIAM Numer. Anal, Ser. B, 1 (1964) 86] 00306 */ 00307 static int LnGammaLanczos(double x, MathResult * result) 00308 { 00309 int k; 00310 double Ag; 00311 double term1, term2; 00312 00313 x -= 1.0; /* Lanczos writes z! instead of Gamma(z) */ 00314 00315 Ag = lanczos_7_c[0]; 00316 for(k=1; k<=8; k++) { Ag += lanczos_7_c[k]/(x+k); } 00317 00318 /* (x+0.5)*log(x+7.5) - (x+7.5) + LogRootTwoPi_ + log(Ag(x)) */ 00319 term1 = (x+0.5)*log((x+7.5)/M_E); 00320 term2 = LogRootTwoPi_ + log(Ag); 00321 result->val = term1 + (term2 - 7.0); 00322 result->err = 2.0 * MATH_DBL_EPSILON * (fabs(term1) + fabs(term2) + 7.0); 00323 result->err += MATH_DBL_EPSILON * fabs(result->val); 00324 00325 return MATH_SUCCESS; 00326 } 00327 00328 /* x = eps near zero 00329 * gives double-precision for |eps| < 0.02 00330 */ 00331 static int LnGammaSgn0(double eps, MathResult * lng, double * sgn) 00332 { 00333 /* calculate series for g(eps) = Gamma(eps) eps - 1/(1+eps) - eps/2 */ 00334 const double c1 = -0.07721566490153286061; 00335 const double c2 = -0.01094400467202744461; 00336 const double c3 = 0.09252092391911371098; 00337 const double c4 = -0.01827191316559981266; 00338 const double c5 = 0.01800493109685479790; 00339 const double c6 = -0.00685088537872380685; 00340 const double c7 = 0.00399823955756846603; 00341 const double c8 = -0.00189430621687107802; 00342 const double c9 = 0.00097473237804513221; 00343 const double c10 = -0.00048434392722255893; 00344 const double g6 = c6+eps*(c7+eps*(c8 + eps*(c9 + eps*c10))); 00345 const double g = eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*g6))))); 00346 00347 /* calculate Gamma(eps) eps, a positive quantity */ 00348 const double gee = g + 1.0/(1.0+eps) + 0.5*eps; 00349 00350 lng->val = log(gee/fabs(eps)); 00351 lng->err = 4.0 * MATH_DBL_EPSILON * fabs(lng->val); 00352 *sgn = MATH_SIGN(eps); 00353 00354 return MATH_SUCCESS; 00355 } 00356 00357 00358 /* x near a negative integer 00359 * Calculates sign as well as log(|gamma(x)|). 00360 * x = -N + eps 00361 * assumes N >= 1 00362 */ 00363 static int LnGammaSgnSing(int N, double eps, MathResult * lng, double * sgn) 00364 { 00365 if(eps == 0.0) 00366 { 00367 lng->val = 0.0; 00368 lng->err = 0.0; 00369 *sgn = 0.0; 00370 MATH_ERROR ("error", MATH_EDOM); 00371 } 00372 else if(N == 1) 00373 { 00374 /* calculate series for 00375 * g = eps gamma(-1+eps) + 1 + eps/2 (1+3eps)/(1-eps^2) 00376 * double-precision for |eps| < 0.02 00377 */ 00378 const double c0 = 0.07721566490153286061; 00379 const double c1 = 0.08815966957356030521; 00380 const double c2 = -0.00436125434555340577; 00381 const double c3 = 0.01391065882004640689; 00382 const double c4 = -0.00409427227680839100; 00383 const double c5 = 0.00275661310191541584; 00384 const double c6 = -0.00124162645565305019; 00385 const double c7 = 0.00065267976121802783; 00386 const double c8 = -0.00032205261682710437; 00387 const double c9 = 0.00016229131039545456; 00388 const double g5 = c5 + eps*(c6 + eps*(c7 + eps*(c8 + eps*c9))); 00389 const double g = eps*(c0 + eps*(c1 + eps*(c2 + eps*(c3 + eps*(c4 + eps*g5))))); 00390 00391 /* calculate eps gamma(-1+eps), a negative quantity */ 00392 const double gam_e = g - 1.0 - 0.5*eps*(1.0+3.0*eps)/(1.0 - eps*eps); 00393 00394 lng->val = log(fabs(gam_e)/fabs(eps)); 00395 lng->err = 2.0 * MATH_DBL_EPSILON * fabs(lng->val); 00396 *sgn = ( eps > 0.0 ? -1.0 : 1.0 ); 00397 return MATH_SUCCESS; 00398 } 00399 else 00400 { 00401 double g; 00402 00403 /* series for sin(Pi(N+1-eps))/(Pi eps) modulo the sign 00404 * double-precision for |eps| < 0.02 00405 */ 00406 const double cs1 = -1.6449340668482264365; 00407 const double cs2 = 0.8117424252833536436; 00408 const double cs3 = -0.1907518241220842137; 00409 const double cs4 = 0.0261478478176548005; 00410 const double cs5 = -0.0023460810354558236; 00411 const double e2 = eps*eps; 00412 const double sin_ser = 1.0 + e2*(cs1+e2*(cs2+e2*(cs3+e2*(cs4+e2*cs5)))); 00413 00414 /* calculate series for ln(gamma(1+N-eps)) 00415 * double-precision for |eps| < 0.02 00416 */ 00417 double aeps = fabs(eps); 00418 double c1, c2, c3, c4, c5, c6, c7; 00419 double lng_ser; 00420 MathResult c0; 00421 MathResult psi_0; 00422 MathResult psi_1; 00423 MathResult psi_2; 00424 MathResult psi_3; 00425 MathResult psi_4; 00426 MathResult psi_5; 00427 MathResult psi_6; 00428 psi_2.val = 0.0; 00429 psi_3.val = 0.0; 00430 psi_4.val = 0.0; 00431 psi_5.val = 0.0; 00432 psi_6.val = 0.0; 00433 LnFact(N, &c0); 00434 PsiInt(N+1, &psi_0); 00435 Psi1Int(N+1, &psi_1); 00436 if(aeps > 0.00001) PsiN(2, N+1.0, &psi_2); 00437 if(aeps > 0.0002) PsiN(3, N+1.0, &psi_3); 00438 if(aeps > 0.001) PsiN(4, N+1.0, &psi_4); 00439 if(aeps > 0.005) PsiN(5, N+1.0, &psi_5); 00440 if(aeps > 0.01) PsiN(6, N+1.0, &psi_6); 00441 c1 = psi_0.val; 00442 c2 = psi_1.val/2.0; 00443 c3 = psi_2.val/6.0; 00444 c4 = psi_3.val/24.0; 00445 c5 = psi_4.val/120.0; 00446 c6 = psi_5.val/720.0; 00447 c7 = psi_6.val/5040.0; 00448 lng_ser = c0.val-eps*(c1-eps*(c2-eps*(c3-eps*(c4-eps*(c5-eps*(c6-eps*c7)))))); 00449 00450 /* calculate 00451 * g = ln(|eps gamma(-N+eps)|) 00452 * = -ln(gamma(1+N-eps)) + ln(|eps Pi/sin(Pi(N+1+eps))|) 00453 */ 00454 g = -lng_ser - log(sin_ser); 00455 00456 lng->val = g - log(fabs(eps)); 00457 lng->err = c0.err + 2.0 * MATH_DBL_EPSILON * (fabs(g) + fabs(lng->val)); 00458 00459 *sgn = ( MATH_IS_ODD(N) ? -1.0 : 1.0 ) * ( eps > 0.0 ? 1.0 : -1.0 ); 00460 00461 return MATH_SUCCESS; 00462 } 00463 } 00464 00465 inline static int LnGamma1Pade(const double eps, MathResult * result) 00466 { 00467 /* Use (2,2) Pade for Log[Gamma[1+eps]]/eps 00468 * plus a correction series. 00469 */ 00470 const double n1 = -1.0017419282349508699871138440; 00471 const double n2 = 1.7364839209922879823280541733; 00472 const double d1 = 1.2433006018858751556055436011; 00473 const double d2 = 5.0456274100274010152489597514; 00474 const double num = (eps + n1) * (eps + n2); 00475 const double den = (eps + d1) * (eps + d2); 00476 const double pade = 2.0816265188662692474880210318 * num / den; 00477 const double c0 = 0.004785324257581753; 00478 const double c1 = -0.01192457083645441; 00479 const double c2 = 0.01931961413960498; 00480 const double c3 = -0.02594027398725020; 00481 const double c4 = 0.03141928755021455; 00482 const double eps5 = eps*eps*eps*eps*eps; 00483 const double corr = eps5 * (c0 + eps*(c1 + eps*(c2 + eps*(c3 + c4*eps)))); 00484 result->val = eps * (pade + corr); 00485 result->err = 2.0 * MATH_DBL_EPSILON * fabs(result->val); 00486 return MATH_SUCCESS; 00487 } 00488 00489 inline static int LnGamma2Pade(const double eps, MathResult * result) 00490 { 00491 /* Use (2,2) Pade for Log[Gamma[2+eps]]/eps 00492 * plus a correction series. 00493 */ 00494 const double n1 = 1.000895834786669227164446568; 00495 const double n2 = 4.209376735287755081642901277; 00496 const double d1 = 2.618851904903217274682578255; 00497 const double d2 = 10.85766559900983515322922936; 00498 const double num = (eps + n1) * (eps + n2); 00499 const double den = (eps + d1) * (eps + d2); 00500 const double pade = 2.85337998765781918463568869 * num/den; 00501 const double c0 = 0.0001139406357036744; 00502 const double c1 = -0.0001365435269792533; 00503 const double c2 = 0.0001067287169183665; 00504 const double c3 = -0.0000693271800931282; 00505 const double c4 = 0.0000407220927867950; 00506 const double eps5 = eps*eps*eps*eps*eps; 00507 const double corr = eps5 * (c0 + eps*(c1 + eps*(c2 + eps*(c3 + c4*eps)))); 00508 result->val = eps * (pade + corr); 00509 result->err = 2.0 * MATH_DBL_EPSILON * fabs(result->val); 00510 return MATH_SUCCESS; 00511 } 00512 00513 00514 /* series for gammastar(x) 00515 * double-precision for x > 10.0 00516 */ 00517 static int GammaStarSeries(const double x, MathResult * result) 00518 { 00519 /* Use the Stirling series for the correction to Log(Gamma(x)), 00520 * which is better behaved and easier to compute than the 00521 * regular Stirling series for Gamma(x). 00522 */ 00523 const double y = 1.0/(x*x); 00524 const double c0 = 1.0/12.0; 00525 const double c1 = -1.0/360.0; 00526 const double c2 = 1.0/1260.0; 00527 const double c3 = -1.0/1680.0; 00528 const double c4 = 1.0/1188.0; 00529 const double c5 = -691.0/360360.0; 00530 const double c6 = 1.0/156.0; 00531 const double c7 = -3617.0/122400.0; 00532 const double ser = c0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*(c5 + y*(c6 + y*c7)))))); 00533 result->val = exp(ser/x); 00534 result->err = 2.0 * MATH_DBL_EPSILON * result->val * MATH_MAX_DBL(1.0, ser/x); 00535 return MATH_SUCCESS; 00536 } 00537 00538 00539 /* Chebyshev expansion for log(gamma(x)/gamma(8)) 00540 * 5 < x < 10 00541 * -1 < t < 1 00542 */ 00543 static double gamma_5_10_data[24] = { 00544 -1.5285594096661578881275075214, 00545 4.8259152300595906319768555035, 00546 0.2277712320977614992970601978, 00547 -0.0138867665685617873604917300, 00548 0.0012704876495201082588139723, 00549 -0.0001393841240254993658962470, 00550 0.0000169709242992322702260663, 00551 -2.2108528820210580075775889168e-06, 00552 3.0196602854202309805163918716e-07, 00553 -4.2705675000079118380587357358e-08, 00554 6.2026423818051402794663551945e-09, 00555 -9.1993973208880910416311405656e-10, 00556 1.3875551258028145778301211638e-10, 00557 -2.1218861491906788718519522978e-11, 00558 3.2821736040381439555133562600e-12, 00559 -5.1260001009953791220611135264e-13, 00560 8.0713532554874636696982146610e-14, 00561 -1.2798522376569209083811628061e-14, 00562 2.0417711600852502310258808643e-15, 00563 -3.2745239502992355776882614137e-16, 00564 5.2759418422036579482120897453e-17, 00565 -8.5354147151695233960425725513e-18, 00566 1.3858639703888078291599886143e-18, 00567 -2.2574398807738626571560124396e-19 00568 }; 00569 00570 static const ChebSeries gamma_5_10_cs = { 00571 gamma_5_10_data, 00572 23, 00573 -1, 1, 00574 11 00575 }; 00576 00577 00578 /* gamma(x) for x >= 1/2 00579 * assumes x >= 1/2 00580 */ 00581 static int GammaXgtHalf(const double x, MathResult * result) 00582 { 00583 /* CHECK_POINTER(result) */ 00584 00585 if(x == 0.5) 00586 { 00587 result->val = 1.77245385090551602729817; 00588 result->err = MATH_DBL_EPSILON * result->val; 00589 return MATH_SUCCESS; 00590 } 00591 else if (x <= (MATH_SF_FACT_NMAX + 1.0) && x == floor(x)) 00592 { 00593 int n = (int) floor (x); 00594 result->val = fact_table[n - 1].f; 00595 result->err = MATH_DBL_EPSILON * result->val; 00596 return MATH_SUCCESS; 00597 } 00598 else if(fabs(x - 1.0) < 0.01) 00599 { 00600 /* Use series for Gamma[1+eps] - 1/(1+eps). 00601 */ 00602 const double eps = x - 1.0; 00603 const double c1 = 0.4227843350984671394; 00604 const double c2 = -0.01094400467202744461; 00605 const double c3 = 0.09252092391911371098; 00606 const double c4 = -0.018271913165599812664; 00607 const double c5 = 0.018004931096854797895; 00608 const double c6 = -0.006850885378723806846; 00609 const double c7 = 0.003998239557568466030; 00610 result->val = 1.0/x + eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*(c6+eps*c7)))))); 00611 result->err = MATH_DBL_EPSILON; 00612 return MATH_SUCCESS; 00613 } 00614 else if(fabs(x - 2.0) < 0.01) 00615 { 00616 /* Use series for Gamma[1 + eps]. 00617 */ 00618 const double eps = x - 2.0; 00619 const double c1 = 0.4227843350984671394; 00620 const double c2 = 0.4118403304264396948; 00621 const double c3 = 0.08157691924708626638; 00622 const double c4 = 0.07424901075351389832; 00623 const double c5 = -0.00026698206874501476832; 00624 const double c6 = 0.011154045718130991049; 00625 const double c7 = -0.002852645821155340816; 00626 const double c8 = 0.0021039333406973880085; 00627 result->val = 1.0 + eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*(c6+eps*(c7+eps*c8))))))); 00628 result->err = MATH_DBL_EPSILON; 00629 return MATH_SUCCESS; 00630 } 00631 else if(x < 5.0) 00632 { 00633 /* Exponentiating the logarithm is fine, as 00634 * long as the exponential is not so large 00635 * that it greatly amplifies the error. 00636 */ 00637 MathResult lg; 00638 LnGammaLanczos(x, &lg); 00639 result->val = exp(lg.val); 00640 result->err = result->val * (lg.err + 2.0 * MATH_DBL_EPSILON); 00641 return MATH_SUCCESS; 00642 } 00643 else if(x < 10.0) 00644 { 00645 /* This is a sticky area. The logarithm 00646 * is too large and the gammastar series 00647 * is not good. 00648 */ 00649 const double gamma_8 = 5040.0; 00650 const double t = (2.0*x - 15.0)/5.0; 00651 MathResult c; 00652 ChebEval(&gamma_5_10_cs, t, &c); 00653 result->val = exp(c.val) * gamma_8; 00654 result->err = result->val * c.err; 00655 result->err += 2.0 * MATH_DBL_EPSILON * result->val; 00656 return MATH_SUCCESS; 00657 } 00658 else if(x < MATH_SF_GAMMA_XMAX) 00659 { 00660 /* We do not want to exponentiate the logarithm 00661 * if x is large because of the inevitable 00662 * inflation of the error. So we carefully 00663 * use pow() and exp() with exact quantities. 00664 */ 00665 double p = pow(x, 0.5*x); 00666 double e = exp(-x); 00667 double q = (p * e) * p; 00668 double pre = M_SQRT2 * M_SQRTPI * q/sqrt(x); 00669 MathResult gstar; 00670 int stat_gs = GammaStarSeries(x, &gstar); 00671 result->val = pre * gstar.val; 00672 result->err = (x + 2.5) * MATH_DBL_EPSILON * result->val; 00673 return stat_gs; 00674 } 00675 else 00676 { 00677 OVERFLOW_ERROR(result); 00678 } 00679 } 00680 00681 int LnGammaE(double x, MathResult * result) 00682 { 00683 /* CHECK_POINTER(result) */ 00684 00685 if(fabs(x - 1.0) < 0.01) 00686 { 00687 /* Note that we must amplify the errors 00688 * from the Pade evaluations because of 00689 * the way we must pass the argument, i.e. 00690 * writing (1-x) is a loss of precision 00691 * when x is near 1. 00692 */ 00693 int stat = LnGamma1Pade(x - 1.0, result); 00694 result->err *= 1.0/(MATH_DBL_EPSILON + fabs(x - 1.0)); 00695 return stat; 00696 } 00697 else if(fabs(x - 2.0) < 0.01) 00698 { 00699 int stat = LnGamma2Pade(x - 2.0, result); 00700 result->err *= 1.0/(MATH_DBL_EPSILON + fabs(x - 2.0)); 00701 return stat; 00702 } 00703 else if(x >= 0.5) 00704 { 00705 return LnGammaLanczos(x, result); 00706 } 00707 else if(x == 0.0) 00708 { 00709 DOMAIN_ERROR(result); 00710 } 00711 else if(fabs(x) < 0.02) 00712 { 00713 double sgn; 00714 return LnGammaSgn0(x, result, &sgn); 00715 } 00716 else if(x > -0.5/(MATH_DBL_EPSILON*M_PI)) 00717 { 00718 /* Try to extract a fractional 00719 * part from x. 00720 */ 00721 double z = 1.0 - x; 00722 double s = sin(M_PI*z); 00723 double as = fabs(s); 00724 if(s == 0.0) 00725 { 00726 DOMAIN_ERROR(result); 00727 } 00728 else if(as < M_PI*0.015) 00729 { 00730 /* x is near a negative integer, -N */ 00731 if(x < INT_MIN + 2.0) 00732 { 00733 result->val = 0.0; 00734 result->err = 0.0; 00735 MATH_ERROR ("error", MATH_EROUND); 00736 } 00737 else 00738 { 00739 int N = -(int)(x - 0.5); 00740 double eps = x + N; 00741 double sgn; 00742 return LnGammaSgnSing(N, eps, result, &sgn); 00743 } 00744 } 00745 else 00746 { 00747 MathResult lg_z; 00748 LnGammaLanczos(z, &lg_z); 00749 result->val = M_LNPI - (log(as) + lg_z.val); 00750 result->err = 2.0 * MATH_DBL_EPSILON * fabs(result->val) + lg_z.err; 00751 return MATH_SUCCESS; 00752 } 00753 } 00754 else 00755 { 00756 /* |x| was too large to extract any fractional part */ 00757 result->val = 0.0; 00758 result->err = 0.0; 00759 MATH_ERROR ("error", MATH_EROUND); 00760 } 00761 } 00762 00763 00764 int LnGammaSgnE(double x, MathResult * result_lg, double * sgn) 00765 { 00766 if(fabs(x - 1.0) < 0.01) 00767 { 00768 int stat = LnGamma1Pade(x - 1.0, result_lg); 00769 result_lg->err *= 1.0/(MATH_DBL_EPSILON + fabs(x - 1.0)); 00770 *sgn = 1.0; 00771 return stat; 00772 } 00773 else if(fabs(x - 2.0) < 0.01) 00774 { 00775 int stat = LnGamma2Pade(x - 2.0, result_lg); 00776 result_lg->err *= 1.0/(MATH_DBL_EPSILON + fabs(x - 2.0)); 00777 *sgn = 1.0; 00778 return stat; 00779 } 00780 else if(x >= 0.5) 00781 { 00782 *sgn = 1.0; 00783 return LnGammaLanczos(x, result_lg); 00784 } 00785 else if(x == 0.0) 00786 { 00787 *sgn = 0.0; 00788 DOMAIN_ERROR(result_lg); 00789 } 00790 else if(fabs(x) < 0.02) 00791 { 00792 return LnGammaSgn0(x, result_lg, sgn); 00793 } 00794 else if(x > -0.5/(MATH_DBL_EPSILON*M_PI)) 00795 { 00796 /* Try to extract a fractional 00797 * part from x. 00798 */ 00799 double z = 1.0 - x; 00800 double s = sin(M_PI*x); 00801 double as = fabs(s); 00802 if(s == 0.0) 00803 { 00804 *sgn = 0.0; 00805 DOMAIN_ERROR(result_lg); 00806 } 00807 else if(as < M_PI*0.015) 00808 { 00809 /* x is near a negative integer, -N */ 00810 if(x < INT_MIN + 2.0) 00811 { 00812 result_lg->val = 0.0; 00813 result_lg->err = 0.0; 00814 *sgn = 0.0; 00815 MATH_ERROR ("error", MATH_EROUND); 00816 } 00817 else 00818 { 00819 int N = -(int)(x - 0.5); 00820 double eps = x + N; 00821 return LnGammaSgnSing(N, eps, result_lg, sgn); 00822 } 00823 } 00824 else 00825 { 00826 MathResult lg_z; 00827 LnGammaLanczos(z, &lg_z); 00828 *sgn = (s > 0.0 ? 1.0 : -1.0); 00829 result_lg->val = M_LNPI - (log(as) + lg_z.val); 00830 result_lg->err = 2.0 * MATH_DBL_EPSILON * fabs(result_lg->val) + lg_z.err; 00831 return MATH_SUCCESS; 00832 } 00833 } 00834 else 00835 { 00836 /* |x| was too large to extract any fractional part */ 00837 result_lg->val = 0.0; 00838 result_lg->err = 0.0; 00839 *sgn = 0.0; 00840 MATH_ERROR ("x too large to extract fraction part", MATH_EROUND); 00841 } 00842 } 00843 00844 00845 int GammaE(const double x, MathResult * result) 00846 { 00847 if(x < 0.5) 00848 { 00849 int rint_x = (int)floor(x+0.5); 00850 double f_x = x - rint_x; 00851 double sgn_gamma = ( MATH_IS_EVEN(rint_x) ? 1.0 : -1.0 ); 00852 double sin_term = sgn_gamma * sin(M_PI * f_x) / M_PI; 00853 00854 if(sin_term == 0.0) 00855 { 00856 DOMAIN_ERROR(result); 00857 } 00858 else if(x > -169.0) 00859 { 00860 MathResult g; 00861 GammaXgtHalf(1.0-x, &g); 00862 if(fabs(sin_term) * g.val * MATH_DBL_MIN < 1.0) 00863 { 00864 result->val = 1.0/(sin_term * g.val); 00865 result->err = fabs(g.err/g.val) * fabs(result->val); 00866 result->err += 2.0 * MATH_DBL_EPSILON * fabs(result->val); 00867 return MATH_SUCCESS; 00868 } 00869 else 00870 { 00871 UNDERFLOW_ERROR(result); 00872 } 00873 } 00874 else 00875 { 00876 /* It is hard to control it here. 00877 * We can only exponentiate the 00878 * logarithm and eat the loss of 00879 * precision. 00880 */ 00881 MathResult lng; 00882 double sgn; 00883 int stat_lng = LnGammaSgnE(x, &lng, &sgn); 00884 int stat_e = ExpMultErr(lng.val, lng.err, sgn, 0.0, result); 00885 return MATH_ERROR_SELECT_2(stat_e, stat_lng); 00886 } 00887 } 00888 else 00889 { 00890 return GammaXgtHalf(x, result); 00891 } 00892 } 00893 00894 int LnFact(const unsigned int n, MathResult * result) 00895 { 00896 /* CHECK_POINTER(result) */ 00897 00898 if(n <= MATH_SF_FACT_NMAX){ 00899 result->val = log(fact_table[n].f); 00900 result->err = 2.0 * MATH_DBL_EPSILON * fabs(result->val); 00901 return MATH_SUCCESS; 00902 } 00903 else { 00904 LnGammaE(n+1.0, result); 00905 return MATH_SUCCESS; 00906 } 00907 } 00908 00909 double LnGamma(const double x) 00910 { 00911 EVAL_RESULT(LnGammaE(x, &result)); 00912 } 00913 00914 double Gamma(const double x) 00915 { 00916 EVAL_RESULT(GammaE(x, &result)); 00917 }
Last updated: February 8, 2011