[gnumeric] Stats: accuracy improvements.
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gnumeric] Stats: accuracy improvements.
- Date: Mon, 13 Jan 2014 20:58:20 +0000 (UTC)
commit 204591361d1cd979d45ccab8be64e37207bd7370
Author: Morten Welinder <terra gnome org>
Date: Mon Jan 13 15:57:58 2014 -0500
Stats: accuracy improvements.
ChangeLog | 14 ++
NEWS | 2 +
plugins/fn-stat/ChangeLog | 4 +
plugins/fn-stat/functions.c | 9 +-
src/mathfunc.c | 386 +++++++++++++++++++++++++++++++++++++------
src/mathfunc.h | 1 +
src/sf-dpq.c | 18 +--
7 files changed, 358 insertions(+), 76 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index d9d8642..ddfea90 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,17 @@
+2014-01-12 Morten Welinder <terra gnome org>
+
+ * src/mathfunc.c (expmx2h): Split out from dnorm.
+ (dhyper): Round the chosen "p" to float to improve accuracy of
+ "q".
+
+2014-01-11 Morten Welinder <terra gnome org>
+
+ * src/mathfunc.c (pow1p): Improve accuracy.
+ (R_D_nonint): Fix rounding problem.
+ (ebd0): New extended version of bd0.
+ (dpois_raw, dbinom_raw): Use ebd0.
+ (dbinom_raw): Handle x==n and x==0 much better.
+
2014-01-08 Morten Welinder <terra gnome org>
* src/wbc-gtk-actions.c (cb_file_sendto): Simplify using
diff --git a/NEWS b/NEWS
index 91266ac..da03491 100644
--- a/NEWS
+++ b/NEWS
@@ -24,6 +24,8 @@ Morten:
* Fix solver problem with constraints.
* Improve accuracy of R.PNORM and R.PLNORM for large arguments.
* Improve accuracy of R.QCAUCHY.
+ * Improve accuracy of R.DHYPER, R.DBINOM, R.DNBINOM, R.DPOIS,
+ RAYLEIGH.
--------------------------------------------------------------------------
Gnumeric 1.12.9
diff --git a/plugins/fn-stat/ChangeLog b/plugins/fn-stat/ChangeLog
index ce410ff..ac3d061 100644
--- a/plugins/fn-stat/ChangeLog
+++ b/plugins/fn-stat/ChangeLog
@@ -1,3 +1,7 @@
+2014-01-12 Morten Welinder <terra gnome org>
+
+ * functions.c (random_rayleigh_pdf): Use expmx2h.
+
2013-11-28 Morten Welinder <terra gnome org>
* Release 1.12.9
diff --git a/plugins/fn-stat/functions.c b/plugins/fn-stat/functions.c
index c261f8f..1b4e532 100644
--- a/plugins/fn-stat/functions.c
+++ b/plugins/fn-stat/functions.c
@@ -4541,12 +4541,14 @@ static GnmFuncHelp const help_rayleigh[] = {
static gnm_float
random_rayleigh_pdf (gnm_float x, gnm_float sigma)
{
+ if (sigma <= 0)
+ return gnm_nan;
+
if (x < 0)
return 0;
else {
gnm_float u = x / sigma;
-
- return (u / sigma) * gnm_exp (-u * u / 2.0);
+ return (u / sigma) * expmx2h (u);
}
}
@@ -4556,9 +4558,6 @@ gnumeric_rayleigh (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
gnm_float x = value_get_as_float (argv[0]);
gnm_float sigma = value_get_as_float (argv[1]);
- if (sigma <= 0)
- return value_new_error_NUM (ei->pos);
-
return value_new_float (random_rayleigh_pdf (x, sigma));
}
diff --git a/src/mathfunc.c b/src/mathfunc.c
index d820242..1b5ca50 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -77,12 +77,238 @@ mathfunc_init (void)
/* Nothing, for the time being. */
}
+/*
+ * A table of logs for scaling purposes. Each value has four parts with
+ * 23 bits in each. That means each part can be multiplied by a double
+ * with at most 30 bits set and not have any rounding error. Note, that
+ * the first entry is log(2).
+ *
+ * Entry i is associated with the value r = 0.5 + i / 256.0. The
+ * argument to log is p/q where q=1024 and p=floor(q / r + 0.5).
+ * Thus r*p/q is close to 1.
+ */
+static const float bd0_scale[128 + 1][4] = {
+ { +0x1.62e430p-1, -0x1.05c610p-29, -0x1.950d88p-54, +0x1.d9cc02p-79 }, /* 128: log(2048/1024.) */
+ { +0x1.5ee02cp-1, -0x1.6dbe98p-25, -0x1.51e540p-50, +0x1.2bfa48p-74 }, /* 129: log(2032/1024.) */
+ { +0x1.5ad404p-1, +0x1.86b3e4p-26, +0x1.9f6534p-50, +0x1.54be04p-74 }, /* 130: log(2016/1024.) */
+ { +0x1.570124p-1, -0x1.9ed750p-25, -0x1.f37dd0p-51, +0x1.10b770p-77 }, /* 131: log(2001/1024.) */
+ { +0x1.5326e4p-1, -0x1.9b9874p-25, -0x1.378194p-49, +0x1.56feb2p-74 }, /* 132: log(1986/1024.) */
+ { +0x1.4f4528p-1, +0x1.aca70cp-28, +0x1.103e74p-53, +0x1.9c410ap-81 }, /* 133: log(1971/1024.) */
+ { +0x1.4b5bd8p-1, -0x1.6a91d8p-25, -0x1.8e43d0p-50, -0x1.afba9ep-77 }, /* 134: log(1956/1024.) */
+ { +0x1.47ae54p-1, -0x1.abb51cp-25, +0x1.19b798p-51, +0x1.45e09cp-76 }, /* 135: log(1942/1024.) */
+ { +0x1.43fa00p-1, -0x1.d06318p-25, -0x1.8858d8p-49, -0x1.1927c4p-75 }, /* 136: log(1928/1024.) */
+ { +0x1.3ffa40p-1, +0x1.1a427cp-25, +0x1.151640p-53, -0x1.4f5606p-77 }, /* 137: log(1913/1024.) */
+ { +0x1.3c7c80p-1, -0x1.19bf48p-34, +0x1.05fc94p-58, -0x1.c096fcp-82 }, /* 138: log(1900/1024.) */
+ { +0x1.38b320p-1, +0x1.6b5778p-25, +0x1.be38d0p-50, -0x1.075e96p-74 }, /* 139: log(1886/1024.) */
+ { +0x1.34e288p-1, +0x1.d9ce1cp-25, +0x1.316eb8p-49, +0x1.2d885cp-73 }, /* 140: log(1872/1024.) */
+ { +0x1.315124p-1, +0x1.c2fc60p-29, -0x1.4396fcp-53, +0x1.acf376p-78 }, /* 141: log(1859/1024.) */
+ { +0x1.2db954p-1, +0x1.720de4p-25, -0x1.d39b04p-49, -0x1.f11176p-76 }, /* 142: log(1846/1024.) */
+ { +0x1.2a1b08p-1, -0x1.562494p-25, +0x1.a7863cp-49, +0x1.85dd64p-73 }, /* 143: log(1833/1024.) */
+ { +0x1.267620p-1, +0x1.3430e0p-29, -0x1.96a958p-56, +0x1.f8e636p-82 }, /* 144: log(1820/1024.) */
+ { +0x1.23130cp-1, +0x1.7bebf4p-25, +0x1.416f1cp-52, -0x1.78dd36p-77 }, /* 145: log(1808/1024.) */
+ { +0x1.1faa34p-1, +0x1.70e128p-26, +0x1.81817cp-50, -0x1.c2179cp-76 }, /* 146: log(1796/1024.) */
+ { +0x1.1bf204p-1, +0x1.3a9620p-28, +0x1.2f94c0p-52, +0x1.9096c0p-76 }, /* 147: log(1783/1024.) */
+ { +0x1.187ce4p-1, -0x1.077870p-27, +0x1.655a80p-51, +0x1.eaafd6p-78 }, /* 148: log(1771/1024.) */
+ { +0x1.1501c0p-1, -0x1.406cacp-25, -0x1.e72290p-49, +0x1.5dd800p-73 }, /* 149: log(1759/1024.) */
+ { +0x1.11cb80p-1, +0x1.787cd0p-25, -0x1.efdc78p-51, -0x1.5380cep-77 }, /* 150: log(1748/1024.) */
+ { +0x1.0e4498p-1, +0x1.747324p-27, -0x1.024548p-51, +0x1.77a5a6p-75 }, /* 151: log(1736/1024.) */
+ { +0x1.0b036cp-1, +0x1.690c74p-25, +0x1.5d0cc4p-50, -0x1.c0e23cp-76 }, /* 152: log(1725/1024.) */
+ { +0x1.077070p-1, -0x1.a769bcp-27, +0x1.452234p-52, +0x1.6ba668p-76 }, /* 153: log(1713/1024.) */
+ { +0x1.04240cp-1, -0x1.a686acp-27, -0x1.ef46b0p-52, -0x1.5ce10cp-76 }, /* 154: log(1702/1024.) */
+ { +0x1.00d22cp-1, +0x1.fc0e10p-25, +0x1.6ee034p-50, -0x1.19a2ccp-74 }, /* 155: log(1691/1024.) */
+ { +0x1.faf588p-2, +0x1.ef1e64p-27, -0x1.26504cp-54, -0x1.b15792p-82 }, /* 156: log(1680/1024.) */
+ { +0x1.f4d87cp-2, +0x1.d7b980p-26, -0x1.a114d8p-50, +0x1.9758c6p-75 }, /* 157: log(1670/1024.) */
+ { +0x1.ee1414p-2, +0x1.2ec060p-26, +0x1.dc00fcp-52, +0x1.f8833cp-76 }, /* 158: log(1659/1024.) */
+ { +0x1.e7e32cp-2, -0x1.ac796cp-27, -0x1.a68818p-54, +0x1.235d02p-78 }, /* 159: log(1649/1024.) */
+ { +0x1.e108a0p-2, -0x1.768ba4p-28, -0x1.f050a8p-52, +0x1.00d632p-82 }, /* 160: log(1638/1024.) */
+ { +0x1.dac354p-2, -0x1.d3a6acp-30, +0x1.18734cp-57, -0x1.f97902p-83 }, /* 161: log(1628/1024.) */
+ { +0x1.d47424p-2, +0x1.7dbbacp-31, -0x1.d5ada4p-56, +0x1.56fcaap-81 }, /* 162: log(1618/1024.) */
+ { +0x1.ce1af0p-2, +0x1.70be7cp-27, +0x1.6f6fa4p-51, +0x1.7955a2p-75 }, /* 163: log(1608/1024.) */
+ { +0x1.c7b798p-2, +0x1.ec36ecp-26, -0x1.07e294p-50, -0x1.ca183cp-75 }, /* 164: log(1598/1024.) */
+ { +0x1.c1ef04p-2, +0x1.c1dfd4p-26, +0x1.888eecp-50, -0x1.fd6b86p-75 }, /* 165: log(1589/1024.) */
+ { +0x1.bb7810p-2, +0x1.478bfcp-26, +0x1.245b8cp-50, +0x1.ea9d52p-74 }, /* 166: log(1579/1024.) */
+ { +0x1.b59da0p-2, -0x1.882b08p-27, +0x1.31573cp-53, -0x1.8c249ap-77 }, /* 167: log(1570/1024.) */
+ { +0x1.af1294p-2, -0x1.b710f4p-27, +0x1.622670p-51, +0x1.128578p-76 }, /* 168: log(1560/1024.) */
+ { +0x1.a925d4p-2, -0x1.0ae750p-27, +0x1.574ed4p-51, +0x1.084996p-75 }, /* 169: log(1551/1024.) */
+ { +0x1.a33040p-2, +0x1.027d30p-29, +0x1.b9a550p-53, -0x1.b2e38ap-78 }, /* 170: log(1542/1024.) */
+ { +0x1.9d31c0p-2, -0x1.5ec12cp-26, -0x1.5245e0p-52, +0x1.2522d0p-79 }, /* 171: log(1533/1024.) */
+ { +0x1.972a34p-2, +0x1.135158p-30, +0x1.a5c09cp-56, +0x1.24b70ep-80 }, /* 172: log(1524/1024.) */
+ { +0x1.911984p-2, +0x1.0995d4p-26, +0x1.3bfb5cp-50, +0x1.2c9dd6p-75 }, /* 173: log(1515/1024.) */
+ { +0x1.8bad98p-2, -0x1.1d6144p-29, +0x1.5b9208p-53, +0x1.1ec158p-77 }, /* 174: log(1507/1024.) */
+ { +0x1.858b58p-2, -0x1.1b4678p-27, +0x1.56cab4p-53, -0x1.2fdc0cp-78 }, /* 175: log(1498/1024.) */
+ { +0x1.7f5fa0p-2, +0x1.3aaf48p-27, +0x1.461964p-51, +0x1.4ae476p-75 }, /* 176: log(1489/1024.) */
+ { +0x1.79db68p-2, -0x1.7e5054p-26, +0x1.673750p-51, -0x1.a11f7ap-76 }, /* 177: log(1481/1024.) */
+ { +0x1.744f88p-2, -0x1.cc0e18p-26, -0x1.1e9d18p-50, -0x1.6c06bcp-78 }, /* 178: log(1473/1024.) */
+ { +0x1.6e08ecp-2, -0x1.5d45e0p-26, -0x1.c73ec8p-50, +0x1.318d72p-74 }, /* 179: log(1464/1024.) */
+ { +0x1.686c80p-2, +0x1.e9b14cp-26, -0x1.13bbd4p-50, -0x1.efeb1cp-78 }, /* 180: log(1456/1024.) */
+ { +0x1.62c830p-2, -0x1.a8c70cp-27, -0x1.5a1214p-51, -0x1.bab3fcp-79 }, /* 181: log(1448/1024.) */
+ { +0x1.5d1bdcp-2, -0x1.4fec6cp-31, +0x1.423638p-56, +0x1.ee3feep-83 }, /* 182: log(1440/1024.) */
+ { +0x1.576770p-2, +0x1.7455a8p-26, -0x1.3ab654p-50, -0x1.26be4cp-75 }, /* 183: log(1432/1024.) */
+ { +0x1.5262e0p-2, -0x1.146778p-26, -0x1.b9f708p-52, -0x1.294018p-77 }, /* 184: log(1425/1024.) */
+ { +0x1.4c9f08p-2, +0x1.e152c4p-26, -0x1.dde710p-53, +0x1.fd2208p-77 }, /* 185: log(1417/1024.) */
+ { +0x1.46d2d8p-2, +0x1.c28058p-26, -0x1.936284p-50, +0x1.9fdd68p-74 }, /* 186: log(1409/1024.) */
+ { +0x1.41b940p-2, +0x1.cce0c0p-26, -0x1.1a4050p-50, +0x1.bc0376p-76 }, /* 187: log(1402/1024.) */
+ { +0x1.3bdd24p-2, +0x1.d6296cp-27, +0x1.425b48p-51, -0x1.cddb2cp-77 }, /* 188: log(1394/1024.) */
+ { +0x1.36b578p-2, -0x1.287ddcp-27, -0x1.2d0f4cp-51, +0x1.38447ep-75 }, /* 189: log(1387/1024.) */
+ { +0x1.31871cp-2, +0x1.2a8830p-27, +0x1.3eae54p-52, -0x1.898136p-77 }, /* 190: log(1380/1024.) */
+ { +0x1.2b9304p-2, -0x1.51d8b8p-28, +0x1.27694cp-52, -0x1.fd852ap-76 }, /* 191: log(1372/1024.) */
+ { +0x1.265620p-2, -0x1.d98f3cp-27, +0x1.a44338p-51, -0x1.56e85ep-78 }, /* 192: log(1365/1024.) */
+ { +0x1.211254p-2, +0x1.986160p-26, +0x1.73c5d0p-51, +0x1.4a861ep-75 }, /* 193: log(1358/1024.) */
+ { +0x1.1bc794p-2, +0x1.fa3918p-27, +0x1.879c5cp-51, +0x1.16107cp-78 }, /* 194: log(1351/1024.) */
+ { +0x1.1675ccp-2, -0x1.4545a0p-26, +0x1.c07398p-51, +0x1.f55c42p-76 }, /* 195: log(1344/1024.) */
+ { +0x1.111ce4p-2, +0x1.f72670p-37, -0x1.b84b5cp-61, +0x1.a4a4dcp-85 }, /* 196: log(1337/1024.) */
+ { +0x1.0c81d4p-2, +0x1.0c150cp-27, +0x1.218600p-51, -0x1.d17312p-76 }, /* 197: log(1331/1024.) */
+ { +0x1.071b84p-2, +0x1.fcd590p-26, +0x1.a3a2e0p-51, +0x1.fe5ef8p-76 }, /* 198: log(1324/1024.) */
+ { +0x1.01ade4p-2, -0x1.bb1844p-28, +0x1.db3cccp-52, +0x1.1f56fcp-77 }, /* 199: log(1317/1024.) */
+ { +0x1.fa01c4p-3, -0x1.12a0d0p-29, -0x1.f71fb0p-54, +0x1.e287a4p-78 }, /* 200: log(1311/1024.) */
+ { +0x1.ef0adcp-3, +0x1.7b8b28p-28, -0x1.35bce4p-52, -0x1.abc8f8p-79 }, /* 201: log(1304/1024.) */
+ { +0x1.e598ecp-3, +0x1.5a87e4p-27, -0x1.134bd0p-51, +0x1.c2cebep-76 }, /* 202: log(1298/1024.) */
+ { +0x1.da85d8p-3, -0x1.df31b0p-27, +0x1.94c16cp-57, +0x1.8fd7eap-82 }, /* 203: log(1291/1024.) */
+ { +0x1.d0fb80p-3, -0x1.bb5434p-28, -0x1.ea5640p-52, -0x1.8ceca4p-77 }, /* 204: log(1285/1024.) */
+ { +0x1.c765b8p-3, +0x1.e4d68cp-27, +0x1.5b59b4p-51, +0x1.76f6c4p-76 }, /* 205: log(1279/1024.) */
+ { +0x1.bdc46cp-3, -0x1.1cbb50p-27, +0x1.2da010p-51, +0x1.eb282cp-75 }, /* 206: log(1273/1024.) */
+ { +0x1.b27980p-3, -0x1.1b9ce0p-27, +0x1.7756f8p-52, +0x1.2ff572p-76 }, /* 207: log(1266/1024.) */
+ { +0x1.a8bed0p-3, -0x1.bbe874p-30, +0x1.85cf20p-56, +0x1.b9cf18p-80 }, /* 208: log(1260/1024.) */
+ { +0x1.9ef83cp-3, +0x1.2769a4p-27, -0x1.85bda0p-52, +0x1.8c8018p-79 }, /* 209: log(1254/1024.) */
+ { +0x1.9525a8p-3, +0x1.cf456cp-27, -0x1.7137d8p-52, -0x1.f158e8p-76 }, /* 210: log(1248/1024.) */
+ { +0x1.8b46f8p-3, +0x1.11b12cp-30, +0x1.9f2104p-54, -0x1.22836ep-78 }, /* 211: log(1242/1024.) */
+ { +0x1.83040cp-3, +0x1.2379e4p-28, +0x1.b71c70p-52, -0x1.990cdep-76 }, /* 212: log(1237/1024.) */
+ { +0x1.790ed4p-3, +0x1.dc4c68p-28, -0x1.910ac8p-52, +0x1.dd1bd6p-76 }, /* 213: log(1231/1024.) */
+ { +0x1.6f0d28p-3, +0x1.5cad68p-28, +0x1.737c94p-52, -0x1.9184bap-77 }, /* 214: log(1225/1024.) */
+ { +0x1.64fee8p-3, +0x1.04bf88p-28, +0x1.6fca28p-52, +0x1.8884a8p-76 }, /* 215: log(1219/1024.) */
+ { +0x1.5c9400p-3, +0x1.d65cb0p-29, -0x1.b2919cp-53, +0x1.b99bcep-77 }, /* 216: log(1214/1024.) */
+ { +0x1.526e60p-3, -0x1.c5e4bcp-27, -0x1.0ba380p-52, +0x1.d6e3ccp-79 }, /* 217: log(1208/1024.) */
+ { +0x1.483bccp-3, +0x1.9cdc7cp-28, -0x1.5ad8dcp-54, -0x1.392d3cp-83 }, /* 218: log(1202/1024.) */
+ { +0x1.3fb25cp-3, -0x1.a6ad74p-27, +0x1.5be6b4p-52, -0x1.4e0114p-77 }, /* 219: log(1197/1024.) */
+ { +0x1.371fc4p-3, -0x1.fe1708p-27, -0x1.78864cp-52, -0x1.27543ap-76 }, /* 220: log(1192/1024.) */
+ { +0x1.2cca10p-3, -0x1.4141b4p-28, -0x1.ef191cp-52, +0x1.00ee08p-76 }, /* 221: log(1186/1024.) */
+ { +0x1.242310p-3, +0x1.3ba510p-27, -0x1.d003c8p-51, +0x1.162640p-76 }, /* 222: log(1181/1024.) */
+ { +0x1.1b72acp-3, +0x1.52f67cp-27, -0x1.fd6fa0p-51, +0x1.1a3966p-77 }, /* 223: log(1176/1024.) */
+ { +0x1.10f8e4p-3, +0x1.129cd8p-30, +0x1.31ef30p-55, +0x1.a73e38p-79 }, /* 224: log(1170/1024.) */
+ { +0x1.08338cp-3, -0x1.005d7cp-27, -0x1.661a9cp-51, +0x1.1f138ap-79 }, /* 225: log(1165/1024.) */
+ { +0x1.fec914p-4, -0x1.c482a8p-29, -0x1.55746cp-54, +0x1.99f932p-80 }, /* 226: log(1160/1024.) */
+ { +0x1.ed1794p-4, +0x1.d06f00p-29, +0x1.75e45cp-53, -0x1.d0483ep-78 }, /* 227: log(1155/1024.) */
+ { +0x1.db5270p-4, +0x1.87d928p-32, -0x1.0f52a4p-57, +0x1.81f4a6p-84 }, /* 228: log(1150/1024.) */
+ { +0x1.c97978p-4, +0x1.af1d24p-29, -0x1.0977d0p-60, -0x1.8839d0p-84 }, /* 229: log(1145/1024.) */
+ { +0x1.b78c84p-4, -0x1.44f124p-28, -0x1.ef7bc4p-52, +0x1.9e0650p-78 }, /* 230: log(1140/1024.) */
+ { +0x1.a58b60p-4, +0x1.856464p-29, +0x1.c651d0p-55, +0x1.b06b0cp-79 }, /* 231: log(1135/1024.) */
+ { +0x1.9375e4p-4, +0x1.5595ecp-28, +0x1.dc3738p-52, +0x1.86c89ap-81 }, /* 232: log(1130/1024.) */
+ { +0x1.814be4p-4, -0x1.c073fcp-28, -0x1.371f88p-53, -0x1.5f4080p-77 }, /* 233: log(1125/1024.) */
+ { +0x1.6f0d28p-4, +0x1.5cad68p-29, +0x1.737c94p-53, -0x1.9184bap-78 }, /* 234: log(1120/1024.) */
+ { +0x1.60658cp-4, -0x1.6c8af4p-28, +0x1.d8ef74p-55, +0x1.c4f792p-80 }, /* 235: log(1116/1024.) */
+ { +0x1.4e0110p-4, +0x1.146b5cp-29, +0x1.73f7ccp-54, -0x1.d28db8p-79 }, /* 236: log(1111/1024.) */
+ { +0x1.3b8758p-4, +0x1.8b1b70p-28, -0x1.20aca4p-52, -0x1.651894p-76 }, /* 237: log(1106/1024.) */
+ { +0x1.28f834p-4, +0x1.43b6a4p-30, -0x1.452af8p-55, +0x1.976892p-80 }, /* 238: log(1101/1024.) */
+ { +0x1.1a0fbcp-4, -0x1.e4075cp-28, +0x1.1fe618p-52, +0x1.9d6dc2p-77 }, /* 239: log(1097/1024.) */
+ { +0x1.075984p-4, -0x1.4ce370p-29, -0x1.d9fc98p-53, +0x1.4ccf12p-77 }, /* 240: log(1092/1024.) */
+ { +0x1.f0a30cp-5, +0x1.162a68p-37, -0x1.e83368p-61, -0x1.d222a6p-86 }, /* 241: log(1088/1024.) */
+ { +0x1.cae730p-5, -0x1.1a8f7cp-31, -0x1.5f9014p-55, +0x1.2720c0p-79 }, /* 242: log(1083/1024.) */
+ { +0x1.ac9724p-5, -0x1.e8ee08p-29, +0x1.a7de04p-54, -0x1.9bba74p-78 }, /* 243: log(1079/1024.) */
+ { +0x1.868a84p-5, -0x1.ef8128p-30, +0x1.dc5eccp-54, -0x1.58d250p-79 }, /* 244: log(1074/1024.) */
+ { +0x1.67f950p-5, -0x1.ed684cp-30, -0x1.f060c0p-55, -0x1.b1294cp-80 }, /* 245: log(1070/1024.) */
+ { +0x1.494accp-5, +0x1.a6c890p-32, -0x1.c3ad48p-56, -0x1.6dc66cp-84 }, /* 246: log(1066/1024.) */
+ { +0x1.22c71cp-5, -0x1.8abe2cp-32, -0x1.7e7078p-56, -0x1.ddc3dcp-86 }, /* 247: log(1061/1024.) */
+ { +0x1.03d5d8p-5, +0x1.79cfbcp-31, -0x1.da7c4cp-58, +0x1.4e7582p-83 }, /* 248: log(1057/1024.) */
+ { +0x1.c98d18p-6, +0x1.a01904p-31, -0x1.854164p-55, +0x1.883c36p-79 }, /* 249: log(1053/1024.) */
+ { +0x1.8b31fcp-6, -0x1.356500p-30, +0x1.c3ab48p-55, +0x1.b69bdap-80 }, /* 250: log(1049/1024.) */
+ { +0x1.3cea44p-6, +0x1.a352bcp-33, -0x1.8865acp-57, -0x1.48159cp-81 }, /* 251: log(1044/1024.) */
+ { +0x1.fc0a8cp-7, -0x1.e07f84p-32, +0x1.e7cf6cp-58, +0x1.3a69c0p-82 }, /* 252: log(1040/1024.) */
+ { +0x1.7dc474p-7, +0x1.f810a8p-31, -0x1.245b5cp-56, -0x1.a1f4f8p-80 }, /* 253: log(1036/1024.) */
+ { +0x1.fe02a8p-8, -0x1.4ef988p-32, +0x1.1f86ecp-57, +0x1.20723cp-81 }, /* 254: log(1032/1024.) */
+ { +0x1.ff00acp-9, -0x1.d4ef44p-33, +0x1.2821acp-63, +0x1.5a6d32p-87 }, /* 255: log(1028/1024.) */
+ { 0, 0, 0, 0 }
+};
+
+#define PAIR_ADD(d, H, L) do { \
+ gnm_float d_ = (d); \
+ gnm_float dh_ = gnm_floor (d_ * 65536 + 0.5) / 65536; \
+ gnm_float dl_ = d_ - dh_; \
+ if (0) g_printerr ("Adding %.50g (%a)\n", d_, d_); \
+ L += dl_; \
+ H += dh_; \
+} while (0)
+
+#define ADD1(d) PAIR_ADD(d,*yh,*yl)
+
+/*
+ * Compute x * gnm_log (x / M) + (M - x)
+ * aka -x * log1pmx ((M - x) / x)
+ *
+ * Deliver the result back in two parts, *yh and *yl.
+ */
+
+static void
+ebd0(gnm_float x, gnm_float M, gnm_float *yh, gnm_float *yl)
+{
+ gnm_float r, f, fg, M1;
+ int e;
+ int i, j;
+ const int Sb = 10;
+ const double S = 1u << Sb;
+ const int N = G_N_ELEMENTS(bd0_scale) - 1;
+
+ *yl = *yh = 0;
+
+ if (x == M) return;
+ if (x == 0) { PAIR_ADD(M, *yh, *yl); return; }
+ if (M == 0) { *yh = gnm_pinf; return; }
+
+#ifdef DEBUG_EBD0
+ g_printerr ("x=%.20g M=%.20g\n", x, M);
+#endif
+
+ /* FIXME: handle overflow/underflow in division below */
+ r = gnm_frexp (M / x, &e);
+ i = gnm_floor ((r - 0.5) * (2 * N) + 0.5);
+ g_assert (i >= 0 && i <= N);
+ f = gnm_floor (S / (0.5 + i / (2.0 * N)) + 0.5);
+ fg = gnm_ldexp (f, -(e + Sb));
+
+ /* We now have (M * fg / x) close to 1. */
+
+ /*
+ * We need to compute this:
+ * (x/M)^x * exp(M-x) =
+ * (M/x)^-x * exp(M-x) =
+ * (M*fg/x)^-x * (fg)^x * exp(M-x) =
+ * (M*fg/x)^-x * (fg)^x * exp(M*fg-x) * exp(M-M*fg)
+ *
+ * In log terms:
+ * log((x/M)^x * exp(M-x)) =
+ * log((M*fg/x)^-x * (fg)^x * exp(M*fg-x) * exp(M-M*fg)) =
+ * log((M*fg/x)^-x * exp(M*fg-x)) + x*log(fg) + (M-M*fg) =
+ * -x*log1pmx((M*fg-x)/x) + x*log(fg) + M - M*fg =
+ *
+ * Note, that fg has at most 10 bits. If M and x are suitably
+ * "nice" -- such as being integers or half-integers -- then
+ * we can compute M*fg as well as x * bd0_scale[.][.] without
+ * rounding errors.
+ */
+
+ for (j = G_N_ELEMENTS(bd0_scale[i]) - 1; j >= 0; j--) {
+ PAIR_ADD(x * bd0_scale[i][j], *yh, *yl); /* Handles x*log(fg*2^e) */
+ PAIR_ADD(-x * e * bd0_scale[0][j], *yh, *yl); /* Handles x*log(1/2^e) */
+ }
+
+ PAIR_ADD(M, *yh, *yl);
+ M1 = gnm_floor (M + 0.5);
+ PAIR_ADD(-M1 * fg, *yh, *yl);
+ PAIR_ADD(-(M-M1) * fg, *yh, *yl);
+
+ PAIR_ADD(-x * log1pmx ((M * fg - x) / x), *yh, *yl);
+
+#ifdef DEBUG_EBD0
+ g_printerr ("res=%.20g + %.20g\n", *yh, *yl);
+#endif
+}
+
+/* Legacy function. */
static gnm_float
bd0(gnm_float x, gnm_float M)
{
- return (gnm_abs (M - x) < x / 4)
- ? -x * log1pmx ((M - x) / x)
- : x * gnm_log (x / M) + (M - x);
+ gnm_float yh, yl;
+ ebd0 (x, M, &yh, &yl);
+ return yh + yl;
}
@@ -148,7 +374,7 @@ bd0(gnm_float x, gnm_float M)
/* additions for density functions (C.Loader) */
#define R_D_fexp(f,x) (give_log ? -0.5*gnm_log(f)+(x) : gnm_exp(x)/gnm_sqrt(f))
#define R_D_forceint(x) gnm_floor((x) + 0.5)
-#define R_D_nonint(x) (gnm_abs((x) - gnm_floor((x)+0.5)) > 1e-7)
+#define R_D_nonint(x) (gnm_abs((x) - gnm_floor((x)+0.25)) > 1e-7)
/* [neg]ative or [non int]eger : */
#define R_D_negInonint(x) (x < 0. || R_D_nonint(x))
@@ -719,8 +945,7 @@ gnm_float ppois(gnm_float x, gnm_float lambda, gboolean lower_tail, gboolean log
static gnm_float dpois_raw(gnm_float x, gnm_float lambda, gboolean give_log)
{
- GnmQuad mfx;
- int efx;
+ gnm_float yh, yl, f2;
/* x >= 0 ; integer for dpois(), but not e.g. for pgamma()!
lambda >= 0
@@ -731,39 +956,13 @@ static gnm_float dpois_raw(gnm_float x, gnm_float lambda, gboolean give_log)
if (x <= lambda * GNM_MIN) return(R_D_exp(-lambda) );
if (lambda < x * GNM_MIN) return(R_D_exp(-lambda + x*gnm_log(lambda) -lgamma1p (x)));
- if (!qfactf (x, &mfx, &efx)) {
- void *state = gnm_quad_start ();
- gnm_float r, elx, eel;
- GnmQuad qr, qx, ql, mlx, mel;
-
- gnm_quad_init (&ql, lambda);
- gnm_quad_init (&qx, x);
- gnm_quad_pow (&mlx, &elx, &ql, &qx);
- gnm_quad_exp (&mel, &eel, &ql);
- gnm_quad_mul (&qr, &mfx, &mel);
- gnm_quad_div (&qr, &mlx, &qr);
- r = gnm_quad_value (&qr);
- if (gnm_finite (r) && r > 0) {
- gnm_float e = elx - eel - efx;
- if (give_log) {
- GnmQuad qt;
- gnm_quad_init (&qt, e);
- gnm_quad_mul (&qt, &qt, &gnm_quad_ln2);
- gnm_quad_log (&qr, &qr);
- gnm_quad_add (&qr, &qr, &qt);
- r = gnm_quad_value (&qr);
- } else {
- r = gnm_ldexp (r, CLAMP (e, G_MININT, G_MAXINT));
- }
- } else
- r = gnm_ninf;
- gnm_quad_end (state);
-
- if (gnm_finite (r))
- return r;
- }
+ ebd0 (x, lambda, &yh, &yl);
+ PAIR_ADD (stirlerr(x), yh, yl);
+ f2 = M_2PIgnum * x;
- return(R_D_fexp( M_2PIgnum*x, -stirlerr(x)-bd0(x,lambda) ));
+ return give_log
+ ? -yl - yh - 0.5 * gnm_log (f2)
+ : gnm_exp (-yl) * gnm_exp (-yh) / gnm_sqrt (f2);
}
gnm_float dpois(gnm_float x, gnm_float lambda, gboolean give_log)
@@ -2098,26 +2297,49 @@ gnm_float pbinom(gnm_float x, gnm_float n, gnm_float p, gboolean lower_tail, gbo
static gnm_float dbinom_raw(gnm_float x, gnm_float n, gnm_float p, gnm_float q, gboolean give_log)
{
- gnm_float f, lc;
+ gnm_float f2, xh, xl, yh, yl;
+
+ /*
+ * We (ought to) have p+q = 1 with the assumption that the smaller is
+ * more accurate that 1-other, i.e., that the bigger may already have
+ * suffered some cancellation.
+ */
if (p == 0) return((x == 0) ? R_D__1 : R_D__0);
if (q == 0) return((x == n) ? R_D__1 : R_D__0);
if (x == 0) {
- if(n == 0) return R_D__1;
- lc = (p < 0.1) ? -bd0(n,n*q) - n*p : n*gnm_log(q);
- return( R_D_exp(lc) );
+ /* Switch p and q to reduce code duplication. */
+ gnm_float t = p;
+ p = q;
+ q = t;
+ x = n;
}
+
if (x == n) {
- lc = (q < 0.1) ? -bd0(n,n*p) - n*q : n*gnm_log(p);
- return( R_D_exp(lc) );
+ /* Probability is p^n. */
+ if (p <= 0.5)
+ return give_log ? n * gnm_log (p) : gnm_pow (p, n);
+ else
+ return give_log ? n * gnm_log1p (-q) : pow1p (-q, n);
}
if (x < 0 || x > n) return( R_D__0 );
- lc = stirlerr(n) - stirlerr(x) - stirlerr(n-x) - bd0(x,n*p) - bd0(n-x,n*q);
- f = (M_2PIgnum*x*(n-x))/n;
+ ebd0 (x, n*p, &xh, &xl);
+ PAIR_ADD(stirlerr(x), xh, xl);
+
+ ebd0 (n-x, n*q, &yh, &yl);
+ PAIR_ADD(stirlerr(n-x), yh, yl);
+
+ PAIR_ADD(yl, xh, xl);
+ PAIR_ADD(yh, xh, xl);
+ PAIR_ADD(-stirlerr(n), xh, xl);
+
+ f2 = (M_2PIgnum*x*(n-x))/n;
- return R_D_fexp(f,lc);
+ return give_log
+ ? -xl - xh - 0.5 * gnm_log (f2)
+ : gnm_exp (-xl) * gnm_exp (-xh) / gnm_sqrt (f2);
}
gnm_float dbinom(gnm_float x, gnm_float n, gnm_float p, gboolean give_log)
@@ -2627,8 +2849,12 @@ gnm_float dhyper(gnm_float x, gnm_float r, gnm_float b, gnm_float n, gboolean gi
if (n < x || r < x || n - x > b) return(R_D__0);
if (n == 0) return((x == 0) ? R_D__1 : R_D__0);
- p = ((gnm_float)n)/((gnm_float)(r+b));
- q = ((gnm_float)(r+b-n))/((gnm_float)(r+b));
+ /*
+ * Round to float to make both p and q numbers with few (<26) bits set.
+ * Unless p<2^-27 that also ensures that p+q=1 without rounding errors.
+ */
+ p = (float)(n / (gnm_float)(r+b));
+ q = 1 - p;
p1 = dbinom_raw(x, r, p,q,give_log);
p2 = dbinom_raw(n-x,b, p,q,give_log);
@@ -3426,9 +3652,7 @@ static gnm_float ptukey_wprob(gnm_float w, gnm_float rr, gnm_float cc)
v = C + binc * 0.5 * xx;
rinsum = pnorm2 (v - w, v);
- elsum += gnm_pow (rinsum, cc - 1) *
- aa *
- gnm_exp(-0.5 * v * v);
+ elsum += gnm_pow (rinsum, cc - 1) * aa * expmx2h(v);
}
elsum *= (binc * cc * M_1_SQRT_2PI);
pr_w += elsum;
@@ -4736,6 +4960,38 @@ pbinom2 (gnm_float x0, gnm_float x1, gnm_float n, gnm_float p)
}
/* ------------------------------------------------------------------------- */
+/**
+ * expmx2h:
+ * @x: a number
+ *
+ * Returns: The result of exp(-0 5* x*@x) with higher accuracy than the
+ * naive formula.
+ */
+gnm_float
+expmx2h (gnm_float x)
+{
+ x = gnm_abs (x);
+
+ if (x < 5)
+ return gnm_exp (-0.5 * x * x);
+ else if (x >= GNM_MAX_EXP * M_LN2gnum + 10)
+ return 0;
+ else {
+ /*
+ * Split x into two parts, x=x1+x2, such that |x2|<=2^-16.
+ * Assuming that we are using IEEE doubles, that means that
+ * x1*x1 is error free for x<1024 (above which we will underflow
+ * anyway). If we are not using IEEE doubles then this is
+ * still an improvement over the naive formula.
+ */
+ gnm_float x1 = gnm_floor (x * 65536 + 0.5) / 65536;
+ gnm_float x2 = x - x1;
+ return (gnm_exp (-0.5 * x1 * x1) *
+ gnm_exp ((-0.5 * x2 - x1) * x2));
+ }
+}
+
+/* ------------------------------------------------------------------------- */
/**
* pow1p:
@@ -4748,10 +5004,30 @@ pbinom2 (gnm_float x0, gnm_float x1, gnm_float n, gnm_float p)
gnm_float
pow1p (gnm_float x, gnm_float y)
{
- if (gnm_abs (x) > 0.5)
+ /*
+ * We defer to the naive algorithm in two cases: (1) when 1+x
+ * is exact (let us hope the compiler does not mess this up),
+ * and (2) when |x|>1/2 and we have no better algorithm.
+ */
+
+ if ((x + 1) - x == 1 || gnm_abs (x) > 0.5)
return gnm_pow (1 + x, y);
- else
- return gnm_exp (y * gnm_log1p (x));
+ else if (y < 0)
+ return 1 / pow1p (x, -y);
+ else {
+ gnm_float x1 = gnm_floor (x * 65536 + 0.5) / 65536;
+ gnm_float x2 = x - x1;
+ gnm_float h, l;
+ ebd0 (y, y*(x+1), &h, &l);
+ PAIR_ADD (-y*x1, h, l);
+ PAIR_ADD (-y*x2, h, l);
+
+#if 0
+ g_printerr ("pow1p (%.20g, %.20g)\n", x, y);
+#endif
+
+ return gnm_exp (-l) * gnm_exp (-h);
+ }
}
/**
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 241d37d..e24fa05 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -36,6 +36,7 @@ gnm_float logspace_add (gnm_float logx, gnm_float logy);
gnm_float logspace_sub (gnm_float logx, gnm_float logy);
gnm_float gnm_owent (gnm_float h, gnm_float a);
gnm_float gnm_logcf (gnm_float x, gnm_float i, gnm_float d);
+gnm_float expmx2h (gnm_float x);
/* "d": density. */
/* "p": distribution function. */
diff --git a/src/sf-dpq.c b/src/sf-dpq.c
index 9dcdba3..c7374df 100644
--- a/src/sf-dpq.c
+++ b/src/sf-dpq.c
@@ -316,22 +316,8 @@ dnorm (gnm_float x, gnm_float mu, gnm_float sigma, gboolean give_log)
return R_D__0;
if (give_log)
return -(M_LN_SQRT_2PI + 0.5 * x * x + gnm_log (sigma));
- else if (x < 5)
- return M_1_SQRT_2PI * gnm_exp (-0.5 * x * x) / sigma;
- else {
- /*
- * Split x into two parts, x=x1+x2, such that |x2|<=2^-16.
- * Assuming that we are using IEEE doubles, that means that
- * x1*x1 is error free for x<1024 (above which we will underflow
- * anyway). If we are not using IEEE doubles then this is
- * still an improvement over the naive formula.
- */
- gnm_float x1 = gnm_floor (x * 65536 + 0.5) / 65536;
- gnm_float x2 = x - x1;
- return M_1_SQRT_2PI / sigma *
- (gnm_exp (-0.5 * x1 * x1) *
- gnm_exp ((-0.5 * x2 - x1) * x2));
- }
+ else
+ return M_1_SQRT_2PI * expmx2h (x) / sigma;
}
gnm_float
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
