STK_Funct_gamma.cpp

Go to the documentation of this file.

/*--------------------------------------------------------------------*/
/*     Copyright (C) 2004-2007  Serge Iovleff

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU Lesser 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 Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this program; if not, write to the
    Free Software Foundation, Inc.,
    59 Temple Place,
    Suite 330,
    Boston, MA 02111-1307
    USA

    Contact : Serge.Iovleff@stkpp.org
*/

/*
 * Project:  Analysis
 * Purpose:  implementation of functions around the gamma function
 * Author:   Serge Iovleff, serge.iovleff@stkpp.org
 **/

#include "../../STKernel/include/STK_Integer.h"
#include "../../STKernel/include/STK_Real.h"
//#include "../../STKernel/include/STK_Algebra.h" // for Vector

#include "../include/STK_Const_Math.h"
#include "../include/STK_Funct_gamma.h"

namespace STK
{

namespace Funct
{
static const Real factorialArray[51] =
{
                                                                1.0, //  0
                                                                1.0, //  1
                                                                2.0, //  2
                                                                6.0, //  3
                                                               24.0, //  4
                                                              120.0, //  5
                                                              720.0, //  6
                                                             5040.0, //  7
                                                            40320.0, //  8
                                                           362880.0, //  9
                                                          3628800.0, // 10
                                                         39916800.0, // 11
                                                        479001600.0, // 12
                                                       6227020800.0, // 13
                                                      87178291200.0, // 14
                                                    1307674368000.0, // 15
                                                   20922789888000.0, // 16
                                                  355687428096000.0, // 17
                                                 6402373705728000.0, // 18
                                               121645100408832000.0, // 19
                                              2432902008176640000.0, // 20
                                             51090942171709440000.0, // 21
                                           1124000727777607680000.0, // 22
                                          25852016738884976640000.0, // 23
                                         620448401733239439360000.0, // 24
                                       15511210043330985984000000.0, // 25
                                      403291461126605635584000000.0, // 26
                                    10888869450418352160768000000.0, // 27
                                   304888344611713860501504000000.0, // 28
                                  8841761993739701954543616000000.0, // 29
                                265252859812191058636308480000000.0, // 30
                               8222838654177922817725562880000000.0, // 31
                             263130836933693530167218012160000000.0, // 32
                            8683317618811886495518194401280000000.0, // 33
                          295232799039604140847618609643520000000.0, // 34
                        10333147966386144929666651337523200000000.0, // 35
                       371993326789901217467999448150835200000000.0, // 36
                     13763753091226345046315979581580902400000000.0, // 37
                    523022617466601111760007224100074291200000000.0, // 38
                  20397882081197443358640281739902897356800000000.0, // 39
                 815915283247897734345611269596115894272000000000.0, // 40
               33452526613163807108170062053440751665152000000000.0, // 41
             1405006117752879898543142606244511569936384000000000.0, // 42
            60415263063373835637355132068513997507264512000000000.0, // 43
          2658271574788448768043625811014615890319638528000000000.0, // 44
        119622220865480194561963161495657715064383733760000000000.0, // 45
       5502622159812088949850305428800254892961651752960000000000.0, // 46
     258623241511168180642964355153611979969197632389120000000000.0, // 47
   12413915592536072670862289047373375038521486354677760000000000.0, // 48
  608281864034267560872252163321295376887552831379210240000000000.0, // 49
30414093201713378043612608166064768844377641568960512000000000000.0  // 50
};

static const Real factorialLnArray[51] =
{
  0.0000000000000000000000000000000000000000000000000000000000000, //  0
  0.0000000000000000000000000000000000000000000000000000000000000, //  1
  0.6931471805599453094172321214581765680755001343602552541206800, //  2
  1.7917594692280550008124773583807022727229906921830047058553743, //  3
  3.1780538303479456196469416012970554088739909609035152140967343, //  4
  4.7874917427820459942477009345232430483995923151720329360093822, //  5
  6.5792512120101009950601782929039453211225830073550376418647565, //  6
  8.5251613610654143001655310363471250507596677369368988303241467, //  7
 10.6046029027452502284172274007216547549861681400176645926861867, //  8
 12.8018274800814696112077178745667061642811492556631634961555754, //  9
 15.1044125730755152952257093292510703718822507442919364721889033, // 10
 17.5023078458738858392876529072161996717039575982293536474074710, // 11
 19.9872144956618861495173623870550785125024484247726136073835254, // 12
 22.5521638531234228855708498286203971173077163695328207238025709, // 13
 25.1912211827386815000934346935217534150203012334749371663826410, // 14
 27.8992713838408915660894392636704667591933931455662043400299833, // 15
 30.6718601060806728037583677495031730314953936830072253565127033, // 16
 33.5050734501368888840079023673762995670835966955929701438099410, // 17
 36.3954452080330535762156249626795275444540779455987243014000097, // 18
 39.3398841871994940362246523945673810816914572068978531199379699, // 19
 42.3356164607534850296598759707099218573680588298868813500919778, // 20
 45.3801388984769080261604739510756272916526341172914919902860623, // 21
 48.4711813518352238796396496504989331595498411055891644196253101, // 22
 51.6066755677643735704464024823091292779922214204296001616239454, // 23
 54.7847293981123191900933440836061846868662123813331153757206798, // 24
 58.0036052229805199392948627500585599659174150898701508195459756, // 25
 61.2617017610020019847655823130820551387981831689906131900857011, // 26
 64.5575386270063310589513180238496322527406548424588615452897841, // 27
 67.8897431371815349828911350102091651185287398407612332419905343, // 28
 71.2570389671680090100744070425710767240232527546839773209122046, // 29
 74.6582363488301643854876437341779666362718448011354997486802268, // 30
 78.0922235533153106314168080587203238467217837316160917204369449, // 31
 81.5579594561150371785029686660112066870992844034173679910403450, // 32
 85.0544670175815174139601574808988616915684818151775346179936070, // 33
 88.5808275421976788036269242202301647952321849621235346594115248, // 34
 92.1361756036870924833330362968995321643948710459739135697835628, // 35
 95.7196945421432024849579910136609367098408524303399229814943115, // 36
 99.3306124547874269293260866846923838737409300175075591451392242, // 37
102.9681986145138126987523462380384139790538094131669432177978644, // 38
106.6317602606434591262010789165262582885065679157498997859516043, // 39
110.3206397147573954290535346141269756322586696730991832702262922, // 40
114.0342117814617032329202979871643832206350801424984849038455551, // 41
117.7718813997450715388381280889882652229951555642633507981603196, // 42
121.5330815154386339623109706023341122585542917491449062133520048, // 43
125.3172711493568951252073784232155946945269988718028338968119326, // 44
129.1239336391272148825986282302868337433475813417168505221939691, // 45
132.9525750356163098828226131835552064298654617909175415183132845, // 46
136.8027226373263684696435638533273801387615122929377655515137929, // 47
140.6739236482342593987077375760826121157110033882015360197312072, // 48
144.5657439463448860089184430629689715749851728473652583966499875, // 49
148.4777669517730320675371938508795234221118756902625490945959633  // 50
};

static const Real factorialHalvesArray[50] =
{
1.772453850905516027298167483341145182797549456122387128213807789, //  0.5
0.886226925452758013649083741670572591398774728061193564106903894, //  1.5
1.329340388179137020473625612505858887098162092091790346160355841, //  2.5
3.323350970447842551184064031264647217745405230229475865400889604, //  3.5
11.63172839656744892914422410942626526210891830580316552890311361, //  4.5
52.34277778455352018114900849241819367949013237611424488006401126, //  5.5
287.8852778150443609963195467083000652371957280686283468403520619, //  6.5
1871.254305797788346476077053603950424041772232446084254462288402, //  7.5
14034.40729348341259857057790202962818031329174334563190846716302, //  8.5
119292.4619946090070878499121672518395326629798184378712219708856, //  9.5
1133278.388948785567334574165588892475560298308275159776608723413, // 10.5
11899423.08396224845701302873868337099338313223688917765439159584, // 11.5
136843365.4655658572556498304948587664239060207242255430255033522, // 12.5
1710542068.319573215695622881185734580298825259052819287818791902, // 13.5
23092317922.31423841189090889600741683403414099721306038555369069, // 14.5
334838609873.5564569724181789921075440934950444595893755905285150, // 15.5
5189998453040.125083072481774377666933449173189123635321653191982, // 16.5
85634974475162.06387069594927723150440191135762053998280727766771, // 17.5
1498612053315336.117737179112351551327033448758359449699127359184, // 18.5
27724322986333718.17813781357850369955011880202964981943385614492, // 19.5
540624298233507504.4736873647808221412273166395781714789601948259, // 20.5
11082798113786903841.71059097800685389515999111135251531868399393, // 21.5
238280159446418432596.7777060271473587459398088940790793517058695, // 22.5
5361303587544414733427.498385610815571783645700116779285413382064, // 23.5
125990634307293746235546.2120618541659369156739527443132072144785, // 24.5
3086770540528696782770882.195515427065454434011842235673576754723, // 25.5
78712648783481767960657495.98564339016908806730197700967620724545, // 26.5
2085885192762266850957423643.619549839480833783502390756419492004, // 27.5
57361842800962338401329150199.53762058572292904631574580153603012, // 28.5
1634812519827426644437880780686.822186693103477819998755343776858, // 29.5
48226969334909086010917483030261.25450744655259568996328264141733, // 30.5
1470922564714727123332983232422968.262477119854168543880120563228, // 31.5
46334060788513904384988971821323500.26802927540630913222379774170, // 32.5
1505856975626701892512141584193013758.710951450705046797273426605, // 33.5
50446208683494513399156743070465960916.81687359861906770865979127, // 34.5
1740394199580560712270907635931075651630.182139152357835948762799, // 35.5
61783994085109905285617221075553185632871.46593990870317618107936, // 36.5
2255115784106511542925028569257691275599808.506806667665930609396, // 37.5
84566841903994182859688571347163422834992819.00525003747239785238, // 38.5
3255823413303776040098009996865791779147223531.702126442687317316, // 39.5
128605024825499153583871394876198775276315329502.2339944861490340, // 40.5
5208503505432715720146791492486050398690770844840.476776689035877, // 41.5
216152895475457702386091846938171091545666990060879.7862325949889, // 42.5
9186498057706952351408903494872271390690847077587390.914885287028, // 43.5
399612665510252427286287302026943805495051847875051504.7975099857, // 44.5
17782763615206233014239784940198999344529807230439791963.48919436, // 45.5
809115744491883602147910214779054470176106228985010534338.7583436, // 46.5
37623882118872587499877824987226032863188939647802989846752.26297, // 47.5
1787134400646447906244196686893236561001474633270642017720732.491, // 48.5
86676018431352723452843539314321973208571519713626137859455525.83, // 49.5
};

static const Real factorialLnHalvesArray[50] =
{
 0.5723649429247000870717136756765293558236474064576557857568114,  //  0.5
-0.1207822376352452223455184457816472122518527279025994683638685,  //  1.5
 0.2846828704729191596324946696827019243201376955598947292501457,  //  2.5
 1.2009736023470742248160218814507129957702389154681571970421136,  //  3.5
 2.4537365708424422205041425034357161573318235106897631313808238,  //  4.5
 3.9578139676187162938774008558225909985513044919750067807295324,  //  5.5
 5.6625620598571415285221123123295437302975112115521687018274201,  //  6.5
 7.5343642367587329551583676324366857670272790219521205641257856,  //  7.5
 9.5492672573009977117371400811272225431248707996831324836524478,  //  8.5
11.6893334207972684825694425775421725106375736779086220168290056,  //  9.5
13.9406252194037636331612378879718494797994528048474955812462858,  // 10.5
16.2920004765672413202446037468793783460085279578918509673196903,  // 11.5
18.7343475119364457016341244572313978963754081383720314551976456,  // 12.5
21.2600761562447011414184110022255966073511107125488116449022614,  // 13.5
23.8627658416890849061869145915349971532180822516568047459856644,  // 14.5
26.5369144911156136239529545024387321906370950312192935707866548,  // 15.5
29.2777545150408145604648867055229128330115338273396302884226928,  // 16.5
32.0811148959473494865048433989523912694052311047395416612552748,  // 17.5
34.9433157768768178567937233541635820704924170542296653175066329,  // 18.5
37.8610865089610969917445869037368526663169945070370462270308656,  // 19.5
40.8315009745307981097760874607665204076942528752597475410639254,  // 20.5
43.8519258606751606042256187123457514279951632102987939205625083,  // 21.5
46.9199787958087777182812291042334218954787992608200940816335135,  // 22.5
50.0334941050191521662552467898464843762238815963738554528948701,  // 23.5
53.1904945261692654436589653381604815170444319640338242319746984,  // 24.5
56.3891676437199467444524387035886644082431012888372913547727987,  // 25.5
59.6278460958843272066799864369261400804032947248855303396840508,  // 26.5
62.9049908288765037314072234544970212826877723434364926886115002,  // 27.5
66.2191768335490293406526942442301616539595804172821723316220358,  // 28.5
69.5690809208236341826397347915823643277689501020437953477740104,  // 29.5
72.9534711841694083238385530438438853837567967570139946231653671,  // 30.5
76.3711978677827742631727100258113235619969783657350004319820707,  // 31.5
79.8211854136143616416513211216413781328535440766021052697901695,  // 32.5
83.3024255029500534428883357749747078091089132412705748540011829,  // 33.5
86.8139709417810741931411756498802539916003460101322149370182870,  // 34.5
90.3549302658183882659259415971547992466147167484351448766309367,  // 35.5
93.9244629622997583778381640082096567752987931407854557306316825,  // 36.5
97.5217752228882041975130407441900227781280458512967981411067014,  // 37.5
101.1461161558645693286925725261067471937512389832963277825460115, // 38.5
104.7967743971583078684426367260568796551345304324553508921032893, // 39.5
108.4730750690653840531983501460801140092325715225177052633947804, // 40.5
112.1743770431778775093620989723120402597470336194484478162128778, // 41.5
115.9000704704145301234203390741452341447008465417514289322801412, // 42.5
119.6495745463449012688534009037863717517391507742454361873693468, // 43.5
123.4223354844395396780146860516126324938056541116306744639050315, // 44.5
127.2178246736117342069152694708243051451348143644033785738910491, // 45.5
131.0355369995686389386568775343746269115016669043851916246488047, // 46.5
134.8749893121619495665640549743813332585235962583282777940195371, // 47.5
138.7357190232025450917566096180371978672110767395356690803494653, // 48.5
142.6172828211459826044560991182829830129444915078341951822050196  // 49.5
};

static const Real gammaLnStirlingErrorArray[100] =
{
.0,    // .0 - place holder only
.081061466795327258219670263594382360138602526362216587182848460, //  1.0
.041340695955409294093822081407117508025352324821833706001828447, //  2.0
.027677925684998339148789292746244666595376266165598211966792637, //  3.0
.020790672103765093111522771767848656333092278023434514193462740, //  4.0
.016644691189821192163194865373593391147387393057402052672667313, //  5.0
.013876128823070747998745727023762908561746034527723640987671827, //  6.0
.011896709945891770095055724117659438620134791435156504061569088, //  7.0
.010411265261972096497478567132534629199517240193372199791694996, //  8.0
.009255462182712732917728636633100136117431183393140500379231566, //  9.0
.008330563433362871256469318659628552209287640052036811021808852, // 10.0
.007573675487951840794972024211595083892931304311272719576790859, // 11.0
.006942840107209529865664152663475362659915619344080694865694464, // 12.0
.006408994188004207068439631082978312575201641632241239328304990, // 13.0
.005951370112758847735624416046469458326423232676465336154472229, // 14.0
.005554733551962801371038689959792284649071034513779735679027273, // 15.0
.005207655919609640440717996857901898650987341592595171730671187, // 16.0
.004901395948434737860716818190967554428646501704652953291129138, // 17.0
.004629153749334028592427213164192323238777346854488973166587706, // 18.0
.004385560249232324268287736348619465701164137927057745630593846, // 19.0
.004166319691996922457462923382218316136328084974019219117664204, // 20.0
.003967954218640859617287636807342814672867964454579813296094458, // 21.0
.003787618068444434577866677068933492001286395253751346675084906, // 22.0
.003622960224683094707381198363902854734886548041576811838862580, // 23.0
.003472021382978766962945115422709529592036365559209217533536432, // 24.0
.003333155636728092875807019117372710250348548537963589183781642, // 25.0
.003204970228055038011184156553815417596431598660576957965823324, // 26.0
.003086278682608777063256241335643979461286348444248364360350104, // 27.0
.002976063983550408826021162556860803706919915505856473202002985, // 28.0
.002873449362352466387552351489066722073724320325243579905777771, // 29.0
.002777674929752693603594903762206672828389910726282288938397223, // 30.0
.002688078828531142867802099230454077687309947839661197283174164, // 31.0
.002604081919251656422419192651896734969115096238105783612691953, // 32.0
.002525175249756784364002445756801401988971047574169202242177115, // 33.0
.002450909735438118343141976894590078973028918848749817676210999, // 34.0
.002380887608234111985727838731372919258117595647678838759060856, // 35.0
.002314755290514683866814115464053161201134427342796041234832959, // 36.0
.002252197424337523742169284607497587626623025083419595637847142, // 37.0
.002192931843287834061023697810637284646555206642873007623063492, // 38.0
.002136705317752500195808379850790425261209590085331931061728931, // 39.0
.002083289938302421748749124892305570437151026815418744284279969, // 40.0
.002032480028256630669288269706350663152648188789683695828992834, // 41.0
.001984089497245795550019025067662482830552665618766876965677141, // 42.0
.001937949563995799461991535384148571870470233159462239696544373, // 43.0
.001893906789600634536900657584008653880134439887268570027997190, // 44.0
.001851821372993179516432272139839382149681486991310654494154255, // 45.0
.001811565668719630626983120170258866922623429447626785947967582, // 46.0
.001773022893912853868077302743511326337715973339340561672494943, // 47.0
.001736085996876597314982286048243593799287794270884898369458756, // 48.0
.001700656664196061708191725688560701052387145123217326353211249, // 49.0
.001666644446983365509949324991037502351984650311585435507034346, // 50.0
.001633965989907857695481447293467650347952910995806610666242257, // 51.0
.001602544349176320672132218463463363238358792776302250799859171, // 52.0
.001572308387716159704965854725696819543607972675053385444352268, // 53.0
.001543192237554763675285768702465545840926034600252049644457878, // 54.0
.001515134820843602919723965079268021392092214030246864192791724, // 55.0
.001488079422197137383510052314484575738722906683837884385882645, // 56.0
.001461973306045267715378856228640422536138121736138335391704498, // 57.0
.001436767373567118932068468170376850640402444971677226134196516, // 58.0
.001412415854510449022614604043223996780378593695548232657545860, // 59.0
.001388876029827013264717462993777049266027800556820666245945813, // 60.0
.001366107981587896049039896180606956774859283253076112800720992, // 61.0
.001344074367099040456088327823086475874783995038527091989833241, // 62.0
.001322740214528256115665122342001450460610220881472963778537242, // 63.0
.001302072737691049536271144579406335881240201561586981276908565, // 64.0
.001282041167932156806185174841585421758722411627618253585847749, // 65.0
.001262616601289715384777593033276423752972227121461953968561174, // 66.0
.001243771859345481473982101843513085274320979873477125723062987, // 67.0
.001225481362352321033918753104920906146070459179179933668143233, // 68.0
.001207721013393527450973008992833538625340451229481907329004317, // 69.0
.001190468092470846417417016052696605317423935681778439551514756, // 70.0
.001173701159542375886638441875640359059407186436576413727618241, // 71.0
.001157399965640261426278583857528660583200257314129537462411392, // 72.0
.001141545371293452806932190482823000296029111777487127052824648, // 73.0
.001126119271564537855340464228357661590400677539896269065627922, // 74.0
.001111104527083370403350203275202408218292587563832173383019685, // 75.0
.001096484900525182775403460362622836646979681244822616323594192, // 76.0
.001082244998038253004347995148792597260457588522436293577817486, // 77.0
.001068370215176956357540908476124371876760948167593252352012952, // 78.0
.001054846686941007759552697948005867222077831034865706038538728, // 79.0
.001041661241561619068541630538108189724254781724001223874392175, // 80.0
.001028801357710777509444775725791195354247445416546206129935072, // 81.0
.001016255124841439713759211079062159605455089228981120834939896, // 82.0
.001004011206394598916929085426692691969752724675423818874895099, // 83.0
.000992058805634328641164086204234888546834019608044835208558260, // 84.0
.000980387633894389604083085646136293236499512995599849768283029, // 85.0
.000968987881040117202834628893480154705367583860232161006273381, // 86.0
.000957850187967355152546746337694541818108892291664654750324331, // 87.0
.000946965620976403329312854005796325025277680931341548879057281, // 88.0
.000936325647873511959772401783962265979175665773627801320106442, // 89.0
.000925922115665562045130688352769690496887231626157673376406747, // 90.0
.000915747229725383337708476995783752559445890490581176873828921, // 91.0
.000905793534315817200223067794124558740392196963249176291113078, // 92.0
.000896053894370256504918044498260962115868508291695650535039653, // 93.0
.000886521478436098067121061427062626885541163079305382450988328, // 94.0
.000877189742695421950545529418272596595408129891999961720819970, // 95.0
.000868052415984352210923475428688927573394487507788091465882714, // 96.0
.000859103485739031402950947590308139924760039653869445385373968, // 97.0
.000850337184802024059392250134312569701177680297232626720844722, // 98.0
.000841747979028312446492568768311415649867638491795533343982122  // 99.0
};

static const Real gammaLnStirlingErrorHalvesArray[100] =
{
.153426409720027345291383939270911715962249932819872372939659995, //  0.5
.054814121051917653896138702348386011314759374997122921204965662, //  1.5
.033162873519936287485110509741062141558537782105586380849058548, //  2.5
.023746163656297495971330279090085871224087656165555981208831773, //  3.5
.018488450532673185230779357482599152592502111911005121168837706, //  4.5
.015134973221917378873513836882209699958873420451413762422942497, //  5.5
.012810465242920226924250655281073870057506875514674115220075686, //  6.5
.011104559758206917326630755197310694482739104197253714621998689, //  7.5
.009799416126158803298390373402005163161850300241432805422834586, //  8.5
.008768700134139385462955047269462148319264061820976524256332105, //  9.5
.007934114564314020547249562490943177847303800766158307695089626, // 10.5
.007244554301320383179546196601545652111448498972082187845029804, // 11.5
.006665247032707682442356180895395724805579774612885765925104933, // 12.5
.006171712263039457647534604797771871219083230507117918000871100, // 13.5
.005746216513010115682026102477088989490505864144177785954648187, // 14.5
.005375599032926834493641719770404915159115615776459393428932388, // 15.5
.005049887331583002045249874245640210850979915303264910284228853, // 16.5
.004761386941714449813554423956530011061672489769656092165036888, // 17.5
.004504066155120685897849725439093705788625430059025533253293369, // 18.5
.004273129932103007365746583767547940287688037167937847585577273, // 19.5
.004064718438875479005132692776282361814648700841036538275115955, // 20.5
.003875689664528467277470746299053970977408675713707125254246501, // 21.5
.003703459975867121072510284340429679225590404999560501066518034, // 22.5
.003545885361874044189390442218932497489825666556790121241666911, // 23.5
.003401171748241482835274831478474488414970695112829872003140039, // 24.5
.003267806405762446983415633746154964376867913993533319174345740, // 25.5
.003144504883064821991494243677611181145479168922725854753213393, // 26.5
.003030169513639539999182395296613618487556875965358914509887857, // 27.5
.002923856655421023235229181962869147435831770318944466549596725, // 28.5
.002824750591511346093675739592636064259999634270232948607515044, // 29.5
.002732142563757402033513848415122396689951021746036945847106575, // 30.5
.002645413798892788555436318674014225582043855216967047113856275, // 31.5
.002564021667551260286522478569210852820333333573295162220589556, // 32.5
.002487488321695942764290166686066147030234395198668701619595500, // 33.5
.002415391307722976678803706140097681250343434194813577371045781, // 34.5
.002347355765761607217827473829168102810645542536481573787682597, // 35.5
.002283047911036127762271776510863033584298088739345540713852618, // 36.5
.002222169558021602090035076865601760208502493676417995035075999, // 37.5
.002164453497832076409803191596096021324766441615666206552301804, // 38.5
.002109659577663921189483608745122205449530446385746998917984704, // 39.5
.002057571360973064868067360397446348792690171650219737853736025, // 40.5
.002007993270447591193925060745473336779306331388438644637924777, // 41.5
.001960748134269273451412492811837009230671303365340805716353806, // 42.5
.001915675070776933137809810851542203018109793042406883517754399, // 43.5
.001872627658307664608685869893418195461205515993910211703107153, // 44.5
.001831472346348536762576874654208018765052591584007876974893374, // 45.5
.001792087071677313146385557660517308631519149365444420777963105, // 46.5
.001754360049287140733656986150079010350616167943103923695764172, // 47.5
.001718188712871740400783869833893232145770387568632778465704605, // 48.5
.001683478783723608350874794629162462064890108286031365591326262, // 49.5
.001650143450247406762664039488740073355159227551597803425876964, // 50.5
.001618102643055092658057584325705222748854647718478756747631518, // 51.5
.001587282392899875354797660326188057289246915997568382981557887, // 52.5
.001557614260611358593597978462839886381856289708771235952187656, // 53.5
.001529034829784896656045678422908813532915491630296862155318513, // 54.5
.001501485254310950424700612460491851735856109836801435394978214, // 55.5
.001474910853950712572427815385758856757170461930743391193481094, // 56.5
.001449260752109393525729544425119958254047478791042248450172758, // 57.5
.001424487550758294210716022733917102582538783579394001775420037, // 58.5
.001400547038135570067058714283678286205256559627296605551991148, // 59.5
.001377397925433403927071037446538221117698617193594124721539348, // 60.5
.001355001609172568501662341636381924907972137013373294795407964, // 61.5
.001333321956387599004960179081345140627983930836244310988249305, // 62.5
.001312325110108196559703272704388133910365080242690689320992664, // 63.5
.001291979312934320038298949371246032616437478436668032097671641, // 64.5
.001272254746771416751441990796291094033025452578721205107641001, // 65.5
.001253123387024833064335156136804425440462360711724892444215557, // 66.5
.001234558869754027032638955733794823668656891374728194445721656, // 67.5
.001216536370462303081373244678951501204430942444380014964787525, // 68.5
.001199032493350213803332672877142461521588387013766571922279510, // 69.5
.001182025169993741884563525536248284455480422866396723643587261, // 70.5
.001165493566524607038870514198105417531684964186189447772438923, // 71.5
.001149417998491846899645598047208909225643160644288785609677267, // 72.5
.001133779852673159899532447929836829412803498593956754811260725, // 73.5
.001118561515183045300193652082422858265371687313303138742416464, // 74.5
.001103746305293960014591823555737678141500083301598576755870739, // 75.5
.001089318414447757486643613794411926607472173791108599977011001, // 76.5
.001075262849988630764454027916160234496048121915426376466064959, // 77.5
.001061565383196552991770894391589130693187124981254720278102308, // 78.5
.001048212501242571885508870091039255740437415226235510761715635, // 79.5
.001035191362724942094212925947649409980698720159845456745577795, // 80.5
.001022489756478551715242865690353097234518477218195413853774395, // 81.5
.001010096063379920287948742256746173475755526743998752766856466, // 82.5
.000997999220896652538330565888998263866277090412496920283546988, // 83.5
.000986188690154005371505629476939383556360990020995750097419973, // 84.5
.000974654425312496505348708284181745201279558962849170710499920, // 85.5
.000963386845069540982855195388740797000921912790444631486117508, // 86.5
.000952376806115199499576328841712860289117767570900540991764240, // 87.5
.000941615578387481510246663134900698304725772551853832331795617, // 88.5
.000931094821986460690961738387983324908175451064875988434315466, // 89.5
.000920806565618901175526963343375677542263352145062493371806017, // 90.5
.000910743186456310188722507266407018763836146605677694240967302, // 91.5
.000900897391299458544421427592926592213225217648487572776568747, // 92.5
.000891262198951561680464610595170617679749184556736415659795466, // 93.5
.000881830923710593597100323580075849732974751023663787779031291, // 94.5
.000872597159898705505258014377097061988892776367102445141351277, // 95.5
.000863554767353521027072301030701774152795229669002347025320781, // 96.5
.000854697857812252154659342926302849688435895237596905019502674, // 97.5
.000846020782125188584738232943998298720583755594186163262272346, // 98.5
.000837518118240214182665925079081837356923152662339645584467363  // 99.5
};

static const Real lanczosCoefArray[21] =
{
   1.5333183020199267370932516012553e0,  // d_1
  -1.1640274608858812982567477805332e1,  // d_2
   4.0053698000222503376927701573076e1,  // d_3
  -8.2667863469173479039227422723581e1,  // d_4
   1.1414465885256804336106748692495e2,  // d_5
  -1.1135645608449754488425056563075e2,  // d_6
   7.9037451549298877731413453151252e1,  // d_7
  -4.1415428804507353801947558814560e1,  // d_8
   1.6094742170165161102085734210327e1,  // d_9
  -4.6223809979028638614212851576524e0,  // d_10
   9.7030884294357827423006360746167e-1, // d_11
  -1.4607332380456449418243363858893e-1, // d_12
   1.5330325530769204955496334450658e-2, // d_13
  -1.0773862404547660506042948153734e-3, // d_14
   4.7911128916072940196391032755132e-5, // d_15
  -1.2437781042887028450811158692678e-6, // d_16
   1.6751019107496606112103160490729e-8, // d_17
  -9.7674656970897286097939311684868e-11,// d_18
   1.8326577220560509759575892664132e-13,// d_19
  -6.4508377189118502115673823719605e-17,// d_20
   1.3382662604773700632782310392171e-21 // d_21
};

static const Real stirlingCoefArray[9] =
{
 +0.0833333333333333333333333333333333333333333333333333333333333333333,
 -0.0027777777777777777777777777777777777777777777777777777777777777777,
 +0.000793650793650793650793650793650793650793650793650793650793650,
 -0.0005952380952380952380952380952380952380952380952380952380952380,
 +0.000841750841750841750841750841750841750841750841750841750841750,
 -0.0019175269175269175269175269175269175269175269175269175269175269,
 +0.006410256410256410256410256410256410256410256410256410256410256410,
 -0.029550653594771241830065359477124183006535947712418300653594771241,
 +0.17964437236883057316493849001588939669435025472177
};

static Real lanczosSerie(Real const& z)
{
  Real sum = 0.0;
  for (Integer k=20; k>=0; k--)
    sum += lanczosCoefArray[k]/(z+k);
  return 2.0240434640140357514731512432760e-10 + sum;
}

static Real gammaLanczos(Real const& z)
{
  // 2 * sqrt(e/pi) = 1.86038273...
  return 1.8603827342052657173362492472666631120594218414085774528900013
       * exp((z-0.5)*(log(z+22.118910)-1.))
       * lanczosSerie(z);
}

static double stirlingSerie(Real const& z)
{
  // use stirling formula
  const Real z2 = z * z;
  if (z <= 50)
    return ( ( stirlingCoefArray[0]
             + ( stirlingCoefArray[1]
               + ( stirlingCoefArray[2]
                 + stirlingCoefArray[3]/z2
                 )/z2
               )/z2
             )/z
           );
    return ( ( stirlingCoefArray[0]
             + ( stirlingCoefArray[1]
                 + stirlingCoefArray[2]/z2
               )/z2
             )/z
           );
}

static Real gammaStirling( Real const& z)
{
  return ( Const::_SQRT2PI_
         * exp((z-0.5)*(log(z)-1.)+stirlingSerie(z)-0.5)
         );
}

static Real gammaLnStirling( Real const& z)
{
  return ( Const::_LNSQRT2PI_ + (z-0.5)*log(z) - z + stirlingSerie(z) );
}
// Last of the Static functions

Real factorial(Integer const& n)
{
  // Check if z is Available and finite
  if (Arithmetic<Integer>::isNA(n)) return n;
  // Negative reals argument not allowed
  if (n < 0)
    throw domain_error("Funct::factorial(n) "
                       "Negative Integer argument");
  // if n is less than 51 return a stored value
  if (n < 51) return factorialArray[n];
  else        return gamma(n + 1.);
}

Real factorial(Real const& z)
{
  // Check if z is Available and finite
  if (Arithmetic<Real>::isNA(z)) return z;
  Real n = floor(z);
  // Negative integers or reals arguments not allowed
  if ((n < 0)||(n != z))
    throw domain_error("Funct::factorial(x) "
                       "negative or not discrete argument");

  // if n is a small number return a stored value
  if (n < 51) return factorialArray[(Integer)n];
  else        return gamma(z + 1.);
}

Real factorialLn(Integer const& n)
{
  // Check if z is Available and finite
  if (Arithmetic<Integer>::isNA(n)) return n;
  // Negative reals argument not allowed
  if (n < 0)
    throw domain_error("Funct::factorial(n) "
                       "Negative Integer argument");
  // if n is a small number return a constant
  if (n < 51) return factorialLnArray[n];
  else        return gammaLn( n + 1.);
}

Real factorialLn(Real const& z)
{
  // Check if z is Available and finite
  if (Arithmetic<Real>::isNA(z)) return z;
  Real n = floor(z);
  // Negative integers or reals arguments not allowed
  if ((n < 0)||(n != z))
    throw domain_error("Funct::factorial(x) "
                       "negative or not discrete argument");
  // if n is a small number return a constant
  if (n < 51) return factorialLnArray[(Integer )n];
  else        return gammaLn(z + 1.);

}

Real gamma(Real const& z)
{
  // Check if z is Available and finite
  if (Arithmetic<Real>::isNA(z)) return z;
  // Negative integer argument not allowed
  if (( z<=0 && z == floor(z))||(Arithmetic<Real>::isInfinite(z)))
    throw domain_error("Funct::gamma(z) "
                       "Null or Negative integer argument");
  // If z is an integer (z integer and z<0 is not possible)
  if ((z == floor(z)))
  {
    // if n is a small number return a constant else stirling is accurate
    if (z < 51) return factorialArray[(Integer )z-1];
    else        return gammaStirling(z);
  }
  // compute the absolute value of x -> y
  Real y = abs(z);
  Real value;
  // compute the sign of the gamma function
  Real signgam = (z<0 && isEven((Integer )floor(y))) ? -1 : 1;
  // if y is an integer and a half, use reflection formula
  if (y == floor(y)+0.5)
  {
    // if n+0.5 is a small number return a constant else stirling is accurate
    if (y<50) value = factorialHalvesArray[(Integer )(y)];
    else      value = gammaStirling(y);
  }
  else // normal case
  {
    if (y <= 8)
    {
      Integer n = (Integer )y; // convert y to Integer value
      Real r = y-n;             // r in [0,1]
      value = gammaLanczos(r);  // use Lanzcos approximation
      // scale result
      for (Real i=0.; i<n; i+=1.) value *= (r+i);
    }
    else
    {
      if (y <=16)
      {
        Real r = y+6.;
        value = gammaStirling(r);
        for (Real i=5.; i>=0.; i-=1.) value /= (y+i);
      }
      else
        value = gammaStirling(y);
    }
  }
  // z >0 terminate
  if (z>0) return value;
  // z <0  -> use reflection formula
  Real den = y*sin(Const::_PI_ *y)*value;
  // overflow
  if (den == 0.0) return signgam*Arithmetic<Real>::infinity();
  // normal case
  return -Const::_PI_/den;
}

Real gammaLn(Real const& z)
{
  // Check if z is Available and finite
  if (Arithmetic<Real>::isNA(z)) return z;
  // Negative integer argument not allowed
  if (( z<=0 && z == floor(z))||(Arithmetic<Real>::isInfinite(z)))
    throw domain_error("Funct::gammaLn(z) "
                       "Null or Negative integer argument");
  // compute the absolute value of x -> y
  Real y = abs(z);
  Real value;
  // If x is an integer
  if (y == floor(z))
  { if (y<=50) value = factorialLnArray[(Integer )y];
    else       value = gammaLnStirling(y);
  }
  // If x is an integer and a half
  if (y == floor(y)+0.5)
  { if (y<=50) value = factorialLnHalvesArray[(Integer )y];
    else       value = gammaLnStirling(y);
  }
  // small values -> use
  if (y <= 16)  return log(abs(gamma(z)));
  else          value = gammaLnStirling(y);
  // z >0 terminate
  if (z>0) return value;
  // z <0  -> use reflectiong formula
  Real sinpiy = abs(sin(Const::_PI_ *y));
  // overflow
  if (sinpiy == 0.0) return -Arithmetic<Real>::infinity();
  // normal case
  return Const::_LNSQRTPI_2_+(z-0.5)*log(y)-z-log(sinpiy)+stirlingSerie(z);
}

Real gammaLnStirlingError(Real const& z)
{
  // Check if z is Available and finite
  if (Arithmetic<Real>::isNA(z)) return z;
  // Negative integer argument not allowed
  if ( z<=0 || Arithmetic<Real>::isInfinite(z))
    throw domain_error("Funct::gammaLnStirlingError(z) "
                       "Null or Negative integer argument");
  // z is a discrete value
  if (z == floor(z))
  {
    if (z<100.0) return gammaLnStirlingErrorArray[(Integer )z];
    else         return stirlingSerie(z);
  }
  // z is a discrete value halves
  if (z== floor(z) + 0.5)
  {
    if (z<100.0) return(gammaLnStirlingErrorHalvesArray[(Integer )z]);
    else         return stirlingSerie(z);
  }
  // other cases
  if (z > 16) return stirlingSerie(z);
  // z <16
  return (gammaLn(z) - (z-0.5)*log(z) + z - Const::_LNSQRT2PI_);
}

Real gammaLnStirlingError(Integer const& z)
{
  // Check if z is Available and finite
  if (Arithmetic<Integer>::isNA(z)) return Arithmetic<Real>::NA();
  // Negative integer argument not allowed
  if ( z<=0 || Arithmetic<Integer>::isInfinite(z))
    throw domain_error("Funct::gammaLnStirlingError(z) "
                       "Null or Negative integer argument");
  // z is a discrete value
  if (z<100.0) return gammaLnStirlingErrorArray[z];
  else         return stirlingSerie(z);
}

void stirlingCoefficients( STK::Vector& A)
{
  Integer shift = A.first()-1;
  Integer n = A.size();

  // initialization
  for (Integer j=1; j<=n; j++)
    A[shift+j] = pow(4., -j)/(2.*j + 1.);

  // for each coefficients
  Real aux = A[shift+1];
  for (Integer k=2; k<=n; k++)
  {
    // Update remaining terms
    Real coef = (k-0.5) * (k-1.) / 6.;
    for (Integer l=k, r=0; l<=n; l++, r++)
    {
      A[shift+l] -= aux* coef;
      coef *= ((l/(r+2.)) * ((l+0.5)/(r+2.5))) / 4.;
    }
    aux = A[shift+k];
    A[shift+k] /= (2.*k-1.);
  }
}

} // namespace Funct

} // namespace STK