/* qlgam
 *
 * Natural logarithm of gamma function
 *
 *
 *
 * SYNOPSIS:
 *
 * int qlgam( x, y );
 * QELT *x, *y;
 *
 * qlgam( x, y );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns the base e (2.718...) logarithm of the absolute
 * value of the gamma function of the argument.
 *
 */

/* qgamma
 *
 * Gamma function
 *
 *
 *
 * SYNOPSIS:
 *
 * int qgamma( x, y );
 * QELT *x, *y;
 *
 * qgamma( x, y );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns gamma function of the argument.
 *
 * qgamma(x) = exp(qlgam(x))
 *
 */
/*
Cephes Math Library Release 2.2:  July, 1992
Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
*/



#include "qhead.h"
extern QELT qpi[], qeul[];


/* constants for Stirling's formula */

#define NG 55
char *gstrs[] = {
 "7645992940484742892248134246724347500528752413412307906683593870759797606269585779977930217515",
  "1518","110","109",
"-3469342247847828789552088659323852541399766785760491146870005891371501266319724897592306597338057",
  "209191710","108","107",
"36373903172617414408151820151593427169231298640581690038930816378281879873386202346572901",
  "642","106","105",
"-319533631363830011287103352796174274671189606078272738327103470162849568365549721224053",
  "1590","104","103",
"3204019410860907078243020782116241775491817197152717450679002501086861530836678158791",
 "4326","102","101",
"-94598037819122125295227433069493721872702841533066936133385696204311395415197247711", 
  "33330","100","99",
"67908260672905495624051117546403605607342195728504487509073961249992947058239",
  "6","98","97",
"-211600449597266513097597728109824233673043954389060234150638733420050668349987259",
 "4501770","96","95",
"1220813806579744469607301679413201203958508415202696621436215105284649447",
  "6","94","93",
"-1295585948207537527989427828538576749659341483719435143023316326829946247",
"1410","92","91",
"1179057279021082799884123351249215083775254949669647116231545215727922535",
  "272118","90","89",
"-1311426488674017507995511424019311843345750275572028644296919890574047",
  "61410","88","87",
"660714619417678653573847847426261496277830686653388931761996983","6",
  "86","85",
"-2024576195935290360231131160111731009989917391198090877281083932477",
  "3404310","84","83",
"1677014149185145836823154509786269900207736027570253414881613","498",
  "82","81",
"-4603784299479457646935574969019046849794257872751288919656867","230010",
  "80","79",
"414846365575400828295179035549542073492199375372400483487","3318","78","77",
"-24655088825935372707687196040585199904365267828865801","30","76","75",
"34152417289221168014330073731472635186688307783087","6","74","73",
"-5827954961669944110438277244641067365282488301844260429","140100870",
  "72","71",
"1505381347333367003803076567377857208511438160235","4686","70","69",
"-78773130858718728141909149208474606244347001","30","68","67",
"1472600022126335654051619428551932342241899101","64722","66","65",
"-106783830147866529886385444979142647942017","510","64","63",
"12300585434086858541953039857403386151","6","62","61",
"-1215233140483755572040304994079820246041491","56786730","60","59",
"84483613348880041862046775994036021","354","58","57",
"-2479392929313226753685415739663229","870","56","55",
"29149963634884862421418123812691","798","54","53",
"-801165718135489957347924991853","1590","52","51",
"495057205241079648212477525","66","50","49",
"-5609403368997817686249127547","46410","48","47",
"596451111593912163277961","282","46","45",
"-27833269579301024235023","690","44","43",
"1520097643918070802691","1806","42","41",
"-261082718496449122051","13530","40","39",
"2929993913841559","6","38","37",
"-26315271553053477373","1919190","36","35",
"2577687858367","6","34","33",
"-7709321041217","510","32","31",
"8615841276005","14322","30","29",
"-23749461029","870","28","27",
"8553103","6","26","25",
"-236364091","2730","24","23",
"854513","138","22","21",
"-174611","330","20","19",
"43867","798","18","17",
"-3617","510","16","15",
"7","6","14","13",
"-691","2730","12","11",
"5","66","10","9",
"-1","30","8","7",
"1","42","6","5",
"-1","30","4","3",
"1","6","2","1",
};


/* This is global so it can be used by qcgamma.c.  */
QELT qgamcof[NG][NQ];

