gamma.c

开放gsl矩阵运算

  { 160, 8.6522164519028464397823346347e142,        0 },
  { 161, 8.7731335619531641891199214170e143,        0 },
  { 162, 1.40165906520826112324473821082e145,       0 },
  { 163, 1.43002077059836576282654719098e146,       0 },
  { 164, 2.29872086694154824212137066574e147,       0 },
  { 165, 2.35953427148730350866380286512e148,       0 },
  { 166, 3.8158766391229700819214753051e149,        0 },
  { 167, 3.9404222333837968594685507847e150,        0 },
  { 168, 6.4106727537265897376280785126e151,        0 },
  { 169, 6.6593135744186166925018508262e152,        0 },
  { 170, 1.08981436813352025539677334714e154,       0 },
  { 171, 1.13874262122558345441781649128e155,       0 },
  { 172, 1.87448071318965483928245015709e156,       0 },
  { 173, 1.97002473472025937614282252992e157,       0 },
  { 174, 3.2615964409499994203514632733e158,        0 },
  { 175, 3.4475432857604539082499394274e159,        0 },
  { 176, 5.7404097360719989798185753611e160,        0 },
  { 177, 6.1021516157960034176023927864e161,        0 },
  { 178, 1.02179293302081581840770641427e163,       0 },
  { 179, 1.09228513922748461175082830877e164,       0 },
  { 180, 1.83922727943746847313387154568e165,       0 },
  { 181, 1.97703610200174714726899923887e166,       0 },
  { 182, 3.3473936485761926211036462131e167,        0 },
  { 183, 3.6179760666631972795022686071e168,        0 },
  { 184, 6.1592043133801944228307090322e169,        0 },
  { 185, 6.6932557233269149670791969232e170,        0 },
  { 186, 1.14561200228871616264651187999e172,       0 },
  { 187, 1.25163882026213309884380982464e173,       0 },
  { 188, 2.15375056430278638577544233437e174,       0 },
  { 189, 2.36559737029543155681480056857e175,       0 },
  { 190, 4.0921260721752941329733404353e176,        0 },
  { 191, 4.5182909772642742735162690860e177,        0 },
  { 192, 7.8568820585765647353088136358e178,        0 },
  { 193, 8.7203015861200493478863993359e179,        0 },
  { 194, 1.52423511936385355864990984535e181,       0 },
  { 195, 1.70045880929340962283784787050e182,       0 },
  { 196, 2.98750083395315297495382329688e183,       0 },
  { 197, 3.3499038543080169569905603049e184,        0 },
  { 198, 5.9152516512272428904085701278e185,        0 },
  { 199, 6.6663086700729537444112150067e186,        0 },
  { 200, 1.18305033024544857808171402556e188,       0 },
  { 201, 1.33992804268466370262665421635e189,       0 },
  { 202, 2.38976166709580612772506233164e190,       0 },
  { 203, 2.72005392664986731633210805920e191,       0 },
  { 204, 4.8751138008754445005591271565e192,        0 },
  { 205, 5.5761105496322279984808215214e193,        0 },
  { 206, 1.00427344298034156711518019425e195,       0 },
  { 207, 1.15425488377387119568553005492e196,       0 },
  { 208, 2.08888876139911045959957480403e197,       0 },
  { 209, 2.41239270708739079898275781478e198,       0 },
  { 210, 4.3866663989381319651591070885e199,        0 },
  { 211, 5.0901486119543945858536189892e200,        0 },
  { 212, 9.2997327657488397661373070276e201,        0 },
  { 213, 1.08420165434628604678682084470e203,       0 },
  { 214, 1.99014281187025170995338370390e204,       0 },
  { 215, 2.33103355684451500059166481610e205,       0 },
  { 216, 4.2987084736397436934993088004e206,        0 },
  { 217, 5.0583428183525975512839126509e207,        0 },
  { 218, 9.3711844725346412518284931849e208,        0 },
  { 219, 1.10777707721921886373117687056e210,       0 },
  { 220, 2.06166058395762107540226850068e211,       0 },
  { 221, 2.44818734065447368884590088393e212,       0 },
  { 222, 4.5768864963859187873930360715e213,        0 },
  { 223, 5.4594577696594763261263589712e214,        0 },
  { 224, 1.02522257519044580837604008002e216,       0 },
  { 225, 1.22837799817338217337843076851e217,       0 },
  { 226, 2.31700301993040752692985058084e218,       0 },
  { 227, 2.78841805585357753356903784452e219,       0 },
  { 228, 5.2827668854413291614000593243e220,        0 },
  { 229, 6.3854773479046925518730966640e221,        0 },
  { 230, 1.21503638365150570712201364459e223,       0 },
  { 231, 1.47504526736598397948268532937e224,       0 },
  { 232, 2.81888441007149324052307165546e225,       0 },
  { 233, 3.4368554729627426721946568174e226,        0 },
  { 234, 6.5961895195672941828239876738e227,        0 },
  { 235, 8.0766103614624452796574435210e228,        0 },
  { 236, 1.55670072661788142714646109101e230,       0 },
  { 237, 1.91415665566659953127881411447e231,       0 },
  { 238, 3.7049477293505577966085773966e232,        0 },
  { 239, 4.5748344070431728797563657336e233,        0 },
  { 240, 8.8918745504413387118605857518e234,        0 },
  { 241, 1.10253509209740466402128414180e236,       0 },
  { 242, 2.15183364120680396827026175195e237,       0 },
  { 243, 2.67916027379669333357172046456e238,       0 },
  { 244, 5.2504740845446016825794386748e239,        0 },
  { 245, 6.5639426708018986672507151382e240,        0 },
  { 246, 1.29161662479797201391454191399e242,       0 },
  { 247, 1.62129383968806897081092663913e243,       0 },
  { 248, 3.2032092294989705945080639467e244,        0 },
  { 249, 4.0370216608232917373192073314e245,        0 },
  { 250, 8.0080230737474264862701598667e246,        0 },
  { 251, 1.01329243686664622606712104019e248,       0 },
  { 252, 2.01802181458435147454008028642e249,       0 },
  { 253, 2.56362986527261495194981623168e250,       0 },
  { 254, 5.1257754090442527453318039275e251,        0 },
  { 255, 6.5372561564451681274720313908e252,        0 },
  { 256, 1.31219850471532870280494180544e254,       0 },
  { 257, 1.68007483220640820876031206743e255,       0 },
  { 258, 3.3854721421655480532367498580e256,        0 },
  { 259, 4.3513938154145972606892082546e257,        0 },
  { 260, 8.8022275696304249384155496309e258,        0 },
  { 261, 1.13571378582320988503988335446e260,       0 },
  { 262, 2.30618362324317133386487400329e261,       0 },
  { 263, 2.98692725671504199765489322224e262,       0 },
  { 264, 6.0883247653619723214032673687e263,        0 },
  { 265, 7.9153572302948612937854670389e264,        0 },
  { 266, 1.61949438758628463749326912007e266,       0 },
  { 267, 2.11340038048872796544071969939e267,       0 },
  { 268, 4.3402449587312428284819612418e268,        0 },
  { 269, 5.6850470235146782270355359914e269,        0 },
  { 270, 1.17186613885743556369012953528e271,       0 },
  { 271, 1.54064774337247779952663025366e272,       0 },
  { 272, 3.1874758976922247332371523360e273,        0 },
  { 273, 4.2059683394068643927077005925e274,        0 },
  { 274, 8.7336839596766957690697974006e275,        0 },
  { 275, 1.15664129333688770799461766294e277,       0 },
  { 276, 2.41049677287076803226326408256e278,       0 },
  { 277, 3.2038963825431789511450909263e279,        0 },
  { 278, 6.7011810285807351296918741495e280,        0 },
  { 279, 8.9388709072954692736948036845e281,        0 },
  { 280, 1.87633068800260583631372476186e283,       0 },
  { 281, 2.51182272495002686590823983534e284,       0 },
  { 282, 5.2912525401673484584047038284e285,        0 },
  { 283, 7.1084583116085760305203187340e286,        0 },
  { 284, 1.50271572140752696218693588728e288,       0 },
  { 285, 2.02591061880844416869829083919e289,       0 },
  { 286, 4.2977669632255271118546366376e290,        0 },
  { 287, 5.8143634759802347641640947085e291,        0 },
  { 288, 1.23775688540895180821413535163e293,       0 },
  { 289, 1.68035104455828784684342337075e294,       0 },
  { 290, 3.5894949676859602438209925197e295,        0 },
  { 291, 4.8898215396646176343143620089e296,        0 },
  { 292, 1.04813253056430039119572981576e298,       0 },
  { 293, 1.43271771112173296685410806860e299,       0 },
  { 294, 3.08150963985904315011544565835e300,       0 },
  { 295, 4.2265172478091122522196188024e301,        0 },
  { 296, 9.1212685339827677243417191487e302,        0 },
  { 297, 1.25527562259930633890922678431e304,       0 },
  /*
  { 298, 2.71813802312686478185383230631e305,       0 },
  { 299, 3.7532741115719259533385880851e306,        0 },
  { 300, 8.1544140693805943455614969189e307,  }
  */
};


