static char rcsid[] = "$Id: pbinom.c 184472 2016-02-18 00:12:16Z twu $"; #ifdef HAVE_CONFIG_H #include #endif #include "pbinom.h" #include #include #include #include /* #include */ /************************************************************************ * Code taken from gsl (GNU scientific library) version 1.8 ************************************************************************/ /* from gsl/gsl_machine.h */ /* magic constants; mostly for the benefit of the implementation */ /* -*-MACHINE CONSTANTS-*- * * PLATFORM: Whiz-O-Matic 9000 * FP_PLATFORM: IEEE-Virtual * HOSTNAME: nnn.lanl.gov * DATE: Fri Nov 20 17:53:26 MST 1998 */ #define GSL_DBL_EPSILON 2.2204460492503131e-16 #define GSL_SQRT_DBL_EPSILON 1.4901161193847656e-08 #define GSL_ROOT3_DBL_EPSILON 6.0554544523933429e-06 #define GSL_ROOT4_DBL_EPSILON 1.2207031250000000e-04 #define GSL_ROOT5_DBL_EPSILON 7.4009597974140505e-04 #define GSL_ROOT6_DBL_EPSILON 2.4607833005759251e-03 #define GSL_LOG_DBL_EPSILON (-3.6043653389117154e+01) #define GSL_DBL_MIN 2.2250738585072014e-308 #define GSL_SQRT_DBL_MIN 1.4916681462400413e-154 #define GSL_ROOT3_DBL_MIN 2.8126442852362996e-103 #define GSL_ROOT4_DBL_MIN 1.2213386697554620e-77 #define GSL_ROOT5_DBL_MIN 2.9476022969691763e-62 #define GSL_ROOT6_DBL_MIN 5.3034368905798218e-52 #define GSL_LOG_DBL_MIN (-7.0839641853226408e+02) #define GSL_DBL_MAX 1.7976931348623157e+308 #define GSL_SQRT_DBL_MAX 1.3407807929942596e+154 #define GSL_ROOT3_DBL_MAX 5.6438030941222897e+102 #define GSL_ROOT4_DBL_MAX 1.1579208923731620e+77 #define GSL_ROOT5_DBL_MAX 4.4765466227572707e+61 #define GSL_ROOT6_DBL_MAX 2.3756689782295612e+51 #define GSL_LOG_DBL_MAX 7.0978271289338397e+02 #define GSL_FLT_EPSILON 1.1920928955078125e-07 #define GSL_SQRT_FLT_EPSILON 3.4526698300124393e-04 #define GSL_ROOT3_FLT_EPSILON 4.9215666011518501e-03 #define GSL_ROOT4_FLT_EPSILON 1.8581361171917516e-02 #define GSL_ROOT5_FLT_EPSILON 4.1234622211652937e-02 #define GSL_ROOT6_FLT_EPSILON 7.0153878019335827e-02 #define GSL_LOG_FLT_EPSILON (-1.5942385152878742e+01) #define GSL_FLT_MIN 1.1754943508222875e-38 #define GSL_SQRT_FLT_MIN 1.0842021724855044e-19 #define GSL_ROOT3_FLT_MIN 2.2737367544323241e-13 #define GSL_ROOT4_FLT_MIN 3.2927225399135965e-10 #define GSL_ROOT5_FLT_MIN 2.5944428542140822e-08 #define GSL_ROOT6_FLT_MIN 4.7683715820312542e-07 #define GSL_LOG_FLT_MIN (-8.7336544750553102e+01) #define GSL_FLT_MAX 3.4028234663852886e+38 #define GSL_SQRT_FLT_MAX 1.8446743523953730e+19 #define GSL_ROOT3_FLT_MAX 6.9814635196223242e+12 #define GSL_ROOT4_FLT_MAX 4.2949672319999986e+09 #define GSL_ROOT5_FLT_MAX 5.0859007855960041e+07 #define GSL_ROOT6_FLT_MAX 2.6422459233807749e+06 #define GSL_LOG_FLT_MAX 8.8722839052068352e+01 #define GSL_SFLT_EPSILON 4.8828125000000000e-04 #define GSL_SQRT_SFLT_EPSILON 2.2097086912079612e-02 #define GSL_ROOT3_SFLT_EPSILON 7.8745065618429588e-02 #define GSL_ROOT4_SFLT_EPSILON 1.4865088937534013e-01 #define GSL_ROOT5_SFLT_EPSILON 2.1763764082403100e-01 #define GSL_ROOT6_SFLT_EPSILON 2.8061551207734325e-01 #define GSL_LOG_SFLT_EPSILON (-7.6246189861593985e+00) /* !MACHINE CONSTANTS! */ #define OVERFLOW_ERROR GSL_DBL_MAX #define UNDERFLOW_ERROR 0.0 /* gsl_math.h */ #define GSL_MAX(a,b) ((a) > (b) ? (a) : (b)) #define GSL_MIN(a,b) ((a) < (b) ? (a) : (b)) #define GSL_IS_ODD(n) ((n) & 1) #define GSL_IS_EVEN(n) (!(GSL_IS_ODD(n))) #define GSL_SIGN(x) ((x) >= 0.0 ? 1 : -1) #ifndef M_E #define M_E 2.71828182845904523536028747135 /* e */ #endif #ifndef M_SQRT2 #define M_SQRT2 1.41421356237309504880168872421 /* sqrt(2) */ #endif #ifndef M_PI #define M_PI 3.14159265358979323846264338328 /* pi */ #endif #ifndef M_SQRTPI #define M_SQRTPI 1.77245385090551602729816748334 /* sqrt(pi) */ #endif #ifndef M_LN2 #define M_LN2 0.69314718055994530941723212146 /* ln(2) */ #endif #ifndef M_LNPI #define M_LNPI 1.14472988584940017414342735135 /* ln(pi) */ #endif #ifndef M_EULER #define M_EULER 0.57721566490153286060651209008 /* Euler constant */ #endif #define LogRootTwoPi_ 0.9189385332046727418 /************************************************************************ * Chebyshev ************************************************************************/ /* data for a Chebyshev series over a given interval */ struct cheb_series_struct { double * c; /* coefficients */ int order; /* order of expansion */ double a; /* lower interval point */ double b; /* upper interval point */ int order_sp; /* effective single precision order */ }; typedef struct cheb_series_struct cheb_series; static inline double cheb_eval (const cheb_series * cs, const double x) { int j; double d = 0.0; double dd = 0.0; double y = (2.0*x - cs->a - cs->b) / (cs->b - cs->a); double y2 = 2.0 * y; double e = 0.0; for (j = cs->order; j>=1; j--) { double temp = d; d = y2*d - dd + cs->c[j]; e += fabs(y2*temp) + fabs(dd) + fabs(cs->c[j]); dd = temp; } { double temp = d; d = y*d - dd + 0.5 * cs->c[0]; e += fabs(y*temp) + fabs(dd) + 0.5 * fabs(cs->c[0]); } return d; } /* The maximum n such that gsl_sf_fact(n) does not give an overflow. */ #define GSL_SF_FACT_NMAX 170 static struct {int n; double f; long i; } fact_table[GSL_SF_FACT_NMAX + 1] = { { 0, 1.0, 1L }, { 1, 1.0, 1L }, { 2, 2.0, 2L }, { 3, 6.0, 6L }, { 4, 24.0, 24L }, { 5, 120.0, 120L }, { 6, 720.0, 720L }, { 7, 5040.0, 5040L }, { 8, 40320.0, 40320L }, { 9, 362880.0, 362880L }, { 10, 3628800.0, 3628800L }, { 11, 39916800.0, 39916800L }, { 12, 479001600.0, 479001600L }, { 13, 6227020800.0, 0 }, { 14, 87178291200.0, 0 }, { 15, 1307674368000.0, 0 }, { 16, 20922789888000.0, 0 }, { 17, 355687428096000.0, 0 }, { 18, 6402373705728000.0, 0 }, { 19, 121645100408832000.0, 0 }, { 20, 2432902008176640000.0, 0 }, { 21, 51090942171709440000.0, 0 }, { 22, 1124000727777607680000.0, 0 }, { 23, 25852016738884976640000.0, 0 }, { 24, 620448401733239439360000.0, 0 }, { 25, 15511210043330985984000000.0, 0 }, { 26, 403291461126605635584000000.0, 0 }, { 27, 10888869450418352160768000000.0, 0 }, { 28, 304888344611713860501504000000.0, 0 }, { 29, 8841761993739701954543616000000.0, 0 }, { 30, 265252859812191058636308480000000.0, 0 }, { 31, 8222838654177922817725562880000000.0, 0 }, { 32, 263130836933693530167218012160000000.0, 0 }, { 33, 8683317618811886495518194401280000000.0, 0 }, { 34, 2.95232799039604140847618609644e38, 0 }, { 35, 1.03331479663861449296666513375e40, 0 }, { 36, 3.71993326789901217467999448151e41, 0 }, { 37, 1.37637530912263450463159795816e43, 0 }, { 38, 5.23022617466601111760007224100e44, 0 }, { 39, 2.03978820811974433586402817399e46, 0 }, { 40, 8.15915283247897734345611269600e47, 0 }, { 41, 3.34525266131638071081700620534e49, 0 }, { 42, 1.40500611775287989854314260624e51, 0 }, { 43, 6.04152630633738356373551320685e52, 0 }, { 44, 2.65827157478844876804362581101e54, 0 }, { 45, 1.19622220865480194561963161496e56, 0 }, { 46, 5.50262215981208894985030542880e57, 0 }, { 47, 2.58623241511168180642964355154e59, 0 }, { 48, 1.24139155925360726708622890474e61, 0 }, { 49, 6.08281864034267560872252163321e62, 0 }, { 50, 3.04140932017133780436126081661e64, 0 }, { 51, 1.55111875328738228022424301647e66, 0 }, { 52, 8.06581751709438785716606368564e67, 0 }, { 53, 4.27488328406002556429801375339e69, 0 }, { 54, 2.30843697339241380472092742683e71, 0 }, { 55, 1.26964033536582759259651008476e73, 0 }, { 56, 7.10998587804863451854045647464e74, 0 }, { 57, 4.05269195048772167556806019054e76, 0 }, { 58, 2.35056133128287857182947491052e78, 0 }, { 59, 1.38683118545689835737939019720e80, 0 }, { 60, 8.32098711274139014427634118320e81, 0 }, { 61, 5.07580213877224798800856812177e83, 0 }, { 62, 3.14699732603879375256531223550e85, 0 }, { 63, 1.982608315404440064116146708360e87, 0 }, { 64, 1.268869321858841641034333893350e89, 0 }, { 65, 8.247650592082470666723170306800e90, 0 }, { 66, 5.443449390774430640037292402480e92, 0 }, { 67, 3.647111091818868528824985909660e94, 0 }, { 68, 2.480035542436830599600990418570e96, 0 }, { 69, 1.711224524281413113724683388810e98, 0 }, { 70, 1.197857166996989179607278372170e100, 0 }, { 71, 8.504785885678623175211676442400e101, 0 }, { 72, 6.123445837688608686152407038530e103, 0 }, { 73, 4.470115461512684340891257138130e105, 0 }, { 74, 3.307885441519386412259530282210e107, 0 }, { 75, 2.480914081139539809194647711660e109, 0 }, { 76, 1.885494701666050254987932260860e111, 0 }, { 77, 1.451830920282858696340707840860e113, 0 }, { 78, 1.132428117820629783145752115870e115, 0 }, { 79, 8.946182130782975286851441715400e116, 0 }, { 80, 7.156945704626380229481153372320e118, 0 }, { 81, 5.797126020747367985879734231580e120, 0 }, { 82, 4.753643337012841748421382069890e122, 0 }, { 83, 3.945523969720658651189747118010e124, 0 }, { 84, 3.314240134565353266999387579130e126, 0 }, { 85, 2.817104114380550276949479442260e128, 0 }, { 86, 2.422709538367273238176552320340e130, 0 }, { 87, 2.107757298379527717213600518700e132, 0 }, { 88, 1.854826422573984391147968456460e134, 0 }, { 89, 1.650795516090846108121691926250e136, 0 }, { 90, 1.485715964481761497309522733620e138, 0 }, { 91, 1.352001527678402962551665687590e140, 0 }, { 92, 1.243841405464130725547532432590e142, 0 }, { 93, 1.156772507081641574759205162310e144, 0 }, { 94, 1.087366156656743080273652852570e146, 0 }, { 95, 1.032997848823905926259970209940e148, 0 }, { 96, 9.916779348709496892095714015400e149, 0 }, { 97, 9.619275968248211985332842594960e151, 0 }, { 98, 9.426890448883247745626185743100e153, 0 }, { 99, 9.332621544394415268169923885600e155, 0 }, { 100, 9.33262154439441526816992388563e157, 0 }, { 101, 9.42594775983835942085162312450e159, 0 }, { 102, 9.61446671503512660926865558700e161, 0 }, { 103, 9.90290071648618040754671525458e163, 0 }, { 104, 1.02990167451456276238485838648e166, 0 }, { 105, 1.08139675824029090050410130580e168, 0 }, { 106, 1.146280563734708354534347384148e170, 0 }, { 107, 1.226520203196137939351751701040e172, 0 }, { 108, 1.324641819451828974499891837120e174, 0 }, { 109, 1.443859583202493582204882102460e176, 0 }, { 110, 1.588245541522742940425370312710e178, 0 }, { 111, 1.762952551090244663872161047110e180, 0 }, { 112, 1.974506857221074023536820372760e182, 0 }, { 113, 2.231192748659813646596607021220e184, 0 }, { 114, 2.543559733472187557120132004190e186, 0 }, { 115, 2.925093693493015690688151804820e188, 0 }, { 116, 3.393108684451898201198256093590e190, 0 }, { 117, 3.96993716080872089540195962950e192, 0 }, { 118, 4.68452584975429065657431236281e194, 0 }, { 119, 5.57458576120760588132343171174e196, 0 }, { 120, 6.68950291344912705758811805409e198, 0 }, { 121, 8.09429852527344373968162284545e200, 0 }, { 122, 9.87504420083360136241157987140e202, 0 }, { 123, 1.21463043670253296757662432419e205, 0 }, { 124, 1.50614174151114087979501416199e207, 0 }, { 125, 1.88267717688892609974376770249e209, 0 }, { 126, 2.37217324288004688567714730514e211, 0 }, { 127, 3.01266001845765954480997707753e213, 0 }, { 128, 3.85620482362580421735677065923e215, 0 }, { 129, 4.97450422247728744039023415041e217, 0 }, { 130, 6.46685548922047367250730439554e219, 0 }, { 131, 8.47158069087882051098456875820e221, 0 }, { 132, 1.11824865119600430744996307608e224, 0 }, { 133, 1.48727070609068572890845089118e226, 0 }, { 134, 1.99294274616151887673732419418e228, 0 }, { 135, 2.69047270731805048359538766215e230, 0 }, { 136, 3.65904288195254865768972722052e232, 0 }, { 137, 5.01288874827499166103492629211e234, 0 }, { 138, 6.91778647261948849222819828311e236, 0 }, { 139, 9.61572319694108900419719561353e238, 0 }, { 140, 1.34620124757175246058760738589e241, 0 }, { 141, 1.89814375907617096942852641411e243, 0 }, { 142, 2.69536413788816277658850750804e245, 0 }, { 143, 3.85437071718007277052156573649e247, 0 }, { 144, 5.55029383273930478955105466055e249, 0 }, { 145, 8.04792605747199194484902925780e251, 0 }, { 146, 1.17499720439091082394795827164e254, 0 }, { 147, 1.72724589045463891120349865931e256, 0 }, { 148, 2.55632391787286558858117801578e258, 0 }, { 149, 3.80892263763056972698595524351e260, 0 }, { 150, 5.71338395644585459047893286526e262, 0 }, { 151, 8.62720977423324043162318862650e264, 0 }, { 152, 1.31133588568345254560672467123e267, 0 }, { 153, 2.00634390509568239477828874699e269, 0 }, { 154, 3.08976961384735088795856467036e271, 0 }, { 155, 4.78914290146339387633577523906e273, 0 }, { 156, 7.47106292628289444708380937294e275, 0 }, { 157, 1.17295687942641442819215807155e278, 0 }, { 158, 1.85327186949373479654360975305e280, 0 }, { 159, 2.94670227249503832650433950735e282, 0 }, { 160, 4.71472363599206132240694321176e284, 0 }, { 161, 7.59070505394721872907517857094e286, 0 }, { 162, 1.22969421873944943411017892849e289, 0 }, { 163, 2.00440157654530257759959165344e291, 0 }, { 164, 3.28721858553429622726333031164e293, 0 }, { 165, 5.42391066613158877498449501421e295, 0 }, { 166, 9.00369170577843736647426172359e297, 0 }, { 167, 1.50361651486499904020120170784e300, 0 }, { 168, 2.52607574497319838753801886917e302, 0 }, { 169, 4.26906800900470527493925188890e304, 0 }, { 170, 7.25741561530799896739672821113e306, 0 }, /* { 171, 1.24101807021766782342484052410e309, 0 }, { 172, 2.13455108077438865629072570146e311, 0 }, { 173, 3.69277336973969237538295546352e313, 0 }, { 174, 6.42542566334706473316634250653e315, 0 }, { 175, 1.12444949108573632830410993864e318, 0 }, { 176, 1.97903110431089593781523349201e320, 0 }, { 177, 3.50288505463028580993296328086e322, 0 }, { 178, 6.23513539724190874168067463993e324, 0 }, { 179, 1.11608923610630166476084076055e327, 0 }, { 180, 2.00896062499134299656951336898e329, 0 }, { 181, 3.63621873123433082379081919786e331, 0 }, { 182, 6.61791809084648209929929094011e333, 0 }, { 183, 1.21107901062490622417177024204e336, 0 }, { 184, 2.22838537954982745247605724535e338, 0 }, { 185, 4.12251295216718078708070590390e340, 0 }, { 186, 7.66787409103095626397011298130e342, 0 }, { 187, 1.43389245502278882136241112750e345, 0 }, { 188, 2.69571781544284298416133291969e347, 0 }, { 189, 5.09490667118697324006491921822e349, 0 }, { 190, 9.68032267525524915612334651460e351, 0 }, { 191, 1.84894163097375258881955918429e354, 0 }, { 192, 3.54996793146960497053355363384e356, 0 }, { 193, 6.85143810773633759312975851330e358, 0 }, { 194, 1.32917899290084949306717315158e361, 0 }, { 195, 2.59189903615665651148098764559e363, 0 }, { 196, 5.08012211086704676250273578535e365, 0 }, { 197, 1.00078405584080821221303894971e368, 0 }, { 198, 1.98155243056480026018181712043e370, 0 }, { 199, 3.94328933682395251776181606966e372, 0 }, { 200, 7.88657867364790503552363213932e374, 0 } */ }; #define PSI_TABLE_NMAX 100 static double psi_table[PSI_TABLE_NMAX+1] = { 0.0, /* Infinity */ /* psi(0) */ -M_EULER, /* psi(1) */ 0.42278433509846713939348790992, /* ... */ 0.92278433509846713939348790992, 1.25611766843180047272682124325, 1.50611766843180047272682124325, 1.70611766843180047272682124325, 1.87278433509846713939348790992, 2.01564147795560999653634505277, 2.14064147795560999653634505277, 2.25175258906672110764745616389, 2.35175258906672110764745616389, 2.44266167997581201673836525479, 2.52599501330914535007169858813, 2.60291809023222227314862166505, 2.67434666166079370172005023648, 2.74101332832746036838671690315, 2.80351332832746036838671690315, 2.86233685773922507426906984432, 2.91789241329478062982462539988, 2.97052399224214905087725697883, 3.02052399224214905087725697883, 3.06814303986119666992487602645, 3.11359758531574212447033057190, 3.15707584618530734186163491973, 3.1987425128519740085283015864, 3.2387425128519740085283015864, 3.2772040513135124700667631249, 3.3142410883505495071038001619, 3.3499553740648352213895144476, 3.3844381326855248765619282407, 3.4177714660188582098952615740, 3.4500295305349872421533260902, 3.4812795305349872421533260902, 3.5115825608380175451836291205, 3.5409943255438998981248055911, 3.5695657541153284695533770196, 3.5973435318931062473311547974, 3.6243705589201332743581818244, 3.6506863483938174848844976139, 3.6763273740348431259101386396, 3.7013273740348431259101386396, 3.7257176179372821503003825420, 3.7495271417468059598241920658, 3.7727829557002943319172153216, 3.7955102284275670591899425943, 3.8177324506497892814121648166, 3.8394715810845718901078169905, 3.8607481768292527411716467777, 3.8815815101625860745049801110, 3.9019896734278921969539597029, 3.9219896734278921969539597029, 3.9415975165651470989147440166, 3.9608282857959163296839747858, 3.9796962103242182164764276160, 3.9982147288427367349949461345, 4.0163965470245549168131279527, 4.0342536898816977739559850956, 4.0517975495308205809735289552, 4.0690389288411654085597358518, 4.0859880813835382899156680552, 4.1026547480502049565823347218, 4.1190481906731557762544658694, 4.1351772229312202923834981274, 4.1510502388042361653993711433, 4.1666752388042361653993711433, 4.1820598541888515500147557587, 4.1972113693403667015299072739, 4.2121367424746950597388624977, 4.2268426248276362362094507330, 4.2413353784508246420065521823, 4.2556210927365389277208378966, 4.2697055997787924488475984600, 4.2835944886676813377364873489, 4.2972931188046676391063503626, 4.3108066323181811526198638761, 4.3241399656515144859531972094, 4.3372978603883565912163551041, 4.3502848733753695782293421171, 4.3631053861958823987421626300, 4.3757636140439836645649474401, 4.3882636140439836645649474401, 4.4006092930563293435772931191, 4.4128044150075488557724150703, 4.4248526077786331931218126607, 4.4367573696833950978837174226, 4.4485220755657480390601880108, 4.4601499825424922251066996387, 4.4716442354160554434975042364, 4.4830078717796918071338678728, 4.4942438268358715824147667492, 4.5053549379469826935258778603, 4.5163439489359936825368668713, 4.5272135141533849868846929582, 4.5379662023254279976373811303, 4.5486045001977684231692960239, 4.5591308159872421073798223397, 4.5695474826539087740464890064, 4.5798567610044242379640147796, 4.5900608426370772991885045755, 4.6001618527380874001986055856 }; #define PSI_1_TABLE_NMAX 100 static double psi_1_table[PSI_1_TABLE_NMAX+1] = { 0.0, /* Infinity */ /* psi(1,0) */ M_PI*M_PI/6.0, /* psi(1,1) */ 0.644934066848226436472415, /* ... */ 0.394934066848226436472415, 0.2838229557371153253613041, 0.2213229557371153253613041, 0.1813229557371153253613041, 0.1535451779593375475835263, 0.1331370146940314251345467, 0.1175120146940314251345467, 0.1051663356816857461222010, 0.0951663356816857461222010, 0.0869018728717683907503002, 0.0799574284273239463058557, 0.0740402686640103368384001, 0.0689382278476838062261552, 0.0644937834032393617817108, 0.0605875334032393617817108, 0.0571273257907826143768665, 0.0540409060376961946237801, 0.0512708229352031198315363, 0.0487708229352031198315363, 0.0465032492390579951149830, 0.0444371335365786562720078, 0.0425467743683366902984728, 0.0408106632572255791873617, 0.0392106632572255791873617, 0.0377313733163971768204978, 0.0363596312039143235969038, 0.0350841209998326909438426, 0.0338950603577399442137594, 0.0327839492466288331026483, 0.0317433665203020901265817, 0.03076680402030209012658168, 0.02984853037475571730748159, 0.02898347847164153045627052, 0.02816715194102928555831133, 0.02739554700275768062003973, 0.02666508681283803124093089, 0.02597256603721476254286995, 0.02531510384129102815759710, 0.02469010384129102815759710, 0.02409521984367056414807896, 0.02352832641963428296894063, 0.02298749353699501850166102, 0.02247096461137518379091722, 0.02197713745088135663042339, 0.02150454765882086513703965, 0.02105185413233829383780923, 0.02061782635456051606003145, 0.02020133322669712580597065, 0.01980133322669712580597065, 0.01941686571420193164987683, 0.01904704322899483105816086, 0.01869104465298913508094477, 0.01834810912486842177504628, 0.01801753061247172756017024, 0.01769865306145131939690494, 0.01739086605006319997554452, 0.01709360088954001329302371, 0.01680632711763538818529605, 0.01652854933985761040751827, 0.01625980437882562975715546, 0.01599965869724394401313881, 0.01574770606433893015574400, 0.01550356543933893015574400, 0.01526687904880638577704578, 0.01503731063741979257227076, 0.01481454387422086185273411, 0.01459828089844231513993134, 0.01438824099085987447620523, 0.01418415935820681325171544, 0.01398578601958352422176106, 0.01379288478501562298719316, 0.01360523231738567365335942, 0.01342261726990576130858221, 0.01324483949212798353080444, 0.01307170929822216635628920, 0.01290304679189732236910755, 0.01273868124291638877278934, 0.01257845051066194236996928, 0.01242220051066194236996928, 0.01226978472038606978956995, 0.01212106372098095378719041, 0.01197590477193174490346273, 0.01183418141592267460867815, 0.01169577311142440471248438, 0.01156056489076458859566448, 0.01142844704164317229232189, 0.01129931481023821361463594, 0.01117306812421372175754719, 0.01104961133409026496742374, 0.01092885297157366069257770, 0.01081070552355853781923177, 0.01069508522063334415522437, 0.01058191183901270133041676, 0.01047110851491297833872701, 0.01036260157046853389428257, 0.01025632035036012704977199, /* ... */ 0.01015219706839427948625679, /* psi(1,99) */ 0.01005016666333357139524567 /* psi(1,100) */ }; inline static double lngamma_1_pade (const double eps) { /* Use (2,2) Pade for Log[Gamma[1+eps]]/eps * plus a correction series. */ const double n1 = -1.0017419282349508699871138440; const double n2 = 1.7364839209922879823280541733; const double d1 = 1.2433006018858751556055436011; const double d2 = 5.0456274100274010152489597514; const double num = (eps + n1) * (eps + n2); const double den = (eps + d1) * (eps + d2); const double pade = 2.0816265188662692474880210318 * num / den; const double c0 = 0.004785324257581753; const double c1 = -0.01192457083645441; const double c2 = 0.01931961413960498; const double c3 = -0.02594027398725020; const double c4 = 0.03141928755021455; const double eps5 = eps*eps*eps*eps*eps; const double corr = eps5 * (c0 + eps*(c1 + eps*(c2 + eps*(c3 + c4*eps)))); return eps * (pade + corr); } inline static double lngamma_2_pade (const double eps) { /* Use (2,2) Pade for Log[Gamma[2+eps]]/eps * plus a correction series. */ const double n1 = 1.000895834786669227164446568; const double n2 = 4.209376735287755081642901277; const double d1 = 2.618851904903217274682578255; const double d2 = 10.85766559900983515322922936; const double num = (eps + n1) * (eps + n2); const double den = (eps + d1) * (eps + d2); const double pade = 2.85337998765781918463568869 * num/den; const double c0 = 0.0001139406357036744; const double c1 = -0.0001365435269792533; const double c2 = 0.0001067287169183665; const double c3 = -0.0000693271800931282; const double c4 = 0.0000407220927867950; const double eps5 = eps*eps*eps*eps*eps; const double corr = eps5 * (c0 + eps*(c1 + eps*(c2 + eps*(c3 + c4*eps)))); return eps * (pade + corr); } /* coefficients for gamma=7, kmax=8 Lanczos method */ static double lanczos_7_c[9] = { 0.99999999999980993227684700473478, 676.520368121885098567009190444019, -1259.13921672240287047156078755283, 771.3234287776530788486528258894, -176.61502916214059906584551354, 12.507343278686904814458936853, -0.13857109526572011689554707, 9.984369578019570859563e-6, 1.50563273514931155834e-7 }; /* Lanczos method for real x > 0; * gamma=7, truncated at 1/(z+8) * [J. SIAM Numer. Anal, Ser. B, 1 (1964) 86] */ static double lngamma_lanczos (double x) { int k; double Ag; double term1, term2; x -= 1.0; /* Lanczos writes z! instead of Gamma(z) */ Ag = lanczos_7_c[0]; for (k=1; k<=8; k++) { Ag += lanczos_7_c[k]/(x+k); } /* (x+0.5)*log(x+7.5) - (x+7.5) + LogRootTwoPi_ + log(Ag(x)) */ term1 = (x+0.5)*log((x+7.5)/M_E); term2 = LogRootTwoPi_ + log(Ag); return term1 + (term2 - 7.0); } /* x = eps near zero * gives double-precision for |eps| < 0.02 */ static double lngamma_sgn_0 (double eps, double * sgn) { /* calculate series for g(eps) = Gamma(eps) eps - 1/(1+eps) - eps/2 */ const double c1 = -0.07721566490153286061; const double c2 = -0.01094400467202744461; const double c3 = 0.09252092391911371098; const double c4 = -0.01827191316559981266; const double c5 = 0.01800493109685479790; const double c6 = -0.00685088537872380685; const double c7 = 0.00399823955756846603; const double c8 = -0.00189430621687107802; const double c9 = 0.00097473237804513221; const double c10 = -0.00048434392722255893; const double g6 = c6+eps*(c7+eps*(c8 + eps*(c9 + eps*c10))); const double g = eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*g6))))); /* calculate Gamma(eps) eps, a positive quantity */ const double gee = g + 1.0/(1.0+eps) + 0.5*eps; *sgn = GSL_SIGN(eps); return log(gee/fabs(eps)); } static double gsl_sf_lnfact (const unsigned int n) { if (n <= GSL_SF_FACT_NMAX) { return log(fact_table[n].f); } else { return gsl_sf_lngamma(n+1.0); } } static double gsl_sf_psi_int (const int n) { if (n <= 0) { fprintf(stderr,"DOMAIN ERROR"); abort(); } else if (n <= PSI_TABLE_NMAX) { return psi_table[n]; } else { /* Abramowitz+Stegun 6.3.18 */ const double c2 = -1.0/12.0; const double c3 = 1.0/120.0; const double c4 = -1.0/252.0; const double c5 = 1.0/240.0; const double ni2 = (1.0/n)*(1.0/n); const double ser = ni2 * (c2 + ni2 * (c3 + ni2 * (c4 + ni2*c5))); return log(n) - 0.5/n + ser; } } static double gsl_sf_psi_1_int (const int n) { if (n <= 0) { fprintf(stderr,"DOMAIN ERROR"); abort(); } else if (n <= PSI_1_TABLE_NMAX) { return psi_1_table[n]; } else { /* Abramowitz+Stegun 6.4.12 * double-precision for n > 100 */ const double c0 = -1.0/30.0; const double c1 = 1.0/42.0; const double c2 = -1.0/30.0; const double ni2 = (1.0/n)*(1.0/n); const double ser = ni2*ni2 * (c0 + ni2*(c1 + c2*ni2)); return (1.0 + 0.5/n + 1.0/(6.0*n*n) + ser) / n; } } static double apsics_data[16] = { -.0204749044678185, -.0101801271534859, .0000559718725387, -.0000012917176570, .0000000572858606, -.0000000038213539, .0000000003397434, -.0000000000374838, .0000000000048990, -.0000000000007344, .0000000000001233, -.0000000000000228, .0000000000000045, -.0000000000000009, .0000000000000002, -.0000000000000000 }; static cheb_series apsi_cs = { apsics_data, 15, -1, 1, 9 }; static double psics_data[23] = { -.038057080835217922, .491415393029387130, -.056815747821244730, .008357821225914313, -.001333232857994342, .000220313287069308, -.000037040238178456, .000006283793654854, -.000001071263908506, .000000183128394654, -.000000031353509361, .000000005372808776, -.000000000921168141, .000000000157981265, -.000000000027098646, .000000000004648722, -.000000000000797527, .000000000000136827, -.000000000000023475, .000000000000004027, -.000000000000000691, .000000000000000118, -.000000000000000020 }; static cheb_series psi_cs = { psics_data, 22, -1, 1, 17 }; /* digamma for x both positive and negative; we do both * cases here because of the way we use even/odd parts * of the function */ static double gsl_sf_psi (const double x) { const double y = fabs(x); double result_c; if (x == 0.0 || x == -1.0 || x == -2.0) { fprintf(stderr,"DOMAIN ERROR"); abort(); } else if (y >= 2.0) { const double t = 8.0/(y*y)-1.0; result_c = cheb_eval(&apsi_cs, t); if (x < 0.0) { const double s = sin(M_PI*x); const double c = cos(M_PI*x); if (fabs(s) < 2.0*GSL_SQRT_DBL_MIN) { fprintf(stderr,"DOMAIN ERROR"); abort(); } else { return log(y) - 0.5/x + result_c - M_PI * c/s; } } else { return log(y) - 0.5/x + result_c; } } else { /* -2 < x < 2 */ if (x < -1.0) { /* x = -2 + v */ const double v = x + 2.0; const double t1 = 1.0/x; const double t2 = 1.0/(x+1.0); const double t3 = 1.0/v; return -(t1 + t2 + t3) + cheb_eval(&psi_cs, 2.0*v-1.0); } else if (x < 0.0) { /* x = -1 + v */ const double v = x + 1.0; const double t1 = 1.0/x; const double t2 = 1.0/v; return -(t1 + t2) + cheb_eval(&psi_cs, 2.0*v-1.0); } else if(x < 1.0) { /* x = v */ const double t1 = 1.0/x; return -t1 + cheb_eval(&psi_cs, 2.0*x-1.0); } else { /* x = 1 + v */ const double v = x - 1.0; return cheb_eval(&psi_cs, 2.0*v-1.0); } } } /* coefficients for Maclaurin summation in hzeta() * B_{2j}/(2j)! */ static double hzeta_c[15] = { 1.00000000000000000000000000000, 0.083333333333333333333333333333, -0.00138888888888888888888888888889, 0.000033068783068783068783068783069, -8.2671957671957671957671957672e-07, 2.0876756987868098979210090321e-08, -5.2841901386874931848476822022e-10, 1.3382536530684678832826980975e-11, -3.3896802963225828668301953912e-13, 8.5860620562778445641359054504e-15, -2.1748686985580618730415164239e-16, 5.5090028283602295152026526089e-18, -1.3954464685812523340707686264e-19, 3.5347070396294674716932299778e-21, -8.9535174270375468504026113181e-23 }; static double gsl_sf_hzeta (const double s, const double q) { if (s <= 1.0 || q <= 0.0) { fprintf(stderr,"DOMAIN_ERROR\n"); abort(); } else { const double max_bits = 54.0; const double ln_term0 = -s * log(q); if (ln_term0 < GSL_LOG_DBL_MIN + 1.0) { return UNDERFLOW_ERROR; } else if (ln_term0 > GSL_LOG_DBL_MAX - 1.0) { return OVERFLOW_ERROR; } else if ((s > max_bits && q < 1.0) || (s > 0.5*max_bits && q < 0.25)) { return pow(q, -s); } else if (s > 0.5*max_bits && q < 1.0) { const double p1 = pow(q, -s); const double p2 = pow(q/(1.0+q), s); const double p3 = pow(q/(2.0+q), s); return p1 * (1.0 + p2 + p3); } else { /* Euler-Maclaurin summation formula * [Moshier, p. 400, with several typo corrections] */ const int jmax = 12; const int kmax = 10; int j, k; const double pmax = pow(kmax + q, -s); double scp = s; double pcp = pmax / (kmax + q); double ans = pmax*((kmax+q)/(s-1.0) + 0.5); for (k=0; k 0.5*GSL_LOG_DBL_MIN) && (ay < 0.8*GSL_SQRT_DBL_MAX && ay > 1.2*GSL_SQRT_DBL_MIN) ) { double ex = exp(x); return y * ex; } else { const double ly = log(ay); const double lnr = x + ly; if (lnr > GSL_LOG_DBL_MAX - 0.01) { return OVERFLOW_ERROR; } else if (lnr < GSL_LOG_DBL_MIN + 0.01) { return UNDERFLOW_ERROR; } else { const double sy = GSL_SIGN(y); const double M = floor(x); const double N = floor(ly); const double a = x - M; const double b = ly - N; const double eMN = exp(M+N); const double eab = exp(a+b); return sy * eMN * eab; } } } /* generic polygamma; assumes n >= 0 and x > 0 */ static double psi_n_xg0 (const int n, const double x) { double result; if (n == 0) { return gsl_sf_psi(x); } else { /* Abramowitz + Stegun 6.4.10 */ double hzeta = gsl_sf_hzeta(n+1.0, x); double ln_nf = gsl_sf_lnfact((unsigned int) n); result = gsl_sf_exp_mult(ln_nf, hzeta); if (GSL_IS_EVEN(n)) { result = -result; } return result; } } static double gsl_sf_psi_1 (const double x) { if (x == 0.0 || x == -1.0 || x == -2.0) { fprintf(stderr,"DOMAIN ERROR\n"); abort(); } else if (x > 0.0) { return psi_n_xg0(1, x); } else if (x > -5.0) { /* Abramowitz + Stegun 6.4.6 */ int M = -floor(x); double fx = x + M; double sum = 0.0; int m; if (fx == 0.0) { fprintf(stderr,"DOMAIN ERROR\n"); abort(); } for (m = 0; m < M; ++m) { sum += 1.0/((x+m)*(x+m)); } return sum + psi_n_xg0(1, fx); } else { /* Abramowitz + Stegun 6.4.7 */ const double sin_px = sin(M_PI * x); const double d = M_PI*M_PI/(sin_px*sin_px); return d - psi_n_xg0(1, 1.0-x); } } static double gsl_sf_psi_n (const int n, const double x) { if (n == 0) { return gsl_sf_psi(x); } else if (n == 1) { return gsl_sf_psi_1(x); } else if (n < 0 || x <= 0.0) { fprintf(stderr,"DOMAIN ERROR"); abort(); } else { double hzeta = gsl_sf_hzeta(n+1.0, x); double ln_nf = gsl_sf_lnfact((unsigned int) n); double result = gsl_sf_exp_mult(ln_nf, hzeta); if (GSL_IS_EVEN(n)) { result = -result; } return result; } } /* x near a negative integer * Calculates sign as well as log(|gamma(x)|). * x = -N + eps * assumes N >= 1 */ static double lngamma_sgn_sing (int N, double eps, double * sgn) { if (eps == 0.0) { fprintf(stderr,"error GSL_EDOM\n"); *sgn = 0.0; return 0.0; } else if (N == 1) { /* calculate series for * g = eps gamma(-1+eps) + 1 + eps/2 (1+3eps)/(1-eps^2) * double-precision for |eps| < 0.02 */ const double c0 = 0.07721566490153286061; const double c1 = 0.08815966957356030521; const double c2 = -0.00436125434555340577; const double c3 = 0.01391065882004640689; const double c4 = -0.00409427227680839100; const double c5 = 0.00275661310191541584; const double c6 = -0.00124162645565305019; const double c7 = 0.00065267976121802783; const double c8 = -0.00032205261682710437; const double c9 = 0.00016229131039545456; const double g5 = c5 + eps*(c6 + eps*(c7 + eps*(c8 + eps*c9))); const double g = eps*(c0 + eps*(c1 + eps*(c2 + eps*(c3 + eps*(c4 + eps*g5))))); /* calculate eps gamma(-1+eps), a negative quantity */ const double gam_e = g - 1.0 - 0.5*eps*(1.0+3.0*eps)/(1.0 - eps*eps); *sgn = ( eps > 0.0 ? -1.0 : 1.0 ); return log(fabs(gam_e)/fabs(eps)); } else { double g; /* series for sin(Pi(N+1-eps))/(Pi eps) modulo the sign * double-precision for |eps| < 0.02 */ const double cs1 = -1.6449340668482264365; const double cs2 = 0.8117424252833536436; const double cs3 = -0.1907518241220842137; const double cs4 = 0.0261478478176548005; const double cs5 = -0.0023460810354558236; const double e2 = eps*eps; const double sin_ser = 1.0 + e2*(cs1+e2*(cs2+e2*(cs3+e2*(cs4+e2*cs5)))); /* calculate series for ln(gamma(1+N-eps)) * double-precision for |eps| < 0.02 */ double aeps = fabs(eps); double c1, c2, c3, c4, c5, c6, c7; double lng_ser; double c0; double psi_0, psi_1; double psi_2 = 0.0, psi_3 = 0.0, psi_4 = 0.0, psi_5 = 0.0, psi_6 = 0.0; c0 = gsl_sf_lnfact(N); psi_0 = gsl_sf_psi_int(N+1); psi_1 = gsl_sf_psi_1_int(N+1); if(aeps > 0.00001) psi_2 = gsl_sf_psi_n(2, N+1.0); if(aeps > 0.0002) psi_3 = gsl_sf_psi_n(3, N+1.0); if(aeps > 0.001) psi_4 = gsl_sf_psi_n(4, N+1.0); if(aeps > 0.005) psi_5 = gsl_sf_psi_n(5, N+1.0); if(aeps > 0.01) psi_6 = gsl_sf_psi_n(6, N+1.0); c1 = psi_0; c2 = psi_1/2.0; c3 = psi_2/6.0; c4 = psi_3/24.0; c5 = psi_4/120.0; c6 = psi_5/720.0; c7 = psi_6/5040.0; lng_ser = c0-eps*(c1-eps*(c2-eps*(c3-eps*(c4-eps*(c5-eps*(c6-eps*c7)))))); /* calculate * g = ln(|eps gamma(-N+eps)|) * = -ln(gamma(1+N-eps)) + ln(|eps Pi/sin(Pi(N+1+eps))|) */ g = -lng_ser - log(sin_ser); *sgn = ( GSL_IS_ODD(N) ? -1.0 : 1.0 ) * ( eps > 0.0 ? 1.0 : -1.0 ); return g - log(fabs(eps)); } } /* Needs to be extern, because calls lngamma_sgn_sing -> gsl_sf_lnfact -> gsl_sf_lngamma */ double gsl_sf_lngamma (double x) { if (fabs(x - 1.0) < 0.01) { /* Note that we must amplify the errors * from the Pade evaluations because of * the way we must pass the argument, i.e. * writing (1-x) is a loss of precision * when x is near 1. */ return lngamma_1_pade(x - 1.0); } else if (fabs(x - 2.0) < 0.01) { return lngamma_2_pade(x - 2.0); } else if (x >= 0.5) { return lngamma_lanczos(x); } else if (x == 0.0) { fprintf(stderr,"x == 0.0\n"); abort(); } else if (fabs(x) < 0.02) { double sgn; return lngamma_sgn_0(x, &sgn); } else if(x > -0.5/(GSL_DBL_EPSILON*M_PI)) { /* Try to extract a fractional * part from x. */ double z = 1.0 - x; double s = sin(M_PI*z); double as = fabs(s); if (s == 0.0) { fprintf(stderr,"s == 0.0\n"); abort(); } else if (as < M_PI*0.015) { /* x is near a negative integer, -N */ if (x < INT_MIN + 2.0) { /* GSL_ERROR ("error", GSL_EROUND); */ return 0.0; } else { int N = -(int)(x - 0.5); double eps = x + N; double sgn; return lngamma_sgn_sing(N, eps, &sgn); } } else { double lg_z; lg_z = lngamma_lanczos(z); return M_LNPI - (log(as) + lg_z); } } else { /* |x| was too large to extract any fractional part */ /* GSL_ERROR ("error", GSL_EROUND); */ return 0.0; } } /* Chebyshev coefficients for Gamma*(3/4(t+1)+1/2), -1 10.0 */ static double gammastar_ser (const double x) { /* Use the Stirling series for the correction to Log(Gamma(x)), * which is better behaved and easier to compute than the * regular Stirling series for Gamma(x). */ const double y = 1.0/(x*x); const double c0 = 1.0/12.0; const double c1 = -1.0/360.0; const double c2 = 1.0/1260.0; const double c3 = -1.0/1680.0; const double c4 = 1.0/1188.0; const double c5 = -691.0/360360.0; const double c6 = 1.0/156.0; const double c7 = -3617.0/122400.0; const double ser = c0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*(c5 + y*(c6 + y*c7)))))); return exp(ser/x); } static double gsl_sf_gammastar (const double x) { if (x <= 0.0) { fprintf(stderr,"Domain error\n"); abort(); } else if (x < 0.5) { double lg = gsl_sf_lngamma(x); const double lx = log(x); const double c = 0.5*(M_LN2+M_LNPI); const double lnr_val = lg - (x-0.5)*lx + x - c; return exp(lnr_val); } else if (x < 2.0) { const double t = 4.0/3.0*(x-0.5) - 1.0; return cheb_eval(&gstar_a_cs, t); } else if (x < 10.0) { const double t = 0.25*(x-2.0) - 1.0; double c = cheb_eval(&gstar_b_cs, t); return c/(x*x) + 1.0 + 1.0/(12.0*x); } else if(x < 1.0/GSL_ROOT4_DBL_EPSILON) { return gammastar_ser(x); } else if(x < 1.0/GSL_DBL_EPSILON) { /* Use Stirling formula for Gamma(x). */ const double xi = 1.0/x; return 1.0 + xi/12.0*(1.0 + xi/24.0*(1.0 - xi*(139.0/180.0 + 571.0/8640.0*xi))); } else { return 1.0; } } /* Chebyshev expansion for log(1 + x(t))/x(t) * * x(t) = (4t-1)/(2(4-t)) * t(x) = (8x+1)/(2(x+2)) * -1/2 < x < 1/2 * -1 < t < 1 */ static double lopx_data[21] = { 2.16647910664395270521272590407, -0.28565398551049742084877469679, 0.01517767255690553732382488171, -0.00200215904941415466274422081, 0.00019211375164056698287947962, -0.00002553258886105542567601400, 2.9004512660400621301999384544e-06, -3.8873813517057343800270917900e-07, 4.7743678729400456026672697926e-08, -6.4501969776090319441714445454e-09, 8.2751976628812389601561347296e-10, -1.1260499376492049411710290413e-10, 1.4844576692270934446023686322e-11, -2.0328515972462118942821556033e-12, 2.7291231220549214896095654769e-13, -3.7581977830387938294437434651e-14, 5.1107345870861673561462339876e-15, -7.0722150011433276578323272272e-16, 9.7089758328248469219003866867e-17, -1.3492637457521938883731579510e-17, 1.8657327910677296608121390705e-18 }; static cheb_series lopx_cs = { lopx_data, 20, -1, 1, 10 }; static double gsl_sf_log_1plusx (const double x) { if (x <= -1.0) { fprintf(stderr,"DOMAIN ERROR\n"); abort(); } else if (fabs(x) < GSL_ROOT6_DBL_EPSILON) { const double c1 = -0.5; const double c2 = 1.0/3.0; const double c3 = -1.0/4.0; const double c4 = 1.0/5.0; const double c5 = -1.0/6.0; const double c6 = 1.0/7.0; const double c7 = -1.0/8.0; const double c8 = 1.0/9.0; const double c9 = -1.0/10.0; const double t = c5 + x*(c6 + x*(c7 + x*(c8 + x*c9))); return x * (1.0 + x*(c1 + x*(c2 + x*(c3 + x*(c4 + x*t))))); } else if (fabs(x) < 0.5) { double t = 0.5*(8.0*x + 1.0)/(x+2.0); return x * cheb_eval(&lopx_cs, t); } else { return log(1.0 + x); } } static double gsl_sf_lnbeta (const double x, const double y) { if (x <= 0.0 || y <= 0.0) { fprintf(stderr,"x <= 0.0 || y <= 0.0"); abort(); } else { const double max = GSL_MAX(x,y); const double min = GSL_MIN(x,y); const double rat = min/max; if (rat < 0.2) { /* min << max, so be careful * with the subtraction */ double lnpre_val; double lnpow_val; double t1, t2, t3; double lnopr; double gsx, gsy, gsxy; gsx = gsl_sf_gammastar(x); gsy = gsl_sf_gammastar(y); gsxy = gsl_sf_gammastar(x+y); lnopr = gsl_sf_log_1plusx(rat); lnpre_val = log(gsx * gsy / gsxy * M_SQRT2*M_SQRTPI); t1 = min * log(rat); t2 = 0.5 * log(min); t3 = (x+y-0.5) * lnopr; lnpow_val = t1 - t2 - t3; return lnpre_val + lnpow_val; } else { double lgx = gsl_sf_lngamma(x); double lgy = gsl_sf_lngamma(y); double lgxy = gsl_sf_lngamma(x+y); return lgx + lgy - lgxy; } } } /* Author: G. Jungman */ /* Modified for cdfs by Brian Gough, June 2003 */ static double beta_cont_frac (const double a, const double b, const double x, const double epsabs) { const unsigned int max_iter = 512; /* control iterations */ const double cutoff = 2.0 * GSL_DBL_MIN; /* control the zero cutoff */ unsigned int iter_count = 0; double cf; /* standard initialization for continued fraction */ double num_term = 1.0; double den_term = 1.0 - (a + b) * x / (a + 1.0); if (fabs(den_term) < cutoff) { /* den_term = GSL_NAN; */ fprintf(stderr,"NAN\n"); abort(); } den_term = 1.0 / den_term; cf = den_term; while (iter_count < max_iter) { const int k = iter_count + 1; double coeff = k * (b - k) * x / (((a - 1.0) + 2 * k) * (a + 2 * k)); double delta_frac; /* first step */ den_term = 1.0 + coeff * den_term; num_term = 1.0 + coeff / num_term; if (fabs(den_term) < cutoff) { /* den_term = GSL_NAN; */ fprintf(stderr,"NAN\n"); abort(); } if (fabs(num_term) < cutoff) { /* num_term = GSL_NAN; */ fprintf(stderr,"NAN\n"); abort(); } den_term = 1.0 / den_term; delta_frac = den_term * num_term; cf *= delta_frac; coeff = -(a + k) * (a + b + k) * x / ((a + 2 * k) * (a + 2 * k + 1.0)); /* second step */ den_term = 1.0 + coeff * den_term; num_term = 1.0 + coeff / num_term; if (fabs(den_term) < cutoff) { /* den_term = GSL_NAN; */ fprintf(stderr,"NAN\n"); abort(); } if (fabs(num_term) < cutoff) { /* num_term = GSL_NAN; */ fprintf(stderr,"NAN\n"); abort(); } den_term = 1.0 / den_term; delta_frac = den_term * num_term; cf *= delta_frac; if (fabs(delta_frac - 1.0) < 2.0 * GSL_DBL_EPSILON) { break; } if (cf * fabs(delta_frac - 1.0) < epsabs) { break; } ++iter_count; } if (iter_count >= max_iter) { /* return GSL_NAN; */ fprintf(stderr,"NAN\n"); abort(); } return cf; } /* The function beta_inc_AXPY(A,Y,a,b,x) computes A * beta_inc(a,b,x) + Y taking account of possible cancellations when using the hypergeometric transformation beta_inc(a,b,x)=1-beta_inc(b,a,1-x). It also adjusts the accuracy of beta_inc() to fit the overall absolute error when A*beta_inc is added to Y. (e.g. if Y >> A*beta_inc then the accuracy of beta_inc can be reduced) */ static double beta_inc_AXPY (const double A, const double Y, const double a, const double b, const double x) { if (x == 0.0) { return A * 0 + Y; } else if (x == 1.0) { return A * 1 + Y; } else { double ln_beta = gsl_sf_lnbeta(a,b); double ln_pre = -ln_beta + a * log (x) + b * log1p(-x); double prefactor = exp(ln_pre); if (x < (a + 1.0) / (a + b + 2.0)) { /* Apply continued fraction directly. */ double epsabs = fabs(Y / (A * prefactor / a)) * GSL_DBL_EPSILON; double cf = beta_cont_frac(a, b, x, epsabs); return A * (prefactor * cf / a) + Y; } else { /* Apply continued fraction after hypergeometric transformation. */ double epsabs = fabs ((A + Y) / (A * prefactor / b)) * GSL_DBL_EPSILON; double cf = beta_cont_frac(b, a, 1.0 - x, epsabs); double term = prefactor * cf / b; if (A == -Y) { return -A * term; } else { return A * (1 - term) + Y; } } } } #if 0 static double gsl_cdf_beta_P (const double x, const double a, const double b) { double P; if (x <= 0.0) { return 0.0; } if (x >= 1.0) { return 1.0; } P = beta_inc_AXPY (1.0, 0.0, a, b, x); return P; } #endif static double gsl_cdf_beta_Q (const double x, const double a, const double b) { double Q; if (x >= 1.0) { return 0.0; } if (x <= 0.0) { return 1.0; } Q = beta_inc_AXPY (-1.0, 1.0, a, b, x); return Q; } static double gsl_cdf_binomial_P (const unsigned int k, const double p, const unsigned int n) { double P; double a; double b; if (p > 1.0 || p < 0.0) { fprintf(stderr,"p < 0 or p > 1\n"); abort(); } if (k >= n) { P = 1.0; } else { a = (double) k + 1.0; b = (double) n - k; P = gsl_cdf_beta_Q (p, a, b); } return P; } #if 0 static double gsl_cdf_binomial_Q (const unsigned int k, const double p, const unsigned int n) { double Q; double a; double b; if (p > 1.0 || p < 0.0) { fprintf(stderr,"p < 0 or p > 1\n"); abort(); } if (k >= n) { Q = 0.0; } else { a = (double) k + 1.0; b = (double) n - k; Q = gsl_cdf_beta_P (p, a, b); } return Q; } #endif /* Same as R function */ double Pbinom (int k, int n, double theta) { return gsl_cdf_binomial_P(k,theta,n); }