[gnumeric] FACT, etc: improve precision.
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gnumeric] FACT, etc: improve precision.
- Date: Fri, 1 Nov 2013 02:22:40 +0000 (UTC)
commit 308b78c786928bfb5538a540eda45c729b036489
Author: Morten Welinder <terra gnome org>
Date: Thu Oct 31 19:31:37 2013 -0400
FACT, etc: improve precision.
This goes slightly insane by tabulating 0! through 50000! in the form
of a quad precision mantissa [0.5;1[ and a power of 2.
That, in turn, makes things like combin fairly easy and spot on.
ChangeLog | 4 +
NEWS | 1 +
src/mathfunc.c | 269 +++++++++++++++++++++++++++++---------------------------
src/mathfunc.h | 1 +
src/numbers.h | 2 +
5 files changed, 149 insertions(+), 128 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 1a52928..2189ed0 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -7,6 +7,10 @@
(qtukey): Ditto.
(qbeta): Ditto.
(J_bessel): Extend xlrg_BESS_IJ to match current R.
+ (permut): Route this into pochhammer.
+ (qfact): New function. Tabulate 0! through 50000! using quad
+ precision arithmetic.
+ (combin, fact, pochhammer): Improve accuracy using qfact.
2013-10-22 Morten Welinder <terra gnome org>
diff --git a/NEWS b/NEWS
index a0e41b9..6c793b7 100644
--- a/NEWS
+++ b/NEWS
@@ -7,6 +7,7 @@ Morten:
* Acquire more special function test cases.
* Improve accuracy of R.QGAMMA and thus R.QCHISQ.
* Improve accuracy of R.QBETA, R.QF, R.QTUKEY, R.QSNORM, and R.QST.
+ * Improve accuracy of COMBIN, PERMUT, POCHHAMMER, FACT, GAMMA.
Xabier RodrÃguez Calvar:
* Fix dialog button order. [#710378]
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 1273ddf..7b4db5a 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -7973,145 +7973,128 @@ pbinom2 (gnm_float x0, gnm_float x1, gnm_float n, gnm_float p)
gnm_float
combin (gnm_float n, gnm_float k)
{
- /* Assume integer arguments. */
+ GnmQuad m1, m2, m3;
+ int e1, e2, e3;
+ gboolean ok;
- if (k < 0 || k > n)
+ if (k < 0 || k > n || n != gnm_floor (n) || k != gnm_floor (k))
return gnm_nan;
- else if (n >= 15)
- return gnm_floor (0.5 + gnm_exp (gnm_lgamma (n + 1) - gnm_lgamma (k + 1) - gnm_lgamma (n - k
+ 1)));
- else
- return fact (n) / fact (k) / fact (n - k);
+
+ k = MIN (k, n - k);
+ if (k == 0)
+ return 1;
+ if (k == 1)
+ return n;
+
+ ok = (n < G_MAXINT &&
+ !qfact ((int)n, &m1, &e1) &&
+ !qfact ((int)k, &m2, &e2) &&
+ !qfact ((int)(n - k), &m3, &e3));
+
+ if (ok) {
+ void *state = gnm_quad_start ();
+ gnm_float c;
+ gnm_quad_mul (&m2, &m2, &m3);
+ gnm_quad_div (&m1, &m1, &m2);
+ c = gnm_ldexp (gnm_quad_value (&m1), e1 - e2 - e3);
+ gnm_quad_end (state);
+ return c;
+ }
+
+ if (k < 30) {
+ void *state = gnm_quad_start ();
+ GnmQuad p, a, b;
+ gnm_float c;
+ int i;
+
+ gnm_quad_init (&p, 1);
+ for (i = 0; i < k; i++) {
+ gnm_quad_init (&a, n - i);
+ gnm_quad_mul (&p, &p, &a);
+
+ gnm_quad_init (&b, i + 1);
+ gnm_quad_div (&p, &p, &b);
+ }
+
+ c = gnm_quad_value (&p);
+ gnm_quad_end (state);
+ return c;
+ }
+
+ {
+ gnm_float lp = pochhammer (n - k + 1, k, TRUE) -
+ gnm_lgamma (k + 1);
+ return gnm_floor (0.5 + gnm_exp (lp));
+ }
}
gnm_float
permut (gnm_float n, gnm_float k)
{
- /* Assume integer arguments. */
- if (k < 0 || k > n)
+ if (k < 0 || k > n || n != gnm_floor (n) || k != gnm_floor (k))
return gnm_nan;
- else if (n >= 15)
- return gnm_floor (0.5 + gnm_exp (gnm_lgamma (n + 1) - gnm_lgamma (n - k + 1)));
- else
- return fact (n) / fact (n - k);
+
+ return pochhammer (n - k + 1, k, FALSE);
}
+gboolean
+qfact (int n, GnmQuad *mant, int *exp2)
+{
+ static GnmQuad mants[50000];
+ static int exp2s[G_N_ELEMENTS (mants)];
+ static int init = 0;
+
+ if (n < 0 || n >= (int)G_N_ELEMENTS(mants)) {
+ *mant = gnm_quad_zero;
+ *exp2 = 0;
+ return TRUE;
+ }
+
+ if (n >= init) {
+ void *state = gnm_quad_start ();
+
+ if (init == 0) {
+ gnm_quad_init (&mants[0], 0.5);
+ exp2s[0] = 1;
+ init++;
+ }
+
+ while (n >= init) {
+ GnmQuad m, s;
+ int e;
+
+ gnm_quad_init (&m, init);
+ gnm_quad_mul (&m, &m, &mants[init - 1]);
+ (void)gnm_frexp (gnm_quad_value (&m), &e);
+ exp2s[init] = exp2s[init - 1] + e;
+ gnm_quad_init (&s, gnm_ldexp (1.0, -e));
+ gnm_quad_mul (&mants[init], &m, &s);
+
+ init++;
+ }
+
+ gnm_quad_end (state);
+ }
+
+ *mant = mants[n];
+ *exp2 = exp2s[n];
+ return FALSE;
+}
+
+
gnm_float
fact (int n)
{
- 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);
+ GnmQuad r;
+ int e;
if (n < 0)
return gnm_nan;
- if (n < N)
- return table[n];
- else {
- gnm_float f = pochhammer (N + 1, n + 1 - N, FALSE);
- return f * table[N - 1];
- }
+ if (qfact (n, &r, &e))
+ return pochhammer (1, n, FALSE);
+
+ return ldexp (gnm_quad_value (&r), e);
}
gnm_float
@@ -8820,11 +8803,41 @@ pochhammer (gnm_float x, gnm_float n, gboolean give_log)
if (x <= 0 || x + n <= 0)
return gnm_nan;
- if (n == gnm_floor (n) && n >= 0 && n <= 40 && !give_log) {
- gnm_float r;
- for (r = 1; n-- > 0; x++)
- r *= x;
- return r;
+ if (n == gnm_floor (n) && n >= 0) {
+ if (x == gnm_floor (x) && x + n < G_MAXINT) {
+ void *state = gnm_quad_start ();
+ GnmQuad m1, m2;
+ int e1, e2;
+ gboolean ok;
+ gnm_float r = 0;
+
+ ok = (!qfact ((int)(x + n - 1), &m1, &e1) &&
+ !qfact ((int)(x - 1), &m2, &e2));
+
+ if (ok) {
+ int de = e1 - e2;
+
+ gnm_quad_div (&m1, &m1, &m2);
+ r = gnm_quad_value (&m1);
+
+ r = give_log
+ ? gnm_log (r) + M_LN2gnum * de
+ : gnm_ldexp (r, de);
+ ok = gnm_finite (r);
+ }
+
+ gnm_quad_end (state);
+
+ if (ok)
+ return r;
+ }
+
+ if (n <= 40 && !give_log) {
+ gnm_float r;
+ for (r = 1; n-- > 0; x++)
+ r *= x;
+ return r;
+ }
}
if (give_log) {
@@ -8857,7 +8870,7 @@ gnm_gamma (gnm_float x)
return x;
if (x > 0) {
- if (x == gnm_floor (x) && x <= 100)
+ if (x == gnm_floor (x) && x < G_MAXINT)
return fact ((int)x - 1);
if (x > 10) {
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 290f436..dd360cd 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -179,6 +179,7 @@ gboolean gnm_matrix_eigen (GnmMatrix const *m, GnmMatrix *EIG, gnm_float *eigenv
gnm_float combin (gnm_float n, gnm_float k);
gnm_float permut (gnm_float n, gnm_float k);
gnm_float fact (int n);
+gboolean qfact (int n, GnmQuad *mant, int *exp2);
gint gnm_float_equal (gnm_float const *a, const gnm_float *b);
guint gnm_float_hash (gnm_float const *d);
diff --git a/src/numbers.h b/src/numbers.h
index c8bd365..c88c7c2 100644
--- a/src/numbers.h
+++ b/src/numbers.h
@@ -121,6 +121,7 @@ gnm_float gnm_erfc (gnm_float x);
#define gnm_quad_mul go_quad_mull
#define gnm_quad_div go_quad_divl
#define gnm_quad_mul12 go_quad_mul12l
+#define gnm_quad_zero go_quad_zerol
#define GnmAccumulator GOAccumulatorl
#define gnm_accumulator_start go_accumulator_startl
#define gnm_accumulator_end go_accumulator_endl
@@ -208,6 +209,7 @@ typedef double gnm_float;
#define gnm_quad_mul go_quad_mul
#define gnm_quad_div go_quad_div
#define gnm_quad_mul12 go_quad_mul12
+#define gnm_quad_zero go_quad_zero
#define GnmAccumulator GOAccumulator
#define gnm_accumulator_start go_accumulator_start
#define gnm_accumulator_end go_accumulator_end
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