#if 0
static QELT ggamcof[NG][NQ] = {
/* -1215233140483755572040304994079820246041491/(56786730*60*59) */
{0xffff,0x4067,0x0000,0x989a,0x1506,0x8967,0x2663,0xf8cc,
0x3b4f,0x4518,0x35e1,0x174b,0x18c9,0xbd60,0xa7d4,0xb5c7,
0x3e81,0xcd5c,0xee1b,0xa1e0,0x619a,0xaa6b,0x15cb,0x7924,},
/* 84483613348880041862046775994036021/(354*58*57) */
{0x0000,0x4060,0x0000,0xe940,0xb372,0x3e6c,0x7d0e,0x7770,
0xe671,0x0431,0x6dcb,0x45c0,0xeb78,0xe86d,0x69d5,0xedfb,
0xfdaa,0x6336,0x8772,0x837d,0x2bb8,0x3de9,0xb8b2,0xd64f,},
/* -2479392929313226753685415739663229/(870*56*55) */
{0xffff,0x405a,0x0000,0xbf58,0x2a43,0x3556,0xfb17,0x24c9,
0x5ab5,0x6cbe,0xc2ef,0x3ba1,0xb5ef,0x128b,0x1478,0x2409,
0xb737,0xc796,0xcbec,0x2ef3,0xba1b,0x5ef1,0x28b1,0x4782,},
/* 29149963634884862421418123812691/(798*54*53) */
{0x0000,0x4054,0x0000,0xa8eb,0xfe48,0xda17,0xdd99,0x9790,
0x760b,0x0ce0,0x256e,0xc758,0x797b,0xf482,0x6900,0x66ba,
0x7710,0xd482,0x31e8,0x3250,0x2fff,0xebcc,0x7550,0x7f8b,},
/* -801165718135489957347924991853/(1590*52*51) */
{0xffff,0x404e,0x0000,0xa0ef,0x80e5,0x7954,0x084c,0xda64,
0x925c,0x6c86,0x491a,0x694d,0xeef0,0x8cb9,0xcebd,0x0717,
0x30f3,0xc513,0xf899,0x9b37,0x6479,0x67b1,0xc1db,0x9e70,},
/* 495057205241079648212477525/(66*50*49) */
{0x0000,0x4048,0x0000,0xa5f7,0xeef9,0xe71a,0xc7c8,0x0326,
0xab4c,0xc8bf,0x3f7c,0x478f,0x4715,0xb086,0x40e9,0x0b3d,
0x95ed,0x5188,0xa0d9,0xd86b,0xa7f9,0xaad3,0x322f,0xcfdf,},
/* -5609403368997817686249127547/(46410*48*47) */
{0xffff,0x4042,0x0000,0xb9e0,0x9405,0x8ad8,0x9016,0xb4f9,
0x2ff9,0x86cd,0xeea2,0x09d8,0xd881,0xad45,0x7156,0x1f50,
0xa7d3,0x0f4b,0x3a8c,0x9209,0x871f,0xa77c,0x3ae3,0x6671,},
/* 596451111593912163277961/(282*46*45) */
{0x0000,0x403c,0x0000,0xe2e1,0x337f,0x5af0,0xbed9,0x0b6b,
0x0a35,0x2d4f,0x335c,0x83da,0x6597,0xd322,0x6a6f,0x46ba,
0x523a,0x04c2,0xcb44,0x4e98,0x6e10,0xdd6c,0x4012,0x1349,},
/* -27833269579301024235023/(690*44*43) */
{0xffff,0x4037,0x0000,0x977d,0x7628,0x7772,0x9bcb,0x4050,
0x9f4f,0xd884,0x644b,0x7203,0x7c5e,0x1516,0x61b4,0xcbd5,
0x6936,0x7a8e,0x3a5c,0x877f,0x7942,0xe745,0x7deb,0xabff,},
/* 1520097643918070802691/(1806*42*41) */
{0x0000,0x4031,0x0000,0xde46,0x6b7c,0x78fb,0xaae3,0xc3a9,
0xe6da,0xeae4,0x6d98,0xeeec,0xac9e,0x8573,0xed1e,0xaac0,
0x952d,0x3b2d,0xcf5d,0x7d42,0x4e93,0x0f91,0x4416,0xc3fc,},
/* -261082718496449122051/(13530*40*39) */
{0xffff,0x402c,0x0000,0xb400,0x5bde,0x03d4,0x642a,0x2435,
0x8171,0x4af6,0x42a2,0x4358,0x1714,0xaf64,0x2a24,0x3581,
0x714a,0xf642,0xa243,0x5817,0x14af,0x642a,0x2435,0x8171,},
/* 2929993913841559/(6*38*37) */
{0x0000,0x4027,0x0000,0xa1bb,0xcde4,0xea01,0x2735,0x0b88,
0x1273,0x50b8,0x8127,0x350b,0x8812,0x7350,0xb881,0x2735,
0x0b88,0x1273,0x50b8,0x8127,0x350b,0x8812,0x7350,0xb881,},
/* -26315271553053477373/(1919190*36*35) */
{0xffff,0x4022,0x0000,0xa228,0x8cec,0xf233,0x76ae,0xa602,
0x4d5c,0x4976,0x1634,0xda88,0xc079,0x3f07,0xa1f8,0xf306,
0x6b7e,0xa521,0xc5cd,0xdbb1,0x324e,0x3cd1,0x0483,0xf9be,},
/* 2577687858367/(6*34*33) */
{0x0000,0x401d,0x0000,0xb694,0xd07b,0x219d,0xbcc4,0x8676,
0xf312,0x19db,0xcc48,0x676f,0x3121,0x9dbc,0xc486,0x76f3,
0x1219,0xdbcc,0x4867,0x6f31,0x219d,0xbcc4,0x8676,0xf312,},
/* -7709321041217/(510*32*31) */
{0xffff,0x4018,0x0000,0xe884,0x4d8a,0x169a,0xbbc4,0x0616,
0x9abb,0xc406,0x169a,0xbbc4,0x0616,0x9abb,0xc406,0x169a,
0xbbc4,0x0616,0x9abb,0xc406,0x169a,0xbbc4,0x0616,0x9abc,},
/* 8615841276005/(14322*30*29) */
{0x0000,0x4014,0x0000,0xa8d1,0x044d,0x3708,0xd1c2,0x19ee,
0x4fdc,0x4469,0xccae,0xdcb0,0x0698,0x234d,0x582b,0x9609,
0x3d5e,0xe562,0xc084,0x6967,0xb29e,0xe369,0xd34f,0x6a58,},
/* -23749461029/(870*28*27) */
{0xffff,0x4010,0x0000,0x8d0c,0xc570,0xe255,0xbf59,0xff6e,
0xec24,0xb48f,0xf1b3,0x94d9,0x2e2f,0xd250,0x842c,0x8d54,
0x1cdb,0xcbc6,0x8b5d,0xebd1,0x01b8,0xd006,0xa04a,0xa199,},
/* 8553103/(6*26*25) */
{0x0000,0x400c,0x0000,0x8911,0xa740,0xda74,0x0da7,0x40da,
0x740d,0xa740,0xda74,0x0da7,0x40da,0x740d,0xa740,0xda74,
0x0da7,0x40da,0x740d,0xa740,0xda74,0x0da7,0x40da,0x740e,},
/* -236364091/(2730*24*23) */
{0xffff,0x4008,0x0000,0x9cd9,0x292e,0x6660,0xd55b,0x3f71,
0x2eb9,0xe07c,0xa39d,0xb44a,0x9292,0xe666,0x0d55,0xb3f7,
0x12eb,0x9e07,0xca39,0xdb44,0xa929,0x2e66,0x60d5,0x5b3f,},
/* 854513/(138*22*21) */
{0x0000,0x4004,0x0000,0xd672,0x2191,0x6700,0x2d3a,0x7a9c,
0x8864,0x59c0,0x0b4e,0x9ea7,0x2219,0x1670,0x02d3,0xa7a9,
0xc886,0x459c,0x00b4,0xe9ea,0x7221,0x9167,0x002d,0x3a7b,},
/* -174611/(330*20*19) */
{0xffff,0x4001,0x0000,0xb23b,0x3808,0xc0f9,0xcf6d,0xedce,
0x7312,0xcc3e,0xa607,0x48b1,0x4c1f,0x4aa7,0x0223,0xb380,
0x8c0f,0x9cf6,0xdedc,0xe731,0x2cc3,0xea60,0x748b,0x14c2,},
/* 43867/(798*18*17) */
{0x0000,0x3ffe,0x0000,0xb7f4,0xb1c0,0xf033,0xffd0,0xc3b7,
0xf4b1,0xc0f0,0x33ff,0xd0c3,0xb7f4,0xb1c0,0xf033,0xffd0,
0xc3b7,0xf4b1,0xc0f0,0x33ff,0xd0c3,0xb7f4,0xb1c0,0xf034,},
/* -3617/(510*16*15) */
{0xffff,0x3ffb,0x0000,0xf214,0x3658,0x7a9c,0xbee1,0x0325,
0x4769,0x8bad,0xcff2,0x1436,0x587a,0x9cbe,0xe103,0x2547,
0x698b,0xadcf,0xf214,0x3658,0x7a9c,0xbee1,0x0325,0x476a,},
/* 7/(6*14*13) */
{0x0000,0x3ff9,0x0000,0xd20d,0x20d2,0x0d20,0xd20d,0x20d2,
0x0d20,0xd20d,0x20d2,0x0d20,0xd20d,0x20d2,0x0d20,0xd20d,
0x20d2,0x0d20,0xd20d,0x20d2,0x0d20,0xd20d,0x20d2,0x0d21,},
/* -691/(2730*12*11) */
{0xffff,0x3ff7,0x0000,0xfb55,0x86cc,0xc9e3,0xe40f,0xb558,
0x6ccc,0x9e3e,0x40fb,0x5586,0xccc9,0xe3e4,0x0fb5,0x586c,
0xcc9e,0x3e40,0xfb55,0x86cc,0xc9e3,0xe40f,0xb558,0x6ccd,},
/* 5/(66*10*9) */
{0x0000,0x3ff6,0x0000,0xdca8,0xf158,0xc7f9,0x1ab8,0x7539,
0xc037,0x2a3c,0x5631,0xfe46,0xae1d,0x4e70,0x0dca,0x8f15,
0x8c7f,0x91ab,0x8753,0x9c03,0x72a3,0xc563,0x1fe4,0x6ae2,},
/* -1/(30*8*7) */
{0xffff,0x3ff6,0x0000,0x9c09,0xc09c,0x09c0,0x9c09,0xc09c,
0x09c0,0x9c09,0xc09c,0x09c0,0x9c09,0xc09c,0x09c0,0x9c09,
0xc09c,0x09c0,0x9c09,0xc09c,0x09c0,0x9c09,0xc09c,0x09c1,},
/* 1/(42*6*5) */
{0x0000,0x3ff6,0x0000,0xd00d,0x00d0,0x0d00,0xd00d,0x00d0,
0x0d00,0xd00d,0x00d0,0x0d00,0xd00d,0x00d0,0x0d00,0xd00d,
0x00d0,0x0d00,0xd00d,0x00d0,0x0d00,0xd00d,0x00d0,0x0d01,},
/* -1/(30*4*3) */
{0xffff,0x3ff8,0x0000,0xb60b,0x60b6,0x0b60,0xb60b,0x60b6,
0x0b60,0xb60b,0x60b6,0x0b60,0xb60b,0x60b6,0x0b60,0xb60b,
0x60b6,0x0b60,0xb60b,0x60b6,0x0b60,0xb60b,0x60b6,0x0b61,},
/* 1/(6*2*1) */
{0x0000,0x3ffd,0x0000,0xaaaa,0xaaaa,0xaaaa,0xaaaa,0xaaaa,
0xaaaa,0xaaaa,0xaaaa,0xaaaa,0xaaaa,0xaaaa,0xaaaa,0xaaaa,
0xaaaa,0xaaaa,0xaaaa,0xaaaa,0xaaaa,0xaaaa,0xaaaa,0xaaab,},
};
#endif