/* Chebyshev coefficients for Gamma*(3/4(t+1)+1/2), -1<t<1
 */
static double gstar_a_data[30] = {
  2.16786447866463034423060819465,
 -0.05533249018745584258035832802,
  0.01800392431460719960888319748,
 -0.00580919269468937714480019814,
  0.00186523689488400339978881560,
 -0.00059746524113955531852595159,
  0.00019125169907783353925426722,
 -0.00006124996546944685735909697,
  0.00001963889633130842586440945,
 -6.3067741254637180272515795142e-06,
  2.0288698405861392526872789863e-06,
 -6.5384896660838465981983750582e-07,
  2.1108698058908865476480734911e-07,
 -6.8260714912274941677892994580e-08,
  2.2108560875880560555583978510e-08,
 -7.1710331930255456643627187187e-09,
  2.3290892983985406754602564745e-09,
 -7.5740371598505586754890405359e-10,
  2.4658267222594334398525312084e-10,
 -8.0362243171659883803428749516e-11,
  2.6215616826341594653521346229e-11,
 -8.5596155025948750540420068109e-12,
  2.7970831499487963614315315444e-12,
 -9.1471771211886202805502562414e-13,
  2.9934720198063397094916415927e-13,
 -9.8026575909753445931073620469e-14,
  3.2116773667767153777571410671e-14,
 -1.0518035333878147029650507254e-14,
  3.4144405720185253938994854173e-15,
 -1.0115153943081187052322643819e-15
};
static cheb_series gstar_a_cs = {
  gstar_a_data,
  29,
  -1, 1,
  17
};


