[gnumeric] FACT: Improve accuracy.

[Morten Welinder <mortenw src gnome org>]

commit ba681f01aed263b2ded18663a65e1c6f7a0c38b9
Author: Morten Welinder <terra gnome org>
Date:   Fri Apr 12 14:16:28 2013 -0400

    FACT: Improve accuracy.
    
    Chained multiply actually will introduce last-bit rounding errors for
    relatively small n.  Discovered for n=29.
    
    For n larger than our table size, rewrite using POCHHAMMER so we avoid
    taking exp(lgamma(...)) which is fairly bad for accuracy.

 ChangeLog      |    1 +
 NEWS           |    1 +
 src/mathfunc.c |  121 ++++++++++++++++++++++++++++++++++++++++++++++++++------
 3 files changed, 111 insertions(+), 12 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 139dbc2..2b3c395 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -4,6 +4,7 @@
        #697850.
        (pnbinom): Allow prob==1 (already fixed in R) and n==0 (now
        reported to R).
+       (fact): Use table for n<=1; rewrite using pochhammer otherwise.
 
 2013-04-11  Morten Welinder  <terra gnome org>
 
diff --git a/NEWS b/NEWS
index fb7bba8..f4bbc13 100644
--- a/NEWS
+++ b/NEWS
@@ -36,6 +36,7 @@ Morten:
        * Fix parsing of 'Sheet1!#REF!' as used by XL.  [#683494]
        * New function NT_OMEGA.
        * Fix minor R.PBINOM problems.
+       * Improve FACT and GAMMA accuracy.
 
 --------------------------------------------------------------------------
 Gnumeric 1.12.1
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 5b52a7e..7f52322 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -6380,23 +6380,120 @@ permut (gnm_float n, gnm_float k)
 gnm_float
 fact (int n)
 {
-       static gnm_float table[100];
-       static gboolean init = FALSE;
+       static const gnm_float table[] = {
+               GNM_const(1.0),
+               GNM_const(1.0),
+               GNM_const(2.0),
+               GNM_const(6.0),
+               GNM_const(24.0),
+               GNM_const(120.0),
+               GNM_const(720.0),
+               GNM_const(5040.0),
+               GNM_const(40320.0),
+               GNM_const(362880.0),
+               GNM_const(3628800.0),
+               GNM_const(39916800.0),
+               GNM_const(479001600.0),
+               GNM_const(6227020800.0),
+               GNM_const(87178291200.0),
+               GNM_const(1307674368000.0),
+               GNM_const(20922789888000.0),
+               GNM_const(355687428096000.0),
+               GNM_const(6402373705728000.0),
+               GNM_const(121645100408832000.0),
+               GNM_const(2432902008176640000.0),
+               GNM_const(51090942171709440000.0),
+               GNM_const(1124000727777607680000.0),
+               GNM_const(25852016738884976640000.0),
+               GNM_const(620448401733239439360000.0),
+               GNM_const(15511210043330985984000000.0),
+               GNM_const(403291461126605635584000000.0),
+               GNM_const(10888869450418352160768000000.0),
+               GNM_const(304888344611713860501504000000.0),
+               GNM_const(8841761993739701954543616000000.0),
+               GNM_const(265252859812191058636308480000000.0),
+               GNM_const(8222838654177922817725562880000000.0),
+               GNM_const(263130836933693530167218012160000000.0),
+               GNM_const(8683317618811886495518194401280000000.0),
+               GNM_const(295232799039604140847618609643520000000.0),
+               GNM_const(10333147966386144929666651337523200000000.0),
+               GNM_const(371993326789901217467999448150835200000000.0),
+               GNM_const(13763753091226345046315979581580902400000000.0),
+               GNM_const(523022617466601111760007224100074291200000000.0),
+               GNM_const(20397882081197443358640281739902897356800000000.0),
+               GNM_const(815915283247897734345611269596115894272000000000.0),
+               GNM_const(33452526613163807108170062053440751665152000000000.0),
+               GNM_const(1405006117752879898543142606244511569936384000000000.0),
+               GNM_const(60415263063373835637355132068513997507264512000000000.0),
+               GNM_const(2658271574788448768043625811014615890319638528000000000.0),
+               GNM_const(119622220865480194561963161495657715064383733760000000000.0),
+               GNM_const(5502622159812088949850305428800254892961651752960000000000.0),
+               GNM_const(258623241511168180642964355153611979969197632389120000000000.0),
+               GNM_const(12413915592536072670862289047373375038521486354677760000000000.0),
+               GNM_const(608281864034267560872252163321295376887552831379210240000000000.0),
+               GNM_const(30414093201713378043612608166064768844377641568960512000000000000.0),
+               GNM_const(1551118753287382280224243016469303211063259720016986112000000000000.0),
+               GNM_const(80658175170943878571660636856403766975289505440883277824000000000000.0),
+               GNM_const(4274883284060025564298013753389399649690343788366813724672000000000000.0),
+               GNM_const(230843697339241380472092742683027581083278564571807941132288000000000000.0),
+               GNM_const(12696403353658275925965100847566516959580321051449436762275840000000000000.0),
+               GNM_const(710998587804863451854045647463724949736497978881168458687447040000000000000.0),
+               GNM_const(40526919504877216755680601905432322134980384796226602145184481280000000000000.0),
+               GNM_const(2350561331282878571829474910515074683828862318181142924420699914240000000000000.0),
+               
GNM_const(138683118545689835737939019720389406345902876772687432540821294940160000000000000.0),
+               
GNM_const(8320987112741390144276341183223364380754172606361245952449277696409600000000000000.0),
+               
GNM_const(507580213877224798800856812176625227226004528988036003099405939480985600000000000000.0),
+               
GNM_const(31469973260387937525653122354950764088012280797258232192163168247821107200000000000000.0),
+               
GNM_const(1982608315404440064116146708361898137544773690227268628106279599612729753600000000000000.0),
+               
GNM_const(126886932185884164103433389335161480802865516174545192198801894375214704230400000000000000.0),
+               
GNM_const(8247650592082470666723170306785496252186258551345437492922123134388955774976000000000000000.0),
+               
GNM_const(544344939077443064003729240247842752644293064388798874532860126869671081148416000000000000000.0),
+               
GNM_const(36471110918188685288249859096605464427167635314049524593701628500267962436943872000000000000000.0),
+               
GNM_const(2480035542436830599600990418569171581047399201355367672371710738018221445712183296000000000000000.0),
+               
GNM_const(171122452428141311372468338881272839092270544893520369393648040923257279754140647424000000000000000.0),
+               
GNM_const(11978571669969891796072783721689098736458938142546425857555362864628009582789845319680000000000000000.0),
+               
GNM_const(850478588567862317521167644239926010288584608120796235886430763388588680378079017697280000000000000000.0),
+               
GNM_const(61234458376886086861524070385274672740778091784697328983823014963978384987221689274204160000000000000000.0),
+               
GNM_const(4470115461512684340891257138125051110076800700282905015819080092370422104067183317016903680000000000000000.0),
+               
GNM_const(330788544151938641225953028221253782145683251820934971170611926835411235700971565459250872320000000000000000.0),
+               
GNM_const(24809140811395398091946477116594033660926243886570122837795894512655842677572867409443815424000000000000000000.0),
+               
GNM_const(1885494701666050254987932260861146558230394535379329335672487982961844043495537923117729972224000000000000000000.0),
+               
GNM_const(145183092028285869634070784086308284983740379224208358846781574688061991349156420080065207861248000000000000000000.0),
+               
GNM_const(11324281178206297831457521158732046228731749579488251990048962825668835325234200766245086213177344000000000000000000.0),
+               
GNM_const(894618213078297528685144171539831652069808216779571907213868063227837990693501860533361810841010176000000000000000000.0),
+               
GNM_const(71569457046263802294811533723186532165584657342365752577109445058227039255480148842668944867280814080000000000000000000.0),
+               
GNM_const(5797126020747367985879734231578109105412357244731625958745865049716390179693892056256184534249745940480000000000000000000.0),
+               
GNM_const(475364333701284174842138206989404946643813294067993328617160934076743994734899148613007131808479167119360000000000000000000.0),
+               
GNM_const(39455239697206586511897471180120610571436503407643446275224357528369751562996629334879591940103770870906880000000000000000000.0),
+               
GNM_const(3314240134565353266999387579130131288000666286242049487118846032383059131291716864129885722968716753156177920000000000000000000.0),
+               
GNM_const(281710411438055027694947944226061159480056634330574206405101912752560026159795933451040286452340924018275123200000000000000000000.0),
+               
GNM_const(24227095383672732381765523203441259715284870552429381750838764496720162249742450276789464634901319465571660595200000000000000000000.0),
+               
GNM_const(2107757298379527717213600518699389595229783738061356212322972511214654115727593174080683423236414793504734471782400000000000000000000.0),
+               
GNM_const(185482642257398439114796845645546284380220968949399346684421580986889562184028199319100141244804501828416633516851200000000000000000000.0),
+               
GNM_const(16507955160908461081216919262453619309839666236496541854913520707833171034378509739399912570787600662729080382999756800000000000000000000.0),
+               
GNM_const(1485715964481761497309522733620825737885569961284688766942216863704985393094065876545992131370884059645617234469978112000000000000000000000.0),
+               
GNM_const(135200152767840296255166568759495142147586866476906677791741734597153670771559994765685283954750449427751168336768008192000000000000000000000.0),
+               
GNM_const(12438414054641307255475324325873553077577991715875414356840239582938137710983519518443046123837041347353107486982656753664000000000000000000000.0),
+               
GNM_const(1156772507081641574759205162306240436214753229576413535186142281213246807121467315215203289516844845303838996289387078090752000000000000000000000.0),
+               
GNM_const(108736615665674308027365285256786601004186803580182872307497374434045199869417927630229109214583415458560865651202385340530688000000000000000000000.0),
+               
GNM_const(10329978488239059262599702099394727095397746340117372869212250571234293987594703124871765375385424468563282236864226607350415360000000000000000000000.0),
+               
GNM_const(991677934870949689209571401541893801158183648651267795444376054838492222809091499987689476037000748982075094738965754305639874560000000000000000000000.0),
+               
GNM_const(96192759682482119853328425949563698712343813919172976158104477319333745612481875498805879175589072651261284189679678167647067832320000000000000000000000.0),
+               
GNM_const(9426890448883247745626185743057242473809693764078951663494238777294707070023223798882976159207729119823605850588608460429412647567360000000000000000000000.0),
+               
GNM_const(933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000.0),
+               
GNM_const(93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000.0)
+       };
+       const int N = G_N_ELEMENTS (table);
 
        if (n < 0)
                return gnm_nan;
 
-       if (n < (int)G_N_ELEMENTS (table)) {
-               if (!init) {
-                       int i;
-                       table[0] = 1;
-                       for (i = 1; i < (int)G_N_ELEMENTS (table); i++)
-                               table[i] = table[i - 1] * i;
-                       init = TRUE;
-               }
+       if (n < N)
                return table[n];
-       } else
-               return gnm_floor (0.5 + gnm_exp (gnm_lgamma (n + 1)));
+       else {
+               gnm_float f = pochhammer (N + 1, n + 1 - N, FALSE);
+               return f * table[N - 1];
+       }
 }
 
 gnm_float