/* log(2*pi)/2 */
#if WORDSIZE ==  32
QELT qgam12[] = {
0,EXPONE-1,0,0xeb3f8e43,0x25f5a534,0x94bc9001,0x44192023,0xcfb08f8d,
0x13458b4d,0xdec6a313,0x3daa155d,0x212f9d7f,0xe00e86bf,0x93eabf90
};
#else
QELT qgam12[] = {
0x0000,EXPONE-1,0x0000,0xeb3f,0x8e43,0x25f5,0xa534,0x94bc,
0x9001,0x4419,0x2023,0xcfb0,0x8f8d,0x1345,0x8b4d,0xdec6,
0xa313,0x3daa,0x155d,0x212f,0x9d7f,0xe00e,0x86bf,0x93eb,
};
#endif

extern QELT q32[], qone[], qhalf[];

/* Initialize Stirling formula coefficients from ASCII strings.  */
int qgamini = 0;
int asctoq(), qfloor(), qinfin();

int initqgam()
{
QELT den[NQ], temp[NQ];
int i, k;

k = 0;
for( i=0; i<NG; i++ )
  {
    asctoq( gstrs[k++], qgamcof[i] );
    asctoq( gstrs[k++], den );
    asctoq( gstrs[k++], temp );
    qmul( den, temp, den );
    asctoq( gstrs[k++], temp );
    qmul( den, temp, den );
    qdiv( den, qgamcof[i], qgamcof[i] );
  }
qgamini = 1;
return 0;
}