/* Chebyshev coefficients for
 * x^2(Gamma*(x) - 1 - 1/(12x)), x = 4(t+1)+2, -1 < t < 1
 */
static double gstar_b_data[] = {
  0.0057502277273114339831606096782,
  0.0004496689534965685038254147807,
 -0.0001672763153188717308905047405,
  0.0000615137014913154794776670946,
 -0.0000223726551711525016380862195,
  8.0507405356647954540694800545e-06,
 -2.8671077107583395569766746448e-06,
  1.0106727053742747568362254106e-06,
 -3.5265558477595061262310873482e-07,
  1.2179216046419401193247254591e-07,
 -4.1619640180795366971160162267e-08,
  1.4066283500795206892487241294e-08,
 -4.6982570380537099016106141654e-09,
  1.5491248664620612686423108936e-09,
 -5.0340936319394885789686867772e-10,
  1.6084448673736032249959475006e-10,
 -5.0349733196835456497619787559e-11,
  1.5357154939762136997591808461e-11,
 -4.5233809655775649997667176224e-12,
  1.2664429179254447281068538964e-12,
 -3.2648287937449326771785041692e-13,
  7.1528272726086133795579071407e-14,
 -9.4831735252566034505739531258e-15,
 -2.3124001991413207293120906691e-15,
  2.8406613277170391482590129474e-15,
 -1.7245370321618816421281770927e-15,
  8.6507923128671112154695006592e-16,
 -3.9506563665427555895391869919e-16,
  1.6779342132074761078792361165e-16,
 -6.0483153034414765129837716260e-17
};
static cheb_series gstar_b_cs = {
  gstar_b_data,
  29,
  -1, 1,
  18
};


