lambda.c

/*
 * Copyright (c) 2002, 2017 Jens Keiner, Stefan Kunis, Daniel Potts
 *
 * This program is free software; you can redistribute it and/or modify it under
 * the terms of the GNU General Public License as published by the Free Software
 * Foundation; either version 2 of the License, or (at your option) any later
 * version.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
 * details.
 *
 * You should have received a copy of the GNU General Public License along with
 * this program; if not, write to the Free Software Foundation, Inc., 51
 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */
 
#include "infft.h"
 
/* Coefficients for Lanzcos's approximation to the Gamma function. Can be
 * regenerated with Mathematica from file lambda.nb. */
 
 #if MANT_DIG == 113
  // IEEE 754 quadruple precision, 128 bits.
  #define N 24
  static const R num[24] =
  {
    K(3.035162425359883494754028782232869726547E21),
    K(3.4967568944064301036001605717507506346E21),
    K(1.9266526566893208886540195401514595829E21),
    K(6.755170664882727663160830237424406199E20),
    K(1.691728531049187527800862627495648317E20),
    K(3.21979351672256057856444116302160246E19),
    K(4.8378495427140832493758744745481812E18),
    K(5.8843103809049324230843820398664955E17),
    K(5.893958514163405862064178891925630E16),
    K(4.919561837722192829918665308020810E15),
    K(3.449165802442404074427531228315120E14),
    K(2.041330296068782505988459692384726E13),
    K(1.022234822943784007524609706893119E12),
    K(4.33137871919821354846952908076307E10),
    K(1.54921950559667418528481770869280E9),
    K(4.6544421199876191938054157935810E7),
    K(1.16527806807504975090675074910053E6),
    K(24024.759267256769471083727721827),
    K(400.96500811342195582435806376976),
    K(5.2829901565447826961703902917085),
    K(0.05289990244125101024092566765994),
    K(0.0003783467106547406854542665695934),
    K(1.7219414217921113919596660801124E-6),
    K(3.747999317071488557713812635427084359354E-9)
  };
  static const R g = K(20.32098218798637390136718750000000000000);
#elif MANT_DIG == 64
// Intel double extended, 80 bits.
  #define N 17
  static const R num[17] =
  {
    K(2.715894658327717377557655133124376674911E12),
    K(3.59017952609791210503852552872112955043E12),
    K(2.22396659973781496931212735323581871017E12),
    K(8.5694083451895624818099258668254858834E11),
    K(2.2988587166874907293359744645339939547E11),
    K(4.552617168754610815813502794395753410E10),
    K(6.884887713165178784550917647709216425E9),
    K(8.11048596140753186476028245385237278E8),
    K(7.52139159654082231449961362311950170E7),
    K(5.50924541722426515169752795795495283E6),
    K(317673.536843541912671493184218236957),
    K(14268.2798984503552014701437332033752),
    K(489.361872040326367021390908360178781),
    K(12.3894133003845444929588321786545861),
    K(0.218362738950461496394157450728168315),
    K(0.00239374952205844918669062799606398310),
    K(0.00001229541408909435212800785616808830746135)
  };
  static const R g = K(12.22522273659706115722656250000000000000);
#elif MANT_DIG == 53
  // IEEE 754 double precision, 64 bits.
  #define N 13
  static const R num[13] =
  {
    K(5.690652191347156388090791033559122686859E7),
    K(1.037940431163445451906271053616070238554E8),
    K(8.63631312881385914554692728897786842234E7),
    K(4.33388893246761383477372374059053331609E7),
    K(1.46055780876850680841416998279135921857E7),
    K(3.48171215498064590882071018964774556468E6),
    K(601859.61716810987866702265336993523025),
    K(75999.293040145426498753034435989091371),
    K(6955.9996025153761403563101155151989875),
    K(449.944556906316811944685860765098840962),
    K(19.5199278824761748284786096623565213621),
    K(0.509841665565667618812517864480469450999),
    K(0.006061842346248906525783753964555936883222)
  };
  static const R g = K(6.024680040776729583740234375000000000000);
#elif MANT_DIG == 24
  // IEEE 754 single precision, 32 bits.
  #define N 6
  static const R num[6] =
  {
    K(14.02614328749964766195705772850038393570),
    K(43.74732405540314316089531289293124360129),
    K(50.59547402616588964511581430025589038612),
    K(26.90456680562548195593733429204228910299),
    K(6.595765571169314946316366571954421695196),
    K(0.6007854010515290065101128585795542383721)
  };
  static const R g = K(1.428456135094165802001953125000000000000);
#else
  // Unknown floating-point type.
  // Assume IEEE 754 double precision, 64 bits.
  #define N 13
  static const R num[13] =
  {
    K(5.690652191347156388090791033559122686859E7),
    K(1.037940431163445451906271053616070238554E8),
    K(8.63631312881385914554692728897786842234E7),
    K(4.33388893246761383477372374059053331609E7),
    K(1.46055780876850680841416998279135921857E7),
    K(3.48171215498064590882071018964774556468E6),
    K(601859.61716810987866702265336993523025),
    K(75999.293040145426498753034435989091371),
    K(6955.9996025153761403563101155151989875),
    K(449.944556906316811944685860765098840962),
    K(19.5199278824761748284786096623565213621),
    K(0.509841665565667618812517864480469450999),
    K(0.006061842346248906525783753964555936883222)
  };
  static const R g = K(6.024680040776729583740234375000000000000);
#endif
 
static inline R evaluate_rational(const R z_)
{
  R z = z_, s1, s2;
  INT i;
 
  if (z <= K(1.0))
  {
    s1 = num[N - 1];
    s2 = K(1.0);
    for (i = N - 2; i >= 0; --i)
    {
      s1 *= z;
      s2 *= z + (R)(i);
      s1 += num[i];
    }
  }
  else
  {
    z = K(1.0)/z;
    s1 = num[0];
    s2 = K(1.0);
    for (i = 1; i < N; ++i)
    {
      s1 *= z;
      s2 *= K(1.0) + (R)(i-1) * z;
      s1 += num[i];
    }
  }
  return s1 / s2;
}
 
R Y(lambda)(const R z, const R eps)
{
  const R d = K(1.0) - eps, zpg = z + g, emh = eps - K(0.5);
  return EXP(-LOG1P(d / (zpg + emh)) * (z + emh)) *
    POW(KE / (zpg + K(0.5)),d) *
    (evaluate_rational(z + eps) / evaluate_rational(z + K(1.0)));
}
 
R Y(lambda2)(const R mu, const R nu)
{
  if (mu == K(0.0))
    return K(1.0);
  else if (nu == K(0.0))
    return K(1.0);
  else
    return
      SQRT(
        POW((mu + nu + g + K(0.5)) / (K(1.0) * (mu + g + K(0.5))), mu) *
        POW((mu + nu + g + K(0.5)) / (K(1.0) * (nu + g + K(0.5))), nu) *
        SQRT(KE * (mu + nu + g + K(0.5)) /
          ((mu + g + K(0.5)) * (nu + g + K(0.5)))) *
        (evaluate_rational(mu + nu + K(1.0)) /
          (evaluate_rational(mu + K(1.0)) * evaluate_rational(nu + K(1.0))))
      );
}