/* natural logarithm of gamma function */

int qlgam( x, y )
QELT *x, *y;
{
QELT v[NQ], w[NQ], g[NQ], xx[NQ], t[NQ];
QELT *p;
long lr;
int i;

if( qgamini == 0 )
  initqgam();

qmov( x, xx );
if( xx[0] != 0 )
 {
/* Test for negative integer.  */
   qfloor( xx, t );
   if( qcmp(xx, t) == 0 )
     {
       qinfin(y);
       y[0] = 0;
       return(0);
     }
 xx[0] = 0;
 qlgam( xx, t );
 qmul( xx, qpi, g ); /* PI/(x * sin(PI*x)) */
 qsin( g, g );
 qmul( xx, g, g );
 qdiv( g, qpi, y );
 y[0] = 0;
 qlog( y, y );
 qsub( t, y, y ); /* ... /gamma(x) */
 return(0);
 }
/* range reduction: transform argument to be greater than 32 */
qmov( qone, v );

if( xx[1] > (QELT) (EXPONE + 5) )
  goto donrecur;

qifrac( xx, &lr, t );

while( lr < 62 )
 {
 qmul( xx, v, v );
 qadd( qone, xx, xx );
 lr += 1;
 }

donrecur:

qdiv( v, qone, v );

/* Avoid overflow of xx * xx.  */
if( (int) (xx[1] - EXPONE) >= (MAXEXP - EXPONE)/2 )
  {
    qclear(g);
    goto asymp;
  }
/*  Asymptotic expansion in 1/x**2 for Stirling's formula */
qmul( xx, xx, w );
qdiv( w, qone, w );

p = &qgamcof[0][0];
qmul( w, p, g );
p += NQ;
qadd( g, p, g );

for( i=0; i<NG-2; i++ )
 {
 qmul( w, g, g );
 p += NQ;
 qadd( g, p, g );
 }

qdiv( xx, g, g );

asymp:

/* g + (x - 0.5)*log(x) - x + qgam12 */
qsub( qhalf, xx, t );
qlog( xx, w );
qmul( t, w, t );
qsub( xx, t, t );
qadd( qgam12, t, t );
qadd( g, t, g );

qlog( v, t );
qadd( t, g, y );
return  0;
}

