1 /* specfunc/gamma.c 2 * 3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman 4 * 5 * This program is free software; you can redistribute it and/or modify 6 * it under the terms of the GNU General Public License as published by 7 * the Free Software Foundation; either version 3 of the License, or (at 8 * your option) any later version. 9 * 10 * This program is distributed in the hope that it will be useful, but 11 * WITHOUT ANY WARRANTY; without even the implied warranty of 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 13 * General Public License for more details. 14 * 15 * You should have received a copy of the GNU General Public License 16 * along with this program; if not, write to the Free Software 17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 18 */ 19 20 /* Author: G. Jungman */ 21 22 #include "gsl__config.h" 23 #include "gsl_math.h" 24 #include "gsl_errno.h" 25 #include "gsl_sf_exp.h" 26 #include "gsl_sf_log.h" 27 #include "gsl_sf_psi.h" 28 #include "gsl_sf_trig.h" 29 #include "gsl_sf_gamma.h" 30 31 #include "gsl_specfunc__error.h" 32 #include "gsl_specfunc__check.h" 33 34 #include "gsl_specfunc__chebyshev.h" 35 #include "gsl_specfunc__cheb_eval.c" 36 37 #define LogRootTwoPi_ 0.9189385332046727418 38 39 40 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/ 41 42 static struct {int n; double f; long i; } fact_table[GSL_SF_FACT_NMAX + 1] = { 43 { 0, 1.0, 1L }, 44 { 1, 1.0, 1L }, 45 { 2, 2.0, 2L }, 46 { 3, 6.0, 6L }, 47 { 4, 24.0, 24L }, 48 { 5, 120.0, 120L }, 49 { 6, 720.0, 720L }, 50 { 7, 5040.0, 5040L }, 51 { 8, 40320.0, 40320L }, 52 53 { 9, 362880.0, 362880L }, 54 { 10, 3628800.0, 3628800L }, 55 { 11, 39916800.0, 39916800L }, 56 { 12, 479001600.0, 479001600L }, 57 58 { 13, 6227020800.0, 0 }, 59 { 14, 87178291200.0, 0 }, 60 { 15, 1307674368000.0, 0 }, 61 { 16, 20922789888000.0, 0 }, 62 { 17, 355687428096000.0, 0 }, 63 { 18, 6402373705728000.0, 0 }, 64 { 19, 121645100408832000.0, 0 }, 65 { 20, 2432902008176640000.0, 0 }, 66 { 21, 51090942171709440000.0, 0 }, 67 { 22, 1124000727777607680000.0, 0 }, 68 { 23, 25852016738884976640000.0, 0 }, 69 { 24, 620448401733239439360000.0, 0 }, 70 { 25, 15511210043330985984000000.0, 0 }, 71 { 26, 403291461126605635584000000.0, 0 }, 72 { 27, 10888869450418352160768000000.0, 0 }, 73 { 28, 304888344611713860501504000000.0, 0 }, 74 { 29, 8841761993739701954543616000000.0, 0 }, 75 { 30, 265252859812191058636308480000000.0, 0 }, 76 { 31, 8222838654177922817725562880000000.0, 0 }, 77 { 32, 263130836933693530167218012160000000.0, 0 }, 78 { 33, 8683317618811886495518194401280000000.0, 0 }, 79 { 34, 2.95232799039604140847618609644e38, 0 }, 80 { 35, 1.03331479663861449296666513375e40, 0 }, 81 { 36, 3.71993326789901217467999448151e41, 0 }, 82 { 37, 1.37637530912263450463159795816e43, 0 }, 83 { 38, 5.23022617466601111760007224100e44, 0 }, 84 { 39, 2.03978820811974433586402817399e46, 0 }, 85 { 40, 8.15915283247897734345611269600e47, 0 }, 86 { 41, 3.34525266131638071081700620534e49, 0 }, 87 { 42, 1.40500611775287989854314260624e51, 0 }, 88 { 43, 6.04152630633738356373551320685e52, 0 }, 89 { 44, 2.65827157478844876804362581101e54, 0 }, 90 { 45, 1.19622220865480194561963161496e56, 0 }, 91 { 46, 5.50262215981208894985030542880e57, 0 }, 92 { 47, 2.58623241511168180642964355154e59, 0 }, 93 { 48, 1.24139155925360726708622890474e61, 0 }, 94 { 49, 6.08281864034267560872252163321e62, 0 }, 95 { 50, 3.04140932017133780436126081661e64, 0 }, 96 { 51, 1.55111875328738228022424301647e66, 0 }, 97 { 52, 8.06581751709438785716606368564e67, 0 }, 98 { 53, 4.27488328406002556429801375339e69, 0 }, 99 { 54, 2.30843697339241380472092742683e71, 0 }, 100 { 55, 1.26964033536582759259651008476e73, 0 }, 101 { 56, 7.10998587804863451854045647464e74, 0 }, 102 { 57, 4.05269195048772167556806019054e76, 0 }, 103 { 58, 2.35056133128287857182947491052e78, 0 }, 104 { 59, 1.38683118545689835737939019720e80, 0 }, 105 { 60, 8.32098711274139014427634118320e81, 0 }, 106 { 61, 5.07580213877224798800856812177e83, 0 }, 107 { 62, 3.14699732603879375256531223550e85, 0 }, 108 { 63, 1.982608315404440064116146708360e87, 0 }, 109 { 64, 1.268869321858841641034333893350e89, 0 }, 110 { 65, 8.247650592082470666723170306800e90, 0 }, 111 { 66, 5.443449390774430640037292402480e92, 0 }, 112 { 67, 3.647111091818868528824985909660e94, 0 }, 113 { 68, 2.480035542436830599600990418570e96, 0 }, 114 { 69, 1.711224524281413113724683388810e98, 0 }, 115 { 70, 1.197857166996989179607278372170e100, 0 }, 116 { 71, 8.504785885678623175211676442400e101, 0 }, 117 { 72, 6.123445837688608686152407038530e103, 0 }, 118 { 73, 4.470115461512684340891257138130e105, 0 }, 119 { 74, 3.307885441519386412259530282210e107, 0 }, 120 { 75, 2.480914081139539809194647711660e109, 0 }, 121 { 76, 1.885494701666050254987932260860e111, 0 }, 122 { 77, 1.451830920282858696340707840860e113, 0 }, 123 { 78, 1.132428117820629783145752115870e115, 0 }, 124 { 79, 8.946182130782975286851441715400e116, 0 }, 125 { 80, 7.156945704626380229481153372320e118, 0 }, 126 { 81, 5.797126020747367985879734231580e120, 0 }, 127 { 82, 4.753643337012841748421382069890e122, 0 }, 128 { 83, 3.945523969720658651189747118010e124, 0 }, 129 { 84, 3.314240134565353266999387579130e126, 0 }, 130 { 85, 2.817104114380550276949479442260e128, 0 }, 131 { 86, 2.422709538367273238176552320340e130, 0 }, 132 { 87, 2.107757298379527717213600518700e132, 0 }, 133 { 88, 1.854826422573984391147968456460e134, 0 }, 134 { 89, 1.650795516090846108121691926250e136, 0 }, 135 { 90, 1.485715964481761497309522733620e138, 0 }, 136 { 91, 1.352001527678402962551665687590e140, 0 }, 137 { 92, 1.243841405464130725547532432590e142, 0 }, 138 { 93, 1.156772507081641574759205162310e144, 0 }, 139 { 94, 1.087366156656743080273652852570e146, 0 }, 140 { 95, 1.032997848823905926259970209940e148, 0 }, 141 { 96, 9.916779348709496892095714015400e149, 0 }, 142 { 97, 9.619275968248211985332842594960e151, 0 }, 143 { 98, 9.426890448883247745626185743100e153, 0 }, 144 { 99, 9.332621544394415268169923885600e155, 0 }, 145 { 100, 9.33262154439441526816992388563e157, 0 }, 146 { 101, 9.42594775983835942085162312450e159, 0 }, 147 { 102, 9.61446671503512660926865558700e161, 0 }, 148 { 103, 9.90290071648618040754671525458e163, 0 }, 149 { 104, 1.02990167451456276238485838648e166, 0 }, 150 { 105, 1.08139675824029090050410130580e168, 0 }, 151 { 106, 1.146280563734708354534347384148e170, 0 }, 152 { 107, 1.226520203196137939351751701040e172, 0 }, 153 { 108, 1.324641819451828974499891837120e174, 0 }, 154 { 109, 1.443859583202493582204882102460e176, 0 }, 155 { 110, 1.588245541522742940425370312710e178, 0 }, 156 { 111, 1.762952551090244663872161047110e180, 0 }, 157 { 112, 1.974506857221074023536820372760e182, 0 }, 158 { 113, 2.231192748659813646596607021220e184, 0 }, 159 { 114, 2.543559733472187557120132004190e186, 0 }, 160 { 115, 2.925093693493015690688151804820e188, 0 }, 161 { 116, 3.393108684451898201198256093590e190, 0 }, 162 { 117, 3.96993716080872089540195962950e192, 0 }, 163 { 118, 4.68452584975429065657431236281e194, 0 }, 164 { 119, 5.57458576120760588132343171174e196, 0 }, 165 { 120, 6.68950291344912705758811805409e198, 0 }, 166 { 121, 8.09429852527344373968162284545e200, 0 }, 167 { 122, 9.87504420083360136241157987140e202, 0 }, 168 { 123, 1.21463043670253296757662432419e205, 0 }, 169 { 124, 1.50614174151114087979501416199e207, 0 }, 170 { 125, 1.88267717688892609974376770249e209, 0 }, 171 { 126, 2.37217324288004688567714730514e211, 0 }, 172 { 127, 3.01266001845765954480997707753e213, 0 }, 173 { 128, 3.85620482362580421735677065923e215, 0 }, 174 { 129, 4.97450422247728744039023415041e217, 0 }, 175 { 130, 6.46685548922047367250730439554e219, 0 }, 176 { 131, 8.47158069087882051098456875820e221, 0 }, 177 { 132, 1.11824865119600430744996307608e224, 0 }, 178 { 133, 1.48727070609068572890845089118e226, 0 }, 179 { 134, 1.99294274616151887673732419418e228, 0 }, 180 { 135, 2.69047270731805048359538766215e230, 0 }, 181 { 136, 3.65904288195254865768972722052e232, 0 }, 182 { 137, 5.01288874827499166103492629211e234, 0 }, 183 { 138, 6.91778647261948849222819828311e236, 0 }, 184 { 139, 9.61572319694108900419719561353e238, 0 }, 185 { 140, 1.34620124757175246058760738589e241, 0 }, 186 { 141, 1.89814375907617096942852641411e243, 0 }, 187 { 142, 2.69536413788816277658850750804e245, 0 }, 188 { 143, 3.85437071718007277052156573649e247, 0 }, 189 { 144, 5.55029383273930478955105466055e249, 0 }, 190 { 145, 8.04792605747199194484902925780e251, 0 }, 191 { 146, 1.17499720439091082394795827164e254, 0 }, 192 { 147, 1.72724589045463891120349865931e256, 0 }, 193 { 148, 2.55632391787286558858117801578e258, 0 }, 194 { 149, 3.80892263763056972698595524351e260, 0 }, 195 { 150, 5.71338395644585459047893286526e262, 0 }, 196 { 151, 8.62720977423324043162318862650e264, 0 }, 197 { 152, 1.31133588568345254560672467123e267, 0 }, 198 { 153, 2.00634390509568239477828874699e269, 0 }, 199 { 154, 3.08976961384735088795856467036e271, 0 }, 200 { 155, 4.78914290146339387633577523906e273, 0 }, 201 { 156, 7.47106292628289444708380937294e275, 0 }, 202 { 157, 1.17295687942641442819215807155e278, 0 }, 203 { 158, 1.85327186949373479654360975305e280, 0 }, 204 { 159, 2.94670227249503832650433950735e282, 0 }, 205 { 160, 4.71472363599206132240694321176e284, 0 }, 206 { 161, 7.59070505394721872907517857094e286, 0 }, 207 { 162, 1.22969421873944943411017892849e289, 0 }, 208 { 163, 2.00440157654530257759959165344e291, 0 }, 209 { 164, 3.28721858553429622726333031164e293, 0 }, 210 { 165, 5.42391066613158877498449501421e295, 0 }, 211 { 166, 9.00369170577843736647426172359e297, 0 }, 212 { 167, 1.50361651486499904020120170784e300, 0 }, 213 { 168, 2.52607574497319838753801886917e302, 0 }, 214 { 169, 4.26906800900470527493925188890e304, 0 }, 215 { 170, 7.25741561530799896739672821113e306, 0 }, 216 217 /* 218 { 171, 1.24101807021766782342484052410e309, 0 }, 219 { 172, 2.13455108077438865629072570146e311, 0 }, 220 { 173, 3.69277336973969237538295546352e313, 0 }, 221 { 174, 6.42542566334706473316634250653e315, 0 }, 222 { 175, 1.12444949108573632830410993864e318, 0 }, 223 { 176, 1.97903110431089593781523349201e320, 0 }, 224 { 177, 3.50288505463028580993296328086e322, 0 }, 225 { 178, 6.23513539724190874168067463993e324, 0 }, 226 { 179, 1.11608923610630166476084076055e327, 0 }, 227 { 180, 2.00896062499134299656951336898e329, 0 }, 228 { 181, 3.63621873123433082379081919786e331, 0 }, 229 { 182, 6.61791809084648209929929094011e333, 0 }, 230 { 183, 1.21107901062490622417177024204e336, 0 }, 231 { 184, 2.22838537954982745247605724535e338, 0 }, 232 { 185, 4.12251295216718078708070590390e340, 0 }, 233 { 186, 7.66787409103095626397011298130e342, 0 }, 234 { 187, 1.43389245502278882136241112750e345, 0 }, 235 { 188, 2.69571781544284298416133291969e347, 0 }, 236 { 189, 5.09490667118697324006491921822e349, 0 }, 237 { 190, 9.68032267525524915612334651460e351, 0 }, 238 { 191, 1.84894163097375258881955918429e354, 0 }, 239 { 192, 3.54996793146960497053355363384e356, 0 }, 240 { 193, 6.85143810773633759312975851330e358, 0 }, 241 { 194, 1.32917899290084949306717315158e361, 0 }, 242 { 195, 2.59189903615665651148098764559e363, 0 }, 243 { 196, 5.08012211086704676250273578535e365, 0 }, 244 { 197, 1.00078405584080821221303894971e368, 0 }, 245 { 198, 1.98155243056480026018181712043e370, 0 }, 246 { 199, 3.94328933682395251776181606966e372, 0 }, 247 { 200, 7.88657867364790503552363213932e374, 0 } 248 */ 249 }; 250 251 static struct {int n; double f; long i; } doub_fact_table[GSL_SF_DOUBLEFACT_NMAX + 1] = { 252 { 0, 1.000000000000000000000000000, 1L }, 253 { 1, 1.000000000000000000000000000, 1L }, 254 { 2, 2.000000000000000000000000000, 2L }, 255 { 3, 3.000000000000000000000000000, 3L }, 256 { 4, 8.000000000000000000000000000, 8L }, 257 { 5, 15.00000000000000000000000000, 15L }, 258 { 6, 48.00000000000000000000000000, 48L }, 259 { 7, 105.0000000000000000000000000, 105L }, 260 { 8, 384.0000000000000000000000000, 384L }, 261 { 9, 945.0000000000000000000000000, 945L }, 262 { 10, 3840.000000000000000000000000, 3840L }, 263 { 11, 10395.00000000000000000000000, 10395L }, 264 { 12, 46080.00000000000000000000000, 46080L }, 265 { 13, 135135.0000000000000000000000, 135135L }, 266 { 14, 645120.00000000000000000000000, 645120L }, 267 { 15, 2.02702500000000000000000000000e6, 2027025L }, 268 { 16, 1.03219200000000000000000000000e7, 10321920L }, 269 { 17, 3.4459425000000000000000000000e7, 34459425L }, 270 { 18, 1.85794560000000000000000000000e8, 185794560L }, 271 { 19, 6.5472907500000000000000000000e8, 0 }, 272 { 20, 3.7158912000000000000000000000e9, 0 }, 273 { 21, 1.37493105750000000000000000000e10, 0 }, 274 { 22, 8.1749606400000000000000000000e10, 0 }, 275 { 23, 3.1623414322500000000000000000e11, 0 }, 276 { 24, 1.96199055360000000000000000000e12, 0 }, 277 { 25, 7.9058535806250000000000000000e12, 0 }, 278 { 26, 5.1011754393600000000000000000e13, 0 }, 279 { 27, 2.13458046676875000000000000000e14, 0 }, 280 { 28, 1.42832912302080000000000000000e15, 0 }, 281 { 29, 6.1902833536293750000000000000e15, 0 }, 282 { 30, 4.2849873690624000000000000000e16, 0 }, 283 { 31, 1.91898783962510625000000000000e17, 0 }, 284 { 32, 1.37119595809996800000000000000e18, 0 }, 285 { 33, 6.3326598707628506250000000000e18, 0 }, 286 { 34, 4.6620662575398912000000000000e19, 0 }, 287 { 35, 2.21643095476699771875000000000e20, 0 }, 288 { 36, 1.67834385271436083200000000000e21, 0 }, 289 { 37, 8.2007945326378915593750000000e21, 0 }, 290 { 38, 6.3777066403145711616000000000e22, 0 }, 291 { 39, 3.1983098677287777081562500000e23, 0 }, 292 { 40, 2.55108265612582846464000000000e24, 0 }, 293 { 41, 1.31130704576879886034406250000e25, 0 }, 294 { 42, 1.07145471557284795514880000000e26, 0 }, 295 { 43, 5.6386202968058350994794687500e26, 0 }, 296 { 44, 4.7144007485205310026547200000e27, 0 }, 297 { 45, 2.53737913356262579476576093750e28, 0 }, 298 { 46, 2.16862434431944426122117120000e29, 0 }, 299 { 47, 1.19256819277443412353990764062e30, 0 }, 300 { 48, 1.04093968527333324538616217600e31, 0 }, 301 { 49, 5.8435841445947272053455474391e31, 0 }, 302 { 50, 5.2046984263666662269308108800e32, 0 }, 303 { 51, 2.98022791374331087472622919392e33, 0 }, 304 { 52, 2.70644318171066643800402165760e34, 0 }, 305 { 53, 1.57952079428395476360490147278e35, 0 }, 306 { 54, 1.46147931812375987652217169510e36, 0 }, 307 { 55, 8.6873643685617511998269581003e36, 0 }, 308 { 56, 8.1842841814930553085241614926e37, 0 }, 309 { 57, 4.9517976900801981839013661172e38, 0 }, 310 { 58, 4.7468848252659720789440136657e39, 0 }, 311 { 59, 2.92156063714731692850180600912e40, 0 }, 312 { 60, 2.84813089515958324736640819942e41, 0 }, 313 { 61, 1.78215198865986332638610166557e42, 0 }, 314 { 62, 1.76584115499894161336717308364e43, 0 }, 315 { 63, 1.12275575285571389562324404931e44, 0 }, 316 { 64, 1.13013833919932263255499077353e45, 0 }, 317 { 65, 7.2979123935621403215510863205e45, 0 }, 318 { 66, 7.4589130387155293748629391053e46, 0 }, 319 { 67, 4.8896013036866340154392278347e47, 0 }, 320 { 68, 5.0720608663265599749067985916e48, 0 }, 321 { 69, 3.3738248995437774706530672060e49, 0 }, 322 { 70, 3.5504426064285919824347590141e50, 0 }, 323 { 71, 2.39541567867608200416367771623e51, 0 }, 324 { 72, 2.55631867662858622735302649017e52, 0 }, 325 { 73, 1.74865344543353986303948473285e53, 0 }, 326 { 74, 1.89167582070515380824123960272e54, 0 }, 327 { 75, 1.31149008407515489727961354964e55, 0 }, 328 { 76, 1.43767362373591689426334209807e56, 0 }, 329 { 77, 1.00984736473786927090530243322e57, 0 }, 330 { 78, 1.12138542651401517752540683649e58, 0 }, 331 { 79, 7.9777941814291672401518892225e58, 0 }, 332 { 80, 8.9710834121121214202032546920e59, 0 }, 333 { 81, 6.4620132869576254645230302702e60, 0 }, 334 { 82, 7.3562883979319395645666688474e61, 0 }, 335 { 83, 5.3634710281748291355541151243e62, 0 }, 336 { 84, 6.1792822542628292342360018318e63, 0 }, 337 { 85, 4.5589503739486047652209978556e64, 0 }, 338 { 86, 5.3141827386660331414429615754e65, 0 }, 339 { 87, 3.9662868253352861457422681344e66, 0 }, 340 { 88, 4.6764808100261091644698061863e67, 0 }, 341 { 89, 3.5299952745484046697106186396e68, 0 }, 342 { 90, 4.2088327290234982480228255677e69, 0 }, 343 { 91, 3.2122956998390482494366629620e70, 0 }, 344 { 92, 3.8721261107016183881809995223e71, 0 }, 345 { 93, 2.98743500085031487197609655470e72, 0 }, 346 { 94, 3.6397985440595212848901395509e73, 0 }, 347 { 95, 2.83806325080779912837729172696e74, 0 }, 348 { 96, 3.4942066022971404334945339689e75, 0 }, 349 { 97, 2.75292135328356515452597297515e76, 0 }, 350 { 98, 3.4243224702511976248246432895e77, 0 }, 351 { 99, 2.72539213975072950298071324540e78, 0 }, 352 { 100, 3.4243224702511976248246432895e79, 0 }, 353 { 101, 2.75264606114823679801052037785e80, 0 }, 354 { 102, 3.4928089196562215773211361553e81, 0 }, 355 { 103, 2.83522544298268390195083598919e82, 0 }, 356 { 104, 3.6325212764424704404139816015e83, 0 }, 357 { 105, 2.97698671513181809704837778865e84, 0 }, 358 { 106, 3.8504725530290186668388204976e85, 0 }, 359 { 107, 3.1853757851910453638417642339e86, 0 }, 360 { 108, 4.1585103572713401601859261374e87, 0 }, 361 { 109, 3.4720596058582394465875230149e88, 0 }, 362 { 110, 4.5743613929984741762045187512e89, 0 }, 363 { 111, 3.8539861625026457857121505465e90, 0 }, 364 { 112, 5.1232847601582910773490610013e91, 0 }, 365 { 113, 4.3550043636279897378547301176e92, 0 }, 366 { 114, 5.8405446265804518281779295415e93, 0 }, 367 { 115, 5.0082550181721881985329396352e94, 0 }, 368 { 116, 6.7750317668333241206863982681e95, 0 }, 369 { 117, 5.8596583712614601922835393732e96, 0 }, 370 { 118, 7.9945374848633224624099499564e97, 0 }, 371 { 119, 6.9729934618011376288174118541e98, 0 }, 372 { 120, 9.5934449818359869548919399477e99, 0 }, 373 { 121, 8.4373220887793765308690683435e100, 0 }, 374 { 122, 1.17040028778399040849681667362e102, 0 }, 375 { 123, 1.03779061691986331329689540625e103, 0 }, 376 { 124, 1.45129635685214810653605267528e104, 0 }, 377 { 125, 1.29723827114982914162111925781e105, 0 }, 378 { 126, 1.82863340963370661423542637086e106, 0 }, 379 { 127, 1.64749260436028300985882145742e107, 0 }, 380 { 128, 2.34065076433114446622134575470e108, 0 }, 381 { 129, 2.12526545962476508271787968008e109, 0 }, 382 { 130, 3.04284599363048780608774948111e110, 0 }, 383 { 131, 2.78409775210844225836042238090e111, 0 }, 384 { 132, 4.0165567115922439040358293151e112, 0 }, 385 { 133, 3.7028500103042282036193617666e113, 0 }, 386 { 134, 5.3821859935336068314080112822e114, 0 }, 387 { 135, 4.9988475139107080748861383849e115, 0 }, 388 { 136, 7.3197729512057052907148953438e116, 0 }, 389 { 137, 6.8484210940576700625940095873e117, 0 }, 390 { 138, 1.01012866726638733011865555744e119, 0 }, 391 { 139, 9.5193053207401613870056733264e119, 0 }, 392 { 140, 1.41418013417294226216611778042e121, 0 }, 393 { 141, 1.34222205022436275556779993902e122, 0 }, 394 { 142, 2.00813579052557801227588724819e123, 0 }, 395 { 143, 1.91937753182083874046195391280e124, 0 }, 396 { 144, 2.89171553835683233767727763739e125, 0 }, 397 { 145, 2.78309742114021617366983317355e126, 0 }, 398 { 146, 4.2219046860009752130088253506e127, 0 }, 399 { 147, 4.0911532090761177752946547651e128, 0 }, 400 { 148, 6.2484189352814433152530615189e129, 0 }, 401 { 149, 6.0958182815234154851890356000e130, 0 }, 402 { 150, 9.3726284029221649728795922783e131, 0 }, 403 { 151, 9.2046856051003573826354437561e132, 0 }, 404 { 152, 1.42463951724416907587769802630e134, 0 }, 405 { 153, 1.40831689758035467954322289468e135, 0 }, 406 { 154, 2.19394485655602037685165496051e136, 0 }, 407 { 155, 2.18289119124954975329199548675e137, 0 }, 408 { 156, 3.4225539762273917878885817384e138, 0 }, 409 { 157, 3.4271391702617931126684329142e139, 0 }, 410 { 158, 5.4076352824392790248639591467e140, 0 }, 411 { 159, 5.4491512807162510491428083336e141, 0 }, 412 { 160, 8.6522164519028464397823346347e142, 0 }, 413 { 161, 8.7731335619531641891199214170e143, 0 }, 414 { 162, 1.40165906520826112324473821082e145, 0 }, 415 { 163, 1.43002077059836576282654719098e146, 0 }, 416 { 164, 2.29872086694154824212137066574e147, 0 }, 417 { 165, 2.35953427148730350866380286512e148, 0 }, 418 { 166, 3.8158766391229700819214753051e149, 0 }, 419 { 167, 3.9404222333837968594685507847e150, 0 }, 420 { 168, 6.4106727537265897376280785126e151, 0 }, 421 { 169, 6.6593135744186166925018508262e152, 0 }, 422 { 170, 1.08981436813352025539677334714e154, 0 }, 423 { 171, 1.13874262122558345441781649128e155, 0 }, 424 { 172, 1.87448071318965483928245015709e156, 0 }, 425 { 173, 1.97002473472025937614282252992e157, 0 }, 426 { 174, 3.2615964409499994203514632733e158, 0 }, 427 { 175, 3.4475432857604539082499394274e159, 0 }, 428 { 176, 5.7404097360719989798185753611e160, 0 }, 429 { 177, 6.1021516157960034176023927864e161, 0 }, 430 { 178, 1.02179293302081581840770641427e163, 0 }, 431 { 179, 1.09228513922748461175082830877e164, 0 }, 432 { 180, 1.83922727943746847313387154568e165, 0 }, 433 { 181, 1.97703610200174714726899923887e166, 0 }, 434 { 182, 3.3473936485761926211036462131e167, 0 }, 435 { 183, 3.6179760666631972795022686071e168, 0 }, 436 { 184, 6.1592043133801944228307090322e169, 0 }, 437 { 185, 6.6932557233269149670791969232e170, 0 }, 438 { 186, 1.14561200228871616264651187999e172, 0 }, 439 { 187, 1.25163882026213309884380982464e173, 0 }, 440 { 188, 2.15375056430278638577544233437e174, 0 }, 441 { 189, 2.36559737029543155681480056857e175, 0 }, 442 { 190, 4.0921260721752941329733404353e176, 0 }, 443 { 191, 4.5182909772642742735162690860e177, 0 }, 444 { 192, 7.8568820585765647353088136358e178, 0 }, 445 { 193, 8.7203015861200493478863993359e179, 0 }, 446 { 194, 1.52423511936385355864990984535e181, 0 }, 447 { 195, 1.70045880929340962283784787050e182, 0 }, 448 { 196, 2.98750083395315297495382329688e183, 0 }, 449 { 197, 3.3499038543080169569905603049e184, 0 }, 450 { 198, 5.9152516512272428904085701278e185, 0 }, 451 { 199, 6.6663086700729537444112150067e186, 0 }, 452 { 200, 1.18305033024544857808171402556e188, 0 }, 453 { 201, 1.33992804268466370262665421635e189, 0 }, 454 { 202, 2.38976166709580612772506233164e190, 0 }, 455 { 203, 2.72005392664986731633210805920e191, 0 }, 456 { 204, 4.8751138008754445005591271565e192, 0 }, 457 { 205, 5.5761105496322279984808215214e193, 0 }, 458 { 206, 1.00427344298034156711518019425e195, 0 }, 459 { 207, 1.15425488377387119568553005492e196, 0 }, 460 { 208, 2.08888876139911045959957480403e197, 0 }, 461 { 209, 2.41239270708739079898275781478e198, 0 }, 462 { 210, 4.3866663989381319651591070885e199, 0 }, 463 { 211, 5.0901486119543945858536189892e200, 0 }, 464 { 212, 9.2997327657488397661373070276e201, 0 }, 465 { 213, 1.08420165434628604678682084470e203, 0 }, 466 { 214, 1.99014281187025170995338370390e204, 0 }, 467 { 215, 2.33103355684451500059166481610e205, 0 }, 468 { 216, 4.2987084736397436934993088004e206, 0 }, 469 { 217, 5.0583428183525975512839126509e207, 0 }, 470 { 218, 9.3711844725346412518284931849e208, 0 }, 471 { 219, 1.10777707721921886373117687056e210, 0 }, 472 { 220, 2.06166058395762107540226850068e211, 0 }, 473 { 221, 2.44818734065447368884590088393e212, 0 }, 474 { 222, 4.5768864963859187873930360715e213, 0 }, 475 { 223, 5.4594577696594763261263589712e214, 0 }, 476 { 224, 1.02522257519044580837604008002e216, 0 }, 477 { 225, 1.22837799817338217337843076851e217, 0 }, 478 { 226, 2.31700301993040752692985058084e218, 0 }, 479 { 227, 2.78841805585357753356903784452e219, 0 }, 480 { 228, 5.2827668854413291614000593243e220, 0 }, 481 { 229, 6.3854773479046925518730966640e221, 0 }, 482 { 230, 1.21503638365150570712201364459e223, 0 }, 483 { 231, 1.47504526736598397948268532937e224, 0 }, 484 { 232, 2.81888441007149324052307165546e225, 0 }, 485 { 233, 3.4368554729627426721946568174e226, 0 }, 486 { 234, 6.5961895195672941828239876738e227, 0 }, 487 { 235, 8.0766103614624452796574435210e228, 0 }, 488 { 236, 1.55670072661788142714646109101e230, 0 }, 489 { 237, 1.91415665566659953127881411447e231, 0 }, 490 { 238, 3.7049477293505577966085773966e232, 0 }, 491 { 239, 4.5748344070431728797563657336e233, 0 }, 492 { 240, 8.8918745504413387118605857518e234, 0 }, 493 { 241, 1.10253509209740466402128414180e236, 0 }, 494 { 242, 2.15183364120680396827026175195e237, 0 }, 495 { 243, 2.67916027379669333357172046456e238, 0 }, 496 { 244, 5.2504740845446016825794386748e239, 0 }, 497 { 245, 6.5639426708018986672507151382e240, 0 }, 498 { 246, 1.29161662479797201391454191399e242, 0 }, 499 { 247, 1.62129383968806897081092663913e243, 0 }, 500 { 248, 3.2032092294989705945080639467e244, 0 }, 501 { 249, 4.0370216608232917373192073314e245, 0 }, 502 { 250, 8.0080230737474264862701598667e246, 0 }, 503 { 251, 1.01329243686664622606712104019e248, 0 }, 504 { 252, 2.01802181458435147454008028642e249, 0 }, 505 { 253, 2.56362986527261495194981623168e250, 0 }, 506 { 254, 5.1257754090442527453318039275e251, 0 }, 507 { 255, 6.5372561564451681274720313908e252, 0 }, 508 { 256, 1.31219850471532870280494180544e254, 0 }, 509 { 257, 1.68007483220640820876031206743e255, 0 }, 510 { 258, 3.3854721421655480532367498580e256, 0 }, 511 { 259, 4.3513938154145972606892082546e257, 0 }, 512 { 260, 8.8022275696304249384155496309e258, 0 }, 513 { 261, 1.13571378582320988503988335446e260, 0 }, 514 { 262, 2.30618362324317133386487400329e261, 0 }, 515 { 263, 2.98692725671504199765489322224e262, 0 }, 516 { 264, 6.0883247653619723214032673687e263, 0 }, 517 { 265, 7.9153572302948612937854670389e264, 0 }, 518 { 266, 1.61949438758628463749326912007e266, 0 }, 519 { 267, 2.11340038048872796544071969939e267, 0 }, 520 { 268, 4.3402449587312428284819612418e268, 0 }, 521 { 269, 5.6850470235146782270355359914e269, 0 }, 522 { 270, 1.17186613885743556369012953528e271, 0 }, 523 { 271, 1.54064774337247779952663025366e272, 0 }, 524 { 272, 3.1874758976922247332371523360e273, 0 }, 525 { 273, 4.2059683394068643927077005925e274, 0 }, 526 { 274, 8.7336839596766957690697974006e275, 0 }, 527 { 275, 1.15664129333688770799461766294e277, 0 }, 528 { 276, 2.41049677287076803226326408256e278, 0 }, 529 { 277, 3.2038963825431789511450909263e279, 0 }, 530 { 278, 6.7011810285807351296918741495e280, 0 }, 531 { 279, 8.9388709072954692736948036845e281, 0 }, 532 { 280, 1.87633068800260583631372476186e283, 0 }, 533 { 281, 2.51182272495002686590823983534e284, 0 }, 534 { 282, 5.2912525401673484584047038284e285, 0 }, 535 { 283, 7.1084583116085760305203187340e286, 0 }, 536 { 284, 1.50271572140752696218693588728e288, 0 }, 537 { 285, 2.02591061880844416869829083919e289, 0 }, 538 { 286, 4.2977669632255271118546366376e290, 0 }, 539 { 287, 5.8143634759802347641640947085e291, 0 }, 540 { 288, 1.23775688540895180821413535163e293, 0 }, 541 { 289, 1.68035104455828784684342337075e294, 0 }, 542 { 290, 3.5894949676859602438209925197e295, 0 }, 543 { 291, 4.8898215396646176343143620089e296, 0 }, 544 { 292, 1.04813253056430039119572981576e298, 0 }, 545 { 293, 1.43271771112173296685410806860e299, 0 }, 546 { 294, 3.08150963985904315011544565835e300, 0 }, 547 { 295, 4.2265172478091122522196188024e301, 0 }, 548 { 296, 9.1212685339827677243417191487e302, 0 }, 549 { 297, 1.25527562259930633890922678431e304, 0 }, 550 /* 551 { 298, 2.71813802312686478185383230631e305, 0 }, 552 { 299, 3.7532741115719259533385880851e306, 0 }, 553 { 300, 8.1544140693805943455614969189e307, } 554 */ 555 }; 556 557 558 /* Chebyshev coefficients for Gamma*(3/4(t+1)+1/2), -1<t<1 559 */ 560 static double gstar_a_data[30] = { 561 2.16786447866463034423060819465, 562 -0.05533249018745584258035832802, 563 0.01800392431460719960888319748, 564 -0.00580919269468937714480019814, 565 0.00186523689488400339978881560, 566 -0.00059746524113955531852595159, 567 0.00019125169907783353925426722, 568 -0.00006124996546944685735909697, 569 0.00001963889633130842586440945, 570 -6.3067741254637180272515795142e-06, 571 2.0288698405861392526872789863e-06, 572 -6.5384896660838465981983750582e-07, 573 2.1108698058908865476480734911e-07, 574 -6.8260714912274941677892994580e-08, 575 2.2108560875880560555583978510e-08, 576 -7.1710331930255456643627187187e-09, 577 2.3290892983985406754602564745e-09, 578 -7.5740371598505586754890405359e-10, 579 2.4658267222594334398525312084e-10, 580 -8.0362243171659883803428749516e-11, 581 2.6215616826341594653521346229e-11, 582 -8.5596155025948750540420068109e-12, 583 2.7970831499487963614315315444e-12, 584 -9.1471771211886202805502562414e-13, 585 2.9934720198063397094916415927e-13, 586 -9.8026575909753445931073620469e-14, 587 3.2116773667767153777571410671e-14, 588 -1.0518035333878147029650507254e-14, 589 3.4144405720185253938994854173e-15, 590 -1.0115153943081187052322643819e-15 591 }; 592 static cheb_series gstar_a_cs = { 593 gstar_a_data, 594 29, 595 -1, 1, 596 17 597 }; 598 599 600 /* Chebyshev coefficients for 601 * x^2(Gamma*(x) - 1 - 1/(12x)), x = 4(t+1)+2, -1 < t < 1 602 */ 603 static double gstar_b_data[] = { 604 0.0057502277273114339831606096782, 605 0.0004496689534965685038254147807, 606 -0.0001672763153188717308905047405, 607 0.0000615137014913154794776670946, 608 -0.0000223726551711525016380862195, 609 8.0507405356647954540694800545e-06, 610 -2.8671077107583395569766746448e-06, 611 1.0106727053742747568362254106e-06, 612 -3.5265558477595061262310873482e-07, 613 1.2179216046419401193247254591e-07, 614 -4.1619640180795366971160162267e-08, 615 1.4066283500795206892487241294e-08, 616 -4.6982570380537099016106141654e-09, 617 1.5491248664620612686423108936e-09, 618 -5.0340936319394885789686867772e-10, 619 1.6084448673736032249959475006e-10, 620 -5.0349733196835456497619787559e-11, 621 1.5357154939762136997591808461e-11, 622 -4.5233809655775649997667176224e-12, 623 1.2664429179254447281068538964e-12, 624 -3.2648287937449326771785041692e-13, 625 7.1528272726086133795579071407e-14, 626 -9.4831735252566034505739531258e-15, 627 -2.3124001991413207293120906691e-15, 628 2.8406613277170391482590129474e-15, 629 -1.7245370321618816421281770927e-15, 630 8.6507923128671112154695006592e-16, 631 -3.9506563665427555895391869919e-16, 632 1.6779342132074761078792361165e-16, 633 -6.0483153034414765129837716260e-17 634 }; 635 static cheb_series gstar_b_cs = { 636 gstar_b_data, 637 29, 638 -1, 1, 639 18 640 }; 641 642 643 /* coefficients for gamma=7, kmax=8 Lanczos method */ 644 static double lanczos_7_c[9] = { 645 0.99999999999980993227684700473478, 646 676.520368121885098567009190444019, 647 -1259.13921672240287047156078755283, 648 771.3234287776530788486528258894, 649 -176.61502916214059906584551354, 650 12.507343278686904814458936853, 651 -0.13857109526572011689554707, 652 9.984369578019570859563e-6, 653 1.50563273514931155834e-7 654 }; 655 656 /* complex version of Lanczos method; this is not safe for export 657 * since it becomes bad in the left half-plane 658 */ 659 static 660 int 661 lngamma_lanczos_complex(double zr, double zi, gsl_sf_result * yr, gsl_sf_result * yi) 662 { 663 int k; 664 gsl_sf_result log1_r, log1_i; 665 gsl_sf_result logAg_r, logAg_i; 666 double Ag_r, Ag_i; 667 double yi_tmp_val, yi_tmp_err; 668 669 zr -= 1.0; /* Lanczos writes z! instead of Gamma(z) */ 670 671 Ag_r = lanczos_7_c[0]; 672 Ag_i = 0.0; 673 for(k=1; k<=8; k++) { 674 double R = zr + k; 675 double I = zi; 676 double a = lanczos_7_c[k] / (R*R + I*I); 677 Ag_r += a * R; 678 Ag_i -= a * I; 679 } 680 681 gsl_sf_complex_log_e(zr + 7.5, zi, &log1_r, &log1_i); 682 gsl_sf_complex_log_e(Ag_r, Ag_i, &logAg_r, &logAg_i); 683 684 /* (z+0.5)*log(z+7.5) - (z+7.5) + LogRootTwoPi_ + log(Ag(z)) */ 685 yr->val = (zr+0.5)*log1_r.val - zi*log1_i.val - (zr+7.5) + LogRootTwoPi_ + logAg_r.val; 686 yi->val = zi*log1_r.val + (zr+0.5)*log1_i.val - zi + logAg_i.val; 687 yr->err = 4.0 * GSL_DBL_EPSILON * fabs(yr->val); 688 yi->err = 4.0 * GSL_DBL_EPSILON * fabs(yi->val); 689 yi_tmp_val = yi->val; 690 yi_tmp_err = yi->err; 691 gsl_sf_angle_restrict_symm_err_e(yi_tmp_val, yi); 692 yi->err += yi_tmp_err; 693 return GSL_SUCCESS; 694 } 695 696 697 /* Lanczos method for real x > 0; 698 * gamma=7, truncated at 1/(z+8) 699 * [J. SIAM Numer. Anal, Ser. B, 1 (1964) 86] 700 */ 701 static 702 int 703 lngamma_lanczos(double x, gsl_sf_result * result) 704 { 705 int k; 706 double Ag; 707 double term1, term2; 708 709 x -= 1.0; /* Lanczos writes z! instead of Gamma(z) */ 710 711 Ag = lanczos_7_c[0]; 712 for(k=1; k<=8; k++) { Ag += lanczos_7_c[k]/(x+k); } 713 714 /* (x+0.5)*log(x+7.5) - (x+7.5) + LogRootTwoPi_ + log(Ag(x)) */ 715 term1 = (x+0.5)*log((x+7.5)/M_E); 716 term2 = LogRootTwoPi_ + log(Ag); 717 result->val = term1 + (term2 - 7.0); 718 result->err = 2.0 * GSL_DBL_EPSILON * (fabs(term1) + fabs(term2) + 7.0); 719 result->err += GSL_DBL_EPSILON * fabs(result->val); 720 721 return GSL_SUCCESS; 722 } 723 724 /* x = eps near zero 725 * gives double-precision for |eps| < 0.02 726 */ 727 static 728 int 729 lngamma_sgn_0(double eps, gsl_sf_result * lng, double * sgn) 730 { 731 /* calculate series for g(eps) = Gamma(eps) eps - 1/(1+eps) - eps/2 */ 732 const double c1 = -0.07721566490153286061; 733 const double c2 = -0.01094400467202744461; 734 const double c3 = 0.09252092391911371098; 735 const double c4 = -0.01827191316559981266; 736 const double c5 = 0.01800493109685479790; 737 const double c6 = -0.00685088537872380685; 738 const double c7 = 0.00399823955756846603; 739 const double c8 = -0.00189430621687107802; 740 const double c9 = 0.00097473237804513221; 741 const double c10 = -0.00048434392722255893; 742 const double g6 = c6+eps*(c7+eps*(c8 + eps*(c9 + eps*c10))); 743 const double g = eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*g6))))); 744 745 /* calculate Gamma(eps) eps, a positive quantity */ 746 const double gee = g + 1.0/(1.0+eps) + 0.5*eps; 747 748 lng->val = log(gee/fabs(eps)); 749 lng->err = 4.0 * GSL_DBL_EPSILON * fabs(lng->val); 750 *sgn = GSL_SIGN(eps); 751 752 return GSL_SUCCESS; 753 } 754 755 756 /* x near a negative integer 757 * Calculates sign as well as log(|gamma(x)|). 758 * x = -N + eps 759 * assumes N >= 1 760 */ 761 static 762 int 763 lngamma_sgn_sing(int N, double eps, gsl_sf_result * lng, double * sgn) 764 { 765 if(eps == 0.0) { 766 lng->val = 0.0; 767 lng->err = 0.0; 768 *sgn = 0.0; 769 GSL_ERROR ("error", GSL_EDOM); 770 } 771 else if(N == 1) { 772 /* calculate series for 773 * g = eps gamma(-1+eps) + 1 + eps/2 (1+3eps)/(1-eps^2) 774 * double-precision for |eps| < 0.02 775 */ 776 const double c0 = 0.07721566490153286061; 777 const double c1 = 0.08815966957356030521; 778 const double c2 = -0.00436125434555340577; 779 const double c3 = 0.01391065882004640689; 780 const double c4 = -0.00409427227680839100; 781 const double c5 = 0.00275661310191541584; 782 const double c6 = -0.00124162645565305019; 783 const double c7 = 0.00065267976121802783; 784 const double c8 = -0.00032205261682710437; 785 const double c9 = 0.00016229131039545456; 786 const double g5 = c5 + eps*(c6 + eps*(c7 + eps*(c8 + eps*c9))); 787 const double g = eps*(c0 + eps*(c1 + eps*(c2 + eps*(c3 + eps*(c4 + eps*g5))))); 788 789 /* calculate eps gamma(-1+eps), a negative quantity */ 790 const double gam_e = g - 1.0 - 0.5*eps*(1.0+3.0*eps)/(1.0 - eps*eps); 791 792 lng->val = log(fabs(gam_e)/fabs(eps)); 793 lng->err = 2.0 * GSL_DBL_EPSILON * fabs(lng->val); 794 *sgn = ( eps > 0.0 ? -1.0 : 1.0 ); 795 return GSL_SUCCESS; 796 } 797 else { 798 double g; 799 800 /* series for sin(Pi(N+1-eps))/(Pi eps) modulo the sign 801 * double-precision for |eps| < 0.02 802 */ 803 const double cs1 = -1.6449340668482264365; 804 const double cs2 = 0.8117424252833536436; 805 const double cs3 = -0.1907518241220842137; 806 const double cs4 = 0.0261478478176548005; 807 const double cs5 = -0.0023460810354558236; 808 const double e2 = eps*eps; 809 const double sin_ser = 1.0 + e2*(cs1+e2*(cs2+e2*(cs3+e2*(cs4+e2*cs5)))); 810 811 /* calculate series for ln(gamma(1+N-eps)) 812 * double-precision for |eps| < 0.02 813 */ 814 double aeps = fabs(eps); 815 double c1, c2, c3, c4, c5, c6, c7; 816 double lng_ser; 817 gsl_sf_result c0; 818 gsl_sf_result psi_0; 819 gsl_sf_result psi_1; 820 gsl_sf_result psi_2; 821 gsl_sf_result psi_3; 822 gsl_sf_result psi_4; 823 gsl_sf_result psi_5; 824 gsl_sf_result psi_6; 825 psi_2.val = 0.0; 826 psi_3.val = 0.0; 827 psi_4.val = 0.0; 828 psi_5.val = 0.0; 829 psi_6.val = 0.0; 830 gsl_sf_lnfact_e(N, &c0); 831 gsl_sf_psi_int_e(N+1, &psi_0); 832 gsl_sf_psi_1_int_e(N+1, &psi_1); 833 if(aeps > 0.00001) gsl_sf_psi_n_e(2, N+1.0, &psi_2); 834 if(aeps > 0.0002) gsl_sf_psi_n_e(3, N+1.0, &psi_3); 835 if(aeps > 0.001) gsl_sf_psi_n_e(4, N+1.0, &psi_4); 836 if(aeps > 0.005) gsl_sf_psi_n_e(5, N+1.0, &psi_5); 837 if(aeps > 0.01) gsl_sf_psi_n_e(6, N+1.0, &psi_6); 838 c1 = psi_0.val; 839 c2 = psi_1.val/2.0; 840 c3 = psi_2.val/6.0; 841 c4 = psi_3.val/24.0; 842 c5 = psi_4.val/120.0; 843 c6 = psi_5.val/720.0; 844 c7 = psi_6.val/5040.0; 845 lng_ser = c0.val-eps*(c1-eps*(c2-eps*(c3-eps*(c4-eps*(c5-eps*(c6-eps*c7)))))); 846 847 /* calculate 848 * g = ln(|eps gamma(-N+eps)|) 849 * = -ln(gamma(1+N-eps)) + ln(|eps Pi/sin(Pi(N+1+eps))|) 850 */ 851 g = -lng_ser - log(sin_ser); 852 853 lng->val = g - log(fabs(eps)); 854 lng->err = c0.err + 2.0 * GSL_DBL_EPSILON * (fabs(g) + fabs(lng->val)); 855 856 *sgn = ( GSL_IS_ODD(N) ? -1.0 : 1.0 ) * ( eps > 0.0 ? 1.0 : -1.0 ); 857 858 return GSL_SUCCESS; 859 } 860 } 861 862 863 /* This gets bad near the negative half axis. However, this 864 * region can be avoided by use of the reflection formula, as usual. 865 * Only the first two terms of the series are kept. 866 */ 867 #if 0 868 static 869 int 870 lngamma_complex_stirling(const double zr, const double zi, double * lg_r, double * arg) 871 { 872 double re_zinv, im_zinv; 873 double re_zinv2, im_zinv2; 874 double re_zinv3, im_zinv3; 875 double re_zhlnz, im_zhlnz; 876 double r, lnr, theta; 877 gsl_sf_complex_log_e(zr, zi, &lnr, &theta); /* z = r e^{i theta} */ 878 r = exp(lnr); 879 re_zinv = (zr/r)/r; 880 im_zinv = -(zi/r)/r; 881 re_zinv2 = re_zinv*re_zinv - im_zinv*im_zinv; 882 re_zinv2 = 2.0*re_zinv*im_zinv; 883 re_zinv3 = re_zinv2*re_zinv - im_zinv2*im_zinv; 884 re_zinv3 = re_zinv2*im_zinv + im_zinv2*re_zinv; 885 re_zhlnz = (zr - 0.5)*lnr - zi*theta; 886 im_zhlnz = zi*lnr + zr*theta; 887 *lg_r = re_zhlnz - zr + 0.5*(M_LN2+M_LNPI) + re_zinv/12.0 - re_zinv3/360.0; 888 *arg = im_zhlnz - zi + 1.0/12.0*im_zinv - im_zinv3/360.0; 889 return GSL_SUCCESS; 890 } 891 #endif /* 0 */ 892 893 894 inline 895 static 896 int 897 lngamma_1_pade(const double eps, gsl_sf_result * result) 898 { 899 /* Use (2,2) Pade for Log[Gamma[1+eps]]/eps 900 * plus a correction series. 901 */ 902 const double n1 = -1.0017419282349508699871138440; 903 const double n2 = 1.7364839209922879823280541733; 904 const double d1 = 1.2433006018858751556055436011; 905 const double d2 = 5.0456274100274010152489597514; 906 const double num = (eps + n1) * (eps + n2); 907 const double den = (eps + d1) * (eps + d2); 908 const double pade = 2.0816265188662692474880210318 * num / den; 909 const double c0 = 0.004785324257581753; 910 const double c1 = -0.01192457083645441; 911 const double c2 = 0.01931961413960498; 912 const double c3 = -0.02594027398725020; 913 const double c4 = 0.03141928755021455; 914 const double eps5 = eps*eps*eps*eps*eps; 915 const double corr = eps5 * (c0 + eps*(c1 + eps*(c2 + eps*(c3 + c4*eps)))); 916 result->val = eps * (pade + corr); 917 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 918 return GSL_SUCCESS; 919 } 920 921 inline 922 static 923 int 924 lngamma_2_pade(const double eps, gsl_sf_result * result) 925 { 926 /* Use (2,2) Pade for Log[Gamma[2+eps]]/eps 927 * plus a correction series. 928 */ 929 const double n1 = 1.000895834786669227164446568; 930 const double n2 = 4.209376735287755081642901277; 931 const double d1 = 2.618851904903217274682578255; 932 const double d2 = 10.85766559900983515322922936; 933 const double num = (eps + n1) * (eps + n2); 934 const double den = (eps + d1) * (eps + d2); 935 const double pade = 2.85337998765781918463568869 * num/den; 936 const double c0 = 0.0001139406357036744; 937 const double c1 = -0.0001365435269792533; 938 const double c2 = 0.0001067287169183665; 939 const double c3 = -0.0000693271800931282; 940 const double c4 = 0.0000407220927867950; 941 const double eps5 = eps*eps*eps*eps*eps; 942 const double corr = eps5 * (c0 + eps*(c1 + eps*(c2 + eps*(c3 + c4*eps)))); 943 result->val = eps * (pade + corr); 944 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 945 return GSL_SUCCESS; 946 } 947 948 949 /* series for gammastar(x) 950 * double-precision for x > 10.0 951 */ 952 static 953 int 954 gammastar_ser(const double x, gsl_sf_result * result) 955 { 956 /* Use the Stirling series for the correction to Log(Gamma(x)), 957 * which is better behaved and easier to compute than the 958 * regular Stirling series for Gamma(x). 959 */ 960 const double y = 1.0/(x*x); 961 const double c0 = 1.0/12.0; 962 const double c1 = -1.0/360.0; 963 const double c2 = 1.0/1260.0; 964 const double c3 = -1.0/1680.0; 965 const double c4 = 1.0/1188.0; 966 const double c5 = -691.0/360360.0; 967 const double c6 = 1.0/156.0; 968 const double c7 = -3617.0/122400.0; 969 const double ser = c0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*(c5 + y*(c6 + y*c7)))))); 970 result->val = exp(ser/x); 971 result->err = 2.0 * GSL_DBL_EPSILON * result->val * GSL_MAX_DBL(1.0, ser/x); 972 return GSL_SUCCESS; 973 } 974 975 976 /* Chebyshev expansion for log(gamma(x)/gamma(8)) 977 * 5 < x < 10 978 * -1 < t < 1 979 */ 980 static double gamma_5_10_data[24] = { 981 -1.5285594096661578881275075214, 982 4.8259152300595906319768555035, 983 0.2277712320977614992970601978, 984 -0.0138867665685617873604917300, 985 0.0012704876495201082588139723, 986 -0.0001393841240254993658962470, 987 0.0000169709242992322702260663, 988 -2.2108528820210580075775889168e-06, 989 3.0196602854202309805163918716e-07, 990 -4.2705675000079118380587357358e-08, 991 6.2026423818051402794663551945e-09, 992 -9.1993973208880910416311405656e-10, 993 1.3875551258028145778301211638e-10, 994 -2.1218861491906788718519522978e-11, 995 3.2821736040381439555133562600e-12, 996 -5.1260001009953791220611135264e-13, 997 8.0713532554874636696982146610e-14, 998 -1.2798522376569209083811628061e-14, 999 2.0417711600852502310258808643e-15, 1000 -3.2745239502992355776882614137e-16, 1001 5.2759418422036579482120897453e-17, 1002 -8.5354147151695233960425725513e-18, 1003 1.3858639703888078291599886143e-18, 1004 -2.2574398807738626571560124396e-19 1005 }; 1006 static const cheb_series gamma_5_10_cs = { 1007 gamma_5_10_data, 1008 23, 1009 -1, 1, 1010 11 1011 }; 1012 1013 1014 /* gamma(x) for x >= 1/2 1015 * assumes x >= 1/2 1016 */ 1017 static 1018 int 1019 gamma_xgthalf(const double x, gsl_sf_result * result) 1020 { 1021 /* CHECK_POINTER(result) */ 1022 1023 if(x == 0.5) { 1024 result->val = 1.77245385090551602729817; 1025 result->err = GSL_DBL_EPSILON * result->val; 1026 return GSL_SUCCESS; 1027 } else if (x <= (GSL_SF_FACT_NMAX + 1.0) && x == floor(x)) { 1028 int n = (int) floor (x); 1029 result->val = fact_table[n - 1].f; 1030 result->err = GSL_DBL_EPSILON * result->val; 1031 return GSL_SUCCESS; 1032 } 1033 else if(fabs(x - 1.0) < 0.01) { 1034 /* Use series for Gamma[1+eps] - 1/(1+eps). 1035 */ 1036 const double eps = x - 1.0; 1037 const double c1 = 0.4227843350984671394; 1038 const double c2 = -0.01094400467202744461; 1039 const double c3 = 0.09252092391911371098; 1040 const double c4 = -0.018271913165599812664; 1041 const double c5 = 0.018004931096854797895; 1042 const double c6 = -0.006850885378723806846; 1043 const double c7 = 0.003998239557568466030; 1044 result->val = 1.0/x + eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*(c6+eps*c7)))))); 1045 result->err = GSL_DBL_EPSILON; 1046 return GSL_SUCCESS; 1047 } 1048 else if(fabs(x - 2.0) < 0.01) { 1049 /* Use series for Gamma[1 + eps]. 1050 */ 1051 const double eps = x - 2.0; 1052 const double c1 = 0.4227843350984671394; 1053 const double c2 = 0.4118403304264396948; 1054 const double c3 = 0.08157691924708626638; 1055 const double c4 = 0.07424901075351389832; 1056 const double c5 = -0.00026698206874501476832; 1057 const double c6 = 0.011154045718130991049; 1058 const double c7 = -0.002852645821155340816; 1059 const double c8 = 0.0021039333406973880085; 1060 result->val = 1.0 + eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*(c6+eps*(c7+eps*c8))))))); 1061 result->err = GSL_DBL_EPSILON; 1062 return GSL_SUCCESS; 1063 } 1064 else if(x < 5.0) { 1065 /* Exponentiating the logarithm is fine, as 1066 * long as the exponential is not so large 1067 * that it greatly amplifies the error. 1068 */ 1069 gsl_sf_result lg; 1070 lngamma_lanczos(x, &lg); 1071 result->val = exp(lg.val); 1072 result->err = result->val * (lg.err + 2.0 * GSL_DBL_EPSILON); 1073 return GSL_SUCCESS; 1074 } 1075 else if(x < 10.0) { 1076 /* This is a sticky area. The logarithm 1077 * is too large and the gammastar series 1078 * is not good. 1079 */ 1080 const double gamma_8 = 5040.0; 1081 const double t = (2.0*x - 15.0)/5.0; 1082 gsl_sf_result c; 1083 cheb_eval_e(&gamma_5_10_cs, t, &c); 1084 result->val = exp(c.val) * gamma_8; 1085 result->err = result->val * c.err; 1086 result->err += 2.0 * GSL_DBL_EPSILON * result->val; 1087 return GSL_SUCCESS; 1088 } 1089 else if(x < GSL_SF_GAMMA_XMAX) { 1090 /* We do not want to exponentiate the logarithm 1091 * if x is large because of the inevitable 1092 * inflation of the error. So we carefully 1093 * use pow() and exp() with exact quantities. 1094 */ 1095 double p = pow(x, 0.5*x); 1096 double e = exp(-x); 1097 double q = (p * e) * p; 1098 double pre = M_SQRT2 * M_SQRTPI * q/sqrt(x); 1099 gsl_sf_result gstar; 1100 int stat_gs = gammastar_ser(x, &gstar); 1101 result->val = pre * gstar.val; 1102 result->err = (x + 2.5) * GSL_DBL_EPSILON * result->val; 1103 return stat_gs; 1104 } 1105 else { 1106 OVERFLOW_ERROR(result); 1107 } 1108 } 1109 1110 1111 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ 1112 1113 1114 int gsl_sf_lngamma_e(double x, gsl_sf_result * result) 1115 { 1116 /* CHECK_POINTER(result) */ 1117 1118 if(fabs(x - 1.0) < 0.01) { 1119 /* Note that we must amplify the errors 1120 * from the Pade evaluations because of 1121 * the way we must pass the argument, i.e. 1122 * writing (1-x) is a loss of precision 1123 * when x is near 1. 1124 */ 1125 int stat = lngamma_1_pade(x - 1.0, result); 1126 result->err *= 1.0/(GSL_DBL_EPSILON + fabs(x - 1.0)); 1127 return stat; 1128 } 1129 else if(fabs(x - 2.0) < 0.01) { 1130 int stat = lngamma_2_pade(x - 2.0, result); 1131 result->err *= 1.0/(GSL_DBL_EPSILON + fabs(x - 2.0)); 1132 return stat; 1133 } 1134 else if(x >= 0.5) { 1135 return lngamma_lanczos(x, result); 1136 } 1137 else if(x == 0.0) { 1138 DOMAIN_ERROR(result); 1139 } 1140 else if(fabs(x) < 0.02) { 1141 double sgn; 1142 return lngamma_sgn_0(x, result, &sgn); 1143 } 1144 else if(x > -0.5/(GSL_DBL_EPSILON*M_PI)) { 1145 /* Try to extract a fractional 1146 * part from x. 1147 */ 1148 double z = 1.0 - x; 1149 double s = sin(M_PI*z); 1150 double as = fabs(s); 1151 if(s == 0.0) { 1152 DOMAIN_ERROR(result); 1153 } 1154 else if(as < M_PI*0.015) { 1155 /* x is near a negative integer, -N */ 1156 if(x < INT_MIN + 2.0) { 1157 result->val = 0.0; 1158 result->err = 0.0; 1159 GSL_ERROR ("error", GSL_EROUND); 1160 } 1161 else { 1162 int N = -(int)(x - 0.5); 1163 double eps = x + N; 1164 double sgn; 1165 return lngamma_sgn_sing(N, eps, result, &sgn); 1166 } 1167 } 1168 else { 1169 gsl_sf_result lg_z; 1170 lngamma_lanczos(z, &lg_z); 1171 result->val = M_LNPI - (log(as) + lg_z.val); 1172 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + lg_z.err; 1173 return GSL_SUCCESS; 1174 } 1175 } 1176 else { 1177 /* |x| was too large to extract any fractional part */ 1178 result->val = 0.0; 1179 result->err = 0.0; 1180 GSL_ERROR ("error", GSL_EROUND); 1181 } 1182 } 1183 1184 1185 int gsl_sf_lngamma_sgn_e(double x, gsl_sf_result * result_lg, double * sgn) 1186 { 1187 if(fabs(x - 1.0) < 0.01) { 1188 int stat = lngamma_1_pade(x - 1.0, result_lg); 1189 result_lg->err *= 1.0/(GSL_DBL_EPSILON + fabs(x - 1.0)); 1190 *sgn = 1.0; 1191 return stat; 1192 } 1193 else if(fabs(x - 2.0) < 0.01) { 1194 int stat = lngamma_2_pade(x - 2.0, result_lg); 1195 result_lg->err *= 1.0/(GSL_DBL_EPSILON + fabs(x - 2.0)); 1196 *sgn = 1.0; 1197 return stat; 1198 } 1199 else if(x >= 0.5) { 1200 *sgn = 1.0; 1201 return lngamma_lanczos(x, result_lg); 1202 } 1203 else if(x == 0.0) { 1204 *sgn = 0.0; 1205 DOMAIN_ERROR(result_lg); 1206 } 1207 else if(fabs(x) < 0.02) { 1208 return lngamma_sgn_0(x, result_lg, sgn); 1209 } 1210 else if(x > -0.5/(GSL_DBL_EPSILON*M_PI)) { 1211 /* Try to extract a fractional 1212 * part from x. 1213 */ 1214 double z = 1.0 - x; 1215 double s = sin(M_PI*x); 1216 double as = fabs(s); 1217 if(s == 0.0) { 1218 *sgn = 0.0; 1219 DOMAIN_ERROR(result_lg); 1220 } 1221 else if(as < M_PI*0.015) { 1222 /* x is near a negative integer, -N */ 1223 if(x < INT_MIN + 2.0) { 1224 result_lg->val = 0.0; 1225 result_lg->err = 0.0; 1226 *sgn = 0.0; 1227 GSL_ERROR ("error", GSL_EROUND); 1228 } 1229 else { 1230 int N = -(int)(x - 0.5); 1231 double eps = x + N; 1232 return lngamma_sgn_sing(N, eps, result_lg, sgn); 1233 } 1234 } 1235 else { 1236 gsl_sf_result lg_z; 1237 lngamma_lanczos(z, &lg_z); 1238 *sgn = (s > 0.0 ? 1.0 : -1.0); 1239 result_lg->val = M_LNPI - (log(as) + lg_z.val); 1240 result_lg->err = 2.0 * GSL_DBL_EPSILON * fabs(result_lg->val) + lg_z.err; 1241 return GSL_SUCCESS; 1242 } 1243 } 1244 else { 1245 /* |x| was too large to extract any fractional part */ 1246 result_lg->val = 0.0; 1247 result_lg->err = 0.0; 1248 *sgn = 0.0; 1249 GSL_ERROR ("x too large to extract fraction part", GSL_EROUND); 1250 } 1251 } 1252 1253 1254 int 1255 gsl_sf_gamma_e(const double x, gsl_sf_result * result) 1256 { 1257 if(x < 0.5) { 1258 int rint_x = (int)floor(x+0.5); 1259 double f_x = x - rint_x; 1260 double sgn_gamma = ( GSL_IS_EVEN(rint_x) ? 1.0 : -1.0 ); 1261 double sin_term = sgn_gamma * sin(M_PI * f_x) / M_PI; 1262 1263 if(sin_term == 0.0) { 1264 DOMAIN_ERROR(result); 1265 } 1266 else if(x > -169.0) { 1267 gsl_sf_result g; 1268 gamma_xgthalf(1.0-x, &g); 1269 if(fabs(sin_term) * g.val * GSL_DBL_MIN < 1.0) { 1270 result->val = 1.0/(sin_term * g.val); 1271 result->err = fabs(g.err/g.val) * fabs(result->val); 1272 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); 1273 return GSL_SUCCESS; 1274 } 1275 else { 1276 UNDERFLOW_ERROR(result); 1277 } 1278 } 1279 else { 1280 /* It is hard to control it here. 1281 * We can only exponentiate the 1282 * logarithm and eat the loss of 1283 * precision. 1284 */ 1285 gsl_sf_result lng; 1286 double sgn; 1287 int stat_lng = gsl_sf_lngamma_sgn_e(x, &lng, &sgn); 1288 int stat_e = gsl_sf_exp_mult_err_e(lng.val, lng.err, sgn, 0.0, result); 1289 return GSL_ERROR_SELECT_2(stat_e, stat_lng); 1290 } 1291 } 1292 else { 1293 return gamma_xgthalf(x, result); 1294 } 1295 } 1296 1297 1298 int 1299 gsl_sf_gammastar_e(const double x, gsl_sf_result * result) 1300 { 1301 /* CHECK_POINTER(result) */ 1302 1303 if(x <= 0.0) { 1304 DOMAIN_ERROR(result); 1305 } 1306 else if(x < 0.5) { 1307 gsl_sf_result lg; 1308 const int stat_lg = gsl_sf_lngamma_e(x, &lg); 1309 const double lx = log(x); 1310 const double c = 0.5*(M_LN2+M_LNPI); 1311 const double lnr_val = lg.val - (x-0.5)*lx + x - c; 1312 const double lnr_err = lg.err + 2.0 * GSL_DBL_EPSILON *((x+0.5)*fabs(lx) + c); 1313 const int stat_e = gsl_sf_exp_err_e(lnr_val, lnr_err, result); 1314 return GSL_ERROR_SELECT_2(stat_lg, stat_e); 1315 } 1316 else if(x < 2.0) { 1317 const double t = 4.0/3.0*(x-0.5) - 1.0; 1318 return cheb_eval_e(&gstar_a_cs, t, result); 1319 } 1320 else if(x < 10.0) { 1321 const double t = 0.25*(x-2.0) - 1.0; 1322 gsl_sf_result c; 1323 cheb_eval_e(&gstar_b_cs, t, &c); 1324 result->val = c.val/(x*x) + 1.0 + 1.0/(12.0*x); 1325 result->err = c.err/(x*x); 1326 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); 1327 return GSL_SUCCESS; 1328 } 1329 else if(x < 1.0/GSL_ROOT4_DBL_EPSILON) { 1330 return gammastar_ser(x, result); 1331 } 1332 else if(x < 1.0/GSL_DBL_EPSILON) { 1333 /* Use Stirling formula for Gamma(x). 1334 */ 1335 const double xi = 1.0/x; 1336 result->val = 1.0 + xi/12.0*(1.0 + xi/24.0*(1.0 - xi*(139.0/180.0 + 571.0/8640.0*xi))); 1337 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 1338 return GSL_SUCCESS; 1339 } 1340 else { 1341 result->val = 1.0; 1342 result->err = 1.0/x; 1343 return GSL_SUCCESS; 1344 } 1345 } 1346 1347 1348 int 1349 gsl_sf_gammainv_e(const double x, gsl_sf_result * result) 1350 { 1351 /* CHECK_POINTER(result) */ 1352 1353 if (x <= 0.0 && x == floor(x)) { /* negative integer */ 1354 result->val = 0.0; 1355 result->err = 0.0; 1356 return GSL_SUCCESS; 1357 } else if(x < 0.5) { 1358 gsl_sf_result lng; 1359 double sgn; 1360 int stat_lng = gsl_sf_lngamma_sgn_e(x, &lng, &sgn); 1361 if(stat_lng == GSL_EDOM) { 1362 result->val = 0.0; 1363 result->err = 0.0; 1364 return GSL_SUCCESS; 1365 } 1366 else if(stat_lng != GSL_SUCCESS) { 1367 result->val = 0.0; 1368 result->err = 0.0; 1369 return stat_lng; 1370 } 1371 else { 1372 return gsl_sf_exp_mult_err_e(-lng.val, lng.err, sgn, 0.0, result); 1373 } 1374 } 1375 else { 1376 gsl_sf_result g; 1377 int stat_g = gamma_xgthalf(x, &g); 1378 if(stat_g == GSL_EOVRFLW) { 1379 UNDERFLOW_ERROR(result); 1380 } 1381 else { 1382 result->val = 1.0/g.val; 1383 result->err = fabs(g.err/g.val) * fabs(result->val); 1384 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); 1385 CHECK_UNDERFLOW(result); 1386 return GSL_SUCCESS; 1387 } 1388 } 1389 } 1390 1391 1392 int 1393 gsl_sf_lngamma_complex_e(double zr, double zi, gsl_sf_result * lnr, gsl_sf_result * arg) 1394 { 1395 if(zr <= 0.5) { 1396 /* Transform to right half plane using reflection; 1397 * in fact we do a little better by stopping at 1/2. 1398 */ 1399 double x = 1.0-zr; 1400 double y = -zi; 1401 gsl_sf_result a, b; 1402 gsl_sf_result lnsin_r, lnsin_i; 1403 1404 int stat_l = lngamma_lanczos_complex(x, y, &a, &b); 1405 int stat_s = gsl_sf_complex_logsin_e(M_PI*zr, M_PI*zi, &lnsin_r, &lnsin_i); 1406 1407 if(stat_s == GSL_SUCCESS) { 1408 int stat_r; 1409 lnr->val = M_LNPI - lnsin_r.val - a.val; 1410 lnr->err = lnsin_r.err + a.err + 2.0 * GSL_DBL_EPSILON * fabs(lnr->val); 1411 arg->val = -lnsin_i.val - b.val; 1412 arg->err = lnsin_i.err + b.err + 2.0 * GSL_DBL_EPSILON * fabs(arg->val); 1413 stat_r = gsl_sf_angle_restrict_symm_e(&(arg->val)); 1414 return GSL_ERROR_SELECT_2(stat_r, stat_l); 1415 } 1416 else { 1417 DOMAIN_ERROR_2(lnr,arg); 1418 } 1419 } 1420 else { 1421 /* otherwise plain vanilla Lanczos */ 1422 return lngamma_lanczos_complex(zr, zi, lnr, arg); 1423 } 1424 } 1425 1426 1427 int gsl_sf_taylorcoeff_e(const int n, const double x, gsl_sf_result * result) 1428 { 1429 /* CHECK_POINTER(result) */ 1430 1431 if(x < 0.0 || n < 0) { 1432 DOMAIN_ERROR(result); 1433 } 1434 else if(n == 0) { 1435 result->val = 1.0; 1436 result->err = 0.0; 1437 return GSL_SUCCESS; 1438 } 1439 else if(n == 1) { 1440 result->val = x; 1441 result->err = 0.0; 1442 return GSL_SUCCESS; 1443 } 1444 else if(x == 0.0) { 1445 result->val = 0.0; 1446 result->err = 0.0; 1447 return GSL_SUCCESS; 1448 } 1449 else { 1450 const double log2pi = M_LNPI + M_LN2; 1451 const double ln_test = n*(log(x)+1.0) + 1.0 - (n+0.5)*log(n+1.0) + 0.5*log2pi; 1452 1453 if(ln_test < GSL_LOG_DBL_MIN+1.0) { 1454 UNDERFLOW_ERROR(result); 1455 } 1456 else if(ln_test > GSL_LOG_DBL_MAX-1.0) { 1457 OVERFLOW_ERROR(result); 1458 } 1459 else { 1460 double product = 1.0; 1461 int k; 1462 for(k=1; k<=n; k++) { 1463 product *= (x/k); 1464 } 1465 result->val = product; 1466 result->err = n * GSL_DBL_EPSILON * product; 1467 CHECK_UNDERFLOW(result); 1468 return GSL_SUCCESS; 1469 } 1470 } 1471 } 1472 1473 1474 int gsl_sf_fact_e(const unsigned int n, gsl_sf_result * result) 1475 { 1476 /* CHECK_POINTER(result) */ 1477 1478 if(n < 18) { 1479 result->val = fact_table[n].f; 1480 result->err = 0.0; 1481 return GSL_SUCCESS; 1482 } 1483 else if(n <= GSL_SF_FACT_NMAX){ 1484 result->val = fact_table[n].f; 1485 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 1486 return GSL_SUCCESS; 1487 } 1488 else { 1489 OVERFLOW_ERROR(result); 1490 } 1491 } 1492 1493 1494 int gsl_sf_doublefact_e(const unsigned int n, gsl_sf_result * result) 1495 { 1496 /* CHECK_POINTER(result) */ 1497 1498 if(n < 26) { 1499 result->val = doub_fact_table[n].f; 1500 result->err = 0.0; 1501 return GSL_SUCCESS; 1502 } 1503 else if(n <= GSL_SF_DOUBLEFACT_NMAX){ 1504 result->val = doub_fact_table[n].f; 1505 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 1506 return GSL_SUCCESS; 1507 } 1508 else { 1509 OVERFLOW_ERROR(result); 1510 } 1511 } 1512 1513 1514 int gsl_sf_lnfact_e(const unsigned int n, gsl_sf_result * result) 1515 { 1516 /* CHECK_POINTER(result) */ 1517 1518 if(n <= GSL_SF_FACT_NMAX){ 1519 result->val = log(fact_table[n].f); 1520 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 1521 return GSL_SUCCESS; 1522 } 1523 else { 1524 gsl_sf_lngamma_e(n+1.0, result); 1525 return GSL_SUCCESS; 1526 } 1527 } 1528 1529 1530 int gsl_sf_lndoublefact_e(const unsigned int n, gsl_sf_result * result) 1531 { 1532 /* CHECK_POINTER(result) */ 1533 1534 if(n <= GSL_SF_DOUBLEFACT_NMAX){ 1535 result->val = log(doub_fact_table[n].f); 1536 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); 1537 return GSL_SUCCESS; 1538 } 1539 else if(GSL_IS_ODD(n)) { 1540 gsl_sf_result lg; 1541 gsl_sf_lngamma_e(0.5*(n+2.0), &lg); 1542 result->val = 0.5*(n+1.0) * M_LN2 - 0.5*M_LNPI + lg.val; 1543 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + lg.err; 1544 return GSL_SUCCESS; 1545 } 1546 else { 1547 gsl_sf_result lg; 1548 gsl_sf_lngamma_e(0.5*n+1.0, &lg); 1549 result->val = 0.5*n*M_LN2 + lg.val; 1550 result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + lg.err; 1551 return GSL_SUCCESS; 1552 } 1553 } 1554 1555 1556 int gsl_sf_lnchoose_e(unsigned int n, unsigned int m, gsl_sf_result * result) 1557 { 1558 /* CHECK_POINTER(result) */ 1559 1560 if(m > n) { 1561 DOMAIN_ERROR(result); 1562 } 1563 else if(m == n || m == 0) { 1564 result->val = 0.0; 1565 result->err = 0.0; 1566 return GSL_SUCCESS; 1567 } 1568 else { 1569 gsl_sf_result nf; 1570 gsl_sf_result mf; 1571 gsl_sf_result nmmf; 1572 if(m*2 > n) m = n-m; 1573 gsl_sf_lnfact_e(n, &nf); 1574 gsl_sf_lnfact_e(m, &mf); 1575 gsl_sf_lnfact_e(n-m, &nmmf); 1576 result->val = nf.val - mf.val - nmmf.val; 1577 result->err = nf.err + mf.err + nmmf.err; 1578 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); 1579 return GSL_SUCCESS; 1580 } 1581 } 1582 1583 1584 int gsl_sf_choose_e(unsigned int n, unsigned int m, gsl_sf_result * result) 1585 { 1586 if(m > n) { 1587 DOMAIN_ERROR(result); 1588 } 1589 else if(m == n || m == 0) { 1590 result->val = 1.0; 1591 result->err = 0.0; 1592 return GSL_SUCCESS; 1593 } 1594 else if (n <= GSL_SF_FACT_NMAX) { 1595 result->val = (fact_table[n].f / fact_table[m].f) / fact_table[n-m].f; 1596 result->err = 6.0 * GSL_DBL_EPSILON * fabs(result->val); 1597 return GSL_SUCCESS; 1598 } else { 1599 if(m*2 < n) m = n-m; 1600 1601 if (n - m < 64) /* compute product for a manageable number of terms */ 1602 { 1603 double prod = 1.0; 1604 unsigned int k; 1605 1606 for(k=n; k>=m+1; k--) { 1607 double tk = (double)k / (double)(k-m); 1608 if(tk > GSL_DBL_MAX/prod) { 1609 OVERFLOW_ERROR(result); 1610 } 1611 prod *= tk; 1612 } 1613 result->val = prod; 1614 result->err = 2.0 * GSL_DBL_EPSILON * prod * fabs(n-m); 1615 return GSL_SUCCESS; 1616 } 1617 else 1618 { 1619 gsl_sf_result lc; 1620 const int stat_lc = gsl_sf_lnchoose_e (n, m, &lc); 1621 const int stat_e = gsl_sf_exp_err_e(lc.val, lc.err, result); 1622 return GSL_ERROR_SELECT_2(stat_lc, stat_e); 1623 } 1624 } 1625 } 1626 1627 1628 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ 1629 1630 #include "gsl_specfunc__eval.h" 1631 1632 double gsl_sf_fact(const unsigned int n) 1633 { 1634 EVAL_RESULT(gsl_sf_fact_e(n, &result)); 1635 } 1636 1637 double gsl_sf_lnfact(const unsigned int n) 1638 { 1639 EVAL_RESULT(gsl_sf_lnfact_e(n, &result)); 1640 } 1641 1642 double gsl_sf_doublefact(const unsigned int n) 1643 { 1644 EVAL_RESULT(gsl_sf_doublefact_e(n, &result)); 1645 } 1646 1647 double gsl_sf_lndoublefact(const unsigned int n) 1648 { 1649 EVAL_RESULT(gsl_sf_lndoublefact_e(n, &result)); 1650 } 1651 1652 double gsl_sf_lngamma(const double x) 1653 { 1654 EVAL_RESULT(gsl_sf_lngamma_e(x, &result)); 1655 } 1656 1657 double gsl_sf_gamma(const double x) 1658 { 1659 EVAL_RESULT(gsl_sf_gamma_e(x, &result)); 1660 } 1661 1662 double gsl_sf_gammastar(const double x) 1663 { 1664 EVAL_RESULT(gsl_sf_gammastar_e(x, &result)); 1665 } 1666 1667 double gsl_sf_gammainv(const double x) 1668 { 1669 EVAL_RESULT(gsl_sf_gammainv_e(x, &result)); 1670 } 1671 1672 double gsl_sf_taylorcoeff(const int n, const double x) 1673 { 1674 EVAL_RESULT(gsl_sf_taylorcoeff_e(n, x, &result)); 1675 } 1676 1677 double gsl_sf_choose(unsigned int n, unsigned int m) 1678 { 1679 EVAL_RESULT(gsl_sf_choose_e(n, m, &result)); 1680 } 1681 1682 double gsl_sf_lnchoose(unsigned int n, unsigned int m) 1683 { 1684 EVAL_RESULT(gsl_sf_lnchoose_e(n, m, &result)); 1685 } 1686