/* 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
};

/* complex version of Lanczos method; this is not safe for export
 * since it becomes bad in the left half-plane
 */
static
int
lngamma_lanczos_complex(double zr, double zi, gsl_sf_result * yr, gsl_sf_result * yi)
{
  int k;
  gsl_sf_result log1_r,    log1_i;
  gsl_sf_result logAg_r,   logAg_i;
  double Ag_r, Ag_i;
  double yi_tmp_val, yi_tmp_err;

  zr -= 1.0; /* Lanczos writes z! instead of Gamma(z) */

  Ag_r = lanczos_7_c[0];
  Ag_i = 0.0;
  for(k=1; k<=8; k++) {
    double R = zr + k;
    double I = zi;
    double a = lanczos_7_c[k] / (R*R + I*I);
    Ag_r +=  a * R;
    Ag_i -=  a * I;
  }

  gsl_sf_complex_log_e(zr + 7.5, zi, &log1_r,  &log1_i);
  gsl_sf_complex_log_e(Ag_r, Ag_i,   &logAg_r, &logAg_i);

  /* (z+0.5)*log(z+7.5) - (z+7.5) + LogRootTwoPi_ + log(Ag(z)) */
  yr->val = (zr+0.5)*log1_r.val - zi*log1_i.val - (zr+7.5) + LogRootTwoPi_ + logAg_r.val;
  yi->val = zi*log1_r.val + (zr+0.5)*log1_i.val - zi + logAg_i.val;
  yr->err = 4.0 * GSL_DBL_EPSILON * fabs(yr->val);
  yi->err = 4.0 * GSL_DBL_EPSILON * fabs(yi->val);
  yi_tmp_val = yi->val;
  yi_tmp_err = yi->err;
  gsl_sf_angle_restrict_symm_err_e(yi_tmp_val, yi);
  yi->err += yi_tmp_err;
  return GSL_SUCCESS;
}


/* Lanczos method for real x > 0;
 * gamma=7, truncated at 1/(z+8) 
 * [J. SIAM Numer. Anal, Ser. B, 1 (1964) 86]
 */
static
int
lngamma_lanczos(double x, gsl_sf_result * result)
{
  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);
  result->val  = term1 + (term2 - 7.0);
  result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(term1) + fabs(term2) + 7.0);
  result->err += GSL_DBL_EPSILON * fabs(result->val);

  return GSL_SUCCESS;
}


/* x = eps near zero
 * gives double-precision for |eps| < 0.02
 */
static
int
lngamma_sgn_0(double eps, gsl_sf_result * lng, 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;

  lng->val = log(gee/fabs(eps));
  lng->err = 4.0 * GSL_DBL_EPSILON * fabs(lng->val);
  *sgn = GSL_SIGN(eps);

  return GSL_SUCCESS;
}


/* x near a negative integer
 * Calculates sign as well as log(|gamma(x)|).
 * x = -N + eps
 * assumes N >= 1
 */
static
int
lngamma_sgn_sing(int N, double eps, gsl_sf_result * lng, double * sgn)
{
  if(eps == 0.0) {
    lng->val = 0.0;
    lng->err = 0.0;
    *sgn = 0.0;
    GSL_ERROR ("error", GSL_EDOM);
  }
  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);

    lng->val = log(fabs(gam_e)/fabs(eps));
    lng->err = 2.0 * GSL_DBL_EPSILON * fabs(lng->val);
    *sgn = ( eps > 0.0 ? -1.0 : 1.0 );
    return GSL_SUCCESS;
  }
  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;
    gsl_sf_result c0;
    gsl_sf_result psi_0;
    gsl_sf_result psi_1;
    gsl_sf_result psi_2;
    gsl_sf_result psi_3;
    gsl_sf_result psi_4;
    gsl_sf_result psi_5;
    gsl_sf_result psi_6;
    psi_2.val = 0.0;