/* qgamma() */
/* gamma function check routine */


int qgamma( xx, y )
QELT *xx, *y;
{
QELT t[NQ], s[NQ], x[NQ];

if( qgamini == 0 )
  initqgam();

qmov( xx, x );
if( x[0] != 0 )
 {
/* Test for negative integer.  */
   qfloor( xx, t );
   if( qcmp(xx, t) == 0 )
     {
       qinfin(y);
       y[0] = 0;
       return 0;
     }
 x[0] = 0;
 qgamma( x, t );
 qmul( x, qpi, s );
 qsin( s, s );
 qmul( x, s, s );
 qmul( s, t, t );
 qdiv( t, qpi, y );
 y[0] = ~y[0];
 return(0);
 }
qlgam( x, y );
qexp( y, y );
return 0;
}


/*  Asymptotic expansion in 1/x**2 for Stirling's formula */
int qlstir( xx, y )
QELT *xx, *y;
{
QELT w[NQ], g[NQ];
QELT *p;
int i;

if( qgamini == 0 )
  initqgam();

qmul( xx, xx, w );
qdiv( w, qone, w );

p = &qgamcof[0][0];
qmul( w, p, g );
p += NQ;
qadd( g, p, g );
for( i=0; i<NG-2; i++ )
 {
 qmul( w, g, g );
 p += NQ;
 qadd( g, p, g );
 }
qdiv( xx, g, y );
return 0;
}