1 /*
2 * mathfunc.c: Mathematical functions.
3 *
4 * Authors:
5 * Ross Ihaka. (See note 1.)
6 * The R Development Core Team. (See note 1.)
7 * Morten Welinder <terra@gnome.org>
8 * Miguel de Icaza (miguel@gnu.org)
9 * Jukka-Pekka Iivonen (iivonen@iki.fi)
10 * Ian Smith (iandjmsmith@aol.com). (See note 2.)
11 */
12
13 /*
14 * NOTE 1: most of this file comes from the "R" package, notably version 2
15 * or newer (we re-sync from time to time).
16 * "R" is distributed under GPL licence, see file COPYING.
17 * The relevant parts are copyright (C) 1998 Ross Ihaka and
18 * 2000-2004 The R Development Core Team.
19 *
20 * Thank you!
21 */
22
23 /*
24 * NOTE 2: the pbeta (and support) code comes from Ian Smith. (Translated
25 * into C, adapted to Gnumeric naming conventions, and R's API conventions
26 * by Morten Welinder. Blame me for problems.)
27 *
28 * Copyright © Ian Smith 2002-2003
29 * Version 1.0.24
30 * Thanks to Jerry W. Lewis for help with testing of and improvements to the code.
31 *
32 * Thank you!
33 */
34
35
36 #include <gnumeric-config.h>
37 #include <gnumeric.h>
38 #include <mathfunc.h>
39 #include <sf-dpq.h>
40 #include <sf-gamma.h>
41 #include <glib/gi18n-lib.h>
42
43 #include <math.h>
44 #include <stdlib.h>
45 #include <float.h>
46 #include <unistd.h>
47 #include <string.h>
48 #include <goffice/goffice.h>
49 #include <glib/gstdio.h>
50 #include <value.h>
51 #include <gutils.h>
52
53 /* R code wants this, so provide it. */
54 #ifndef IEEE_754
55 #define IEEE_754
56 #endif
57
58 #define M_SQRT_32 GNM_const(5.656854249492380195206754896838) /* sqrt(32) */
59 #define M_1_SQRT_2PI GNM_const(0.398942280401432677939946059934) /* 1/sqrt(2pi) */
60 #define M_2PIgnum (2 * M_PIgnum)
61
62 #define ML_ERROR(cause) do { } while(0)
63 #define MATHLIB_WARNING g_warning
64 #define REprintf g_warning
65
fmin2(gnm_float x,gnm_float y)66 static inline gnm_float fmin2 (gnm_float x, gnm_float y) { return MIN (x, y); }
fmax2(gnm_float x,gnm_float y)67 static inline gnm_float fmax2 (gnm_float x, gnm_float y) { return MAX (x, y); }
68
69 #define MATHLIB_STANDALONE
70 #define ML_ERR_return_NAN { return gnm_nan; }
71 static void pnorm_both (gnm_float x, gnm_float *cum, gnm_float *ccum, int i_tail, gboolean log_p);
72
73 /* MW ---------------------------------------------------------------------- */
74
75 void
mathfunc_init(void)76 mathfunc_init (void)
77 {
78 /* Nothing, for the time being. */
79 }
80
81 /*
82 * A table of logs for scaling purposes. Each value has four parts with
83 * 23 bits in each. That means each part can be multiplied by a double
84 * with at most 30 bits set and not have any rounding error. Note, that
85 * the first entry is log(2).
86 *
87 * Entry i is associated with the value r = 0.5 + i / 256.0. The
88 * argument to log is p/q where q=1024 and p=floor(q / r + 0.5).
89 * Thus r*p/q is close to 1.
90 */
91 static const float bd0_scale[128 + 1][4] = {
92 { +0x1.62e430p-1, -0x1.05c610p-29, -0x1.950d88p-54, +0x1.d9cc02p-79 }, /* 128: log(2048/1024.) */
93 { +0x1.5ee02cp-1, -0x1.6dbe98p-25, -0x1.51e540p-50, +0x1.2bfa48p-74 }, /* 129: log(2032/1024.) */
94 { +0x1.5ad404p-1, +0x1.86b3e4p-26, +0x1.9f6534p-50, +0x1.54be04p-74 }, /* 130: log(2016/1024.) */
95 { +0x1.570124p-1, -0x1.9ed750p-25, -0x1.f37dd0p-51, +0x1.10b770p-77 }, /* 131: log(2001/1024.) */
96 { +0x1.5326e4p-1, -0x1.9b9874p-25, -0x1.378194p-49, +0x1.56feb2p-74 }, /* 132: log(1986/1024.) */
97 { +0x1.4f4528p-1, +0x1.aca70cp-28, +0x1.103e74p-53, +0x1.9c410ap-81 }, /* 133: log(1971/1024.) */
98 { +0x1.4b5bd8p-1, -0x1.6a91d8p-25, -0x1.8e43d0p-50, -0x1.afba9ep-77 }, /* 134: log(1956/1024.) */
99 { +0x1.47ae54p-1, -0x1.abb51cp-25, +0x1.19b798p-51, +0x1.45e09cp-76 }, /* 135: log(1942/1024.) */
100 { +0x1.43fa00p-1, -0x1.d06318p-25, -0x1.8858d8p-49, -0x1.1927c4p-75 }, /* 136: log(1928/1024.) */
101 { +0x1.3ffa40p-1, +0x1.1a427cp-25, +0x1.151640p-53, -0x1.4f5606p-77 }, /* 137: log(1913/1024.) */
102 { +0x1.3c7c80p-1, -0x1.19bf48p-34, +0x1.05fc94p-58, -0x1.c096fcp-82 }, /* 138: log(1900/1024.) */
103 { +0x1.38b320p-1, +0x1.6b5778p-25, +0x1.be38d0p-50, -0x1.075e96p-74 }, /* 139: log(1886/1024.) */
104 { +0x1.34e288p-1, +0x1.d9ce1cp-25, +0x1.316eb8p-49, +0x1.2d885cp-73 }, /* 140: log(1872/1024.) */
105 { +0x1.315124p-1, +0x1.c2fc60p-29, -0x1.4396fcp-53, +0x1.acf376p-78 }, /* 141: log(1859/1024.) */
106 { +0x1.2db954p-1, +0x1.720de4p-25, -0x1.d39b04p-49, -0x1.f11176p-76 }, /* 142: log(1846/1024.) */
107 { +0x1.2a1b08p-1, -0x1.562494p-25, +0x1.a7863cp-49, +0x1.85dd64p-73 }, /* 143: log(1833/1024.) */
108 { +0x1.267620p-1, +0x1.3430e0p-29, -0x1.96a958p-56, +0x1.f8e636p-82 }, /* 144: log(1820/1024.) */
109 { +0x1.23130cp-1, +0x1.7bebf4p-25, +0x1.416f1cp-52, -0x1.78dd36p-77 }, /* 145: log(1808/1024.) */
110 { +0x1.1faa34p-1, +0x1.70e128p-26, +0x1.81817cp-50, -0x1.c2179cp-76 }, /* 146: log(1796/1024.) */
111 { +0x1.1bf204p-1, +0x1.3a9620p-28, +0x1.2f94c0p-52, +0x1.9096c0p-76 }, /* 147: log(1783/1024.) */
112 { +0x1.187ce4p-1, -0x1.077870p-27, +0x1.655a80p-51, +0x1.eaafd6p-78 }, /* 148: log(1771/1024.) */
113 { +0x1.1501c0p-1, -0x1.406cacp-25, -0x1.e72290p-49, +0x1.5dd800p-73 }, /* 149: log(1759/1024.) */
114 { +0x1.11cb80p-1, +0x1.787cd0p-25, -0x1.efdc78p-51, -0x1.5380cep-77 }, /* 150: log(1748/1024.) */
115 { +0x1.0e4498p-1, +0x1.747324p-27, -0x1.024548p-51, +0x1.77a5a6p-75 }, /* 151: log(1736/1024.) */
116 { +0x1.0b036cp-1, +0x1.690c74p-25, +0x1.5d0cc4p-50, -0x1.c0e23cp-76 }, /* 152: log(1725/1024.) */
117 { +0x1.077070p-1, -0x1.a769bcp-27, +0x1.452234p-52, +0x1.6ba668p-76 }, /* 153: log(1713/1024.) */
118 { +0x1.04240cp-1, -0x1.a686acp-27, -0x1.ef46b0p-52, -0x1.5ce10cp-76 }, /* 154: log(1702/1024.) */
119 { +0x1.00d22cp-1, +0x1.fc0e10p-25, +0x1.6ee034p-50, -0x1.19a2ccp-74 }, /* 155: log(1691/1024.) */
120 { +0x1.faf588p-2, +0x1.ef1e64p-27, -0x1.26504cp-54, -0x1.b15792p-82 }, /* 156: log(1680/1024.) */
121 { +0x1.f4d87cp-2, +0x1.d7b980p-26, -0x1.a114d8p-50, +0x1.9758c6p-75 }, /* 157: log(1670/1024.) */
122 { +0x1.ee1414p-2, +0x1.2ec060p-26, +0x1.dc00fcp-52, +0x1.f8833cp-76 }, /* 158: log(1659/1024.) */
123 { +0x1.e7e32cp-2, -0x1.ac796cp-27, -0x1.a68818p-54, +0x1.235d02p-78 }, /* 159: log(1649/1024.) */
124 { +0x1.e108a0p-2, -0x1.768ba4p-28, -0x1.f050a8p-52, +0x1.00d632p-82 }, /* 160: log(1638/1024.) */
125 { +0x1.dac354p-2, -0x1.d3a6acp-30, +0x1.18734cp-57, -0x1.f97902p-83 }, /* 161: log(1628/1024.) */
126 { +0x1.d47424p-2, +0x1.7dbbacp-31, -0x1.d5ada4p-56, +0x1.56fcaap-81 }, /* 162: log(1618/1024.) */
127 { +0x1.ce1af0p-2, +0x1.70be7cp-27, +0x1.6f6fa4p-51, +0x1.7955a2p-75 }, /* 163: log(1608/1024.) */
128 { +0x1.c7b798p-2, +0x1.ec36ecp-26, -0x1.07e294p-50, -0x1.ca183cp-75 }, /* 164: log(1598/1024.) */
129 { +0x1.c1ef04p-2, +0x1.c1dfd4p-26, +0x1.888eecp-50, -0x1.fd6b86p-75 }, /* 165: log(1589/1024.) */
130 { +0x1.bb7810p-2, +0x1.478bfcp-26, +0x1.245b8cp-50, +0x1.ea9d52p-74 }, /* 166: log(1579/1024.) */
131 { +0x1.b59da0p-2, -0x1.882b08p-27, +0x1.31573cp-53, -0x1.8c249ap-77 }, /* 167: log(1570/1024.) */
132 { +0x1.af1294p-2, -0x1.b710f4p-27, +0x1.622670p-51, +0x1.128578p-76 }, /* 168: log(1560/1024.) */
133 { +0x1.a925d4p-2, -0x1.0ae750p-27, +0x1.574ed4p-51, +0x1.084996p-75 }, /* 169: log(1551/1024.) */
134 { +0x1.a33040p-2, +0x1.027d30p-29, +0x1.b9a550p-53, -0x1.b2e38ap-78 }, /* 170: log(1542/1024.) */
135 { +0x1.9d31c0p-2, -0x1.5ec12cp-26, -0x1.5245e0p-52, +0x1.2522d0p-79 }, /* 171: log(1533/1024.) */
136 { +0x1.972a34p-2, +0x1.135158p-30, +0x1.a5c09cp-56, +0x1.24b70ep-80 }, /* 172: log(1524/1024.) */
137 { +0x1.911984p-2, +0x1.0995d4p-26, +0x1.3bfb5cp-50, +0x1.2c9dd6p-75 }, /* 173: log(1515/1024.) */
138 { +0x1.8bad98p-2, -0x1.1d6144p-29, +0x1.5b9208p-53, +0x1.1ec158p-77 }, /* 174: log(1507/1024.) */
139 { +0x1.858b58p-2, -0x1.1b4678p-27, +0x1.56cab4p-53, -0x1.2fdc0cp-78 }, /* 175: log(1498/1024.) */
140 { +0x1.7f5fa0p-2, +0x1.3aaf48p-27, +0x1.461964p-51, +0x1.4ae476p-75 }, /* 176: log(1489/1024.) */
141 { +0x1.79db68p-2, -0x1.7e5054p-26, +0x1.673750p-51, -0x1.a11f7ap-76 }, /* 177: log(1481/1024.) */
142 { +0x1.744f88p-2, -0x1.cc0e18p-26, -0x1.1e9d18p-50, -0x1.6c06bcp-78 }, /* 178: log(1473/1024.) */
143 { +0x1.6e08ecp-2, -0x1.5d45e0p-26, -0x1.c73ec8p-50, +0x1.318d72p-74 }, /* 179: log(1464/1024.) */
144 { +0x1.686c80p-2, +0x1.e9b14cp-26, -0x1.13bbd4p-50, -0x1.efeb1cp-78 }, /* 180: log(1456/1024.) */
145 { +0x1.62c830p-2, -0x1.a8c70cp-27, -0x1.5a1214p-51, -0x1.bab3fcp-79 }, /* 181: log(1448/1024.) */
146 { +0x1.5d1bdcp-2, -0x1.4fec6cp-31, +0x1.423638p-56, +0x1.ee3feep-83 }, /* 182: log(1440/1024.) */
147 { +0x1.576770p-2, +0x1.7455a8p-26, -0x1.3ab654p-50, -0x1.26be4cp-75 }, /* 183: log(1432/1024.) */
148 { +0x1.5262e0p-2, -0x1.146778p-26, -0x1.b9f708p-52, -0x1.294018p-77 }, /* 184: log(1425/1024.) */
149 { +0x1.4c9f08p-2, +0x1.e152c4p-26, -0x1.dde710p-53, +0x1.fd2208p-77 }, /* 185: log(1417/1024.) */
150 { +0x1.46d2d8p-2, +0x1.c28058p-26, -0x1.936284p-50, +0x1.9fdd68p-74 }, /* 186: log(1409/1024.) */
151 { +0x1.41b940p-2, +0x1.cce0c0p-26, -0x1.1a4050p-50, +0x1.bc0376p-76 }, /* 187: log(1402/1024.) */
152 { +0x1.3bdd24p-2, +0x1.d6296cp-27, +0x1.425b48p-51, -0x1.cddb2cp-77 }, /* 188: log(1394/1024.) */
153 { +0x1.36b578p-2, -0x1.287ddcp-27, -0x1.2d0f4cp-51, +0x1.38447ep-75 }, /* 189: log(1387/1024.) */
154 { +0x1.31871cp-2, +0x1.2a8830p-27, +0x1.3eae54p-52, -0x1.898136p-77 }, /* 190: log(1380/1024.) */
155 { +0x1.2b9304p-2, -0x1.51d8b8p-28, +0x1.27694cp-52, -0x1.fd852ap-76 }, /* 191: log(1372/1024.) */
156 { +0x1.265620p-2, -0x1.d98f3cp-27, +0x1.a44338p-51, -0x1.56e85ep-78 }, /* 192: log(1365/1024.) */
157 { +0x1.211254p-2, +0x1.986160p-26, +0x1.73c5d0p-51, +0x1.4a861ep-75 }, /* 193: log(1358/1024.) */
158 { +0x1.1bc794p-2, +0x1.fa3918p-27, +0x1.879c5cp-51, +0x1.16107cp-78 }, /* 194: log(1351/1024.) */
159 { +0x1.1675ccp-2, -0x1.4545a0p-26, +0x1.c07398p-51, +0x1.f55c42p-76 }, /* 195: log(1344/1024.) */
160 { +0x1.111ce4p-2, +0x1.f72670p-37, -0x1.b84b5cp-61, +0x1.a4a4dcp-85 }, /* 196: log(1337/1024.) */
161 { +0x1.0c81d4p-2, +0x1.0c150cp-27, +0x1.218600p-51, -0x1.d17312p-76 }, /* 197: log(1331/1024.) */
162 { +0x1.071b84p-2, +0x1.fcd590p-26, +0x1.a3a2e0p-51, +0x1.fe5ef8p-76 }, /* 198: log(1324/1024.) */
163 { +0x1.01ade4p-2, -0x1.bb1844p-28, +0x1.db3cccp-52, +0x1.1f56fcp-77 }, /* 199: log(1317/1024.) */
164 { +0x1.fa01c4p-3, -0x1.12a0d0p-29, -0x1.f71fb0p-54, +0x1.e287a4p-78 }, /* 200: log(1311/1024.) */
165 { +0x1.ef0adcp-3, +0x1.7b8b28p-28, -0x1.35bce4p-52, -0x1.abc8f8p-79 }, /* 201: log(1304/1024.) */
166 { +0x1.e598ecp-3, +0x1.5a87e4p-27, -0x1.134bd0p-51, +0x1.c2cebep-76 }, /* 202: log(1298/1024.) */
167 { +0x1.da85d8p-3, -0x1.df31b0p-27, +0x1.94c16cp-57, +0x1.8fd7eap-82 }, /* 203: log(1291/1024.) */
168 { +0x1.d0fb80p-3, -0x1.bb5434p-28, -0x1.ea5640p-52, -0x1.8ceca4p-77 }, /* 204: log(1285/1024.) */
169 { +0x1.c765b8p-3, +0x1.e4d68cp-27, +0x1.5b59b4p-51, +0x1.76f6c4p-76 }, /* 205: log(1279/1024.) */
170 { +0x1.bdc46cp-3, -0x1.1cbb50p-27, +0x1.2da010p-51, +0x1.eb282cp-75 }, /* 206: log(1273/1024.) */
171 { +0x1.b27980p-3, -0x1.1b9ce0p-27, +0x1.7756f8p-52, +0x1.2ff572p-76 }, /* 207: log(1266/1024.) */
172 { +0x1.a8bed0p-3, -0x1.bbe874p-30, +0x1.85cf20p-56, +0x1.b9cf18p-80 }, /* 208: log(1260/1024.) */
173 { +0x1.9ef83cp-3, +0x1.2769a4p-27, -0x1.85bda0p-52, +0x1.8c8018p-79 }, /* 209: log(1254/1024.) */
174 { +0x1.9525a8p-3, +0x1.cf456cp-27, -0x1.7137d8p-52, -0x1.f158e8p-76 }, /* 210: log(1248/1024.) */
175 { +0x1.8b46f8p-3, +0x1.11b12cp-30, +0x1.9f2104p-54, -0x1.22836ep-78 }, /* 211: log(1242/1024.) */
176 { +0x1.83040cp-3, +0x1.2379e4p-28, +0x1.b71c70p-52, -0x1.990cdep-76 }, /* 212: log(1237/1024.) */
177 { +0x1.790ed4p-3, +0x1.dc4c68p-28, -0x1.910ac8p-52, +0x1.dd1bd6p-76 }, /* 213: log(1231/1024.) */
178 { +0x1.6f0d28p-3, +0x1.5cad68p-28, +0x1.737c94p-52, -0x1.9184bap-77 }, /* 214: log(1225/1024.) */
179 { +0x1.64fee8p-3, +0x1.04bf88p-28, +0x1.6fca28p-52, +0x1.8884a8p-76 }, /* 215: log(1219/1024.) */
180 { +0x1.5c9400p-3, +0x1.d65cb0p-29, -0x1.b2919cp-53, +0x1.b99bcep-77 }, /* 216: log(1214/1024.) */
181 { +0x1.526e60p-3, -0x1.c5e4bcp-27, -0x1.0ba380p-52, +0x1.d6e3ccp-79 }, /* 217: log(1208/1024.) */
182 { +0x1.483bccp-3, +0x1.9cdc7cp-28, -0x1.5ad8dcp-54, -0x1.392d3cp-83 }, /* 218: log(1202/1024.) */
183 { +0x1.3fb25cp-3, -0x1.a6ad74p-27, +0x1.5be6b4p-52, -0x1.4e0114p-77 }, /* 219: log(1197/1024.) */
184 { +0x1.371fc4p-3, -0x1.fe1708p-27, -0x1.78864cp-52, -0x1.27543ap-76 }, /* 220: log(1192/1024.) */
185 { +0x1.2cca10p-3, -0x1.4141b4p-28, -0x1.ef191cp-52, +0x1.00ee08p-76 }, /* 221: log(1186/1024.) */
186 { +0x1.242310p-3, +0x1.3ba510p-27, -0x1.d003c8p-51, +0x1.162640p-76 }, /* 222: log(1181/1024.) */
187 { +0x1.1b72acp-3, +0x1.52f67cp-27, -0x1.fd6fa0p-51, +0x1.1a3966p-77 }, /* 223: log(1176/1024.) */
188 { +0x1.10f8e4p-3, +0x1.129cd8p-30, +0x1.31ef30p-55, +0x1.a73e38p-79 }, /* 224: log(1170/1024.) */
189 { +0x1.08338cp-3, -0x1.005d7cp-27, -0x1.661a9cp-51, +0x1.1f138ap-79 }, /* 225: log(1165/1024.) */
190 { +0x1.fec914p-4, -0x1.c482a8p-29, -0x1.55746cp-54, +0x1.99f932p-80 }, /* 226: log(1160/1024.) */
191 { +0x1.ed1794p-4, +0x1.d06f00p-29, +0x1.75e45cp-53, -0x1.d0483ep-78 }, /* 227: log(1155/1024.) */
192 { +0x1.db5270p-4, +0x1.87d928p-32, -0x1.0f52a4p-57, +0x1.81f4a6p-84 }, /* 228: log(1150/1024.) */
193 { +0x1.c97978p-4, +0x1.af1d24p-29, -0x1.0977d0p-60, -0x1.8839d0p-84 }, /* 229: log(1145/1024.) */
194 { +0x1.b78c84p-4, -0x1.44f124p-28, -0x1.ef7bc4p-52, +0x1.9e0650p-78 }, /* 230: log(1140/1024.) */
195 { +0x1.a58b60p-4, +0x1.856464p-29, +0x1.c651d0p-55, +0x1.b06b0cp-79 }, /* 231: log(1135/1024.) */
196 { +0x1.9375e4p-4, +0x1.5595ecp-28, +0x1.dc3738p-52, +0x1.86c89ap-81 }, /* 232: log(1130/1024.) */
197 { +0x1.814be4p-4, -0x1.c073fcp-28, -0x1.371f88p-53, -0x1.5f4080p-77 }, /* 233: log(1125/1024.) */
198 { +0x1.6f0d28p-4, +0x1.5cad68p-29, +0x1.737c94p-53, -0x1.9184bap-78 }, /* 234: log(1120/1024.) */
199 { +0x1.60658cp-4, -0x1.6c8af4p-28, +0x1.d8ef74p-55, +0x1.c4f792p-80 }, /* 235: log(1116/1024.) */
200 { +0x1.4e0110p-4, +0x1.146b5cp-29, +0x1.73f7ccp-54, -0x1.d28db8p-79 }, /* 236: log(1111/1024.) */
201 { +0x1.3b8758p-4, +0x1.8b1b70p-28, -0x1.20aca4p-52, -0x1.651894p-76 }, /* 237: log(1106/1024.) */
202 { +0x1.28f834p-4, +0x1.43b6a4p-30, -0x1.452af8p-55, +0x1.976892p-80 }, /* 238: log(1101/1024.) */
203 { +0x1.1a0fbcp-4, -0x1.e4075cp-28, +0x1.1fe618p-52, +0x1.9d6dc2p-77 }, /* 239: log(1097/1024.) */
204 { +0x1.075984p-4, -0x1.4ce370p-29, -0x1.d9fc98p-53, +0x1.4ccf12p-77 }, /* 240: log(1092/1024.) */
205 { +0x1.f0a30cp-5, +0x1.162a68p-37, -0x1.e83368p-61, -0x1.d222a6p-86 }, /* 241: log(1088/1024.) */
206 { +0x1.cae730p-5, -0x1.1a8f7cp-31, -0x1.5f9014p-55, +0x1.2720c0p-79 }, /* 242: log(1083/1024.) */
207 { +0x1.ac9724p-5, -0x1.e8ee08p-29, +0x1.a7de04p-54, -0x1.9bba74p-78 }, /* 243: log(1079/1024.) */
208 { +0x1.868a84p-5, -0x1.ef8128p-30, +0x1.dc5eccp-54, -0x1.58d250p-79 }, /* 244: log(1074/1024.) */
209 { +0x1.67f950p-5, -0x1.ed684cp-30, -0x1.f060c0p-55, -0x1.b1294cp-80 }, /* 245: log(1070/1024.) */
210 { +0x1.494accp-5, +0x1.a6c890p-32, -0x1.c3ad48p-56, -0x1.6dc66cp-84 }, /* 246: log(1066/1024.) */
211 { +0x1.22c71cp-5, -0x1.8abe2cp-32, -0x1.7e7078p-56, -0x1.ddc3dcp-86 }, /* 247: log(1061/1024.) */
212 { +0x1.03d5d8p-5, +0x1.79cfbcp-31, -0x1.da7c4cp-58, +0x1.4e7582p-83 }, /* 248: log(1057/1024.) */
213 { +0x1.c98d18p-6, +0x1.a01904p-31, -0x1.854164p-55, +0x1.883c36p-79 }, /* 249: log(1053/1024.) */
214 { +0x1.8b31fcp-6, -0x1.356500p-30, +0x1.c3ab48p-55, +0x1.b69bdap-80 }, /* 250: log(1049/1024.) */
215 { +0x1.3cea44p-6, +0x1.a352bcp-33, -0x1.8865acp-57, -0x1.48159cp-81 }, /* 251: log(1044/1024.) */
216 { +0x1.fc0a8cp-7, -0x1.e07f84p-32, +0x1.e7cf6cp-58, +0x1.3a69c0p-82 }, /* 252: log(1040/1024.) */
217 { +0x1.7dc474p-7, +0x1.f810a8p-31, -0x1.245b5cp-56, -0x1.a1f4f8p-80 }, /* 253: log(1036/1024.) */
218 { +0x1.fe02a8p-8, -0x1.4ef988p-32, +0x1.1f86ecp-57, +0x1.20723cp-81 }, /* 254: log(1032/1024.) */
219 { +0x1.ff00acp-9, -0x1.d4ef44p-33, +0x1.2821acp-63, +0x1.5a6d32p-87 }, /* 255: log(1028/1024.) */
220 { 0, 0, 0, 0 }
221 };
222
223 #define PAIR_ADD(d, H, L) do { \
224 gnm_float d_ = (d); \
225 gnm_float dh_ = gnm_floor (d_ * 65536 + 0.5) / 65536; \
226 gnm_float dl_ = d_ - dh_; \
227 if (0) g_printerr ("Adding %.50g (%a)\n", d_, d_); \
228 L += dl_; \
229 H += dh_; \
230 } while (0)
231
232 #define ADD1(d) PAIR_ADD(d,*yh,*yl)
233
234 /*
235 * Compute x * gnm_log (x / M) + (M - x)
236 * aka -x * log1pmx ((M - x) / x)
237 *
238 * Deliver the result back in two parts, *yh and *yl.
239 */
240
241 static void
ebd0(gnm_float x,gnm_float M,gnm_float * yh,gnm_float * yl)242 ebd0(gnm_float x, gnm_float M, gnm_float *yh, gnm_float *yl)
243 {
244 gnm_float r, f, fg, M1;
245 int e;
246 int i, j;
247 const int Sb = 10;
248 const double S = 1u << Sb;
249 const int N = G_N_ELEMENTS(bd0_scale) - 1;
250 const double e4 = GNM_EPSILON * GNM_EPSILON * GNM_EPSILON * GNM_EPSILON;
251
252 if (gnm_isnan (x) || gnm_isnan (M)) {
253 *yl = *yh = x;
254 return;
255 }
256
257 *yl = *yh = 0;
258
259 if (x == M) return;
260 if (x < M * e4) { ADD1(M); return; }
261 if (M == 0) { *yh = gnm_pinf; return; }
262
263 if (M < x * e4) {
264 ADD1(x * (gnm_log(x) - 1));
265 ADD1(-x * gnm_log(M));
266 return;
267 }
268
269 #ifdef DEBUG_EBD0
270 g_printerr ("x=%.20g M=%.20g\n", x, M);
271 #endif
272
273 r = gnm_frexp (M / x, &e);
274 i = gnm_floor ((r - 0.5) * (2 * N) + 0.5);
275 g_assert (i >= 0 && i <= N);
276 f = gnm_floor (S / (0.5 + i / (2.0 * N)) + 0.5);
277 fg = gnm_ldexp (f, -(e + Sb));
278
279 /* We now have (M * fg / x) close to 1. */
280
281 /*
282 * We need to compute this:
283 * (x/M)^x * exp(M-x) =
284 * (M/x)^-x * exp(M-x) =
285 * (M*fg/x)^-x * (fg)^x * exp(M-x) =
286 * (M*fg/x)^-x * (fg)^x * exp(M*fg-x) * exp(M-M*fg)
287 *
288 * In log terms:
289 * log((x/M)^x * exp(M-x)) =
290 * log((M*fg/x)^-x * (fg)^x * exp(M*fg-x) * exp(M-M*fg)) =
291 * log((M*fg/x)^-x * exp(M*fg-x)) + x*log(fg) + (M-M*fg) =
292 * -x*log1pmx((M*fg-x)/x) + x*log(fg) + M - M*fg =
293 *
294 * Note, that fg has at most 10 bits. If M and x are suitably
295 * "nice" -- such as being integers or half-integers -- then
296 * we can compute M*fg as well as x * bd0_scale[.][.] without
297 * rounding errors.
298 */
299
300 for (j = G_N_ELEMENTS(bd0_scale[i]) - 1; j >= 0; j--) {
301 ADD1(x * bd0_scale[i][j]); /* Handles x*log(fg*2^e) */
302 ADD1(-x * e * bd0_scale[0][j]); /* Handles x*log(1/2^e) */
303 }
304
305 ADD1(M);
306 M1 = gnm_floor (M + 0.5);
307 ADD1(-M1 * fg);
308 ADD1(-(M-M1) * fg);
309
310 ADD1(-x * log1pmx ((M * fg - x) / x));
311
312 #ifdef DEBUG_EBD0
313 g_printerr ("res=%.20g + %.20g\n", *yh, *yl);
314 #endif
315 }
316
317 #undef ADD1
318
319 /* Legacy function. */
320 static gnm_float
bd0(gnm_float x,gnm_float M)321 bd0(gnm_float x, gnm_float M)
322 {
323 gnm_float yh, yl;
324 ebd0 (x, M, &yh, &yl);
325 return yh + yl;
326 }
327
328 /* ------------------------------------------------------------------------- */
329 /* --- BEGIN MAGIC R SOURCE MARKER --- */
330
331 /* The following source code was imported from the R project. */
332 /* It was automatically transformed by tools/import-R. */
333
334 /* Imported src/nmath/dpq.h from R. */
335 /* Utilities for `dpq' handling (density/probability/quantile) */
336
337 /* give_log in "d"; log_p in "p" & "q" : */
338 #define give_log log_p
339 /* "DEFAULT" */
340 /* --------- */
341 #define R_D__0 (log_p ? gnm_ninf : 0.) /* 0 */
342 #define R_D__1 (log_p ? 0. : 1.) /* 1 */
343 #define R_DT_0 (lower_tail ? R_D__0 : R_D__1) /* 0 */
344 #define R_DT_1 (lower_tail ? R_D__1 : R_D__0) /* 1 */
345
346 #define R_D_Lval(p) (lower_tail ? (p) : (1 - (p))) /* p */
347 #define R_D_Cval(p) (lower_tail ? (1 - (p)) : (p)) /* 1 - p */
348
349 #define R_D_val(x) (log_p ? gnm_log(x) : (x)) /* x in pF(x,..) */
350 #define R_D_qIv(p) (log_p ? gnm_exp(p) : (p)) /* p in qF(p,..) */
351 #define R_D_exp(x) (log_p ? (x) : gnm_exp(x)) /* gnm_exp(x) */
352 #define R_D_log(p) (log_p ? (p) : gnm_log(p)) /* gnm_log(p) */
353 #define R_D_Clog(p) (log_p ? gnm_log1p(-(p)) : (1 - (p)))/* [log](1-p) */
354
355 /* gnm_log(1-gnm_exp(x)): R_D_LExp(x) == (gnm_log1p(- R_D_qIv(x))) but even more stable:*/
356 #define R_D_LExp(x) (log_p ? swap_log_tail(x) : gnm_log1p(-x))
357
358 /*till 1.8.x:
359 * #define R_DT_val(x) R_D_val(R_D_Lval(x))
360 * #define R_DT_Cval(x) R_D_val(R_D_Cval(x)) */
361 #define R_DT_val(x) (lower_tail ? R_D_val(x) : R_D_Clog(x))
362 #define R_DT_Cval(x) (lower_tail ? R_D_Clog(x) : R_D_val(x))
363
364 /*#define R_DT_qIv(p) R_D_Lval(R_D_qIv(p)) * p in qF ! */
365 #define R_DT_qIv(p) (log_p ? (lower_tail ? gnm_exp(p) : - gnm_expm1(p)) \
366 : R_D_Lval(p))
367
368 /*#define R_DT_CIv(p) R_D_Cval(R_D_qIv(p)) * 1 - p in qF */
369 #define R_DT_CIv(p) (log_p ? (lower_tail ? -gnm_expm1(p) : gnm_exp(p)) \
370 : R_D_Cval(p))
371
372 #define R_DT_exp(x) R_D_exp(R_D_Lval(x)) /* gnm_exp(x) */
373 #define R_DT_Cexp(x) R_D_exp(R_D_Cval(x)) /* gnm_exp(1 - x) */
374
375 #define R_DT_log(p) (lower_tail? R_D_log(p) : R_D_LExp(p))/* gnm_log(p) in qF */
376 #define R_DT_Clog(p) (lower_tail? R_D_LExp(p): R_D_log(p))/* gnm_log1p (-p) in qF*/
377 #define R_DT_Log(p) (lower_tail? (p) : swap_log_tail(p))
378 /* == R_DT_log when we already "know" log_p == TRUE :*/
379
380
381 #define R_Q_P01_check(p) \
382 if ((log_p && p > 0) || \
383 (!log_p && (p < 0 || p > 1)) ) \
384 ML_ERR_return_NAN
385
386
387 /* additions for density functions (C.Loader) */
388 #define R_D_fexp(f,x) (give_log ? -0.5*gnm_log(f)+(x) : gnm_exp(x)/gnm_sqrt(f))
389 #define R_D_forceint(x) gnm_floor((x) + 0.5)
390 #define R_D_nonint(x) (gnm_abs((x) - gnm_floor((x)+0.25)) > 1e-7)
391 /* [neg]ative or [non int]eger : */
392 #define R_D_negInonint(x) (x < 0. || R_D_nonint(x))
393
394 #define R_D_nonint_check(x) \
395 if(R_D_nonint(x)) { \
396 MATHLIB_WARNING("non-integer x = %" GNM_FORMAT_f "", x); \
397 return R_D__0; \
398 }
399
400 /* ------------------------------------------------------------------------ */
401 /* Imported src/nmath/ftrunc.c from R. */
402 /*
403 * Mathlib : A C Library of Special Functions
404 * Copyright (C) 1998 Ross Ihaka
405 *
406 * This program is free software; you can redistribute it and/or modify
407 * it under the terms of the GNU General Public License as published by
408 * the Free Software Foundation; either version 2 of the License, or
409 * (at your option) any later version.
410 *
411 * This program is distributed in the hope that it will be useful,
412 * but WITHOUT ANY WARRANTY; without even the implied warranty of
413 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
414 * GNU General Public License for more details.
415 *
416 * You should have received a copy of the GNU General Public License
417 * along with this program; if not, write to the Free Software
418 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
419 * 02110-1301 USA.
420 *
421 * SYNOPSIS
422 *
423 * #include <Rmath.h>
424 * double ftrunc(double x);
425 *
426 * DESCRIPTION
427 *
428 * Truncation toward zero.
429 */
430
431
gnm_trunc(gnm_float x)432 gnm_float gnm_trunc(gnm_float x)
433 {
434 if(x >= 0) return gnm_floor(x);
435 else return gnm_ceil(x);
436 }
437
438 /* ------------------------------------------------------------------------ */
439 /* Imported src/nmath/pnorm.c from R. */
440 /*
441 * Mathlib : A C Library of Special Functions
442 * Copyright (C) 1998 Ross Ihaka
443 * Copyright (C) 2000-2002 The R Development Core Team
444 * Copyright (C) 2003 The R Foundation
445 *
446 * This program is free software; you can redistribute it and/or modify
447 * it under the terms of the GNU General Public License as published by
448 * the Free Software Foundation; either version 2 of the License, or
449 * (at your option) any later version.
450 *
451 * This program is distributed in the hope that it will be useful,
452 * but WITHOUT ANY WARRANTY; without even the implied warranty of
453 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
454 * GNU General Public License for more details.
455 *
456 * You should have received a copy of the GNU General Public License
457 * along with this program; if not, write to the Free Software
458 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
459 * 02110-1301 USA.
460 *
461 * SYNOPSIS
462 *
463 * #include <Rmath.h>
464 *
465 * double pnorm5(double x, double mu, double sigma, int lower_tail,int log_p);
466 * {pnorm (..) is synonymous and preferred inside R}
467 *
468 * void pnorm_both(double x, double *cum, double *ccum,
469 * int i_tail, int log_p);
470 *
471 * DESCRIPTION
472 *
473 * The main computation evaluates near-minimax approximations derived
474 * from those in "Rational Chebyshev approximations for the error
475 * function" by W. J. Cody, Math. Comp., 1969, 631-637. This
476 * transportable program uses rational functions that theoretically
477 * approximate the normal distribution function to at least 18
478 * significant decimal digits. The accuracy achieved depends on the
479 * arithmetic system, the compiler, the intrinsic functions, and
480 * proper selection of the machine-dependent constants.
481 *
482 * REFERENCE
483 *
484 * Cody, W. D. (1993).
485 * ALGORITHM 715: SPECFUN - A Portable FORTRAN Package of
486 * Special Function Routines and Test Drivers".
487 * ACM Transactions on Mathematical Software. 19, 22-32.
488 *
489 * EXTENSIONS
490 *
491 * The "_both" , lower, upper, and log_p variants were added by
492 * Martin Maechler, Jan.2000;
493 * as well as log1p() and similar improvements later on.
494 *
495 * James M. Rath contributed bug report PR#699 and patches correcting SIXTEN
496 * and if() clauses {with a bug: "|| instead of &&" -> PR #2883) more in line
497 * with the original Cody code.
498 */
499
pnorm(gnm_float x,gnm_float mu,gnm_float sigma,gboolean lower_tail,gboolean log_p)500 gnm_float pnorm(gnm_float x, gnm_float mu, gnm_float sigma, gboolean lower_tail, gboolean log_p)
501 {
502 gnm_float p, cp = 0.;
503
504 /* Note: The structure of these checks has been carefully thought through.
505 * For example, if x == mu and sigma == 0, we get the correct answer 1.
506 */
507 #ifdef IEEE_754
508 if(gnm_isnan(x) || gnm_isnan(mu) || gnm_isnan(sigma))
509 return x + mu + sigma;
510 #endif
511 if(!gnm_finite(x) && mu == x) return gnm_nan;/* x-mu is NaN */
512 if (sigma <= 0) {
513 if(sigma < 0) ML_ERR_return_NAN;
514 /* sigma = 0 : */
515 return (x < mu) ? R_DT_0 : R_DT_1;
516 }
517 p = (x - mu) / sigma;
518 if(!gnm_finite(p))
519 return (x < mu) ? R_DT_0 : R_DT_1;
520 x = p;
521
522 pnorm_both(x, &p, &cp, (lower_tail ? 0 : 1), log_p);
523
524 return(lower_tail ? p : cp);
525 }
526
527 #define SIXTEN 16 /* Cutoff allowing exact "*" and "/" */
528
pnorm_both(gnm_float x,gnm_float * cum,gnm_float * ccum,int i_tail,gboolean log_p)529 void pnorm_both(gnm_float x, gnm_float *cum, gnm_float *ccum, int i_tail, gboolean log_p)
530 {
531 /* i_tail in {0,1,2} means: "lower", "upper", or "both" :
532 if(lower) return *cum := P[X <= x]
533 if(upper) return *ccum := P[X > x] = 1 - P[X <= x]
534 */
535 static const gnm_float a[5] = {
536 GNM_const(2.2352520354606839287),
537 GNM_const(161.02823106855587881),
538 GNM_const(1067.6894854603709582),
539 GNM_const(18154.981253343561249),
540 GNM_const(0.065682337918207449113)
541 };
542 static const gnm_float b[4] = {
543 GNM_const(47.20258190468824187),
544 GNM_const(976.09855173777669322),
545 GNM_const(10260.932208618978205),
546 GNM_const(45507.789335026729956)
547 };
548 static const gnm_float c[9] = {
549 GNM_const(0.39894151208813466764),
550 GNM_const(8.8831497943883759412),
551 GNM_const(93.506656132177855979),
552 GNM_const(597.27027639480026226),
553 GNM_const(2494.5375852903726711),
554 GNM_const(6848.1904505362823326),
555 GNM_const(11602.651437647350124),
556 GNM_const(9842.7148383839780218),
557 GNM_const(1.0765576773720192317e-8)
558 };
559 static const gnm_float d[8] = {
560 GNM_const(22.266688044328115691),
561 GNM_const(235.38790178262499861),
562 GNM_const(1519.377599407554805),
563 GNM_const(6485.558298266760755),
564 GNM_const(18615.571640885098091),
565 GNM_const(34900.952721145977266),
566 GNM_const(38912.003286093271411),
567 GNM_const(19685.429676859990727)
568 };
569 static const gnm_float p[6] = {
570 GNM_const(0.21589853405795699),
571 GNM_const(0.1274011611602473639),
572 GNM_const(0.022235277870649807),
573 GNM_const(0.001421619193227893466),
574 GNM_const(2.9112874951168792e-5),
575 GNM_const(0.02307344176494017303)
576 };
577 static const gnm_float q[5] = {
578 GNM_const(1.28426009614491121),
579 GNM_const(0.468238212480865118),
580 GNM_const(0.0659881378689285515),
581 GNM_const(0.00378239633202758244),
582 GNM_const(7.29751555083966205e-5)
583 };
584
585 gnm_float xden, xnum, temp, del, eps, xsq, y;
586 #ifdef NO_DENORMS
587 gnm_float min = GNM_MIN;
588 #endif
589 int i, lower, upper;
590
591 #ifdef IEEE_754
592 if(gnm_isnan(x)) { *cum = *ccum = x; return; }
593 #endif
594
595 /* Consider changing these : */
596 eps = GNM_EPSILON * 0.5;
597
598 /* i_tail in {0,1,2} =^= {lower, upper, both} */
599 lower = i_tail != 1;
600 upper = i_tail != 0;
601
602 y = gnm_abs(x);
603 if (y <= GNM_const(0.67448975)) { /* qnorm(3/4) = .6744.... -- earlier had GNM_const(0.66291) */
604 if (y > eps) {
605 xsq = x * x;
606 xnum = a[4] * xsq;
607 xden = xsq;
608 for (i = 0; i < 3; ++i) {
609 xnum = (xnum + a[i]) * xsq;
610 xden = (xden + b[i]) * xsq;
611 }
612 } else xnum = xden = 0.0;
613
614 temp = x * (xnum + a[3]) / (xden + b[3]);
615 if(lower) *cum = 0.5 + temp;
616 if(upper) *ccum = 0.5 - temp;
617 if(log_p) {
618 if(lower) *cum = gnm_log(*cum);
619 if(upper) *ccum = gnm_log(*ccum);
620 }
621 }
622 else if (y <= M_SQRT_32) {
623
624 /* Evaluate pnorm for 0.674.. = qnorm(3/4) < |x| <= gnm_sqrt(32) ~= 5.657 */
625
626 xnum = c[8] * y;
627 xden = y;
628 for (i = 0; i < 7; ++i) {
629 xnum = (xnum + c[i]) * y;
630 xden = (xden + d[i]) * y;
631 }
632 temp = (xnum + c[7]) / (xden + d[7]);
633
634 #define do_del(X) \
635 xsq = gnm_trunc(X * SIXTEN) / SIXTEN; \
636 del = (X - xsq) * (X + xsq); \
637 if(log_p) { \
638 *cum = (-xsq * xsq * 0.5) + (-del * 0.5) + gnm_log(temp); \
639 if((lower && x > 0.) || (upper && x <= 0.)) \
640 *ccum = gnm_log1p(-gnm_exp(-xsq * xsq * 0.5) * \
641 gnm_exp(-del * 0.5) * temp); \
642 } \
643 else { \
644 *cum = gnm_exp(-xsq * xsq * 0.5) * gnm_exp(-del * 0.5) * temp; \
645 *ccum = 1.0 - *cum; \
646 }
647
648 #define swap_tail \
649 if (x > 0.) {/* swap ccum <--> cum */ \
650 temp = *cum; if(lower) *cum = *ccum; *ccum = temp; \
651 }
652
653 do_del(y);
654 swap_tail;
655 }
656
657 /* else |x| > gnm_sqrt(32) = 5.657 :
658 * the next two case differentiations were really for lower=T, log=F
659 * Particularly *not* for log_p !
660
661 * Cody had (-37.5193 < x && x < 8.2924) ; R originally had y < 50
662 *
663 * Note that we do want symmetry(0), lower/upper -> hence use y
664 */
665 else if(log_p
666 /* ^^^^^ MM FIXME: can speedup for log_p and much larger |x| !
667 * Then, make use of Abramowitz & Stegun, 26.2.13, something like
668
669 xsq = x*x;
670
671 if(xsq * GNM_EPSILON < 1.)
672 del = (1. - (1. - 5./(xsq+6.)) / (xsq+4.)) / (xsq+2.);
673 else
674 del = 0.;
675 *cum = -.5*xsq - M_LN_SQRT_2PI - gnm_log(x) + gnm_log1p(-del);
676 *ccum = gnm_log1p(-gnm_exp(*cum)); /.* ~ gnm_log(1) = 0 *./
677
678 swap_tail;
679
680 */
681 || (lower && -37.5193 < x && x < 8.2924)
682 || (upper && -8.2924 < x && x < 37.5193)
683 ) {
684
685 /* Evaluate pnorm for x in (-37.5, -5.657) union (5.657, 37.5) */
686 xsq = 1.0 / (x * x);
687 xnum = p[5] * xsq;
688 xden = xsq;
689 for (i = 0; i < 4; ++i) {
690 xnum = (xnum + p[i]) * xsq;
691 xden = (xden + q[i]) * xsq;
692 }
693 temp = xsq * (xnum + p[4]) / (xden + q[4]);
694 temp = (M_1_SQRT_2PI - temp) / y;
695
696 do_del(x);
697 swap_tail;
698 }
699 else { /* no log_p , large x such that probs are 0 or 1 */
700 if(x > 0) { *cum = 1.; *ccum = 0.; }
701 else { *cum = 0.; *ccum = 1.; }
702 }
703
704
705 #ifdef NO_DENORMS
706 /* do not return "denormalized" -- we do in R */
707 if(log_p) {
708 if(*cum > -min) *cum = -0.;
709 if(*ccum > -min)*ccum = -0.;
710 }
711 else {
712 if(*cum < min) *cum = 0.;
713 if(*ccum < min) *ccum = 0.;
714 }
715 #endif
716 return;
717 }
718 /* Cleaning up done by tools/import-R: */
719 #undef SIXTEN
720 #undef do_del
721 #undef swap_tail
722
723 /* ------------------------------------------------------------------------ */
724 /* Imported src/nmath/qnorm.c from R. */
725 /*
726 * Mathlib : A C Library of Special Functions
727 * Copyright (C) 1998 Ross Ihaka
728 * Copyright (C) 2000 The R Development Core Team
729 * based on AS 111 (C) 1977 Royal Statistical Society
730 * and on AS 241 (C) 1988 Royal Statistical Society
731 *
732 * This program is free software; you can redistribute it and/or modify
733 * it under the terms of the GNU General Public License as published by
734 * the Free Software Foundation; either version 2 of the License, or
735 * (at your option) any later version.
736 *
737 * This program is distributed in the hope that it will be useful,
738 * but WITHOUT ANY WARRANTY; without even the implied warranty of
739 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
740 * GNU General Public License for more details.
741 *
742 * You should have received a copy of the GNU General Public License
743 * along with this program; if not, write to the Free Software
744 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
745 * 02110-1301 USA.
746 *
747 * SYNOPSIS
748 *
749 * double qnorm5(double p, double mu, double sigma,
750 * int lower_tail, int log_p)
751 * {qnorm (..) is synonymous and preferred inside R}
752 *
753 * DESCRIPTION
754 *
755 * Compute the quantile function for the normal distribution.
756 *
757 * For small to moderate probabilities, algorithm referenced
758 * below is used to obtain an initial approximation which is
759 * polished with a final Newton step.
760 *
761 * For very large arguments, an algorithm of Wichura is used.
762 *
763 * REFERENCE
764 *
765 * Beasley, J. D. and S. G. Springer (1977).
766 * Algorithm AS 111: The percentage points of the normal distribution,
767 * Applied Statistics, 26, 118-121.
768 *
769 * Wichura, M.J. (1988).
770 * Algorithm AS 241: The Percentage Points of the Normal Distribution.
771 * Applied Statistics, 37, 477-484.
772 */
773
774
qnorm(gnm_float p,gnm_float mu,gnm_float sigma,gboolean lower_tail,gboolean log_p)775 gnm_float qnorm(gnm_float p, gnm_float mu, gnm_float sigma, gboolean lower_tail, gboolean log_p)
776 {
777 gnm_float p_, q, r, val;
778
779 #ifdef IEEE_754
780 if (gnm_isnan(p) || gnm_isnan(mu) || gnm_isnan(sigma))
781 return p + mu + sigma;
782 #endif
783 if (p == R_DT_0) return gnm_ninf;
784 if (p == R_DT_1) return gnm_pinf;
785 R_Q_P01_check(p);
786
787 if(sigma < 0) ML_ERR_return_NAN;
788 if(sigma == 0) return mu;
789
790 p_ = R_DT_qIv(p);/* real lower_tail prob. p */
791 q = p_ - 0.5;
792
793 #ifdef DEBUG_qnorm
794 REprintf("qnorm(p=%10.7" GNM_FORMAT_g ", m=%" GNM_FORMAT_g ", s=%" GNM_FORMAT_g ", l.t.= %d, log= %d): q = %" GNM_FORMAT_g "\n",
795 p,mu,sigma, lower_tail, log_p, q);
796 #endif
797
798
799 /*-- use AS 241 --- */
800 /* gnm_float ppnd16_(gnm_float *p, long *ifault)*/
801 /* ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3
802
803 Produces the normal deviate Z corresponding to a given lower
804 tail area of P; Z is accurate to about 1 part in 10**16.
805
806 (original fortran code used PARAMETER(..) for the coefficients
807 and provided hash codes for checking them...)
808 */
809 if (gnm_abs(q) <= .425) {/* 0.075 <= p <= 0.925 */
810 r = GNM_const(.180625) - q * q;
811 val =
812 q * (((((((r * GNM_const(2509.0809287301226727) +
813 GNM_const(33430.575583588128105)) * r + GNM_const(67265.770927008700853)) * r +
814 GNM_const(45921.953931549871457)) * r + GNM_const(13731.693765509461125)) * r +
815 GNM_const(1971.5909503065514427)) * r + GNM_const(133.14166789178437745)) * r +
816 GNM_const(3.387132872796366608))
817 / (((((((r * GNM_const(5226.495278852854561) +
818 GNM_const(28729.085735721942674)) * r + GNM_const(39307.89580009271061)) * r +
819 GNM_const(21213.794301586595867)) * r + GNM_const(5394.1960214247511077)) * r +
820 GNM_const(687.1870074920579083)) * r + GNM_const(42.313330701600911252)) * r + 1.);
821 }
822 else { /* closer than 0.075 from {0,1} boundary */
823
824 /* r = min(p, 1-p) < 0.075 */
825 if (q > 0)
826 r = R_DT_CIv(p);/* 1-p */
827 else
828 r = p_;/* = R_DT_Iv(p) ^= p */
829
830 r = gnm_sqrt(- ((log_p &&
831 ((lower_tail && q <= 0) || (!lower_tail && q > 0))) ?
832 p : /* else */ gnm_log(r)));
833 /* r = gnm_sqrt(-gnm_log(r)) <==> min(p, 1-p) = gnm_exp( - r^2 ) */
834 #ifdef DEBUG_qnorm
835 REprintf("\t close to 0 or 1: r = %7" GNM_FORMAT_g "\n", r);
836 #endif
837
838 if (r <= 5.) { /* <==> min(p,1-p) >= gnm_exp(-25) ~= 1.3888e-11 */
839 r += -1.6;
840 val = (((((((r * GNM_const(7.7454501427834140764e-4) +
841 GNM_const(.0227238449892691845833)) * r + GNM_const(.24178072517745061177)) *
842 r + GNM_const(1.27045825245236838258)) * r +
843 GNM_const(3.64784832476320460504)) * r + GNM_const(5.7694972214606914055)) *
844 r + GNM_const(4.6303378461565452959)) * r +
845 GNM_const(1.42343711074968357734))
846 / (((((((r *
847 GNM_const(1.05075007164441684324e-9) + GNM_const(5.475938084995344946e-4)) *
848 r + GNM_const(.0151986665636164571966)) * r +
849 GNM_const(.14810397642748007459)) * r + GNM_const(.68976733498510000455)) *
850 r + GNM_const(1.6763848301838038494)) * r +
851 GNM_const(2.05319162663775882187)) * r + 1.);
852 }
853 else { /* very close to 0 or 1 */
854 r += -5.;
855 val = (((((((r * GNM_const(2.01033439929228813265e-7) +
856 GNM_const(2.71155556874348757815e-5)) * r +
857 GNM_const(.0012426609473880784386)) * r + GNM_const(.026532189526576123093)) *
858 r + GNM_const(.29656057182850489123)) * r +
859 GNM_const(1.7848265399172913358)) * r + GNM_const(5.4637849111641143699)) *
860 r + GNM_const(6.6579046435011037772))
861 / (((((((r *
862 GNM_const(2.04426310338993978564e-15) + GNM_const(1.4215117583164458887e-7))*
863 r + GNM_const(1.8463183175100546818e-5)) * r +
864 GNM_const(7.868691311456132591e-4)) * r + GNM_const(.0148753612908506148525))
865 * r + GNM_const(.13692988092273580531)) * r +
866 GNM_const(.59983220655588793769)) * r + 1.);
867 }
868
869 if(q < 0.0)
870 val = -val;
871 /* return (q >= 0.)? r : -r ;*/
872 }
873 return mu + sigma * val;
874 }
875
876
877 /* ------------------------------------------------------------------------ */
878 /* Imported src/nmath/ppois.c from R. */
879 /*
880 * Mathlib : A C Library of Special Functions
881 * Copyright (C) 1998 Ross Ihaka
882 * Copyright (C) 2000 The R Development Core Team
883 *
884 * This program is free software; you can redistribute it and/or modify
885 * it under the terms of the GNU General Public License as published by
886 * the Free Software Foundation; either version 2 of the License, or
887 * (at your option) any later version.
888 *
889 * This program is distributed in the hope that it will be useful,
890 * but WITHOUT ANY WARRANTY; without even the implied warranty of
891 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
892 * GNU General Public License for more details.
893 *
894 * You should have received a copy of the GNU General Public License
895 * along with this program; if not, write to the Free Software
896 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
897 * USA.
898 *
899 * DESCRIPTION
900 *
901 * The distribution function of the Poisson distribution.
902 */
903
904
ppois(gnm_float x,gnm_float lambda,gboolean lower_tail,gboolean log_p)905 gnm_float ppois(gnm_float x, gnm_float lambda, gboolean lower_tail, gboolean log_p)
906 {
907 #ifdef IEEE_754
908 if (gnm_isnan(x) || gnm_isnan(lambda))
909 return x + lambda;
910 #endif
911 if(lambda < 0.) ML_ERR_return_NAN;
912
913 x = gnm_fake_floor(x);
914 if (x < 0) return R_DT_0;
915 if (lambda == 0.) return R_DT_1;
916 if (!gnm_finite(x)) return R_DT_1;
917
918 return pgamma(lambda, x + 1, 1., !lower_tail, log_p);
919 }
920
921 /* ------------------------------------------------------------------------ */
922 /* Imported src/nmath/dpois.c from R. */
923 /*
924 * AUTHOR
925 * Catherine Loader, catherine@research.bell-labs.com.
926 * October 23, 2000.
927 *
928 * Merge in to R:
929 * Copyright (C) 2000, The R Core Development Team
930 *
931 * This program is free software; you can redistribute it and/or modify
932 * it under the terms of the GNU General Public License as published by
933 * the Free Software Foundation; either version 2 of the License, or
934 * (at your option) any later version.
935 *
936 * This program is distributed in the hope that it will be useful,
937 * but WITHOUT ANY WARRANTY; without even the implied warranty of
938 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
939 * GNU General Public License for more details.
940 *
941 * You should have received a copy of the GNU General Public License
942 * along with this program; if not, write to the Free Software
943 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
944 * USA.
945 *
946 *
947 * DESCRIPTION
948 *
949 * dpois() checks argument validity and calls dpois_raw().
950 *
951 * dpois_raw() computes the Poisson probability lb^x exp(-lb) / x!.
952 * This does not check that x is an integer, since dgamma() may
953 * call this with a fractional x argument. Any necessary argument
954 * checks should be done in the calling function.
955 *
956 */
957
958
dpois_raw(gnm_float x,gnm_float lambda,gboolean give_log)959 static gnm_float dpois_raw(gnm_float x, gnm_float lambda, gboolean give_log)
960 {
961 gnm_float yh, yl, f2;
962
963 /* x >= 0 ; integer for dpois(), but not e.g. for pgamma()!
964 lambda >= 0
965 */
966 if (lambda == 0) return( (x == 0) ? R_D__1 : R_D__0 );
967 if (!gnm_finite(lambda)) return R_D__0;
968 if (x < 0) return( R_D__0 );
969 if (x <= lambda * GNM_MIN) return(R_D_exp(-lambda) );
970 if (lambda < x * GNM_MIN) return(R_D_exp(-lambda + x*gnm_log(lambda) -lgamma1p (x)));
971
972 ebd0 (x, lambda, &yh, &yl);
973 PAIR_ADD (stirlerr(x), yh, yl);
974 f2 = M_2PIgnum * x;
975
976 return give_log
977 ? -yl - yh - 0.5 * gnm_log (f2)
978 : gnm_exp (-yl) * gnm_exp (-yh) / gnm_sqrt (f2);
979 }
980
dpois(gnm_float x,gnm_float lambda,gboolean give_log)981 gnm_float dpois(gnm_float x, gnm_float lambda, gboolean give_log)
982 {
983 #ifdef IEEE_754
984 if(gnm_isnan(x) || gnm_isnan(lambda))
985 return x + lambda;
986 #endif
987
988 if (lambda < 0) ML_ERR_return_NAN;
989 R_D_nonint_check(x);
990 if (x < 0 || !gnm_finite(x))
991 return R_D__0;
992
993 x = R_D_forceint(x);
994
995 return( dpois_raw(x,lambda,give_log) );
996 }
997
998 /* ------------------------------------------------------------------------ */
999 /* Imported src/nmath/dgamma.c from R. */
1000 /*
1001 * AUTHOR
1002 * Catherine Loader, catherine@research.bell-labs.com.
1003 * October 23, 2000.
1004 *
1005 * Merge in to R:
1006 * Copyright (C) 2000 The R Core Development Team
1007 * Copyright (C) 2004 The R Foundation
1008 *
1009 * This program is free software; you can redistribute it and/or modify
1010 * it under the terms of the GNU General Public License as published by
1011 * the Free Software Foundation; either version 2 of the License, or
1012 * (at your option) any later version.
1013 *
1014 * This program is distributed in the hope that it will be useful,
1015 * but WITHOUT ANY WARRANTY; without even the implied warranty of
1016 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1017 * GNU General Public License for more details.
1018 *
1019 * You should have received a copy of the GNU General Public License
1020 * along with this program; if not, write to the Free Software
1021 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
1022 * USA.
1023 *
1024 *
1025 * DESCRIPTION
1026 *
1027 * Computes the density of the gamma distribution,
1028 *
1029 * 1/s (x/s)^{a-1} exp(-x/s)
1030 * p(x;a,s) = -----------------------
1031 * (a-1)!
1032 *
1033 * where `s' is the scale (= 1/lambda in other parametrizations)
1034 * and `a' is the shape parameter ( = alpha in other contexts).
1035 *
1036 * The old (R 1.1.1) version of the code is available via `#define D_non_pois'
1037 */
1038
1039
dgamma(gnm_float x,gnm_float shape,gnm_float scale,gboolean give_log)1040 gnm_float dgamma(gnm_float x, gnm_float shape, gnm_float scale, gboolean give_log)
1041 {
1042 gnm_float pr;
1043 #ifdef IEEE_754
1044 if (gnm_isnan(x) || gnm_isnan(shape) || gnm_isnan(scale))
1045 return x + shape + scale;
1046 #endif
1047 if (shape <= 0 || scale <= 0) ML_ERR_return_NAN;
1048 if (x < 0)
1049 return R_D__0;
1050 if (x == 0) {
1051 if (shape < 1) return gnm_pinf;
1052 if (shape > 1) return R_D__0;
1053 /* else */
1054 return give_log ? -gnm_log(scale) : 1 / scale;
1055 }
1056
1057 if (shape < 1) {
1058 pr = dpois_raw(shape, x/scale, give_log);
1059 return give_log ? pr + gnm_log(shape/x) : pr*shape/x;
1060 }
1061 /* else shape >= 1 */
1062 pr = dpois_raw(shape-1, x/scale, give_log);
1063 return give_log ? pr - gnm_log(scale) : pr/scale;
1064 }
1065
1066 /* ------------------------------------------------------------------------ */
1067 /* Imported src/nmath/pgamma.c from R. */
1068 /*
1069 * Mathlib : A C Library of Special Functions
1070 * Copyright (C) 2005 Morten Welinder <terra@gnome.org>
1071 * Copyright (C) 2005 The R Foundation
1072 *
1073 * This program is free software; you can redistribute it and/or modify
1074 * it under the terms of the GNU General Public License as published by
1075 * the Free Software Foundation; either version 2 of the License, or
1076 * (at your option) any later version.
1077 *
1078 * This program is distributed in the hope that it will be useful,
1079 * but WITHOUT ANY WARRANTY; without even the implied warranty of
1080 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1081 * GNU General Public License for more details.
1082 *
1083 * A copy of the GNU General Public License is available via WWW at
1084 * http://www.gnu.org/copyleft/gpl.html. You can also obtain it by
1085 * writing to the Free Software Foundation, Inc., 51 Franklin St, Fifth
1086 * Floor, Boston, MA 02110-1301 USA.
1087 *
1088 * SYNOPSIS
1089 *
1090 * #include <Rmath.h>
1091 *
1092 * double pgamma (double x, double alph, double scale,
1093 * int lower_tail, int log_p)
1094 *
1095 * double log1pmx (double x)
1096 * double lgamma1p (double a)
1097 *
1098 * double logspace_add (double logx, double logy)
1099 * double logspace_sub (double logx, double logy)
1100 *
1101 *
1102 * DESCRIPTION
1103 *
1104 * This function computes the distribution function for the
1105 * gamma distribution with shape parameter alph and scale parameter
1106 * scale. This is also known as the incomplete gamma function.
1107 * See Abramowitz and Stegun (6.5.1) for example.
1108 *
1109 * NOTES
1110 *
1111 * Complete redesign by Morten Welinder, originally for Gnumeric.
1112 * Improvements (e.g. "while NEEDED_SCALE") by Martin Maechler
1113 * The old version can be activated by compiling with -DR_USE_OLD_PGAMMA
1114 *
1115 * REFERENCES
1116 *
1117 */
1118
1119 /*----------- DEBUGGING -------------
1120 * make CFLAGS='-DDEBUG_p -g -I/usr/local/include -I../include'
1121 */
1122
1123
1124 /* Scalefactor:= (2^32)^8 = 2^256 = GNM_const(1.157921e+77) */
1125 #define SQR(x) ((x)*(x))
1126 static const gnm_float scalefactor = SQR(SQR(SQR(4294967296.0)));
1127 #undef SQR
1128
1129 /* If |x| > |k| * M_cutoff, then log[ gnm_exp(-x) * k^x ] =~= -x */
1130 static const gnm_float M_cutoff = M_LN2gnum * GNM_MAX_EXP / GNM_EPSILON;/*=GNM_const(3.196577e18)*/
1131
1132 /* Continued fraction for calculation of
1133 * 1/i + x/(i+d) + x^2/(i+2*d) + x^3/(i+3*d) + ...
1134 *
1135 * auxilary in log1pmx() and lgamma1p()
1136 */
1137 gnm_float
gnm_logcf(gnm_float x,gnm_float i,gnm_float d)1138 gnm_logcf (gnm_float x, gnm_float i, gnm_float d)
1139 {
1140 gnm_float c1 = 2 * d;
1141 gnm_float c2 = i + d;
1142 gnm_float c4 = c2 + d;
1143 gnm_float a1 = c2;
1144 gnm_float b1 = i * (c2 - i * x);
1145 gnm_float b2 = d * d * x;
1146 gnm_float a2 = c4 * c2 - b2;
1147 const gnm_float cfVSmall = 1.0e-14;/* ~= relative tolerance */
1148
1149 #if 0
1150 assert (i > 0);
1151 assert (d >= 0);
1152 #endif
1153
1154 b2 = c4 * b1 - i * b2;
1155
1156 while (gnm_abs (a2 * b1 - a1 * b2) > gnm_abs (cfVSmall * b1 * b2)) {
1157 gnm_float c3 = c2*c2*x;
1158 c2 += d;
1159 c4 += d;
1160 a1 = c4 * a2 - c3 * a1;
1161 b1 = c4 * b2 - c3 * b1;
1162
1163 c3 = c1 * c1 * x;
1164 c1 += d;
1165 c4 += d;
1166 a2 = c4 * a1 - c3 * a2;
1167 b2 = c4 * b1 - c3 * b2;
1168
1169 if (gnm_abs (b2) > scalefactor) {
1170 a1 /= scalefactor;
1171 b1 /= scalefactor;
1172 a2 /= scalefactor;
1173 b2 /= scalefactor;
1174 } else if (gnm_abs (b2) < 1 / scalefactor) {
1175 a1 *= scalefactor;
1176 b1 *= scalefactor;
1177 a2 *= scalefactor;
1178 b2 *= scalefactor;
1179 }
1180 }
1181
1182 return a2 / b2;
1183 }
1184
1185 /* Accurate calculation of gnm_log1p (x)-x, particularly for small x. */
log1pmx(gnm_float x)1186 gnm_float log1pmx (gnm_float x)
1187 {
1188 static const gnm_float minLog1Value = GNM_const(-0.79149064);
1189 static const gnm_float two = 2;
1190
1191 if (x > 1 || x < minLog1Value)
1192 return gnm_log1p(x) - x;
1193 else { /* expand in [x/(2+x)]^2 */
1194 gnm_float term = x / (2 + x);
1195 gnm_float y = term * term;
1196 if (gnm_abs(x) < 1e-2)
1197 return term * ((((two / 9 * y + two / 7) * y + two / 5) * y +
1198 two / 3) * y - x);
1199 else
1200 return term * (2 * y * gnm_logcf (y, 3, 2) - x);
1201 }
1202 }
1203
1204
1205 /*
1206 * Compute the log of a sum from logs of terms, i.e.,
1207 *
1208 * log (exp (logx) + exp (logy))
1209 *
1210 * without causing overflows and without throwing away large handfuls
1211 * of accuracy.
1212 */
logspace_add(gnm_float logx,gnm_float logy)1213 gnm_float logspace_add (gnm_float logx, gnm_float logy)
1214 {
1215 return fmax2 (logx, logy) + gnm_log1p (gnm_exp (-gnm_abs (logx - logy)));
1216 }
1217
1218
1219 /*
1220 * Compute the log of a difference from logs of terms, i.e.,
1221 *
1222 * log (exp (logx) - exp (logy))
1223 *
1224 * without causing overflows and without throwing away large handfuls
1225 * of accuracy.
1226 */
logspace_sub(gnm_float logx,gnm_float logy)1227 gnm_float logspace_sub (gnm_float logx, gnm_float logy)
1228 {
1229 return logx + gnm_log1p (-gnm_exp (logy - logx));
1230 }
1231
1232
1233 #ifndef R_USE_OLD_PGAMMA
1234
1235 /* dpois_wrap (x_P_1, lambda, g_log) ==
1236 * dpois (x_P_1 - 1, lambda, g_log)
1237 */
1238 static gnm_float
dpois_wrap(gnm_float x_plus_1,gnm_float lambda,gboolean give_log)1239 dpois_wrap (gnm_float x_plus_1, gnm_float lambda, gboolean give_log)
1240 {
1241 #ifdef DEBUG_p
1242 REprintf (" dpois_wrap(x+1=%.14" GNM_FORMAT_g ", lambda=%.14" GNM_FORMAT_g ", log=%d)\n",
1243 x_plus_1, lambda, give_log);
1244 #endif
1245 if (!gnm_finite(lambda))
1246 return R_D__0;
1247 if (x_plus_1 > 1)
1248 return dpois_raw (x_plus_1 - 1, lambda, give_log);
1249 if (lambda > gnm_abs(x_plus_1 - 1) * M_cutoff)
1250 return R_D_exp(-lambda - gnm_lgamma(x_plus_1));
1251 else {
1252 gnm_float d = dpois_raw (x_plus_1, lambda, give_log);
1253 #ifdef DEBUG_p
1254 REprintf (" -> d=dpois_raw(..)=%.14" GNM_FORMAT_g "\n", d);
1255 #endif
1256 return give_log
1257 ? d + gnm_log (x_plus_1 / lambda)
1258 : d * (x_plus_1 / lambda);
1259 }
1260 }
1261
1262 /*
1263 * Abramowitz and Stegun 6.5.29 [right]
1264 */
1265 static gnm_float
pgamma_smallx(gnm_float x,gnm_float alph,gboolean lower_tail,gboolean log_p)1266 pgamma_smallx (gnm_float x, gnm_float alph, gboolean lower_tail, gboolean log_p)
1267 {
1268 gnm_float sum = 0, c = alph, n = 0, term;
1269
1270 #ifdef DEBUG_p
1271 REprintf (" pg_smallx(x=%.12" GNM_FORMAT_g ", alph=%.12" GNM_FORMAT_g "): ", x, alph);
1272 #endif
1273
1274 /*
1275 * Relative to 6.5.29 all terms have been multiplied by alph
1276 * and the first, thus being 1, is omitted.
1277 */
1278
1279 do {
1280 n++;
1281 c *= -x / n;
1282 term = c / (alph + n);
1283 sum += term;
1284 } while (gnm_abs (term) > GNM_EPSILON * gnm_abs (sum));
1285
1286 #ifdef DEBUG_p
1287 REprintf (" conv.sum=%" GNM_FORMAT_g ";", sum);
1288 #endif
1289 if (lower_tail) {
1290 gnm_float f1 = log_p ? gnm_log1p (sum) : 1 + sum;
1291 gnm_float f2;
1292 if (alph > 1) {
1293 f2 = dpois_raw (alph, x, log_p);
1294 f2 = log_p ? f2 + x : f2 * gnm_exp (x);
1295 } else if (log_p)
1296 f2 = alph * gnm_log (x) - lgamma1p (alph);
1297 else
1298 f2 = gnm_pow (x, alph) / gnm_exp (lgamma1p (alph));
1299 #ifdef DEBUG_p
1300 REprintf (" (f1,f2)= (%" GNM_FORMAT_g ",%" GNM_FORMAT_g ")\n", f1,f2);
1301 #endif
1302 return log_p ? f1 + f2 : f1 * f2;
1303 } else {
1304 gnm_float lf2 = alph * gnm_log (x) - lgamma1p (alph);
1305 #ifdef DEBUG_p
1306 REprintf (" 1:%.14" GNM_FORMAT_g " 2:%.14" GNM_FORMAT_g "\n", alph * gnm_log (x), lgamma1p (alph));
1307 REprintf (" sum=%.14" GNM_FORMAT_g " gnm_log1p (sum)=%.14" GNM_FORMAT_g " lf2=%.14" GNM_FORMAT_g "\n",
1308 sum, gnm_log1p (sum), lf2);
1309 #endif
1310 if (log_p)
1311 return swap_log_tail (gnm_log1p (sum) + lf2);
1312 else {
1313 gnm_float f1m1 = sum;
1314 gnm_float f2m1 = gnm_expm1 (lf2);
1315 return -(f1m1 + f2m1 + f1m1 * f2m1);
1316 }
1317 }
1318 } /* pgamma_smallx() */
1319
1320 static gnm_float
pd_upper_series(gnm_float x,gnm_float y,gboolean log_p)1321 pd_upper_series (gnm_float x, gnm_float y, gboolean log_p)
1322 {
1323 gnm_float term = x / y;
1324 gnm_float sum = term;
1325
1326 do {
1327 y++;
1328 term *= x / y;
1329 sum += term;
1330 } while (term > sum * GNM_EPSILON);
1331
1332 /* sum = \sum_{n=1}^ oo x^n / (y*(y+1)*...*(y+n-1))
1333 * = \sum_{n=0}^ oo x^(n+1) / (y*(y+1)*...*(y+n))
1334 * = x/y * (1 + \sum_{n=1}^oo x^n / ((y+1)*...*(y+n)))
1335 * ~ x/y + o(x/y) {which happens when alph -> Inf}
1336 */
1337 return log_p ? gnm_log (sum) : sum;
1338 }
1339
1340 /* Continued fraction for calculation of
1341 * ???
1342 * = (i / d) + o(i/d)
1343 */
1344 static gnm_float
pd_lower_cf(gnm_float i,gnm_float d)1345 pd_lower_cf (gnm_float i, gnm_float d)
1346 {
1347 gnm_float f = -42, of, rf = i / d;
1348
1349 gnm_float c1 = 0, c2, c3, c4;
1350 gnm_float a1 = 0, b1 = 1;
1351 gnm_float a2 = i, b2 = d;
1352
1353 #define NEEDED_SCALE \
1354 (b2 > scalefactor) { \
1355 a1 /= scalefactor; \
1356 b1 /= scalefactor; \
1357 a2 /= scalefactor; \
1358 b2 /= scalefactor; \
1359 }
1360
1361 #define max_it 200000
1362
1363 #ifdef DEBUG_p
1364 REprintf("pd_lower_cf(i=%.14" GNM_FORMAT_g ", d=%.14" GNM_FORMAT_g ")\n", i, d);
1365 #endif
1366
1367 while NEEDED_SCALE
1368
1369 if(a2 == 0)
1370 return 0;/* when d >>>> i originally */
1371
1372 c2 = a2;
1373 c4 = b2;
1374
1375 while (c1 < max_it) {
1376 c1++;
1377 c2--;
1378 c3 = c1 * c2;
1379 c4 += 2;
1380 a1 = c4 * a2 + c3 * a1;
1381 b1 = c4 * b2 + c3 * b1;
1382
1383 c1++;
1384 c2--;
1385 c3 = c1 * c2;
1386 c4 += 2;
1387 a2 = c4 * a1 + c3 * a2;
1388 b2 = c4 * b1 + c3 * b2;
1389
1390 if NEEDED_SCALE
1391
1392 if (b2 != 0) {
1393 of = f;
1394 f = a2 / b2;
1395 /* convergence check: relative; absolute for small f : */
1396 if (gnm_abs (f - of) <= GNM_EPSILON * fmax2(rf, gnm_abs(f)))
1397 return f;
1398 }
1399 }
1400
1401 REprintf(" ** NON-convergence in pgamma()'s pd_lower_cf() f= %" GNM_FORMAT_g ".\n", f);
1402 return f;/* should not happen ... */
1403 } /* pd_lower_cf() */
1404 #undef NEEDED_SCALE
1405
1406
1407 static gnm_float
pd_lower_series(gnm_float lambda,gnm_float y)1408 pd_lower_series (gnm_float lambda, gnm_float y)
1409 {
1410 gnm_float term = 1, sum = 0;
1411
1412 #ifdef DEBUG_p
1413 REprintf("pd_lower_series(lam=%.14" GNM_FORMAT_g ", y=%.14" GNM_FORMAT_g ") ...", lambda, y);
1414 #endif
1415 while (y >= 1 && term > sum * GNM_EPSILON) {
1416 term *= y / lambda;
1417 sum += term;
1418 y--;
1419 }
1420 /* sum = \sum_{n=0}^ oo y*(y-1)*...*(y - n) / lambda^(n+1)
1421 * = y/lambda * (1 + \sum_{n=1}^Inf (y-1)*...*(y-n) / lambda^n
1422 * ~ y/lambda + o(y/lambda)
1423 */
1424 #ifdef DEBUG_p
1425 REprintf(" done: term=%" GNM_FORMAT_g ", sum=%" GNM_FORMAT_g ", y= %" GNM_FORMAT_g "\n", term, sum, y);
1426 #endif
1427
1428 if (y != gnm_floor (y)) {
1429 /*
1430 * The series does not converge as the terms start getting
1431 * bigger (besides flipping sign) for y < -lambda.
1432 */
1433 gnm_float f;
1434 #ifdef DEBUG_p
1435 REprintf(" y not int: add another term ");
1436 #endif
1437 f = pd_lower_cf (y, lambda + 1 - y);
1438 #ifdef DEBUG_p
1439 REprintf(" (= %.14" GNM_FORMAT_g ") * term = %.14" GNM_FORMAT_g " to sum %" GNM_FORMAT_g "\n", f, term * f, sum);
1440 #endif
1441 sum += term * f;
1442 }
1443
1444 return sum;
1445 } /* pd_lower_series() */
1446
1447 /*
1448 * Asymptotic expansion to calculate the probability that poisson variate
1449 * has value <= x.
1450 */
1451 static gnm_float
ppois_asymp(gnm_float x,gnm_float lambda,gboolean lower_tail,gboolean log_p)1452 ppois_asymp (gnm_float x, gnm_float lambda, gboolean lower_tail, gboolean log_p)
1453 {
1454 static const gnm_float coef15 = 1 / GNM_const(12.);
1455 static const gnm_float coef25 = 1 / GNM_const(288.);
1456 static const gnm_float coef35 = -139 / GNM_const(51840.);
1457 static const gnm_float coef45 = -571 / GNM_const(2488320.);
1458 static const gnm_float coef55 = 163879 / GNM_const(209018880.);
1459 static const gnm_float coef65 = 5246819 / GNM_const(75246796800.);
1460 static const gnm_float coef75 = -534703531 / GNM_const(902961561600.);
1461 static const gnm_float coef1 = 2 / GNM_const(3.);
1462 static const gnm_float coef2 = -4 / GNM_const(135.);
1463 static const gnm_float coef3 = 8 / GNM_const(2835.);
1464 static const gnm_float coef4 = 16 / GNM_const(8505.);
1465 static const gnm_float coef5 = -8992 / GNM_const(12629925.);
1466 static const gnm_float coef6 = -334144 / GNM_const(492567075.);
1467 static const gnm_float coef7 = 698752 / GNM_const(1477701225.);
1468 static const gnm_float two = 2;
1469
1470 gnm_float dfm, pt_,s2pt,res1,res2,elfb,term;
1471 gnm_float ig2,ig3,ig4,ig5,ig6,ig7,ig25,ig35,ig45,ig55,ig65,ig75;
1472 gnm_float f, np, nd;
1473
1474 dfm = lambda - x;
1475 pt_ = -x * log1pmx (dfm / x);
1476 s2pt = gnm_sqrt (2 * pt_);
1477 if (dfm < 0) s2pt = -s2pt;
1478
1479 ig2 = 1.0 + pt_;
1480 term = pt_ * pt_ * 0.5;
1481 ig3 = ig2 + term;
1482 term *= pt_ / 3;
1483 ig4 = ig3 + term;
1484 term *= pt_ / 4;
1485 ig5 = ig4 + term;
1486 term *= pt_ / 5;
1487 ig6 = ig5 + term;
1488 term *= pt_ / 6;
1489 ig7 = ig6 + term;
1490
1491 term = pt_ * (two / 3);
1492 ig25 = 1.0 + term;
1493 term *= pt_ * (two / 5);
1494 ig35 = ig25 + term;
1495 term *= pt_ * (two / 7);
1496 ig45 = ig35 + term;
1497 term *= pt_ * (two / 9);
1498 ig55 = ig45 + term;
1499 term *= pt_ * (two / 11);
1500 ig65 = ig55 + term;
1501 term *= pt_ * (two / 13);
1502 ig75 = ig65 + term;
1503
1504 elfb = ((((((coef75/x + coef65)/x + coef55)/x + coef45)/x + coef35)/x +
1505 coef25)/x + coef15) + x;
1506 res1 = ((((((ig7*coef7/x + ig6*coef6)/x + ig5*coef5)/x + ig4*coef4)/x +
1507 ig3*coef3)/x + ig2*coef2)/x + coef1)*gnm_sqrt(x);
1508 res2 = ((((((ig75*coef75/x + ig65*coef65)/x + ig55*coef55)/x + ig45*coef45)/
1509 x + ig35*coef35)/x + ig25*coef25)/x + coef15)*s2pt;
1510
1511 if (!lower_tail) elfb = -elfb;
1512 f = (res1 + res2) / elfb;
1513
1514 np = pnorm (s2pt, 0.0, 1.0, !lower_tail, log_p);
1515 nd = dnorm (s2pt, 0.0, 1.0, log_p);
1516
1517 #ifdef DEBUG_p
1518 REprintf ("pp*_asymp(): f=%.14" GNM_FORMAT_g " np=%.14" GNM_FORMAT_g " nd=%.14" GNM_FORMAT_g " f*nd=%.14" GNM_FORMAT_g "\n",
1519 f, np, nd, f * nd);
1520 #endif
1521
1522 if (log_p)
1523 return (f >= 0)
1524 ? logspace_add (np, gnm_log (gnm_abs (f)) + nd)
1525 : logspace_sub (np, gnm_log (gnm_abs (f)) + nd);
1526 else
1527 return np + f * nd;
1528 } /* ppois_asymp() */
1529
1530
1531 static gnm_float
pgamma_raw(gnm_float x,gnm_float alph,gboolean lower_tail,gboolean log_p)1532 pgamma_raw (gnm_float x, gnm_float alph, gboolean lower_tail, gboolean log_p)
1533 {
1534 gnm_float res;
1535
1536 #ifdef DEBUG_p
1537 REprintf("pgamma_raw(x=%.14" GNM_FORMAT_g ", alph=%.14" GNM_FORMAT_g ", low=%d, log=%d)\n",
1538 x, alph, lower_tail, log_p);
1539 #endif
1540 if (x < 1) {
1541 res = pgamma_smallx (x, alph, lower_tail, log_p);
1542 } else if (x <= alph - 1 && x < 0.8 * (alph + 50)) {/* incl. large alph */
1543 gnm_float sum = pd_upper_series (x, alph, log_p);/* = x/alph + o(x/alph) */
1544 gnm_float d = dpois_wrap (alph, x, log_p);
1545 #ifdef DEBUG_p
1546 REprintf(" alph `large': sum=pd_upper*()= %.12" GNM_FORMAT_g ", d=dpois_w(*)= %.12" GNM_FORMAT_g " ",
1547 sum, d);
1548 #endif
1549 if (!lower_tail)
1550 res = log_p
1551 ? swap_log_tail (d + sum)
1552 : 1 - d * sum;
1553 else
1554 res = log_p ? sum + d : sum * d;
1555 } else if (alph - 1 < x && alph < 0.8 * (x + 50)) {/* incl. large x */
1556 gnm_float sum;
1557 gnm_float d = dpois_wrap (alph, x, log_p);
1558 #ifdef DEBUG_p
1559 REprintf(" x `large': d=dpois_w(*)= %.14" GNM_FORMAT_g " ", d);
1560 #endif
1561
1562 if (alph < 1) {
1563 if (x * GNM_EPSILON > 1 - alph)
1564 sum = R_D__1;
1565 else {
1566 gnm_float f = pd_lower_cf (alph, x - (alph - 1)) * x / alph;
1567 /* = [alph/(x - alph+1) + o(alph/(x-alph+1))] * x/alph = 1 + o(1) */
1568 sum = log_p ? gnm_log (f) : f;
1569 }
1570 } else {
1571 sum = pd_lower_series (x, alph - 1);/* = (alph-1)/x + o((alph-1)/x) */
1572 sum = log_p ? gnm_log1p (sum) : 1 + sum;
1573 }
1574 #ifdef DEBUG_p
1575 REprintf(", sum= %.14" GNM_FORMAT_g "\n", sum);
1576 #endif
1577 if (!lower_tail)
1578 res = log_p ? sum + d : sum * d;
1579 else
1580 res = log_p
1581 ? swap_log_tail (d + sum)
1582 : 1 - d * sum;
1583 } else {
1584 #ifdef DEBUG_p
1585 REprintf(" using ppois_asymp()\n");
1586 #endif
1587 res = ppois_asymp (alph - 1, x, !lower_tail, log_p);
1588 }
1589
1590 /*
1591 * We lose a fair amount of accuracy to underflow in the cases
1592 * where the final result is very close to DBL_MIN. In those
1593 * cases, simply redo via log space.
1594 */
1595 if (!log_p && res < GNM_MIN / GNM_EPSILON) {
1596 /* with(.Machine, gnm_float.xmin / gnm_float.eps) #|-> GNM_const(1.002084e-292) */
1597 #ifdef DEBUG_p
1598 REprintf(" very small res=%.14" GNM_FORMAT_g "; -> recompute via log\n", res);
1599 #endif
1600 return gnm_exp (pgamma_raw (x, alph, lower_tail, 1));
1601 } else
1602 return res;
1603 }
1604
1605
pgamma(gnm_float x,gnm_float alph,gnm_float scale,gboolean lower_tail,gboolean log_p)1606 gnm_float pgamma(gnm_float x, gnm_float alph, gnm_float scale, gboolean lower_tail, gboolean log_p)
1607 {
1608 #ifdef IEEE_754
1609 if (gnm_isnan(x) || gnm_isnan(alph) || gnm_isnan(scale))
1610 return x + alph + scale;
1611 #endif
1612 if(alph <= 0. || scale <= 0.)
1613 ML_ERR_return_NAN;
1614 x /= scale;
1615 #ifdef IEEE_754
1616 if (gnm_isnan(x)) /* eg. original x = scale = +Inf */
1617 return x;
1618 #endif
1619 if (x <= 0.) /* also for scale=Inf and finite x */
1620 return R_DT_0;
1621
1622 return pgamma_raw (x, alph, lower_tail, log_p);
1623 }
1624 /* From: terra@gnome.org (Morten Welinder)
1625 * To: R-bugs@biostat.ku.dk
1626 * Cc: maechler@stat.math.ethz.ch
1627 * Subject: Re: [Rd] pgamma discontinuity (PR#7307)
1628 * Date: Tue, 11 Jan 2005 13:57:26 -0500 (EST)
1629
1630 * this version of pgamma appears to be quite good and certainly a vast
1631 * improvement over current R code. (I last looked at 2.0.1) Apart from
1632 * type naming, this is what I have been using for Gnumeric 1.4.1.
1633
1634 * This could be included into R as-is, but you might want to benefit from
1635 * making gnm_logcf, log1pmx, lgamma1p, and possibly logspace_add/logspace_sub
1636 * available to other parts of R.
1637
1638 * MM: I've not (yet?) taken gnm_logcf(), but the other four
1639 */
1640
1641
1642 #else
1643 /* R_USE_OLD_PGAMMA */
1644 /*
1645 * Copyright (C) 1998 Ross Ihaka
1646 * Copyright (C) 1999-2000 The R Development Core Team
1647 * Copyright (C) 2003-2004 The R Foundation
1648 * based on AS 239 (C) 1988 Royal Statistical Society
1649 *
1650 * ................
1651 *
1652 * NOTES
1653 *
1654 * This function is an adaptation of Algorithm 239 from the
1655 * Applied Statistics Series. The algorithm is faster than
1656 * those by W. Fullerton in the FNLIB library and also the
1657 * TOMS 542 alorithm of W. Gautschi. It provides comparable
1658 * accuracy to those algorithms and is considerably simpler.
1659 *
1660 * REFERENCES
1661 *
1662 * Algorithm AS 239, Incomplete Gamma Function
1663 * Applied Statistics 37, 1988.
1664 */
1665
pgamma(gnm_float x,gnm_float alph,gnm_float scale,gboolean lower_tail,gboolean log_p)1666 gnm_float pgamma(gnm_float x, gnm_float alph, gnm_float scale, gboolean lower_tail, gboolean log_p)
1667 {
1668 const gnm_float
1669 xbig = 1.0e+8,
1670 xlarge = 1.0e+37,
1671
1672 /* normal approx. for alph > alphlimit */
1673 alphlimit = 1e5;/* was 1000. till R.1.8.x */
1674
1675 gnm_float pn1, pn2, pn3, pn4, pn5, pn6, arg, a, b, c, an, osum, sum;
1676 long n;
1677 int pearson;
1678
1679 /* check that we have valid values for x and alph */
1680
1681 #ifdef IEEE_754
1682 if (gnm_isnan(x) || gnm_isnan(alph) || gnm_isnan(scale))
1683 return x + alph + scale;
1684 #endif
1685 #ifdef DEBUG_p
1686 REprintf("pgamma(x=%4" GNM_FORMAT_g ", alph=%4" GNM_FORMAT_g ", scale=%4" GNM_FORMAT_g "): ",x,alph,scale);
1687 #endif
1688 if(alph <= 0. || scale <= 0.)
1689 ML_ERR_return_NAN;
1690
1691 x /= scale;
1692 #ifdef DEBUG_p
1693 REprintf("-> x=%4" GNM_FORMAT_g "; ",x);
1694 #endif
1695 #ifdef IEEE_754
1696 if (gnm_isnan(x)) /* eg. original x = scale = Inf */
1697 return x;
1698 #endif
1699 if (x <= 0.)
1700 return R_DT_0;
1701
1702 #define USE_PNORM \
1703 pn1 = gnm_sqrt(alph) * 3. * (gnm_pow(x/alph, GNM_const(1.)/3.) + 1. / (9. * alph) - 1.); \
1704 return pnorm(pn1, 0., 1., lower_tail, log_p);
1705
1706 if (alph > alphlimit) { /* use a normal approximation */
1707 USE_PNORM;
1708 }
1709
1710 if (x > xbig * alph) {
1711 if (x > GNM_MAX * alph)
1712 /* if x is extremely large __compared to alph__ then return 1 */
1713 return R_DT_1;
1714 else { /* this only "helps" when log_p = TRUE */
1715 USE_PNORM;
1716 }
1717 }
1718
1719 if (x <= 1. || x < alph) {
1720
1721 pearson = 1;/* use pearson's series expansion. */
1722
1723 arg = alph * gnm_log(x) - x - gnm_lgamma(alph + 1.);
1724 #ifdef DEBUG_p
1725 REprintf("Pearson arg=%" GNM_FORMAT_g " ", arg);
1726 #endif
1727 c = 1.;
1728 sum = 1.;
1729 a = alph;
1730 do {
1731 a += 1.;
1732 c *= x / a;
1733 sum += c;
1734 } while (c > GNM_EPSILON * sum);
1735 }
1736 else { /* x >= max( 1, alph) */
1737
1738 pearson = 0;/* use a continued fraction expansion */
1739
1740 arg = alph * gnm_log(x) - x - gnm_lgamma(alph);
1741 #ifdef DEBUG_p
1742 REprintf("Cont.Fract. arg=%" GNM_FORMAT_g " ", arg);
1743 #endif
1744 a = 1. - alph;
1745 b = a + x + 1.;
1746 pn1 = 1.;
1747 pn2 = x;
1748 pn3 = x + 1.;
1749 pn4 = x * b;
1750 sum = pn3 / pn4;
1751 for (n = 1; ; n++) {
1752 a += 1.;/* = n+1 -alph */
1753 b += 2.;/* = 2(n+1)-alph+x */
1754 an = a * n;
1755 pn5 = b * pn3 - an * pn1;
1756 pn6 = b * pn4 - an * pn2;
1757 if (gnm_abs(pn6) > 0.) {
1758 osum = sum;
1759 sum = pn5 / pn6;
1760 if (gnm_abs(osum - sum) <= GNM_EPSILON * fmin2(1., sum))
1761 break;
1762 }
1763 pn1 = pn3;
1764 pn2 = pn4;
1765 pn3 = pn5;
1766 pn4 = pn6;
1767 if (gnm_abs(pn5) >= xlarge) {
1768 /* re-scale the terms in continued fraction if they are large */
1769 #ifdef DEBUG_p
1770 REprintf(" [r] ");
1771 #endif
1772 pn1 /= xlarge;
1773 pn2 /= xlarge;
1774 pn3 /= xlarge;
1775 pn4 /= xlarge;
1776 }
1777 }
1778 }
1779
1780 arg += gnm_log(sum);
1781
1782 lower_tail = (lower_tail == pearson);
1783
1784 if (log_p && lower_tail)
1785 return(arg);
1786 /* else */
1787 /* sum = gnm_exp(arg); and return if(lower_tail) sum else 1-sum : */
1788 return (lower_tail) ? gnm_exp(arg) : (log_p ? swap_log_tail(arg) : -gnm_expm1(arg));
1789 }
1790
1791 #endif
1792 /* R_USE_OLD_PGAMMA */
1793 /* Cleaning up done by tools/import-R: */
1794 #undef USE_PNORM
1795 #undef max_it
1796
1797 /* ------------------------------------------------------------------------ */
1798 /* Imported src/nmath/dt.c from R. */
1799 /*
1800 * AUTHOR
1801 * Catherine Loader, catherine@research.bell-labs.com.
1802 * October 23, 2000.
1803 *
1804 * Merge in to R:
1805 * Copyright (C) 2000, The R Core Development Team
1806 *
1807 * This program is free software; you can redistribute it and/or modify
1808 * it under the terms of the GNU General Public License as published by
1809 * the Free Software Foundation; either version 2 of the License, or
1810 * (at your option) any later version.
1811 *
1812 * This program is distributed in the hope that it will be useful,
1813 * but WITHOUT ANY WARRANTY; without even the implied warranty of
1814 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1815 * GNU General Public License for more details.
1816 *
1817 * You should have received a copy of the GNU General Public License
1818 * along with this program; if not, write to the Free Software
1819 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
1820 * USA.
1821 *
1822 *
1823 * DESCRIPTION
1824 *
1825 * The t density is evaluated as
1826 * sqrt(n/2) / ((n+1)/2) * Gamma((n+3)/2) / Gamma((n+2)/2).
1827 * * (1+x^2/n)^(-n/2)
1828 * / sqrt( 2 pi (1+x^2/n) )
1829 *
1830 * This form leads to a stable computation for all
1831 * values of n, including n -> 0 and n -> infinity.
1832 */
1833
1834
dt(gnm_float x,gnm_float n,gboolean give_log)1835 gnm_float dt(gnm_float x, gnm_float n, gboolean give_log)
1836 {
1837 gnm_float t, u;
1838 #ifdef IEEE_754
1839 if (gnm_isnan(x) || gnm_isnan(n))
1840 return x + n;
1841 #endif
1842
1843 if (n <= 0) ML_ERR_return_NAN;
1844 if(!gnm_finite(x))
1845 return R_D__0;
1846 if(!gnm_finite(n))
1847 return dnorm(x, 0., 1., give_log);
1848
1849 t = -bd0(n/2.,(n+1)/2.) + stirlerr((n+1)/2.) - stirlerr(n/2.);
1850 if ( x*x > 0.2*n )
1851 u = gnm_log1p (x*x/n ) * n/2;
1852 else
1853 u = -bd0(n/2.,(n+x*x)/2.) + x*x/2.;
1854
1855 return R_D_fexp(M_2PIgnum*(1+x*x/n), t-u);
1856 }
1857
1858 /* ------------------------------------------------------------------------ */
1859 /* Imported src/nmath/pt.c from R. */
1860 /*
1861 * R : A Computer Language for Statistical Data Analysis
1862 * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka
1863 * Copyright (C) 2000 The R Development Core Team
1864 *
1865 * This program is free software; you can redistribute it and/or modify
1866 * it under the terms of the GNU General Public License as published by
1867 * the Free Software Foundation; either version 2 of the License, or
1868 * (at your option) any later version.
1869 *
1870 * This program is distributed in the hope that it will be useful,
1871 * but WITHOUT ANY WARRANTY; without even the implied warranty of
1872 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1873 * GNU General Public License for more details.
1874 *
1875 * You should have received a copy of the GNU General Public License
1876 * along with this program; if not, write to the Free Software
1877 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
1878 * USA.
1879 */
1880
1881
pt(gnm_float x,gnm_float n,gboolean lower_tail,gboolean log_p)1882 gnm_float pt(gnm_float x, gnm_float n, gboolean lower_tail, gboolean log_p)
1883 {
1884 /* return P[ T <= x ] where
1885 * T ~ t_{n} (t distrib. with n degrees of freedom).
1886
1887 * --> ./pnt.c for NON-central
1888 */
1889 gnm_float val;
1890 #ifdef IEEE_754
1891 if (gnm_isnan(x) || gnm_isnan(n))
1892 return x + n;
1893 #endif
1894 if (n <= 0.0) ML_ERR_return_NAN;
1895
1896 if(!gnm_finite(x))
1897 return (x < 0) ? R_DT_0 : R_DT_1;
1898 if(!gnm_finite(n))
1899 return pnorm(x, 0.0, 1.0, lower_tail, log_p);
1900 if (0 && n > 4e5) { /*-- Fixme(?): test should depend on `n' AND `x' ! */
1901 /* Approx. from Abramowitz & Stegun 26.7.8 (p.949) */
1902 val = 1./(4.*n);
1903 return pnorm(x*(1. - val)/gnm_sqrt(1. + x*x*2.*val), 0.0, 1.0,
1904 lower_tail, log_p);
1905 }
1906
1907 val = (n > x * x)
1908 ? pbeta (x * x / (n + x * x), 0.5, n / 2, /*lower_tail*/0, log_p)
1909 : pbeta (n / (n + x * x), n / 2.0, 0.5, /*lower_tail*/1, log_p);
1910
1911 /* Use "1 - v" if lower_tail and x > 0 (but not both):*/
1912 if(x <= 0.)
1913 lower_tail = !lower_tail;
1914
1915 if(log_p) {
1916 if(lower_tail) return gnm_log1p(-0.5*gnm_exp(val));
1917 else return val - M_LN2gnum; /* = gnm_log(.5* pbeta(....)) */
1918 }
1919 else {
1920 val /= 2.;
1921 return R_D_Cval(val);
1922 }
1923 }
1924
1925 /* ------------------------------------------------------------------------ */
1926 /* Imported src/nmath/qt.c from R. */
1927 /*
1928 * Mathlib : A C Library of Special Functions
1929 * Copyright (C) 1998 Ross Ihaka
1930 * Copyright (C) 2000-2002 The R Development Core Team
1931 * Copyright (C) 2003 The R Foundation
1932 *
1933 * This program is free software; you can redistribute it and/or modify
1934 * it under the terms of the GNU General Public License as published by
1935 * the Free Software Foundation; either version 2 of the License, or
1936 * (at your option) any later version.
1937 *
1938 * This program is distributed in the hope that it will be useful,
1939 * but WITHOUT ANY WARRANTY; without even the implied warranty of
1940 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1941 * GNU General Public License for more details.
1942 *
1943 * You should have received a copy of the GNU General Public License
1944 * along with this program; if not, write to the Free Software
1945 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
1946 * USA.
1947 *
1948 * DESCRIPTION
1949 *
1950 * The "Student" t distribution quantile function.
1951 *
1952 * NOTES
1953 *
1954 * This is a C translation of the Fortran routine given in:
1955 * Hill, G.W (1970) "Algorithm 396: Student's t-quantiles"
1956 * CACM 13(10), 619-620.
1957 *
1958 * ADDITIONS:
1959 * - lower_tail, log_p
1960 * - using expm1() : takes care of Lozy (1979) "Remark on Algo.", TOMS
1961 * - Apply 2-term Taylor expansion as in
1962 * Hill, G.W (1981) "Remark on Algo.396", ACM TOMS 7, 250-1
1963 * - Improve the formula decision for 1 < df < 2
1964 */
1965
1966
qt(gnm_float p,gnm_float ndf,gboolean lower_tail,gboolean log_p)1967 gnm_float qt(gnm_float p, gnm_float ndf, gboolean lower_tail, gboolean log_p)
1968 {
1969 const gnm_float eps = 1.e-12;
1970
1971 gnm_float a, b, c, d, p_, P, q, x, y;
1972 gboolean neg;
1973
1974 #ifdef IEEE_754
1975 if (gnm_isnan(p) || gnm_isnan(ndf))
1976 return p + ndf;
1977 #endif
1978 if (p == R_DT_0) return gnm_ninf;
1979 if (p == R_DT_1) return gnm_pinf;
1980 R_Q_P01_check(p);
1981
1982 if (ndf < 1) /* FIXME: not yet treated here */
1983 ML_ERR_return_NAN;
1984
1985 /* FIXME: This test should depend on ndf AND p !!
1986 * ----- and in fact should be replaced by
1987 * something like Abramowitz & Stegun 26.7.5 (p.949)
1988 */
1989 if (ndf > 1e20) return qnorm(p, 0., 1., lower_tail, log_p);
1990
1991 p_ = R_D_qIv(p); /* note: gnm_exp(p) may underflow to 0; fix later */
1992
1993 if((lower_tail && p_ > 0.5) || (!lower_tail && p_ < 0.5)) {
1994 neg = FALSE; P = 2 * R_D_Cval(p_);
1995 } else {
1996 neg = TRUE; P = 2 * R_D_Lval(p_);
1997 } /* 0 <= P <= 1 ; P = 2*min(p_, 1 - p_) in all cases */
1998
1999 if (gnm_abs(ndf - 2) < eps) { /* df ~= 2 */
2000 if(P > 0)
2001 q = gnm_sqrt(2 / (P * (2 - P)) - 2);
2002 else { /* P = 0, but maybe = gnm_exp(p) ! */
2003 if(log_p) q = M_SQRT2gnum * gnm_exp(- .5 * R_D_Lval(p));
2004 else q = gnm_pinf;
2005 }
2006 }
2007 else if (ndf < 1 + eps) { /* df ~= 1 (df < 1 excluded above): Cauchy */
2008 if(P > 0)
2009 q = gnm_cotpi(P / 2);
2010
2011 else { /* P = 0, but maybe p_ = gnm_exp(p) ! */
2012 if(log_p) q = M_1_PI * gnm_exp(-R_D_Lval(p));/* cot(e) ~ 1/e */
2013 else q = gnm_pinf;
2014 }
2015 }
2016 else { /*-- usual case; including, e.g., df = 1.1 */
2017 a = 1 / (ndf - 0.5);
2018 b = 48 / (a * a);
2019 c = ((20700 * a / b - 98) * a - 16) * a + 96.36;
2020 d = ((94.5 / (b + c) - 3) / b + 1) * gnm_sqrt(a * M_PI_2gnum) * ndf;
2021 if(P > 0 || !log_p)
2022 y = gnm_pow(d * P, 2 / ndf);
2023 else /* P = 0 && log_p; P = 2*gnm_exp(p) */
2024 y = gnm_exp(2 / ndf * (gnm_log(d) + M_LN2gnum + R_D_Lval(p)));
2025
2026 if ((ndf < 2.1 && P > 0.5) || y > 0.05 + a) { /* P > P0(df) */
2027 /* Asymptotic inverse expansion about normal */
2028 if(P > 0 || !log_p)
2029 x = qnorm(0.5 * P, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE);
2030 else /* P = 0 && log_p; P = 2*gnm_exp(p') */
2031 x = qnorm( p, 0., 1., lower_tail, /*log_p*/TRUE);
2032
2033 y = x * x;
2034 if (ndf < 5)
2035 c += 0.3 * (ndf - 4.5) * (x + 0.6);
2036 c = (((0.05 * d * x - 5) * x - 7) * x - 2) * x + b + c;
2037 y = (((((0.4 * y + 6.3) * y + 36) * y + 94.5) / c
2038 - y - 3) / b + 1) * x;
2039 y = gnm_expm1(a * y * y);
2040 } else {
2041 y = ((1 / (((ndf + 6) / (ndf * y) - 0.089 * d - 0.822)
2042 * (ndf + 2) * 3) + 0.5 / (ndf + 4))
2043 * y - 1) * (ndf + 1) / (ndf + 2) + 1 / y;
2044 }
2045 q = gnm_sqrt(ndf * y);
2046
2047 /* Now apply 2-term Taylor expansion improvement (1-term = Newton):
2048 * as by Hill (1981) [ref.above] */
2049
2050 /* FIXME: This is can be far from optimal when log_p = TRUE !
2051 * and probably also improvable when lower_tail = FALSE */
2052 x = (pt(q, ndf, /*lower_tail = */FALSE, /*log_p = */FALSE) - P/2) /
2053 dt(q, ndf, /* give_log = */FALSE);
2054 /* Newton (=Taylor 1 term):
2055 * q += x;
2056 * Taylor 2-term : */
2057 q += x * (1. + x * q * (ndf + 1) / (2 * (q * q + ndf)));
2058 }
2059 if(neg) q = -q;
2060 return q;
2061 }
2062
2063 /* ------------------------------------------------------------------------ */
2064 /* Imported src/nmath/pchisq.c from R. */
2065 /*
2066 * Mathlib : A C Library of Special Functions
2067 * Copyright (C) 1998 Ross Ihaka
2068 * Copyright (C) 2000 The R Development Core Team
2069 *
2070 * This program is free software; you can redistribute it and/or modify
2071 * it under the terms of the GNU General Public License as published by
2072 * the Free Software Foundation; either version 2 of the License, or
2073 * (at your option) any later version.
2074 *
2075 * This program is distributed in the hope that it will be useful, but
2076 * WITHOUT ANY WARRANTY; without even the implied warranty of
2077 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2078 * General Public License for more details.
2079 *
2080 * You should have received a copy of the GNU General Public License
2081 * along with this program; if not, write to the Free Software
2082 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2083 * USA.
2084 *
2085 * DESCRIPTION
2086 *
2087 * The distribution function of the chi-squared distribution.
2088 */
2089
2090
pchisq(gnm_float x,gnm_float df,gboolean lower_tail,gboolean log_p)2091 gnm_float pchisq(gnm_float x, gnm_float df, gboolean lower_tail, gboolean log_p)
2092 {
2093 return pgamma(x, df/2., 2., lower_tail, log_p);
2094 }
2095
2096 /* ------------------------------------------------------------------------ */
2097 /* Imported src/nmath/qchisq.c from R. */
2098 /*
2099 * Mathlib : A C Library of Special Functions
2100 * Copyright (C) 1998 Ross Ihaka
2101 * Copyright (C) 2000 The R Development Core Team
2102 *
2103 * This program is free software; you can redistribute it and/or modify
2104 * it under the terms of the GNU General Public License as published by
2105 * the Free Software Foundation; either version 2 of the License, or
2106 * (at your option) any later version.
2107 *
2108 * This program is distributed in the hope that it will be useful,
2109 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2110 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2111 * GNU General Public License for more details.
2112 *
2113 * You should have received a copy of the GNU General Public License
2114 * along with this program; if not, write to the Free Software
2115 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2116 * USA.
2117 *
2118 * DESCRIPTION
2119 *
2120 * The quantile function of the chi-squared distribution.
2121 */
2122
2123
qchisq(gnm_float p,gnm_float df,gboolean lower_tail,gboolean log_p)2124 gnm_float qchisq(gnm_float p, gnm_float df, gboolean lower_tail, gboolean log_p)
2125 {
2126 return qgamma(p, 0.5 * df, 2.0, lower_tail, log_p);
2127 }
2128
2129 /* ------------------------------------------------------------------------ */
2130 /* Imported src/nmath/dweibull.c from R. */
2131 /*
2132 * Mathlib : A C Library of Special Functions
2133 * Copyright (C) 1998 Ross Ihaka
2134 * Copyright (C) 2000 The R Development Core Team
2135 *
2136 * This program is free software; you can redistribute it and/or modify
2137 * it under the terms of the GNU General Public License as published by
2138 * the Free Software Foundation; either version 2 of the License, or
2139 * (at your option) any later version.
2140 *
2141 * This program is distributed in the hope that it will be useful,
2142 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2143 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2144 * GNU General Public License for more details.
2145 *
2146 * You should have received a copy of the GNU General Public License
2147 * along with this program; if not, write to the Free Software
2148 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2149 * USA.
2150 *
2151 * DESCRIPTION
2152 *
2153 * The density function of the Weibull distribution.
2154 */
2155
2156
dweibull(gnm_float x,gnm_float shape,gnm_float scale,gboolean give_log)2157 gnm_float dweibull(gnm_float x, gnm_float shape, gnm_float scale, gboolean give_log)
2158 {
2159 gnm_float tmp1, tmp2;
2160 #ifdef IEEE_754
2161 if (gnm_isnan(x) || gnm_isnan(shape) || gnm_isnan(scale))
2162 return x + shape + scale;
2163 #endif
2164 if (shape <= 0 || scale <= 0) ML_ERR_return_NAN;
2165
2166 if (x < 0) return R_D__0;
2167 if (!gnm_finite(x)) return R_D__0;
2168 tmp1 = gnm_pow(x / scale, shape - 1);
2169 tmp2 = tmp1 * (x / scale);
2170 return give_log ?
2171 -tmp2 + gnm_log(shape * tmp1 / scale) :
2172 shape * tmp1 * gnm_exp(-tmp2) / scale;
2173 }
2174
2175 /* ------------------------------------------------------------------------ */
2176 /* Imported src/nmath/pweibull.c from R. */
2177 /*
2178 * Mathlib : A C Library of Special Functions
2179 * Copyright (C) 1998 Ross Ihaka
2180 * Copyright (C) 2000-2002 The R Development Core Team
2181 *
2182 * This program is free software; you can redistribute it and/or modify
2183 * it under the terms of the GNU General Public License as published by
2184 * the Free Software Foundation; either version 2 of the License, or
2185 * (at your option) any later version.
2186 *
2187 * This program is distributed in the hope that it will be useful,
2188 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2189 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2190 * GNU General Public License for more details.
2191 *
2192 * You should have received a copy of the GNU General Public License
2193 * along with this program; if not, write to the Free Software
2194 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2195 * USA.
2196 *
2197 * DESCRIPTION
2198 *
2199 * The distribution function of the Weibull distribution.
2200 */
2201
2202
pweibull(gnm_float x,gnm_float shape,gnm_float scale,gboolean lower_tail,gboolean log_p)2203 gnm_float pweibull(gnm_float x, gnm_float shape, gnm_float scale, gboolean lower_tail, gboolean log_p)
2204 {
2205 #ifdef IEEE_754
2206 if (gnm_isnan(x) || gnm_isnan(shape) || gnm_isnan(scale))
2207 return x + shape + scale;
2208 #endif
2209 if(shape <= 0 || scale <= 0) ML_ERR_return_NAN;
2210
2211 if (x <= 0)
2212 return R_DT_0;
2213 x = -gnm_pow(x / scale, shape);
2214 if (lower_tail)
2215 return (log_p
2216 /* gnm_log(1 - gnm_exp(x)) for x < 0 : */
2217 ? (x > -M_LN2gnum ? gnm_log(-gnm_expm1(x)) : gnm_log1p(-gnm_exp(x)))
2218 : -gnm_expm1(x));
2219 /* else: !lower_tail */
2220 return R_D_exp(x);
2221 }
2222
2223 /* ------------------------------------------------------------------------ */
2224 /* Imported src/nmath/pbinom.c from R. */
2225 /*
2226 * Mathlib : A C Library of Special Functions
2227 * Copyright (C) 1998 Ross Ihaka
2228 * Copyright (C) 2000, 2002 The R Development Core Team
2229 * Copyright (C) 2004 The R Foundation
2230 *
2231 * This program is free software; you can redistribute it and/or modify
2232 * it under the terms of the GNU General Public License as published by
2233 * the Free Software Foundation; either version 2 of the License, or
2234 * (at your option) any later version.
2235 *
2236 * This program is distributed in the hope that it will be useful,
2237 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2238 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2239 * GNU General Public License for more details.
2240 *
2241 * You should have received a copy of the GNU General Public License
2242 * along with this program; if not, write to the Free Software
2243 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2244 * USA.
2245 *
2246 * DESCRIPTION
2247 *
2248 * The distribution function of the binomial distribution.
2249 */
2250
pbinom(gnm_float x,gnm_float n,gnm_float p,gboolean lower_tail,gboolean log_p)2251 gnm_float pbinom(gnm_float x, gnm_float n, gnm_float p, gboolean lower_tail, gboolean log_p)
2252 {
2253 #ifdef IEEE_754
2254 if (gnm_isnan(x) || gnm_isnan(n) || gnm_isnan(p))
2255 return x + n + p;
2256 if (!gnm_finite(n) || !gnm_finite(p)) ML_ERR_return_NAN;
2257
2258 #endif
2259 if(R_D_nonint(n)) ML_ERR_return_NAN;
2260 n = R_D_forceint(n);
2261 if(n <= 0 || p < 0 || p > 1) ML_ERR_return_NAN;
2262
2263 x = gnm_fake_floor(x);
2264 if (x < 0.0) return R_DT_0;
2265 if (n <= x) return R_DT_1;
2266 return pbeta(p, x + 1, n - x, !lower_tail, log_p);
2267 }
2268
2269 /* ------------------------------------------------------------------------ */
2270 /* Imported src/nmath/dbinom.c from R. */
2271 /*
2272 * AUTHOR
2273 * Catherine Loader, catherine@research.bell-labs.com.
2274 * October 23, 2000.
2275 *
2276 * Merge in to R:
2277 * Copyright (C) 2000, The R Core Development Team
2278 *
2279 * This program is free software; you can redistribute it and/or modify
2280 * it under the terms of the GNU General Public License as published by
2281 * the Free Software Foundation; either version 2 of the License, or
2282 * (at your option) any later version.
2283 *
2284 * This program is distributed in the hope that it will be useful,
2285 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2286 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2287 * GNU General Public License for more details.
2288 *
2289 * You should have received a copy of the GNU General Public License
2290 * along with this program; if not, write to the Free Software
2291 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2292 * USA.
2293 *
2294 *
2295 * DESCRIPTION
2296 *
2297 * To compute the binomial probability, call dbinom(x,n,p).
2298 * This checks for argument validity, and calls dbinom_raw().
2299 *
2300 * dbinom_raw() does the actual computation; note this is called by
2301 * other functions in addition to dbinom()).
2302 * (1) dbinom_raw() has both p and q arguments, when one may be represented
2303 * more accurately than the other (in particular, in df()).
2304 * (2) dbinom_raw() does NOT check that inputs x and n are integers. This
2305 * should be done in the calling function, where necessary.
2306 * (3) Also does not check for 0 <= p <= 1 and 0 <= q <= 1 or NaN's.
2307 * Do this in the calling function.
2308 */
2309
2310
dbinom_raw(gnm_float x,gnm_float n,gnm_float p,gnm_float q,gboolean give_log)2311 static gnm_float dbinom_raw(gnm_float x, gnm_float n, gnm_float p, gnm_float q, gboolean give_log)
2312 {
2313 gnm_float f2, xh, xl, yh, yl;
2314
2315 /*
2316 * We (ought to) have p+q = 1 with the assumption that the smaller is
2317 * more accurate that 1-other, i.e., that the bigger may already have
2318 * suffered some cancellation.
2319 */
2320
2321 if (p == 0) return((x == 0) ? R_D__1 : R_D__0);
2322 if (q == 0) return((x == n) ? R_D__1 : R_D__0);
2323
2324 if (x == 0) {
2325 /* Switch p and q to reduce code duplication. */
2326 gnm_float t = p;
2327 p = q;
2328 q = t;
2329 x = n;
2330 }
2331
2332 if (x == n) {
2333 /* Probability is p^n. */
2334 if (p <= 0.5)
2335 return give_log ? n * gnm_log (p) : gnm_pow (p, n);
2336 else
2337 return give_log ? n * gnm_log1p (-q) : pow1p (-q, n);
2338 }
2339 if (x < 0 || x > n) return( R_D__0 );
2340
2341 ebd0 (x, n*p, &xh, &xl);
2342 PAIR_ADD(stirlerr(x), xh, xl);
2343
2344 ebd0 (n-x, n*q, &yh, &yl);
2345 PAIR_ADD(stirlerr(n-x), yh, yl);
2346
2347 PAIR_ADD(yl, xh, xl);
2348 PAIR_ADD(yh, xh, xl);
2349 PAIR_ADD(-stirlerr(n), xh, xl);
2350
2351 f2 = (M_2PIgnum*x*(n-x))/n;
2352
2353 return give_log
2354 ? -xl - xh - 0.5 * gnm_log (f2)
2355 : gnm_exp (-xl) * gnm_exp (-xh) / gnm_sqrt (f2);
2356 }
2357
dbinom(gnm_float x,gnm_float n,gnm_float p,gboolean give_log)2358 gnm_float dbinom(gnm_float x, gnm_float n, gnm_float p, gboolean give_log)
2359 {
2360 #ifdef IEEE_754
2361 /* NaNs propagated correctly */
2362 if (gnm_isnan(x) || gnm_isnan(n) || gnm_isnan(p)) return x + n + p;
2363 #endif
2364
2365 if (p < 0 || p > 1 || R_D_negInonint(n))
2366 ML_ERR_return_NAN;
2367 R_D_nonint_check(x);
2368
2369 n = R_D_forceint(n);
2370 x = R_D_forceint(x);
2371
2372 return dbinom_raw(x,n,p,1-p,give_log);
2373 }
2374
2375 /* ------------------------------------------------------------------------ */
2376 /* Imported src/nmath/qbinom.c from R. */
2377 /*
2378 * Mathlib : A C Library of Special Functions
2379 * Copyright (C) 1998 Ross Ihaka
2380 * Copyright (C) 2000, 2002 The R Development Core Team
2381 * Copyright (C) 2003--2004 The R Foundation
2382 *
2383 * This program is free software; you can redistribute it and/or modify
2384 * it under the terms of the GNU General Public License as published by
2385 * the Free Software Foundation; either version 2 of the License, or
2386 * (at your option) any later version.
2387 *
2388 * This program is distributed in the hope that it will be useful,
2389 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2390 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2391 * GNU General Public License for more details.
2392 *
2393 * You should have received a copy of the GNU General Public License
2394 * along with this program; if not, write to the Free Software
2395 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2396 * USA.
2397 *
2398 * DESCRIPTION
2399 *
2400 * The quantile function of the binomial distribution.
2401 *
2402 * METHOD
2403 *
2404 * Uses the Cornish-Fisher Expansion to include a skewness
2405 * correction to a normal approximation. This gives an
2406 * initial value which never seems to be off by more than
2407 * 1 or 2. A search is then conducted of values close to
2408 * this initial start point.
2409 */
2410
2411
qbinom(gnm_float p,gnm_float n,gnm_float pr,gboolean lower_tail,gboolean log_p)2412 gnm_float qbinom(gnm_float p, gnm_float n, gnm_float pr, gboolean lower_tail, gboolean log_p)
2413 {
2414 gnm_float q, mu, sigma, gamma, z, y;
2415
2416 #ifdef IEEE_754
2417 if (gnm_isnan(p) || gnm_isnan(n) || gnm_isnan(pr))
2418 return p + n + pr;
2419 #endif
2420 if(!gnm_finite(p) || !gnm_finite(n) || !gnm_finite(pr))
2421 ML_ERR_return_NAN;
2422 R_Q_P01_check(p);
2423
2424 if(n != gnm_floor(n + 0.5)) ML_ERR_return_NAN;
2425 if (pr < 0 || pr > 1 || n < 0)
2426 ML_ERR_return_NAN;
2427
2428 if (pr == 0. || n == 0) return 0.;
2429 if (p == R_DT_0) return 0.;
2430 if (p == R_DT_1) return n;
2431
2432 q = 1 - pr;
2433 if(q == 0.) return n; /* covers the full range of the distribution */
2434 mu = n * pr;
2435 sigma = gnm_sqrt(n * pr * q);
2436 gamma = (q - pr) / sigma;
2437
2438 #ifdef DEBUG_qbinom
2439 REprintf("qbinom(p=%7" GNM_FORMAT_g ", n=%" GNM_FORMAT_g ", pr=%7" GNM_FORMAT_g ", l.t.=%d, log=%d): sigm=%" GNM_FORMAT_g ", gam=%" GNM_FORMAT_g "\n",
2440 p,n,pr, lower_tail, log_p, sigma, gamma);
2441 #endif
2442 /* Note : "same" code in qpois.c, qbinom.c, qnbinom.c --
2443 * FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */
2444 if(!lower_tail || log_p) {
2445 p = R_DT_qIv(p); /* need check again (cancellation!): */
2446 if (p == 0.) return 0.;
2447 if (p == 1.) return n;
2448 }
2449 /* temporary hack --- FIXME --- */
2450 if (p + 1.01*GNM_EPSILON >= 1.) return n;
2451
2452 /* y := approx.value (Cornish-Fisher expansion) : */
2453 z = qnorm(p, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE);
2454 y = gnm_floor(mu + sigma * (z + gamma * (z*z - 1) / 6) + 0.5);
2455 if(y > n) /* way off */ y = n;
2456
2457 #ifdef DEBUG_qbinom
2458 REprintf(" new (p,1-p)=(%7" GNM_FORMAT_g ",%7" GNM_FORMAT_g "), z=qnorm(..)=%7" GNM_FORMAT_g ", y=%5" GNM_FORMAT_g "\n", p, 1-p, z, y);
2459 #endif
2460 z = pbinom(y, n, pr, /*lower_tail*/TRUE, /*log_p*/FALSE);
2461
2462 /* fuzz to ensure left continuity: */
2463 p *= 1 - 64*GNM_EPSILON;
2464
2465 /*-- Fixme, here y can be way off --
2466 should use interval search instead of primitive stepping down or up */
2467
2468 #ifdef maybe_future
2469 if((lower_tail && z >= p) || (!lower_tail && z <= p)) {
2470 #else
2471 if(z >= p) {
2472 #endif
2473 /* search to the left */
2474 #ifdef DEBUG_qbinom
2475 REprintf("\tnew z=%7" GNM_FORMAT_g " >= p = %7" GNM_FORMAT_g " --> search to left (y--) ..\n", z,p);
2476 #endif
2477 for(;;) {
2478 if(y == 0 ||
2479 (z = pbinom(y - 1, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) < p)
2480 return y;
2481 y = y - 1;
2482 }
2483 }
2484 else { /* search to the right */
2485 #ifdef DEBUG_qbinom
2486 REprintf("\tnew z=%7" GNM_FORMAT_g " < p = %7" GNM_FORMAT_g " --> search to right (y++) ..\n", z,p);
2487 #endif
2488 for(;;) {
2489 y = y + 1;
2490 if(y == n ||
2491 (z = pbinom(y, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) >= p)
2492 return y;
2493 }
2494 }
2495 }
2496
2497 #if 0
2498 }
2499 #endif
2500
2501 /* ------------------------------------------------------------------------ */
2502 /* Imported src/nmath/dnbinom.c from R. */
2503 /*
2504 * AUTHOR
2505 * Catherine Loader, catherine@research.bell-labs.com.
2506 * October 23, 2000 and Feb, 2001.
2507 *
2508 * Merge in to R:
2509 * Copyright (C) 2000--2001, The R Core Development Team
2510 *
2511 * This program is free software; you can redistribute it and/or modify
2512 * it under the terms of the GNU General Public License as published by
2513 * the Free Software Foundation; either version 2 of the License, or
2514 * (at your option) any later version.
2515 *
2516 * This program is distributed in the hope that it will be useful,
2517 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2518 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2519 * GNU General Public License for more details.
2520 *
2521 * You should have received a copy of the GNU General Public License
2522 * along with this program; if not, write to the Free Software
2523 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2524 * USA.
2525 *
2526 *
2527 * DESCRIPTION
2528 *
2529 * Computes the negative binomial distribution. For integer n,
2530 * this is probability of x failures before the nth success in a
2531 * sequence of Bernoulli trials. We do not enforce integer n, since
2532 * the distribution is well defined for non-integers,
2533 * and this can be useful for e.g. overdispersed discrete survival times.
2534 */
2535
2536
dnbinom(gnm_float x,gnm_float n,gnm_float p,gboolean give_log)2537 gnm_float dnbinom(gnm_float x, gnm_float n, gnm_float p, gboolean give_log)
2538 {
2539 gnm_float prob;
2540
2541 #ifdef IEEE_754
2542 if (gnm_isnan(x) || gnm_isnan(n) || gnm_isnan(p))
2543 return x + n + p;
2544 #endif
2545
2546 if (p < 0 || p > 1 || n <= 0) ML_ERR_return_NAN;
2547 R_D_nonint_check(x);
2548 if (x < 0 || !gnm_finite(x)) return R_D__0;
2549 x = R_D_forceint(x);
2550
2551 prob = dbinom_raw(n, x+n, p, 1-p, give_log);
2552 p = ((gnm_float)n)/(n+x);
2553 return((give_log) ? gnm_log(p) + prob : p * prob);
2554 }
2555
2556 /* ------------------------------------------------------------------------ */
2557 /* Imported src/nmath/pnbinom.c from R. */
2558 /*
2559 * Mathlib : A C Library of Special Functions
2560 * Copyright (C) 1998 Ross Ihaka
2561 * Copyright (C) 2000 The R Development Core Team
2562 *
2563 * This program is free software; you can redistribute it and/or modify
2564 * it under the terms of the GNU General Public License as published by
2565 * the Free Software Foundation; either version 2 of the License, or
2566 * (at your option) any later version.
2567 *
2568 * This program is distributed in the hope that it will be useful,
2569 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2570 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2571 * GNU General Public License for more details.
2572 *
2573 * You should have received a copy of the GNU General Public License
2574 * along with this program; if not, write to the Free Software
2575 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2576 * USA.
2577 *
2578 * DESCRIPTION
2579 *
2580 * The distribution function of the negative binomial distribution.
2581 *
2582 * NOTES
2583 *
2584 * x = the number of failures before the n-th success
2585 */
2586
2587
pnbinom(gnm_float x,gnm_float n,gnm_float p,gboolean lower_tail,gboolean log_p)2588 gnm_float pnbinom(gnm_float x, gnm_float n, gnm_float p, gboolean lower_tail, gboolean log_p)
2589 {
2590 #ifdef IEEE_754
2591 if (gnm_isnan(x) || gnm_isnan(n) || gnm_isnan(p))
2592 return x + n + p;
2593 if(!gnm_finite(n) || !gnm_finite(p)) ML_ERR_return_NAN;
2594 #endif
2595 if (n < 0 || p <= 0 || p > 1) ML_ERR_return_NAN;
2596
2597 x = gnm_fake_floor(x);
2598 if (x < 0) return R_DT_0;
2599 if (!gnm_finite(x)) return R_DT_1;
2600 return pbeta(p, n, x + 1, lower_tail, log_p);
2601 }
2602
2603 /* ------------------------------------------------------------------------ */
2604 /* Imported src/nmath/qnbinom.c from R. */
2605 /*
2606 * Mathlib : A C Library of Special Functions
2607 * Copyright (C) 1998 Ross Ihaka
2608 * Copyright (C) 2000 The R Development Core Team
2609 *
2610 * This program is free software; you can redistribute it and/or modify
2611 * it under the terms of the GNU General Public License as published by
2612 * the Free Software Foundation; either version 2 of the License, or
2613 * (at your option) any later version.
2614 *
2615 * This program is distributed in the hope that it will be useful,
2616 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2617 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2618 * GNU General Public License for more details.
2619 *
2620 * You should have received a copy of the GNU General Public License
2621 * along with this program; if not, write to the Free Software
2622 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2623 * USA.
2624 *
2625 * SYNOPSIS
2626 *
2627 * #include <Rmath.h>
2628 * double qnbinom(double p, double n, double pr, int lower_tail, int log_p)
2629 *
2630 * DESCRIPTION
2631 *
2632 * The quantile function of the negative binomial distribution.
2633 *
2634 * NOTES
2635 *
2636 * x = the number of failures before the n-th success
2637 *
2638 * METHOD
2639 *
2640 * Uses the Cornish-Fisher Expansion to include a skewness
2641 * correction to a normal approximation. This gives an
2642 * initial value which never seems to be off by more than
2643 * 1 or 2. A search is then conducted of values close to
2644 * this initial start point.
2645 */
2646
2647
qnbinom(gnm_float p,gnm_float n,gnm_float pr,gboolean lower_tail,gboolean log_p)2648 gnm_float qnbinom(gnm_float p, gnm_float n, gnm_float pr, gboolean lower_tail, gboolean log_p)
2649 {
2650 gnm_float P, Q, mu, sigma, gamma, z, y;
2651
2652 #ifdef IEEE_754
2653 if (gnm_isnan(p) || gnm_isnan(n) || gnm_isnan(pr))
2654 return p + n + pr;
2655 #endif
2656 R_Q_P01_check(p);
2657 if (pr <= 0 || pr >= 1 || n <= 0) ML_ERR_return_NAN;
2658
2659 if (p == R_DT_0) return 0;
2660 if (p == R_DT_1) return gnm_pinf;
2661 Q = 1.0 / pr;
2662 P = (1.0 - pr) * Q;
2663 mu = n * P;
2664 sigma = gnm_sqrt(n * P * Q);
2665 gamma = (Q + P)/sigma;
2666
2667 /* Note : "same" code in qpois.c, qbinom.c, qnbinom.c --
2668 * FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */
2669 if(!lower_tail || log_p) {
2670 p = R_DT_qIv(p); /* need check again (cancellation!): */
2671 if (p == R_DT_0) return 0;
2672 if (p == R_DT_1) return gnm_pinf;
2673 }
2674 /* temporary hack --- FIXME --- */
2675 if (p + 1.01*GNM_EPSILON >= 1.) return gnm_pinf;
2676
2677 /* y := approx.value (Cornish-Fisher expansion) : */
2678 z = qnorm(p, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE);
2679 y = gnm_floor(mu + sigma * (z + gamma * (z*z - 1) / 6) + 0.5);
2680
2681 z = pnbinom(y, n, pr, /*lower_tail*/TRUE, /*log_p*/FALSE);
2682
2683 /* fuzz to ensure left continuity: */
2684 p *= 1 - 64*GNM_EPSILON;
2685
2686 /*-- Fixme, here y can be way off --
2687 should use interval search instead of primitive stepping down or up */
2688
2689 #ifdef maybe_future
2690 if((lower_tail && z >= p) || (!lower_tail && z <= p)) {
2691 #else
2692 if(z >= p) {
2693 #endif
2694 /* search to the left */
2695 for(;;) {
2696 if(y == 0 ||
2697 (z = pnbinom(y - 1, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) < p)
2698 return y;
2699 y = y - 1;
2700 }
2701 }
2702 else { /* search to the right */
2703
2704 for(;;) {
2705 y = y + 1;
2706 if((z = pnbinom(y, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) >= p)
2707 return y;
2708 }
2709 }
2710 }
2711
2712 #if 0
2713 }
2714 #endif
2715 /* ------------------------------------------------------------------------ */
2716 /* Imported src/nmath/dbeta.c from R. */
2717 /*
2718 * AUTHOR
2719 * Catherine Loader, catherine@research.bell-labs.com.
2720 * October 23, 2000.
2721 *
2722 * Merge in to R:
2723 * Copyright (C) 2000, The R Core Development Team
2724 *
2725 * This program is free software; you can redistribute it and/or modify
2726 * it under the terms of the GNU General Public License as published by
2727 * the Free Software Foundation; either version 2 of the License, or
2728 * (at your option) any later version.
2729 *
2730 * This program is distributed in the hope that it will be useful,
2731 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2732 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2733 * GNU General Public License for more details.
2734 *
2735 * You should have received a copy of the GNU General Public License
2736 * along with this program; if not, write to the Free Software
2737 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2738 * USA.
2739 *
2740 *
2741 * DESCRIPTION
2742 * Beta density,
2743 * (a+b-1)! a-1 b-1
2744 * p(x;a,b) = ------------ x (1-x)
2745 * (a-1)!(b-1)!
2746 *
2747 * = (a+b-1) dbinom(a-1; a+b-2,x)
2748 *
2749 * We must modify this when a<1 or b<1, to avoid passing negative
2750 * arguments to dbinom_raw. Note that the modifications require
2751 * division by x and/or 1-x, so cannot be used except where necessary.
2752 */
2753
2754
dbeta(gnm_float x,gnm_float a,gnm_float b,gboolean give_log)2755 gnm_float dbeta(gnm_float x, gnm_float a, gnm_float b, gboolean give_log)
2756 {
2757 gnm_float f, p;
2758 volatile gnm_float am1, bm1; /* prevent roundoff trouble on some
2759 platforms */
2760
2761 #ifdef IEEE_754
2762 /* NaNs propagated correctly */
2763 if (gnm_isnan(x) || gnm_isnan(a) || gnm_isnan(b)) return x + a + b;
2764 #endif
2765
2766 if (a <= 0 || b <= 0) ML_ERR_return_NAN;
2767 if (x < 0 || x > 1) return(R_D__0);
2768 if (x == 0) {
2769 if(a > 1) return(R_D__0);
2770 if(a < 1) return(gnm_pinf);
2771 /* a == 1 : */ return(R_D_val(b));
2772 }
2773 if (x == 1) {
2774 if(b > 1) return(R_D__0);
2775 if(b < 1) return(gnm_pinf);
2776 /* b == 1 : */ return(R_D_val(a));
2777 }
2778 if (a < 1) {
2779 if (b < 1) { /* a,b < 1 */
2780 f = a*b/((a+b)*x*(1-x));
2781 p = dbinom_raw(a,a+b, x,1-x, give_log);
2782 }
2783 else { /* a < 1 <= b */
2784 f = a/x;
2785 bm1 = b - 1;
2786 p = dbinom_raw(a,a+bm1, x,1-x, give_log);
2787 }
2788 }
2789 else {
2790 if (b < 1) { /* a >= 1 > b */
2791 f = b/(1-x);
2792 am1 = a - 1;
2793 p = dbinom_raw(am1,am1+b, x,1-x, give_log);
2794 }
2795 else { /* a,b >= 1 */
2796 f = a+b-1;
2797 am1 = a - 1;
2798 bm1 = b - 1;
2799 p = dbinom_raw(am1,am1+bm1, x,1-x, give_log);
2800 }
2801 }
2802 return( (give_log) ? p + gnm_log(f) : p*f );
2803 }
2804
2805 /* ------------------------------------------------------------------------ */
2806 /* Imported src/nmath/dhyper.c from R. */
2807 /*
2808 * AUTHOR
2809 * Catherine Loader, catherine@research.bell-labs.com.
2810 * October 23, 2000.
2811 *
2812 * Merge in to R:
2813 * Copyright (C) 2000, 2001 The R Core Development Team
2814 *
2815 * This program is free software; you can redistribute it and/or modify
2816 * it under the terms of the GNU General Public License as published by
2817 * the Free Software Foundation; either version 2 of the License, or
2818 * (at your option) any later version.
2819 *
2820 * This program is distributed in the hope that it will be useful,
2821 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2822 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2823 * GNU General Public License for more details.
2824 *
2825 * You should have received a copy of the GNU General Public License
2826 * along with this program; if not, write to the Free Software
2827 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2828 * USA.
2829 *
2830 *
2831 * DESCRIPTION
2832 *
2833 * Given a sequence of r successes and b failures, we sample n (\le b+r)
2834 * items without replacement. The hypergeometric probability is the
2835 * probability of x successes:
2836 *
2837 * choose(r, x) * choose(b, n-x)
2838 * p(x; r,b,n) = ----------------------------- =
2839 * choose(r+b, n)
2840 *
2841 * dbinom(x,r,p) * dbinom(n-x,b,p)
2842 * = --------------------------------
2843 * dbinom(n,r+b,p)
2844 *
2845 * for any p. For numerical stability, we take p=n/(r+b); with this choice,
2846 * the denominator is not exponentially small.
2847 */
2848
2849
dhyper(gnm_float x,gnm_float r,gnm_float b,gnm_float n,gboolean give_log)2850 gnm_float dhyper(gnm_float x, gnm_float r, gnm_float b, gnm_float n, gboolean give_log)
2851 {
2852 gnm_float p, q, p1, p2, p3;
2853
2854 #ifdef IEEE_754
2855 if (gnm_isnan(x) || gnm_isnan(r) || gnm_isnan(b) || gnm_isnan(n))
2856 return x + r + b + n;
2857 #endif
2858
2859 if (R_D_negInonint(r) || R_D_negInonint(b) || R_D_negInonint(n) || n > r+b)
2860 ML_ERR_return_NAN;
2861 if (R_D_negInonint(x))
2862 return(R_D__0);
2863
2864 x = R_D_forceint(x);
2865 r = R_D_forceint(r);
2866 b = R_D_forceint(b);
2867 n = R_D_forceint(n);
2868
2869 if (n < x || r < x || n - x > b) return(R_D__0);
2870 if (n == 0) return((x == 0) ? R_D__1 : R_D__0);
2871
2872 /*
2873 * Round to float to make both p and q numbers with few (<26) bits set.
2874 * Unless p<2^-27 that also ensures that p+q=1 without rounding errors.
2875 */
2876 p = (float)(n / (gnm_float)(r+b));
2877 q = 1 - p;
2878
2879 p1 = dbinom_raw(x, r, p,q,give_log);
2880 p2 = dbinom_raw(n-x,b, p,q,give_log);
2881 p3 = dbinom_raw(n,r+b, p,q,give_log);
2882
2883 return( (give_log) ? p1 + p2 - p3 : p1*p2/p3 );
2884 }
2885
2886
2887 /* ------------------------------------------------------------------------ */
2888 /* Imported src/nmath/phyper.c from R. */
2889 /*
2890 * Mathlib : A C Library of Special Functions
2891 * Copyright (C) 1998 Ross Ihaka
2892 * Copyright (C) 1999-2000 The R Development Core Team
2893 * Copyright (C) 2004 Morten Welinder
2894 * Copyright (C) 2004 The R Foundation
2895 *
2896 * This program is free software; you can redistribute it and/or modify
2897 * it under the terms of the GNU General Public License as published by
2898 * the Free Software Foundation; either version 2 of the License, or
2899 * (at your option) any later version.
2900 *
2901 * This program is distributed in the hope that it will be useful,
2902 * but WITHOUT ANY WARRANTY; without even the implied warranty of
2903 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2904 * GNU General Public License for more details.
2905 *
2906 * You should have received a copy of the GNU General Public License
2907 * along with this program; if not, write to the Free Software
2908 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
2909 * USA.
2910 *
2911 * DESCRIPTION
2912 *
2913 * The distribution function of the hypergeometric distribution.
2914 *
2915 * Current implementation based on posting
2916 * From: Morten Welinder <terra@gnome.org>
2917 * Cc: R-bugs@biostat.ku.dk
2918 * Subject: [Rd] phyper accuracy and efficiency (PR#6772)
2919 * Date: Thu, 15 Apr 2004 18:06:37 +0200 (CEST)
2920 ......
2921
2922 The current version has very serious cancellation issues. For example,
2923 if you ask for a small right-tail you are likely to get total cancellation.
2924 For example, phyper(59, 150, 150, 60, FALSE, FALSE) gives 6.372680161e-14.
2925 The right answer is dhyper(0, 150, 150, 60, FALSE) which is 5.111204798e-22.
2926
2927 phyper is also really slow for large arguments.
2928
2929 Therefore, I suggest using the code below. This is a sniplet from Gnumeric ...
2930 The code isn't perfect. In fact, if x*(NR+NB) is close to n*NR,
2931 then this code can take a while. Not longer than the old code, though.
2932
2933 -- Thanks to Ian Smith for ideas.
2934 */
2935
2936
pdhyper(gnm_float x,gnm_float NR,gnm_float NB,gnm_float n,gboolean log_p)2937 static gnm_float pdhyper (gnm_float x, gnm_float NR, gnm_float NB, gnm_float n, gboolean log_p)
2938 {
2939 /*
2940 * Calculate
2941 *
2942 * phyper (x, NR, NB, n, TRUE, FALSE)
2943 * [log] ----------------------------------
2944 * dhyper (x, NR, NB, n, FALSE)
2945 *
2946 * without actually calling phyper. This assumes that
2947 *
2948 * x * (NR + NB) <= n * NR
2949 *
2950 */
2951 gnm_float sum = 0;
2952 gnm_float term = 1;
2953
2954 while (x > 0 && term > GNM_EPSILON * sum) {
2955 term *= x * (NB - n + x) / (n + 1 - x) / (NR + 1 - x);
2956 sum += term;
2957 x--;
2958 }
2959
2960 return log_p ? gnm_log1p(sum) : 1 + sum;
2961 }
2962
2963
2964 /* FIXME: The old phyper() code was basically used in ./qhyper.c as well
2965 * ----- We need to sync this again!
2966 */
phyper(gnm_float x,gnm_float NR,gnm_float NB,gnm_float n,gboolean lower_tail,gboolean log_p)2967 gnm_float phyper (gnm_float x, gnm_float NR, gnm_float NB, gnm_float n,
2968 gboolean lower_tail, gboolean log_p)
2969 {
2970 /* Sample of n balls from NR red and NB black ones; x are red */
2971
2972 gnm_float d, pd;
2973
2974 #ifdef IEEE_754
2975 if(gnm_isnan(x) || gnm_isnan(NR) || gnm_isnan(NB) || gnm_isnan(n))
2976 return x + NR + NB + n;
2977 #endif
2978
2979 x = gnm_fake_floor(x);
2980 NR = R_D_forceint(NR);
2981 NB = R_D_forceint(NB);
2982 n = R_D_forceint(n);
2983
2984 if (NR < 0 || NB < 0 || !gnm_finite(NR + NB) || n < 0 || n > NR + NB)
2985 ML_ERR_return_NAN;
2986
2987 if (x * (NR + NB) > n * NR) {
2988 /* Swap tails. */
2989 gnm_float oldNB = NB;
2990 NB = NR;
2991 NR = oldNB;
2992 x = n - x - 1;
2993 lower_tail = !lower_tail;
2994 }
2995
2996 if (x < 0)
2997 return R_DT_0;
2998 /* Warning: the following line is not in R: */
2999 if (x >= NR)
3000 return R_DT_1;
3001
3002 d = dhyper (x, NR, NB, n, log_p);
3003 pd = pdhyper(x, NR, NB, n, log_p);
3004
3005 return log_p ? R_DT_Log(d + pd) : R_D_Lval(d * pd);
3006 }
3007
3008 /* ------------------------------------------------------------------------ */
3009 /* Imported src/nmath/dexp.c from R. */
3010 /*
3011 * Mathlib : A C Library of Special Functions
3012 * Copyright (C) 1998 Ross Ihaka
3013 * Copyright (C) 2000 The R Development Core Team
3014 *
3015 * This program is free software; you can redistribute it and/or modify
3016 * it under the terms of the GNU General Public License as published by
3017 * the Free Software Foundation; either version 2 of the License, or
3018 * (at your option) any later version.
3019 *
3020 * This program is distributed in the hope that it will be useful,
3021 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3022 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3023 * GNU General Public License for more details.
3024 *
3025 * You should have received a copy of the GNU General Public License
3026 * along with this program; if not, write to the Free Software
3027 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3028 * USA.
3029 *
3030 * DESCRIPTION
3031 *
3032 * The density of the exponential distribution.
3033 */
3034
3035
dexp(gnm_float x,gnm_float scale,gboolean give_log)3036 gnm_float dexp(gnm_float x, gnm_float scale, gboolean give_log)
3037 {
3038 #ifdef IEEE_754
3039 /* NaNs propagated correctly */
3040 if (gnm_isnan(x) || gnm_isnan(scale)) return x + scale;
3041 #endif
3042 if (scale <= 0.0) ML_ERR_return_NAN;
3043
3044 if (x < 0.)
3045 return R_D__0;
3046 return (give_log ?
3047 (-x / scale) - gnm_log(scale) :
3048 gnm_exp(-x / scale) / scale);
3049 }
3050
3051 /* ------------------------------------------------------------------------ */
3052 /* Imported src/nmath/pexp.c from R. */
3053 /*
3054 * Mathlib : A C Library of Special Functions
3055 * Copyright (C) 1998 Ross Ihaka
3056 * Copyright (C) 2000-2002 The R Development Core Team
3057 *
3058 * This program is free software; you can redistribute it and/or modify
3059 * it under the terms of the GNU General Public License as published by
3060 * the Free Software Foundation; either version 2 of the License, or
3061 * (at your option) any later version.
3062 *
3063 * This program is distributed in the hope that it will be useful,
3064 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3065 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3066 * GNU General Public License for more details.
3067 *
3068 * You should have received a copy of the GNU General Public License
3069 * along with this program; if not, write to the Free Software
3070 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3071 * USA.
3072 *
3073 * DESCRIPTION
3074 *
3075 * The distribution function of the exponential distribution.
3076 */
3077
pexp(gnm_float x,gnm_float scale,gboolean lower_tail,gboolean log_p)3078 gnm_float pexp(gnm_float x, gnm_float scale, gboolean lower_tail, gboolean log_p)
3079 {
3080 #ifdef IEEE_754
3081 if (gnm_isnan(x) || gnm_isnan(scale))
3082 return x + scale;
3083 if (scale < 0) ML_ERR_return_NAN;
3084 #else
3085 if (scale <= 0) ML_ERR_return_NAN;
3086 #endif
3087
3088 if (x <= 0.)
3089 return R_DT_0;
3090 /* same as weibull( shape = 1): */
3091 x = -(x / scale);
3092 if (lower_tail)
3093 return (log_p
3094 /* gnm_log(1 - gnm_exp(x)) for x < 0 : */
3095 ? (x > -M_LN2gnum ? gnm_log(-gnm_expm1(x)) : gnm_log1p(-gnm_exp(x)))
3096 : -gnm_expm1(x));
3097 /* else: !lower_tail */
3098 return R_D_exp(x);
3099 }
3100
3101 /* ------------------------------------------------------------------------ */
3102 /* Imported src/nmath/dgeom.c from R. */
3103 /*
3104 * AUTHOR
3105 * Catherine Loader, catherine@research.bell-labs.com.
3106 * October 23, 2000.
3107 *
3108 * Merge in to R:
3109 * Copyright (C) 2000, 2001 The R Core Development Team
3110 *
3111 * This program is free software; you can redistribute it and/or modify
3112 * it under the terms of the GNU General Public License as published by
3113 * the Free Software Foundation; either version 2 of the License, or
3114 * (at your option) any later version.
3115 *
3116 * This program is distributed in the hope that it will be useful,
3117 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3118 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3119 * GNU General Public License for more details.
3120 *
3121 * You should have received a copy of the GNU General Public License
3122 * along with this program; if not, write to the Free Software
3123 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3124 * USA.
3125 *
3126 *
3127 * DESCRIPTION
3128 *
3129 * Computes the geometric probabilities, Pr(X=x) = p(1-p)^x.
3130 */
3131
3132
dgeom(gnm_float x,gnm_float p,gboolean give_log)3133 gnm_float dgeom(gnm_float x, gnm_float p, gboolean give_log)
3134 {
3135 gnm_float prob;
3136
3137 #ifdef IEEE_754
3138 if (gnm_isnan(x) || gnm_isnan(p)) return x + p;
3139 #endif
3140
3141 if (p < 0 || p > 1) ML_ERR_return_NAN;
3142
3143 R_D_nonint_check(x);
3144 if (x < 0 || !gnm_finite(x) || p == 0) return R_D__0;
3145 x = R_D_forceint(x);
3146
3147 /* prob = (1-p)^x, stable for small p */
3148 prob = dbinom_raw(0.,x, p,1-p, give_log);
3149
3150 return((give_log) ? gnm_log(p) + prob : p*prob);
3151 }
3152
3153 /* ------------------------------------------------------------------------ */
3154 /* Imported src/nmath/pgeom.c from R. */
3155 /*
3156 * Mathlib : A C Library of Special Functions
3157 * Copyright (C) 1998 Ross Ihaka
3158 * Copyright (C) 2000, 2001 The R Development Core Team
3159 * Copyright (C) 2004 The R Foundation
3160 *
3161 * This program is free software; you can redistribute it and/or modify
3162 * it under the terms of the GNU General Public License as published by
3163 * the Free Software Foundation; either version 2 of the License, or
3164 * (at your option) any later version.
3165 *
3166 * This program is distributed in the hope that it will be useful,
3167 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3168 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3169 * GNU General Public License for more details.
3170 *
3171 * You should have received a copy of the GNU General Public License
3172 * along with this program; if not, write to the Free Software
3173 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3174 * USA.
3175 *
3176 * DESCRIPTION
3177 *
3178 * The distribution function of the geometric distribution.
3179 */
3180
3181
pgeom(gnm_float x,gnm_float p,gboolean lower_tail,gboolean log_p)3182 gnm_float pgeom(gnm_float x, gnm_float p, gboolean lower_tail, gboolean log_p)
3183 {
3184 #ifdef IEEE_754
3185 if (gnm_isnan(x) || gnm_isnan(p))
3186 return x + p;
3187 #endif
3188 x = gnm_fake_floor(x);
3189 if(p < 0 || p > 1) ML_ERR_return_NAN;
3190
3191 if (x < 0. || p == 0.) return R_DT_0;
3192 if (!gnm_finite(x)) return R_DT_1;
3193
3194 if(p == 1.) { /* we cannot assume IEEE */
3195 x = lower_tail ? 1: 0;
3196 return log_p ? gnm_log(x) : x;
3197 }
3198 x = gnm_log1p(-p) * (x + 1);
3199 if (log_p)
3200 return R_DT_Clog(x);
3201 else
3202 return lower_tail ? -gnm_expm1(x) : gnm_exp(x);
3203 }
3204
3205 /* ------------------------------------------------------------------------ */
3206 /* Imported src/nmath/dcauchy.c from R. */
3207 /*
3208 * Mathlib : A C Library of Special Functions
3209 * Copyright (C) 1998 Ross Ihaka
3210 * Copyright (C) 2000 The R Development Core Team
3211 *
3212 * This program is free software; you can redistribute it and/or modify
3213 * it under the terms of the GNU General Public License as published by
3214 * the Free Software Foundation; either version 2 of the License, or
3215 * (at your option) any later version.
3216 *
3217 * This program is distributed in the hope that it will be useful,
3218 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3219 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3220 * GNU General Public License for more details.
3221 *
3222 * You should have received a copy of the GNU General Public License
3223 * along with this program; if not, write to the Free Software
3224 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
3225 * 02110-1301 USA.
3226 *
3227 * DESCRIPTION
3228 *
3229 * The density of the Cauchy distribution.
3230 */
3231
3232
dcauchy(gnm_float x,gnm_float location,gnm_float scale,gboolean give_log)3233 gnm_float dcauchy(gnm_float x, gnm_float location, gnm_float scale, gboolean give_log)
3234 {
3235 gnm_float y;
3236 #ifdef IEEE_754
3237 /* NaNs propagated correctly */
3238 if (gnm_isnan(x) || gnm_isnan(location) || gnm_isnan(scale))
3239 return x + location + scale;
3240 #endif
3241 if (scale <= 0) ML_ERR_return_NAN;
3242
3243 y = (x - location) / scale;
3244 return give_log ?
3245 - gnm_log(M_PIgnum * scale * (1. + y * y)) :
3246 1. / (M_PIgnum * scale * (1. + y * y));
3247 }
3248
3249 /* ------------------------------------------------------------------------ */
3250 /* Imported src/nmath/pcauchy.c from R. */
3251 /*
3252 * Mathlib : A C Library of Special Functions
3253 * Copyright (C) 1998 Ross Ihaka
3254 * Copyright (C) 2000 The R Development Core Team
3255 * Copyright (C) 2004 The R Foundation
3256 *
3257 * This program is free software; you can redistribute it and/or modify
3258 * it under the terms of the GNU General Public License as published by
3259 * the Free Software Foundation; either version 2 of the License, or
3260 * (at your option) any later version.
3261 *
3262 * This program is distributed in the hope that it will be useful,
3263 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3264 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3265 * GNU General Public License for more details.
3266 *
3267 * You should have received a copy of the GNU General Public License
3268 * along with this program; if not, write to the Free Software
3269 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3270 * USA.
3271 *
3272 * DESCRIPTION
3273 *
3274 * The distribution function of the Cauchy distribution.
3275 */
3276
3277
pcauchy(gnm_float x,gnm_float location,gnm_float scale,gboolean lower_tail,gboolean log_p)3278 gnm_float pcauchy(gnm_float x, gnm_float location, gnm_float scale,
3279 gboolean lower_tail, gboolean log_p)
3280 {
3281 #ifdef IEEE_754
3282 if (gnm_isnan(x) || gnm_isnan(location) || gnm_isnan(scale))
3283 return x + location + scale;
3284 #endif
3285 if (scale <= 0) ML_ERR_return_NAN;
3286
3287 x = (x - location) / scale;
3288 if (gnm_isnan(x)) ML_ERR_return_NAN;
3289 #ifdef IEEE_754
3290 if(!gnm_finite(x)) {
3291 if(x < 0) return R_DT_0;
3292 else return R_DT_1;
3293 }
3294 #endif
3295 if (!lower_tail)
3296 x = -x;
3297
3298 if (log_p && x > 0)
3299 return R_D_Clog(gnm_atan2pi (1, x));
3300 else
3301 return R_D_val(gnm_atan2pi (1, -x));
3302 }
3303
3304 /* ------------------------------------------------------------------------ */
3305 /* Imported src/nmath/df.c from R. */
3306 /*
3307 * AUTHOR
3308 * Catherine Loader, catherine@research.bell-labs.com.
3309 * October 23, 2000.
3310 *
3311 * Merge in to R:
3312 * Copyright (C) 2000, The R Core Development Team
3313 *
3314 * This program is free software; you can redistribute it and/or modify
3315 * it under the terms of the GNU General Public License as published by
3316 * the Free Software Foundation; either version 2 of the License, or
3317 * (at your option) any later version.
3318 *
3319 * This program is distributed in the hope that it will be useful,
3320 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3321 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3322 * GNU General Public License for more details.
3323 *
3324 * You should have received a copy of the GNU General Public License
3325 * along with this program; if not, write to the Free Software
3326 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3327 * USA.
3328 *
3329 *
3330 * DESCRIPTION
3331 *
3332 * The density function of the F distribution.
3333 * To evaluate it, write it as a Binomial probability with p = x*m/(n+x*m).
3334 * For m >= 2, we use the simplest conversion.
3335 * For m < 2, (m-2)/2 < 0 so the conversion will not work, and we must use
3336 * a second conversion.
3337 * Note the division by p; this seems unavoidable
3338 * for m < 2, since the F density has a singularity as x (or p) -> 0.
3339 */
3340
3341
df(gnm_float x,gnm_float m,gnm_float n,gboolean give_log)3342 gnm_float df(gnm_float x, gnm_float m, gnm_float n, gboolean give_log)
3343 {
3344 gnm_float p, q, f, dens;
3345
3346 #ifdef IEEE_754
3347 if (gnm_isnan(x) || gnm_isnan(m) || gnm_isnan(n))
3348 return x + m + n;
3349 #endif
3350 if (m <= 0 || n <= 0) ML_ERR_return_NAN;
3351 if (x <= 0.) return(R_D__0);
3352
3353 f = 1./(n+x*m);
3354 q = n*f;
3355 p = x*m*f;
3356
3357 if (m >= 2) {
3358 f = m*q/2;
3359 dens = dbinom_raw((m-2)/2, (m+n-2)/2, p, q, give_log);
3360 }
3361 else {
3362 f = m*m*q / (2*p*(m+n));
3363 dens = dbinom_raw(m/2, (m+n)/2, p, q, give_log);
3364 }
3365 return(give_log ? gnm_log(f)+dens : f*dens);
3366 }
3367
3368 /* ------------------------------------------------------------------------ */
3369 /* Imported src/nmath/dchisq.c from R. */
3370 /*
3371 * Mathlib : A C Library of Special Functions
3372 * Copyright (C) 1998 Ross Ihaka
3373 * Copyright (C) 2000 The R Development Core Team
3374 *
3375 * This program is free software; you can redistribute it and/or modify
3376 * it under the terms of the GNU General Public License as published by
3377 * the Free Software Foundation; either version 2 of the License, or
3378 * (at your option) any later version.
3379 *
3380 * This program is distributed in the hope that it will be useful, but
3381 * WITHOUT ANY WARRANTY; without even the implied warranty of
3382 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3383 * General Public License for more details.
3384 *
3385 * You should have received a copy of the GNU General Public License
3386 * along with this program; if not, write to the Free Software
3387 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3388 * USA
3389 *
3390 * DESCRIPTION
3391 *
3392 * The density of the chi-squared distribution.
3393 */
3394
3395
dchisq(gnm_float x,gnm_float df,gboolean give_log)3396 gnm_float dchisq(gnm_float x, gnm_float df, gboolean give_log)
3397 {
3398 return dgamma(x, df / 2., 2., give_log);
3399 }
3400
3401 /* ------------------------------------------------------------------------ */
3402 /* Imported src/nmath/qweibull.c from R. */
3403 /*
3404 * Mathlib : A C Library of Special Functions
3405 * Copyright (C) 1998 Ross Ihaka
3406 * Copyright (C) 2000 The R Development Core Team
3407 *
3408 * This program is free software; you can redistribute it and/or modify
3409 * it under the terms of the GNU General Public License as published by
3410 * the Free Software Foundation; either version 2 of the License, or
3411 * (at your option) any later version.
3412 *
3413 * This program is distributed in the hope that it will be useful,
3414 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3415 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3416 * GNU General Public License for more details.
3417 *
3418 * You should have received a copy of the GNU General Public License
3419 * along with this program; if not, write to the Free Software
3420 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3421 * USA.
3422 *
3423 * DESCRIPTION
3424 *
3425 * The quantile function of the Weibull distribution.
3426 */
3427
3428
qweibull(gnm_float p,gnm_float shape,gnm_float scale,gboolean lower_tail,gboolean log_p)3429 gnm_float qweibull(gnm_float p, gnm_float shape, gnm_float scale, gboolean lower_tail, gboolean log_p)
3430 {
3431 #ifdef IEEE_754
3432 if (gnm_isnan(p) || gnm_isnan(shape) || gnm_isnan(scale))
3433 return p + shape + scale;
3434 #endif
3435 R_Q_P01_check(p);
3436 if (shape <= 0 || scale <= 0) ML_ERR_return_NAN;
3437
3438 if (p == R_D__0) return 0;
3439 if (p == R_D__1) return gnm_pinf;
3440 return scale * gnm_pow(- R_DT_Clog(p), 1./shape) ;
3441 }
3442
3443 /* ------------------------------------------------------------------------ */
3444 /* Imported src/nmath/qexp.c from R. */
3445 /*
3446 * Mathlib : A C Library of Special Functions
3447 * Copyright (C) 1998 Ross Ihaka
3448 * Copyright (C) 2000 The R Development Core Team
3449 *
3450 * This program is free software; you can redistribute it and/or modify
3451 * it under the terms of the GNU General Public License as published by
3452 * the Free Software Foundation; either version 2 of the License, or
3453 * (at your option) any later version.
3454 *
3455 * This program is distributed in the hope that it will be useful,
3456 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3457 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3458 * GNU General Public License for more details.
3459 *
3460 * You should have received a copy of the GNU General Public License
3461 * along with this program; if not, write to the Free Software
3462 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3463 * USA.
3464 *
3465 * DESCRIPTION
3466 *
3467 * The quantile function of the exponential distribution.
3468 *
3469 */
3470
3471
qexp(gnm_float p,gnm_float scale,gboolean lower_tail,gboolean log_p)3472 gnm_float qexp(gnm_float p, gnm_float scale, gboolean lower_tail, gboolean log_p)
3473 {
3474 #ifdef IEEE_754
3475 if (gnm_isnan(p) || gnm_isnan(scale))
3476 return p + scale;
3477 #endif
3478 R_Q_P01_check(p);
3479 if (scale < 0) ML_ERR_return_NAN;
3480
3481 if (p == R_DT_0)
3482 return 0;
3483
3484 return - scale * R_DT_Clog(p);
3485 }
3486
3487 /* ------------------------------------------------------------------------ */
3488 /* Imported src/nmath/qgeom.c from R. */
3489 /*
3490 * Mathlib : A C Library of Special Functions
3491 * Copyright (C) 1998 Ross Ihaka
3492 * Copyright (C) 2000 The R Development Core Team
3493 * Copyright (C) 2004 The R Foundation
3494 *
3495 * This program is free software; you can redistribute it and/or modify
3496 * it under the terms of the GNU General Public License as published by
3497 * the Free Software Foundation; either version 2 of the License, or
3498 * (at your option) any later version.
3499 *
3500 * This program is distributed in the hope that it will be useful,
3501 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3502 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3503 * GNU General Public License for more details.
3504 *
3505 * You should have received a copy of the GNU General Public License
3506 * along with this program; if not, write to the Free Software
3507 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
3508 * USA.
3509 *
3510 * DESCRIPTION
3511 *
3512 * The quantile function of the geometric distribution.
3513 */
3514
3515
qgeom(gnm_float p,gnm_float prob,gboolean lower_tail,gboolean log_p)3516 gnm_float qgeom(gnm_float p, gnm_float prob, gboolean lower_tail, gboolean log_p)
3517 {
3518 gnm_float q;
3519
3520 R_Q_P01_check(p);
3521 if (prob <= 0 || prob > 1) ML_ERR_return_NAN;
3522
3523 #ifdef IEEE_754
3524 if (gnm_isnan(p) || gnm_isnan(prob))
3525 return p + prob;
3526 if (p == R_DT_1) return gnm_pinf;
3527 #else
3528 if (p == R_DT_1) ML_ERR_return_NAN;
3529 #endif
3530 if (p == R_DT_0) return 0;
3531
3532 /* add a fuzz to ensure left continuity */
3533 q = gnm_ceil(R_DT_Clog(p) / gnm_log1p(- prob) - 1 - 1e-12);
3534 return MAX (q, 0.0);
3535 }
3536
3537 /* ------------------------------------------------------------------------ */
3538 /*
3539 * Based on code imported from R by hand. Heavily modified to enhance
3540 * accuracy. See bug 700132.
3541 */
3542 /*
3543 * Mathlib : A C Library of Special Functions
3544 * Copyright (C) 1998 Ross Ihaka
3545 * Copyright (C) 2000--2007 The R Core Team
3546 *
3547 * This program is free software; you can redistribute it and/or modify
3548 * it under the terms of the GNU General Public License as published by
3549 * the Free Software Foundation; either version 2 of the License, or
3550 * (at your option) any later version.
3551 *
3552 * This program is distributed in the hope that it will be useful,
3553 * but WITHOUT ANY WARRANTY; without even the implied warranty of
3554 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3555 * GNU General Public License for more details.
3556 *
3557 * You should have received a copy of the GNU General Public License
3558 * along with this program; if not, a copy is available at
3559 * http://www.r-project.org/Licenses/
3560 *
3561 * SYNOPSIS
3562 *
3563 * #include <Rmath.h>
3564 * double ptukey(q, rr, cc, df, lower_tail, log_p);
3565 *
3566 * DESCRIPTION
3567 *
3568 * Computes the probability that the maximum of rr studentized
3569 * ranges, each based on cc means and with df degrees of freedom
3570 * for the standard error, is less than q.
3571 *
3572 * The algorithm is based on that of the reference.
3573 *
3574 * REFERENCE
3575 *
3576 * Copenhaver, Margaret Diponzio & Holland, Burt S.
3577 * Multiple comparisons of simple effects in
3578 * the two-way analysis of variance with fixed effects.
3579 * Journal of Statistical Computation and Simulation,
3580 * Vol.30, pp.1-15, 1988.
3581 */
3582
ptukey_wprob(gnm_float w,gnm_float rr,gnm_float cc)3583 static gnm_float ptukey_wprob(gnm_float w, gnm_float rr, gnm_float cc)
3584 {
3585 /* wprob() :
3586
3587 This function calculates probability integral of Hartley's
3588 form of the range.
3589
3590 w = value of range
3591 rr = no. of rows or groups
3592 cc = no. of columns or treatments
3593 pr_w = returned probability integral
3594
3595 program will not terminate if ir is raised.
3596
3597 bb = upper limit of legendre integration
3598 nleg = order of legendre quadrature
3599 wlar = value of range above which wincr1 intervals are used to
3600 calculate second part of integral,
3601 else wincr2 intervals are used.
3602 or modifying a calculation.
3603
3604 M_1_SQRT_2PI = 1 / sqrt(2 * pi); from abramowitz & stegun, p. 3.
3605 M_SQRT2gnum = sqrt(2)
3606 xleg = legendre 12-point nodes
3607 aleg = legendre 12-point coefficients
3608 */
3609
3610 static const gboolean debug = FALSE;
3611 static const gnm_float xleg[] = {
3612 GNM_const(0.981560634246719250690549090149),
3613 GNM_const(0.904117256370474856678465866119),
3614 GNM_const(0.769902674194304687036893833213),
3615 GNM_const(0.587317954286617447296702418941),
3616 GNM_const(0.367831498998180193752691536644),
3617 GNM_const(0.125233408511468915472441369464)
3618 };
3619 static const gnm_float aleg[G_N_ELEMENTS (xleg)] = {
3620 GNM_const(0.047175336386511827194615961485),
3621 GNM_const(0.106939325995318430960254718194),
3622 GNM_const(0.160078328543346226334652529543),
3623 GNM_const(0.203167426723065921749064455810),
3624 GNM_const(0.233492536538354808760849898925),
3625 GNM_const(0.249147045813402785000562436043)
3626 };
3627 const int nleg = G_N_ELEMENTS (xleg) * 2;
3628 gnm_float pr_w, binc, qsqz, blb;
3629 int i, jj;
3630
3631 qsqz = w * 0.5;
3632
3633 /* find (2F(w/2) - 1) ^ cc */
3634 /* (first term in integral of hartley's form). */
3635
3636 /*
3637 * Use different formulas for different size of qsqz to avoid
3638 * cancellation for pnorm or squeezing erf's result against 1.
3639 */
3640 pr_w = qsqz > 1
3641 ? pow1p (-2.0 * pnorm (qsqz, 0, 1, FALSE, FALSE), cc)
3642 : gnm_pow (gnm_erf (qsqz / M_SQRT2gnum), cc);
3643 if (pr_w >= 1.0)
3644 return 1.0;
3645
3646 /* find the integral of second term of hartley's form */
3647 /* Limits of integration are (w/2, Inf). Large cc means that */
3648 /* that we need smaller step size. The formula for binc is */
3649 /* a heuristic. */
3650
3651 /* blb and bub are lower and upper limits of integration. */
3652
3653 blb = qsqz;
3654 binc = 3.0 / gnm_log1p (cc);
3655
3656 /* integrate over each interval */
3657
3658 for (i = 1; TRUE; i++) {
3659 gnm_float C = blb + binc * 0.5;
3660 gnm_float elsum = 0.0;
3661
3662 /* legendre quadrature with order = nleg */
3663
3664 for (jj = 0; jj < nleg; jj++) {
3665 gnm_float xx, aa, v, rinsum;
3666
3667 if (nleg / 2 <= jj) {
3668 xx = xleg[nleg - 1 - jj];
3669 aa = aleg[nleg - 1 - jj];
3670 } else {
3671 xx = -xleg[jj];
3672 aa = aleg[jj];
3673 }
3674 v = C + binc * 0.5 * xx;
3675
3676 rinsum = pnorm2 (v - w, v);
3677 elsum += gnm_pow (rinsum, cc - 1) * aa * expmx2h(v);
3678 }
3679 elsum *= (binc * cc * M_1_SQRT_2PI);
3680 pr_w += elsum;
3681
3682 if (pr_w >= 1) {
3683 /* Nothing more will contribute. */
3684 pr_w = 1;
3685 break;
3686 }
3687
3688 if (elsum <= pr_w * (GNM_EPSILON / 2)) {
3689 /* This interval contributed nothing. We're done. */
3690 if (debug) {
3691 g_printerr ("End at i=%d for w=%g cc=%g dP=%g P=%g\n",
3692 i, w, cc, elsum, pr_w);
3693 }
3694 break;
3695 }
3696
3697 blb += binc;
3698 }
3699
3700 if (0) g_printerr ("w=%g pr_w=%.20g\n", w, pr_w);
3701
3702 return gnm_pow (pr_w, rr);
3703 }
3704
3705 static gnm_float
ptukey_otsum(gnm_float u0,gnm_float u1,gnm_float f2,gnm_float f2lf,gnm_float q,gnm_float rr,gnm_float cc)3706 ptukey_otsum (gnm_float u0, gnm_float u1, gnm_float f2, gnm_float f2lf,
3707 gnm_float q, gnm_float rr, gnm_float cc)
3708 {
3709 gboolean debug = FALSE;
3710 static const gnm_float xlegq[] = {
3711 GNM_const(0.989400934991649932596154173450),
3712 GNM_const(0.944575023073232576077988415535),
3713 GNM_const(0.865631202387831743880467897712),
3714 GNM_const(0.755404408355003033895101194847),
3715 GNM_const(0.617876244402643748446671764049),
3716 GNM_const(0.458016777657227386342419442984),
3717 GNM_const(0.281603550779258913230460501460),
3718 GNM_const(0.950125098376374401853193354250e-1)
3719 };
3720 static const gnm_float alegq[G_N_ELEMENTS (xlegq)] = {
3721 GNM_const(0.271524594117540948517805724560e-1),
3722 GNM_const(0.622535239386478928628438369944e-1),
3723 GNM_const(0.951585116824927848099251076022e-1),
3724 GNM_const(0.124628971255533872052476282192),
3725 GNM_const(0.149595988816576732081501730547),
3726 GNM_const(0.169156519395002538189312079030),
3727 GNM_const(0.182603415044923588866763667969),
3728 GNM_const(0.189450610455068496285396723208)
3729 };
3730 const int nlegq = G_N_ELEMENTS (xlegq) * 2;
3731 int jj;
3732 gnm_float C = 0.5 * (u0 + u1);
3733 gnm_float L = u1 - u0;
3734 gnm_float otsum = 0.0;
3735
3736 for (jj = 0; jj < nlegq; jj++) {
3737 gnm_float xx, aa, u, t1, wprb;
3738
3739 if (nlegq / 2 <= jj) {
3740 xx = xlegq[nlegq - 1 - jj];
3741 aa = alegq[nlegq - 1 - jj];
3742 } else {
3743 xx = -xlegq[jj];
3744 aa = alegq[jj];
3745 }
3746
3747 u = xx * (0.5 * L) + C;
3748
3749 t1 = f2lf + (f2 - 1) * gnm_log(u) - u * f2;
3750
3751 wprb = ptukey_wprob(q * gnm_sqrt(u), rr, cc);
3752 otsum += aa * (wprb * (0.5 * L) * gnm_exp(t1));
3753 }
3754
3755 if (debug)
3756 g_printerr ("Integral over [%g;%g] was %g\n",
3757 u0, u1, otsum);
3758
3759
3760 return otsum;
3761 }
3762
3763
3764 static gnm_float
R_ptukey(gnm_float q,gnm_float rr,gnm_float cc,gnm_float df,gboolean lower_tail,gboolean log_p)3765 R_ptukey(gnm_float q, gnm_float rr, gnm_float cc, gnm_float df,
3766 gboolean lower_tail, gboolean log_p)
3767 {
3768 /* function ptukey() [was qprob() ]:
3769
3770 q = value of studentized range
3771 rr = no. of rows or groups
3772 cc = no. of columns or treatments
3773 df = degrees of freedom of error term
3774 ir[0] = error flag = 1 if wprob probability > 1
3775 ir[1] = error flag = 1 if qprob probability > 1
3776
3777 qprob = returned probability integral over [0, q]
3778
3779 The program will not terminate if ir[0] or ir[1] are raised.
3780
3781 All references in wprob to Abramowitz and Stegun
3782 are from the following reference:
3783
3784 Abramowitz, Milton and Stegun, Irene A.
3785 Handbook of Mathematical Functions.
3786 New York: Dover publications, Inc. (1970).
3787
3788 All constants taken from this text are
3789 given to 25 significant digits.
3790
3791 nlegq = order of legendre quadrature
3792 eps = max. allowable value of integral
3793 eps1 & eps2 = values below which there is
3794 no contribution to integral.
3795
3796 d.f. <= dhaf: integral is divided into ulen1 length intervals. else
3797 d.f. <= dquar: integral is divided into ulen2 length intervals. else
3798 d.f. <= deigh: integral is divided into ulen3 length intervals. else
3799 d.f. <= dlarg: integral is divided into ulen4 length intervals.
3800
3801 d.f. > dlarg: the range is used to calculate integral.
3802
3803 xlegq = legendre 16-point nodes
3804 alegq = legendre 16-point coefficients
3805
3806 The coefficients and nodes for the legendre quadrature used in
3807 qprob and wprob were calculated using the algorithms found in:
3808
3809 Stroud, A. H. and Secrest, D.
3810 Gaussian Quadrature Formulas.
3811 Englewood Cliffs,
3812 New Jersey: Prentice-Hall, Inc, 1966.
3813
3814 All values matched the tables (provided in same reference)
3815 to 30 significant digits.
3816
3817 F(x) = .5 + erf(x / sqrt(2)) / 2 for x > 0
3818
3819 F(x) = erfc( -x / sqrt(2)) / 2 for x < 0
3820
3821 where F(x) is standard normal c. d. f.
3822
3823 if degrees of freedom large, approximate integral
3824 with range distribution.
3825 */
3826
3827 static const gnm_float dhaf = 100.0;
3828 static const gnm_float dquar = 800.0;
3829 static const gnm_float deigh = 5000.0;
3830 static const gnm_float dlarg = 25000.0;
3831 static const gnm_float ulen1 = 1.0;
3832 static const gnm_float ulen2 = 0.5;
3833 static const gnm_float ulen3 = 0.25;
3834 static const gnm_float ulen4 = 0.125;
3835 gnm_float ans, f2, f2lf, ulen, u0, u1;
3836 int i;
3837 gboolean fail = FALSE;
3838 gboolean debug = TRUE;
3839
3840 #ifdef IEEE_754
3841 if (gnm_isnan(q) || gnm_isnan(rr) || gnm_isnan(cc) || gnm_isnan(df))
3842 ML_ERR_return_NAN;
3843 #endif
3844
3845 if (q <= 0)
3846 return R_DT_0;
3847
3848 /* FIXME: Special case for cc==2&&rr=1: we have explicit formula */
3849
3850
3851 /* df must be > 1 */
3852 /* there must be at least two values */
3853
3854 if (df < 2 || rr < 1 || cc < 2) ML_ERR_return_NAN;
3855
3856 if(!gnm_finite(q))
3857 return R_DT_1;
3858
3859 if (df > dlarg)
3860 return R_DT_val(ptukey_wprob(q, rr, cc));
3861
3862 /* calculate leading constant */
3863
3864 f2 = df * 0.5;
3865 /* gnm_lgamma(u) = log(gamma(u)) */
3866 f2lf = f2 * gnm_log(f2) - gnm_lgamma(f2);
3867
3868 /* integral is divided into unit, half-unit, quarter-unit, or */
3869 /* eighth-unit length intervals depending on the value of the */
3870 /* degrees of freedom. */
3871
3872 if (df <= dhaf) ulen = ulen1;
3873 else if (df <= dquar) ulen = ulen2;
3874 else if (df <= deigh) ulen = ulen3;
3875 else ulen = ulen4;
3876
3877 /* integrate over each subinterval */
3878 ans = 0.0;
3879
3880 /*
3881 * Integration for [0;ulen/2].
3882 *
3883 * We use a sequence of intervals each covering an ever increasing
3884 * fraction of what is left. Breaking the interval takes care of
3885 * some very un-polynomial behaviour of the integrand:
3886 * (1) Infinite derivative at 0 for low df.
3887 * (2) Infinite roots (due to underflow) in the left part of an
3888 * interval with meaningful contributions from its right part
3889 */
3890 u1 = ulen / 2;
3891 for (i = 1; TRUE; i++) {
3892 gnm_float otsum;
3893
3894 u0 = u1 / (i + 1);
3895 otsum = ptukey_otsum (u0, u1, f2, f2lf, q, rr, cc);
3896
3897 ans += otsum;
3898 if (otsum <= ans * (GNM_EPSILON / 2)) {
3899 /* This interval contributed nothing. We're done. */
3900 break;
3901 }
3902
3903 if (i == 20) {
3904 if (debug)
3905 g_printerr ("PTUKEY FAIL LEFT: %d q=%g cc=%g df=%g otsum=%g ans=%g\n",
3906 i, q, cc, df, otsum, ans);
3907 fail = TRUE;
3908 break;
3909 }
3910
3911 u1 = u0;
3912 }
3913
3914 /*
3915 * Integration for [ulen/2;Infinity].
3916 *
3917 * We use a sequence of intervals starting with length ulen, but
3918 * when we think we're in the tail we ramp up the length.
3919 */
3920 u0 = ulen / 2;
3921 for (i = 1; TRUE; i++) {
3922 gnm_float otsum;
3923
3924 u1 = u0 + ulen;
3925 otsum = ptukey_otsum (u0, u1, f2, f2lf, q, rr, cc);
3926
3927 ans += otsum;
3928 if (otsum < ans * GNM_EPSILON) {
3929 /*
3930 * This interval contributed nothing. This can either be
3931 * because the integrand is so flat that we haven't seen
3932 * anyting yet or because we're done. Or both.
3933 *
3934 * The maximum of the integrand falls not far from
3935 * u=(f-2)/f,
3936 * but let's go a little further.
3937 */
3938 if (ans > 0 || u0 > 2) {
3939 break;
3940 }
3941 }
3942
3943 if (i == 150) {
3944 if (debug)
3945 g_printerr ("PTUKEY FAIL RIGHT: %i %g %g\n",
3946 i, otsum, ans);
3947 fail = TRUE;
3948 break;
3949 }
3950
3951 /*
3952 * When we see a contribution that added less than 0.1%
3953 * to the result, start ramping up the interval size.
3954 * Note that this will not trigger unless ans>0.
3955 */
3956 if (otsum < ans / 1000)
3957 ulen *= 2;
3958
3959 u0 = u1;
3960 }
3961
3962 if (fail)
3963 ML_ERROR(ME_PRECISION);
3964 ans = MIN (ans, 1.0);
3965 return R_DT_val(ans);
3966 }
3967
3968 /* ------------------------------------------------------------------------ */
3969 /* --- END MAGIC R SOURCE MARKER --- */
3970
3971 /* Silly order-of-arguments wrapper. */
3972 gnm_float
ptukey(gnm_float x,gnm_float nmeans,gnm_float df,gnm_float nranges,gboolean lower_tail,gboolean log_p)3973 ptukey(gnm_float x, gnm_float nmeans, gnm_float df, gnm_float nranges, gboolean lower_tail, gboolean log_p)
3974 {
3975 return R_ptukey (x, nranges, nmeans, df, lower_tail, log_p);
3976 }
3977
3978 static gnm_float
ptukey1(gnm_float x,const gnm_float shape[],gboolean lower_tail,gboolean log_p)3979 ptukey1 (gnm_float x, const gnm_float shape[],
3980 gboolean lower_tail, gboolean log_p)
3981 {
3982 return ptukey (x, shape[0], shape[1], shape[2], lower_tail, log_p);
3983 }
3984
3985 gnm_float
qtukey(gnm_float p,gnm_float nmeans,gnm_float df,gnm_float nranges,gboolean lower_tail,gboolean log_p)3986 qtukey (gnm_float p, gnm_float nmeans, gnm_float df, gnm_float nranges, gboolean lower_tail, gboolean log_p)
3987 {
3988 gnm_float x0, shape[3];
3989
3990 if (!log_p && p > 0.9) {
3991 /* We're far into the tail. Flip. */
3992 p = 1 - p;
3993 lower_tail = !lower_tail;
3994 }
3995
3996 /* This is accurate for nmeans==2 and nranges==1. */
3997 x0 = M_SQRT2gnum * qt ((1 + p) / 2, df, lower_tail, log_p);
3998
3999 shape[0] = nmeans;
4000 shape[1] = df;
4001 shape[2] = nranges;
4002
4003 return pfuncinverter (p, shape, lower_tail, log_p, 0, gnm_pinf, x0, ptukey1, NULL);
4004 }
4005
4006 static gnm_float
logspace_signed_add(gnm_float logx,gnm_float logabsy,gboolean ypos)4007 logspace_signed_add (gnm_float logx, gnm_float logabsy, gboolean ypos)
4008 {
4009 return ypos
4010 ? logspace_add (logx, logabsy)
4011 : logspace_sub (logx, logabsy);
4012 }
4013
4014 /* Accurate gnm_log (1 - gnm_exp (p)) for p <= 0. */
4015 gnm_float
swap_log_tail(gnm_float lp)4016 swap_log_tail (gnm_float lp)
4017 {
4018 if (lp > -1 / gnm_log (2))
4019 return gnm_log (-gnm_expm1 (lp)); /* Good formula for lp near zero. */
4020 else
4021 return gnm_log1p (-gnm_exp (lp)); /* Good formula for small lp. */
4022 }
4023
4024
4025 /* --- BEGIN IANDJMSMITH SOURCE MARKER --- */
4026
4027 /* Calculation of logfbit(x)-logfbit(1+x). y2 must be < 1. */
4028 static gnm_float
logfbitdif(gnm_float x)4029 logfbitdif (gnm_float x)
4030 {
4031 gnm_float y = 1 / (2 * x + 3);
4032 gnm_float y2 = y * y;
4033 return y2 * gnm_logcf (y2, 3, 2);
4034 }
4035
4036 /*
4037 * lfbc{1-7} from this Mathematica program:
4038 *
4039 * Table[Numerator[BernoulliB[2n]/(2n(2n - 1))], {n, 1, 22}]
4040 * Table[Denominator[BernoulliB[2n]/(2n(2n - 1))], {n, 1, 22}]
4041 */
4042 static const gnm_float lfbc1 = GNM_const (1.0) / 12;
4043 static const gnm_float lfbc2 = GNM_const (1.0) / 30;
4044 static const gnm_float lfbc3 = GNM_const (1.0) / 105;
4045 static const gnm_float lfbc4 = GNM_const (1.0) / 140;
4046 static const gnm_float lfbc5 = GNM_const (1.0) / 99;
4047 static const gnm_float lfbc6 = GNM_const (691.0) / 30030;
4048 static const gnm_float lfbc7 = GNM_const (1.0) / 13;
4049 /* lfbc{8,9} to make logfbit(6) and logfbit(7) exact. */
4050 static const gnm_float lfbc8 = GNM_const (3.5068606896459316479e-01);
4051 static const gnm_float lfbc9 = GNM_const (1.6769998201671114808);
4052
4053 /* This is also stirlerr(x+1). */
4054 static gnm_float
logfbit(gnm_float x)4055 logfbit (gnm_float x)
4056 {
4057 /*
4058 * Error part of Stirling's formula where
4059 * log(x!) = log(sqrt(twopi))+(x+0.5)*log(x+1)-(x+1)+logfbit(x).
4060 *
4061 * Are we ever concerned about the relative error involved in this
4062 * function? I don't think so.
4063 */
4064 if (x >= 1e10) return 1 / (12 * (x + 1));
4065 else if (x >= 6) {
4066 gnm_float x1 = x + 1;
4067 gnm_float x2 = 1 / (x1 * x1);
4068 gnm_float x3 =
4069 x2 * (lfbc2 - x2 *
4070 (lfbc3 - x2 *
4071 (lfbc4 - x2 *
4072 (lfbc5 - x2 *
4073 (lfbc6 - x2 *
4074 (lfbc7 - x2 *
4075 (lfbc8 - x2 *
4076 lfbc9)))))));
4077 return lfbc1 * (1 - x3) / x1;
4078 }
4079 else if (x == 5) return GNM_const (0.13876128823070747998745727023762908562e-1);
4080 else if (x == 4) return GNM_const (0.16644691189821192163194865373593391145e-1);
4081 else if (x == 3) return GNM_const (0.20790672103765093111522771767848656333e-1);
4082 else if (x == 2) return GNM_const (0.27677925684998339148789292746244666596e-1);
4083 else if (x == 1) return GNM_const (0.41340695955409294093822081407117508025e-1);
4084 else if (x == 0) return GNM_const (0.81061466795327258219670263594382360138e-1);
4085 else if (x > -1) {
4086 gnm_float x1 = x;
4087 gnm_float x2 = 0;
4088 while (x1 < 6) {
4089 x2 += logfbitdif (x1);
4090 x1++;
4091 }
4092 return x2 + logfbit (x1);
4093 }
4094 else return gnm_pinf;
4095 }
4096
4097 /* Calculation of logfbit1(x)-logfbit1(1+x). */
4098 static gnm_float
logfbit1dif(gnm_float x)4099 logfbit1dif (gnm_float x)
4100 {
4101 return (logfbitdif (x) - 1 / (4 * (x + 1) * (x + 2))) / (x + 1.5);
4102 }
4103
4104 /* Derivative logfbit. */
4105 static gnm_float
logfbit1(gnm_float x)4106 logfbit1 (gnm_float x)
4107 {
4108 if (x >= 1e10) return -lfbc1 * gnm_pow (x + 1, -2);
4109 else if (x >= 6) {
4110 gnm_float x1 = x + 1;
4111 gnm_float x2 = 1 / (x1 * x1);
4112 gnm_float x3 =
4113 x2 * (3 * lfbc2 - x2*
4114 (5 * lfbc3 - x2 *
4115 (7 * lfbc4 - x2 *
4116 (9 * lfbc5 - x2 *
4117 (11 * lfbc6 - x2 *
4118 (13 * lfbc7 - x2 *
4119 (15 * lfbc8 - x2 *
4120 17 * lfbc9)))))));
4121 return -lfbc1 * (1 - x3) * x2;
4122 }
4123 else if (x > -1) {
4124 gnm_float x1 = x;
4125 gnm_float x2 = 0;
4126 while (x1 < 6) {
4127 x2 += logfbit1dif (x1);
4128 x1++;
4129 }
4130 return x2 + logfbit1 (x1);
4131 }
4132 else return gnm_ninf;
4133 }
4134
4135
4136 /* Calculation of logfbit3(x)-logfbit3(1+x). */
4137 static gnm_float
logfbit3dif(gnm_float x)4138 logfbit3dif (gnm_float x)
4139 {
4140 return -(2 * x + 3) * gnm_pow ((x + 1) * (x + 2), -3);
4141 }
4142
4143
4144 /* Third derivative logfbit. */
4145 static gnm_float
logfbit3(gnm_float x)4146 logfbit3 (gnm_float x)
4147 {
4148 if (x >= 1e10) return -0.5 * gnm_pow (x + 1, -4);
4149 else if (x >= 6) {
4150 gnm_float x1 = x + 1;
4151 gnm_float x2 = 1 / (x1 * x1);
4152 gnm_float x3 =
4153 x2 * (60 * lfbc2 - x2 *
4154 (210 * lfbc3 - x2 *
4155 (504 * lfbc4 - x2 *
4156 (990 * lfbc5 - x2 *
4157 (1716 * lfbc6 - x2 *
4158 (2730 * lfbc7 - x2 *
4159 (4080 * lfbc8 - x2 *
4160 5814 * lfbc9)))))));
4161 return -lfbc1 * (6 - x3) * x2 * x2;
4162 }
4163 else if (x > -1) {
4164 gnm_float x1 = x;
4165 gnm_float x2 = 0;
4166 while (x1 < 6) {
4167 x2 += logfbit3dif (x1);
4168 x1++;
4169 }
4170 return x2 + logfbit3 (x1);
4171 }
4172 else return gnm_ninf;
4173 }
4174
4175 /* Calculation of logfbit5(x)-logfbit5(1+x). */
4176 static gnm_float
logfbit5dif(gnm_float x)4177 logfbit5dif (gnm_float x)
4178 {
4179 return -6 * (2 * x + 3) * ((5 * x + 15) * x + 12) *
4180 gnm_pow ((x + 1) * (x + 2), -5);
4181 }
4182
4183 /* Fifth derivative logfbit. */
4184 static gnm_float
logfbit5(gnm_float x)4185 logfbit5 (gnm_float x)
4186 {
4187 if (x >= 1e10) return -10 * gnm_pow (x + 1, -6);
4188 else if (x >= 6) {
4189 gnm_float x1 = x + 1;
4190 gnm_float x2 = 1 / (x1 * x1);
4191 gnm_float x3 =
4192 x2 * (2520 * lfbc2 - x2 *
4193 (15120 * lfbc3 - x2 *
4194 (55440 * lfbc4 - x2 *
4195 (154440 * lfbc5 - x2 *
4196 (360360 * lfbc6 - x2 *
4197 (742560 * lfbc7 - x2 *
4198 (1395360 * lfbc8 - x2 *
4199 2441880 * lfbc9)))))));
4200 return -lfbc1 * (120 - x3) * x2 * x2 * x2;
4201 } else if (x > -1) {
4202 gnm_float x1 = x;
4203 gnm_float x2 = 0;
4204 while (x1 < 6) {
4205 x2 += logfbit5dif (x1);
4206 x1++;
4207 }
4208 return x2 + logfbit5 (x1);
4209 }
4210 else return gnm_ninf;
4211 }
4212
4213 /* Calculation of logfbit7(x)-logfbit7(1+x). */
4214 static gnm_float
logfbit7dif(gnm_float x)4215 logfbit7dif (gnm_float x)
4216 {
4217 return -120 *
4218 (2 * x + 3) *
4219 ((((14 * x + 84) * x + 196) * x + 210) * x + 87) *
4220 gnm_pow ((x + 1) * (x + 2), -7);
4221 }
4222
4223 /* Seventh derivative logfbit. */
4224 static gnm_float
logfbit7(gnm_float x)4225 logfbit7 (gnm_float x)
4226 {
4227 if (x >= 1e10) return -420 * gnm_pow (x + 1, -8);
4228 else if (x >= 6) {
4229 gnm_float x1 = x + 1;
4230 gnm_float x2 = 1 / (x1 * x1);
4231 gnm_float x3 =
4232 x2 * (181440 * lfbc2 - x2 *
4233 (1663200 * lfbc3 - x2 *
4234 (8648640 * lfbc4 - x2 *
4235 (32432400 * lfbc5 - x2 *
4236 (98017920 * lfbc6 - x2 *
4237 (253955520 * lfbc7 - x2 *
4238 (586051200 * lfbc8 - x2 *
4239 1235591280 * lfbc9)))))));
4240 return -lfbc1 * (5040 - x3) * x2 * x2 * x2 * x2;
4241 } else if (x > -1) {
4242 gnm_float x1 = x;
4243 gnm_float x2 = 0;
4244 while (x1 < 6) {
4245 x2 += logfbit7dif (x1);
4246 x1++;
4247 }
4248 return x2 + logfbit7 (x1);
4249 }
4250 else return gnm_ninf;
4251 }
4252
4253
4254 static gnm_float
lfbaccdif(gnm_float a,gnm_float b)4255 lfbaccdif (gnm_float a, gnm_float b)
4256 {
4257 if (a > 0.03 * (a + b))
4258 return logfbit (a + b) - logfbit (b);
4259 else {
4260 gnm_float a2 = a * a;
4261 gnm_float ab2 = a / 2 + b;
4262 return a * (logfbit1 (ab2) + a2 / 24 *
4263 (logfbit3 (ab2) + a2 / 80 *
4264 (logfbit5 (ab2) + a2 / 168 *
4265 logfbit7 (ab2))));
4266 }
4267 }
4268
4269 static gnm_float
compbfunc(gnm_float x,gnm_float a,gnm_float b)4270 compbfunc (gnm_float x, gnm_float a, gnm_float b)
4271 {
4272 const gnm_float sumAcc = 5E-16;
4273 gnm_float term = x;
4274 gnm_float d = 2;
4275 gnm_float sum = term / (a + 1);
4276 while (gnm_abs (term) > gnm_abs (sum * sumAcc)) {
4277 term *= (d - b) * x / d;
4278 sum += term / (a + d);
4279 d++;
4280 }
4281 return a * (b - 1) * sum;
4282 }
4283
4284 static gnm_float
pbeta_smalla(gnm_float x,gnm_float a,gnm_float b,gboolean lower_tail,gboolean log_p)4285 pbeta_smalla (gnm_float x, gnm_float a, gnm_float b, gboolean lower_tail, gboolean log_p)
4286 {
4287 gnm_float r;
4288
4289 #if 0
4290 assert (a >= 0 && b >= 0);
4291 assert (a < 1);
4292 assert (b < 1 || (1 + b) * x <= 1);
4293 #endif
4294
4295 if (x > 0.5) {
4296 gnm_float olda = a;
4297 a = b;
4298 b = olda;
4299 x = 1 - x;
4300 lower_tail = !lower_tail;
4301 }
4302
4303 r = (a + b + 0.5) * log1pmx (a / (1 + b)) +
4304 a * (a - 0.5) / (1 + b) +
4305 lfbaccdif (a, b);
4306 r += a * gnm_log ((1 + b) * x) - lgamma1p (a);
4307 if (lower_tail) {
4308 if (log_p)
4309 return r + gnm_log1p (-compbfunc (x, a, b)) + gnm_log (b / (a + b));
4310 else
4311 return gnm_exp (r) * (1 - compbfunc (x, a, b)) * (b / (a + b));
4312 } else {
4313 /* x=0.500000001 a=0.5 b=0.000001 ends up here [swapped]
4314 * with r=-7.94418987455065e-08 and cbf=-3.16694087508444e-07.
4315 *
4316 * x=0.0000001 a=0.999999 b=0.02 end up here with
4317 * r=-16.098276918385 and cbf=-4.89999787339858e-08.
4318 */
4319 if (log_p) {
4320 return swap_log_tail (r + gnm_log1p (-compbfunc (x, a, b)) + gnm_log (b / (a + b)));
4321 } else {
4322 r = -gnm_expm1 (r);
4323 r += compbfunc (x, a, b) * (1 - r);
4324 r += (a / (a + b)) * (1 - r);
4325 return r;
4326 }
4327 }
4328 }
4329
4330 /* Cumulative Students t-distribution, with odd parameterisation for
4331 * use with binApprox.
4332 * p is x*x/(k+x*x)
4333 * q is 1-p
4334 * logqk2 is LN(q)*k/2
4335 * approxtdistDens returns with approx density function, for use in
4336 * binApprox
4337 */
4338 static gnm_float
tdistexp(gnm_float p,gnm_float q,gnm_float logqk2,gnm_float k,gboolean log_p,gnm_float * approxtdistDens)4339 tdistexp (gnm_float p, gnm_float q, gnm_float logqk2, gnm_float k,
4340 gboolean log_p, gnm_float *approxtdistDens)
4341 {
4342 const gnm_float sumAcc = 5E-16;
4343 const gnm_float cfVSmall = 1.0e-14;
4344 const gnm_float lstpi = gnm_log (2 * M_PIgnum) / 2;
4345
4346 if (gnm_floor (k / 2) * 2 == k) {
4347 gnm_float ldens = logqk2 + logfbit (k - 1) - 2 * logfbit (k * 0.5 - 1) - lstpi;
4348 *approxtdistDens = R_D_exp (ldens);
4349 } else {
4350 gnm_float ldens = logqk2 + k * log1pmx (1 / k) + 2 * logfbit ((k - 1) * 0.5) - logfbit (k - 1) - lstpi;
4351 *approxtdistDens = R_D_exp (ldens);
4352 }
4353
4354 if (k * p < 4 * q) {
4355 gnm_float sum = 0;
4356 gnm_float aki = k + 1;
4357 gnm_float ai = 3;
4358 gnm_float term = aki * p / ai;
4359
4360 while (term > sumAcc * sum) {
4361 sum += term;
4362 ai += 2;
4363 aki += 2;
4364 term *= aki * p / ai;
4365 }
4366 sum += term;
4367
4368 return log_p
4369 ? logspace_sub (-M_LN2gnum, *approxtdistDens + gnm_log1p (sum) + gnm_log (k * p) / 2)
4370 : 0.5 - *approxtdistDens * (sum + 1) * gnm_sqrt (k * p);
4371 } else {
4372 gnm_float q1 = 2 * (1 + q);
4373 gnm_float q8 = 8 * q;
4374 gnm_float a1 = 0;
4375 gnm_float b1 = 1;
4376 gnm_float c1 = -6 * q;
4377 gnm_float a2 = 1;
4378 gnm_float b2 = (k - 1) * p + 3;
4379 gnm_float cadd = -14 * q;
4380 gnm_float c2 = b2 + q1;
4381
4382 while (gnm_abs (a2 * b1 - a1 * b2) > gnm_abs (cfVSmall * b1 * b2)) {
4383 a1 = c2 * a2 + c1 * a1;
4384 b1 = c2 * b2 + c1 * b1;
4385 c1 += cadd;
4386 cadd -= q8;
4387 c2 += q1;
4388 a2 = c2 * a1 + c1 * a2;
4389 b2 = c2 * b1 + c1 * b2;
4390 c1 += cadd;
4391 cadd -= q8;
4392 c2 += q1;
4393
4394 if (gnm_abs (b2) > scalefactor) {
4395 a1 *= 1 / scalefactor;
4396 b1 *= 1 / scalefactor;
4397 a2 *= 1 / scalefactor;
4398 b2 *= 1 / scalefactor;
4399 } else if (gnm_abs (b2) < 1 / scalefactor) {
4400 a1 *= scalefactor;
4401 b1 *= scalefactor;
4402 a2 *= scalefactor;
4403 b2 *= scalefactor;
4404 }
4405 }
4406
4407 return log_p
4408 ? *approxtdistDens + gnm_log1p (-q * a2 / b2) - gnm_log (k * p) / 2
4409 : *approxtdistDens * (1 - q * a2 / b2) / gnm_sqrt (k * p);
4410 }
4411 }
4412
4413
4414 /* Asymptotic expansion to calculate the probability that binomial variate
4415 * has value <= a.
4416 * (diffFromMean = (a+b)*p-a).
4417 */
4418 static gnm_float
binApprox(gnm_float a,gnm_float b,gnm_float diffFromMean,gboolean lower_tail,gboolean log_p)4419 binApprox (gnm_float a, gnm_float b, gnm_float diffFromMean,
4420 gboolean lower_tail, gboolean log_p)
4421 {
4422 gnm_float pq1, res, t;
4423 gnm_float ib05, ib15, ib25, ib35, ib3;
4424 gnm_float elfb, coef15, coef25, coef35;
4425 gnm_float approxtdistDens;
4426
4427 gnm_float n = a + b;
4428 gnm_float n1 = n + 1;
4429 gnm_float lvv = b - n * diffFromMean;
4430 gnm_float lval = (a * log1pmx (lvv / (a * n1)) +
4431 b * log1pmx (-lvv / (b * n1))) / n;
4432 gnm_float tp = -gnm_expm1 (lval);
4433 gnm_float mfac = n1 * tp;
4434 gnm_float ib2 = 1 + mfac;
4435 gnm_float t1 = (n + 2) * tp;
4436
4437 mfac = 2 * mfac;
4438
4439 ib3 = ib2 + mfac*t1;
4440 ib05 = tdistexp (tp, 1 - tp, n1 * lval, 2 * n1, log_p, &approxtdistDens);
4441
4442 ib15 = gnm_sqrt (mfac);
4443 mfac = t1 * (GNM_const (2.0) / 3);
4444 ib25 = 1 + mfac;
4445 ib35 = ib25 + mfac * (n + 3) * tp * (GNM_const (2.0) / 5);
4446
4447 pq1 = (n * n) / (a * b);
4448
4449 res = (ib2 * (1 + 2 * pq1) / 135 - 2 * ib3 * ((2 * pq1 - 43) * pq1 - 22) / (2835 * (n + 3))) / (n + 2);
4450 res = (GNM_const (1.0) / 3 - res) * 2 * gnm_sqrt (pq1 / n1) * (a - b) / n;
4451 if (lvv > 0) {
4452 res = -res;
4453 lower_tail = !lower_tail;
4454 }
4455
4456 n1 = (n + 1.5) * (n + 2.5);
4457 coef15 = (-17 + 2 * pq1) / (24 * (n + 1.5));
4458 coef25 = (-503 + 4 * pq1 * (19 + pq1)) / (1152 * n1);
4459 coef35 = (-315733 + pq1 * (53310 + pq1 * (8196 - 1112 * pq1))) /
4460 (414720 * n1 * (n + 3.5));
4461 elfb = (coef35 + coef25) + coef15;
4462 res += ib15 * ((coef35 * ib35 + coef25 * ib25) + coef15);
4463
4464 t = log_p
4465 ? logspace_signed_add (ib05,
4466 gnm_log (gnm_abs (res)) + approxtdistDens - gnm_log1p (elfb),
4467 res >= 0)
4468 : ib05 + res * approxtdistDens / (1 + elfb);
4469
4470 return lower_tail
4471 ? t
4472 : log_p ? swap_log_tail (t) : (1 - t);
4473 }
4474
4475 /* Probability that binomial variate with sample size i+j and
4476 * event prob p (=1-q) has value i (diffFromMean = (i+j)*p-i)
4477 */
4478 static gnm_float
binomialTerm(gnm_float i,gnm_float j,gnm_float p,gnm_float q,gnm_float diffFromMean,gboolean log_p)4479 binomialTerm (gnm_float i, gnm_float j, gnm_float p, gnm_float q,
4480 gnm_float diffFromMean, gboolean log_p)
4481 {
4482 const gnm_float minLog1Value = -0.79149064;
4483 gnm_float c1,c2,c3;
4484 gnm_float c4,c5,c6,ps,logbinomialTerm,dfm;
4485 gnm_float t;
4486
4487 if (i == 0 && j <= 0)
4488 return R_D__1;
4489
4490 if (i <= -1 || j < 0)
4491 return R_D__0;
4492
4493 c1 = (i + 1) + j;
4494 if (p < q) {
4495 c2 = i;
4496 c3 = j;
4497 ps = p;
4498 dfm = diffFromMean;
4499 } else {
4500 c3 = i;
4501 c2 = j;
4502 ps = q;
4503 dfm = -diffFromMean;
4504 }
4505
4506 c5 = (dfm - (1 - ps)) / (c2 + 1);
4507 c6 = -(dfm + ps) / (c3 + 1);
4508
4509 if (c5 < minLog1Value) {
4510 if (c2 == 0) {
4511 logbinomialTerm = c3 * gnm_log1p (-ps);
4512 return log_p ? logbinomialTerm : gnm_exp (logbinomialTerm);
4513 } else if (ps == 0 && c2 > 0) {
4514 return R_D__0;
4515 } else {
4516 t = gnm_log ((ps * c1) / (c2 + 1)) - c5;
4517 }
4518 } else {
4519 t = log1pmx (c5);
4520 }
4521
4522 c4 = logfbit (i + j) - logfbit (i) - logfbit (j);
4523 logbinomialTerm = c4 + c2 * t - c5 + (c3 * log1pmx (c6) - c6);
4524
4525 return log_p
4526 ? logbinomialTerm + 0.5 * gnm_log (c1 / ((c2 + 1) * (c3 + 1) * 2 * M_PIgnum))
4527 : gnm_exp (logbinomialTerm) * gnm_sqrt (c1 / ((c2 + 1) * (c3 + 1) * 2 * M_PIgnum));
4528 }
4529
4530
4531 /*
4532 * Probability that binomial variate with sample size ii+jj
4533 * and event prob pp (=1-qq) has value <=i.
4534 * (diffFromMean = (ii+jj)*pp-ii).
4535 */
4536 static gnm_float
binomialcf(gnm_float ii,gnm_float jj,gnm_float pp,gnm_float qq,gnm_float diffFromMean,gboolean lower_tail,gboolean log_p)4537 binomialcf (gnm_float ii, gnm_float jj, gnm_float pp, gnm_float qq,
4538 gnm_float diffFromMean, gboolean lower_tail, gboolean log_p)
4539 {
4540 const gnm_float sumAlways = 0;
4541 const gnm_float sumFactor = 6;
4542 const gnm_float cfSmall = 1.0e-15;
4543
4544 gnm_float prob,p,q,a1,a2,b1,b2,c1,c2,c3,c4,n1,q1,dfm;
4545 gnm_float i,j,ni,nj,numb,ip1;
4546 gboolean swapped;
4547
4548 ip1 = ii + 1;
4549 if (ii > -1 && (jj <= 0 || pp == 0)) {
4550 return R_DT_1;
4551 } else if (ii > -1 && ii < 0) {
4552 gnm_float f;
4553 ii = -ii;
4554 ip1 = ii;
4555 f = ii / ((ii + jj) * pp);
4556 prob = binomialTerm (ii, jj, pp, qq, (ii + jj) * pp - ii, log_p);
4557 prob = log_p ? prob + gnm_log (f) : prob * f;
4558 ii--;
4559 diffFromMean = (ii + jj) * pp - ii;
4560 } else
4561 prob = binomialTerm (ii, jj, pp, qq, diffFromMean, log_p);
4562
4563 n1 = (ii + 3) + jj;
4564 if (ii < 0)
4565 swapped = FALSE;
4566 else if (pp > qq)
4567 swapped = (n1 * qq >= jj + 1);
4568 else
4569 swapped = (n1 * pp <= ii + 2);
4570
4571 if (prob == R_D__0) {
4572 if (swapped == !lower_tail)
4573 return R_D__0;
4574 else
4575 return R_D__1;
4576 }
4577
4578 if (swapped) {
4579 j = ip1;
4580 ip1 = jj;
4581 i = jj - 1;
4582 p = qq;
4583 q = pp;
4584 dfm = 1 - diffFromMean;
4585 } else {
4586 i = ii;
4587 j = jj;
4588 p = pp;
4589 q = qq;
4590 dfm = diffFromMean;
4591 }
4592
4593 if (i > sumAlways) {
4594 numb = gnm_floor (sumFactor * gnm_sqrt (p + 0.5) * gnm_exp (gnm_log (n1 * p * q) / 3));
4595 numb = gnm_floor (numb - dfm);
4596 if (numb > i) numb = gnm_floor (i);
4597 } else
4598 numb = gnm_floor (i);
4599 if (numb < 0) numb = 0;
4600
4601 a1 = 0;
4602 b1 = 1;
4603 q1 = q + 1;
4604 a2 = (i - numb) * q;
4605 b2 = dfm + numb + 1;
4606 c1 = 0;
4607
4608 c2 = a2;
4609 c4 = b2;
4610
4611 while (gnm_abs (a2 * b1 - a1 * b2) > gnm_abs (cfSmall * b1 * b2)) {
4612 c1++;
4613 c2 -= q;
4614 c3 = c1 * c2;
4615 c4 += q1;
4616 a1 = c4 * a2 + c3 * a1;
4617 b1 = c4 * b2 + c3 * b1;
4618 c1++;
4619 c2 -= q;
4620 c3 = c1 * c2;
4621 c4 += q1;
4622 a2 = c4 * a1 + c3 * a2;
4623 b2 = c4 * b1 + c3 * b2;
4624
4625 if (gnm_abs (b2) > scalefactor) {
4626 a1 *= 1 / scalefactor;
4627 b1 *= 1 / scalefactor;
4628 a2 *= 1 / scalefactor;
4629 b2 *= 1 / scalefactor;
4630 } else if (gnm_abs (b2) < 1 / scalefactor) {
4631 a1 *= scalefactor;
4632 b1 *= scalefactor;
4633 a2 *= scalefactor;
4634 b2 *= scalefactor;
4635 }
4636 }
4637
4638 a1 = a2 / b2;
4639
4640 ni = (i - numb + 1) * q;
4641 nj = (j + numb) * p;
4642 while (numb > 0) {
4643 a1 = (1 + a1) * (ni / nj);
4644 ni = ni + q;
4645 nj = nj - p;
4646 numb--;
4647 }
4648
4649 prob = log_p ? prob + gnm_log1p (a1) : prob * (1 + a1);
4650
4651 if (swapped) {
4652 if (log_p)
4653 prob += gnm_log (ip1 * q / nj);
4654 else
4655 prob *= ip1 * q / nj;
4656 }
4657
4658 if (swapped == !lower_tail)
4659 return prob;
4660 else
4661 return log_p ? swap_log_tail (prob) : 1 - prob;
4662 }
4663
4664 static gnm_float
binomial(gnm_float ii,gnm_float jj,gnm_float pp,gnm_float qq,gnm_float diffFromMean,gboolean lower_tail,gboolean log_p)4665 binomial (gnm_float ii, gnm_float jj, gnm_float pp, gnm_float qq,
4666 gnm_float diffFromMean, gboolean lower_tail, gboolean log_p)
4667 {
4668 gnm_float mij = fmin2 (ii, jj);
4669
4670 if (mij > 500 && gnm_abs (diffFromMean) < 0.01 * mij)
4671 return binApprox (jj - 1, ii, diffFromMean, lower_tail, log_p);
4672
4673 return binomialcf (ii, jj, pp, qq, diffFromMean, lower_tail, log_p);
4674 }
4675
4676 gnm_float
pbeta(gnm_float x,gnm_float a,gnm_float b,gboolean lower_tail,gboolean log_p)4677 pbeta (gnm_float x, gnm_float a, gnm_float b, gboolean lower_tail, gboolean log_p)
4678 {
4679 gnm_float am1;
4680
4681 if (gnm_isnan (x) || gnm_isnan (a) || gnm_isnan (b))
4682 return x + a + b;
4683
4684 if (x <= 0) return R_DT_0;
4685 if (x >= 1) return R_DT_1;
4686
4687 if (a < 1 && (b < 1 || (1 + b) * x <= 1))
4688 return pbeta_smalla (x, a, b, lower_tail, log_p);
4689
4690 if (b < 1 && (1 + a) * (1 - x) <= 1)
4691 return pbeta_smalla (1 - x, b, a, !lower_tail, log_p);
4692
4693 if (a < 1)
4694 return binomial (-a, b, x, 1 - x, 0, !lower_tail, log_p);
4695
4696 if (b < 1)
4697 return binomial (-b, a, 1 - x, x, 0, lower_tail, log_p);
4698
4699 am1 = a - 1;
4700 return binomial (am1, b, x, 1 - x, (am1 + b) * x - am1,
4701 !lower_tail, log_p);
4702 }
4703
4704 /* --- END IANDJMSMITH SOURCE MARKER --- */
4705 /* ------------------------------------------------------------------------ */
4706
4707 gnm_float
pf(gnm_float x,gnm_float n1,gnm_float n2,gboolean lower_tail,gboolean log_p)4708 pf (gnm_float x, gnm_float n1, gnm_float n2, gboolean lower_tail, gboolean log_p)
4709 {
4710 #ifdef IEEE_754
4711 if (gnm_isnan (x) || gnm_isnan (n1) || gnm_isnan (n2))
4712 return x + n2 + n1;
4713 #endif
4714 if (n1 <= 0 || n2 <= 0) ML_ERR_return_NAN;
4715
4716 if (x <= 0)
4717 return R_DT_0;
4718
4719 /* Avoid squeezing pbeta's first parameter against 1. */
4720 if (n1 * x > n2)
4721 return pbeta (n2 / (n2 + n1 * x), n2 / 2, n1 / 2,
4722 !lower_tail, log_p);
4723 else
4724 return pbeta (n1 * x / (n2 + n1 * x), n1 / 2, n2 / 2,
4725 lower_tail, log_p);
4726 }
4727
4728 /* ------------------------------------------------------------------------ */
4729
4730 static gnm_float
ppois1(gnm_float x,const gnm_float * plambda,gboolean lower_tail,gboolean log_p)4731 ppois1 (gnm_float x, const gnm_float *plambda,
4732 gboolean lower_tail, gboolean log_p)
4733 {
4734 return ppois (x, *plambda, lower_tail, log_p);
4735 }
4736
4737 gnm_float
qpois(gnm_float p,gnm_float lambda,gboolean lower_tail,gboolean log_p)4738 qpois (gnm_float p, gnm_float lambda, gboolean lower_tail, gboolean log_p)
4739 {
4740 gnm_float mu, sigma, gamma, y, z;
4741
4742 if (!(lambda >= 0))
4743 return gnm_nan;
4744
4745 mu = lambda;
4746 sigma = gnm_sqrt (lambda);
4747 gamma = 1 / sigma;
4748
4749 /* Cornish-Fisher expansion: */
4750 z = qnorm (p, 0., 1., lower_tail, log_p);
4751 y = mu + sigma * (z + gamma * (z * z - 1) / 6);
4752
4753 return discpfuncinverter (p, &lambda, lower_tail, log_p,
4754 0, gnm_pinf, y,
4755 ppois1);
4756 }
4757
4758 /* ------------------------------------------------------------------------ */
4759
4760 static gnm_float
dgamma1(gnm_float x,const gnm_float * palpha,gboolean log_p)4761 dgamma1 (gnm_float x, const gnm_float *palpha, gboolean log_p)
4762 {
4763 return dgamma (x, *palpha, 1, log_p);
4764 }
4765
4766 static gnm_float
pgamma1(gnm_float x,const gnm_float * palpha,gboolean lower_tail,gboolean log_p)4767 pgamma1 (gnm_float x, const gnm_float *palpha,
4768 gboolean lower_tail, gboolean log_p)
4769 {
4770 return pgamma (x, *palpha, 1, lower_tail, log_p);
4771 }
4772
4773
4774 gnm_float
qgamma(gnm_float p,gnm_float alpha,gnm_float scale,gboolean lower_tail,gboolean log_p)4775 qgamma (gnm_float p, gnm_float alpha, gnm_float scale,
4776 gboolean lower_tail, gboolean log_p)
4777 {
4778 gnm_float res1, x0, v;
4779
4780 #ifdef IEEE_754
4781 if (gnm_isnan(p) || gnm_isnan(alpha) || gnm_isnan(scale))
4782 return p + alpha + scale;
4783 #endif
4784 R_Q_P01_check(p);
4785 if (alpha <= 0) ML_ERR_return_NAN;
4786
4787 if (!log_p && p > 0.9) {
4788 /* We're far into the tail. Flip. */
4789 p = 1 - p;
4790 lower_tail = !lower_tail;
4791 }
4792
4793 /* Make an initial guess, x0, assuming scale==1. */
4794 v = 2 * alpha;
4795 if (v < -1.24 * R_DT_log (p))
4796 x0 = gnm_pow (R_DT_qIv (p) * alpha * gnm_exp (gnm_lgamma (alpha) + alpha * M_LN2gnum),
4797 1 / alpha) / 2;
4798 else {
4799 gnm_float x1 = qnorm (p, 0, 1, lower_tail, log_p);
4800 gnm_float p1 = 0.222222 / v;
4801 x0 = v * gnm_pow (x1 * gnm_sqrt (p1) + 1 - p1, 3) / 2;
4802 }
4803
4804 res1 = pfuncinverter (p, &alpha, lower_tail, log_p, 0, gnm_pinf, x0,
4805 pgamma1, dgamma1);
4806
4807 return res1 * scale;
4808 }
4809
4810 /* ------------------------------------------------------------------------ */
4811
4812 static gnm_float
dbeta1(gnm_float x,const gnm_float shape[],gboolean log_p)4813 dbeta1 (gnm_float x, const gnm_float shape[], gboolean log_p)
4814 {
4815 return dbeta (x, shape[0], shape[1], log_p);
4816 }
4817
4818 static gnm_float
pbeta1(gnm_float x,const gnm_float shape[],gboolean lower_tail,gboolean log_p)4819 pbeta1 (gnm_float x, const gnm_float shape[],
4820 gboolean lower_tail, gboolean log_p)
4821 {
4822 return pbeta (x, shape[0], shape[1], lower_tail, log_p);
4823 }
4824
4825 static gnm_float
abramowitz_stegun_26_5_22(gnm_float p,gnm_float a,gnm_float b,gboolean lower_tail,gboolean log_p)4826 abramowitz_stegun_26_5_22 (gnm_float p, gnm_float a, gnm_float b,
4827 gboolean lower_tail, gboolean log_p)
4828 {
4829 gnm_float yp = qnorm (p, 0, 1, !lower_tail, log_p);
4830 gnm_float ta = 1 / (2 * a - 1);
4831 gnm_float tb = 1 / (2 * b - 1);
4832 gnm_float h = 2 / (ta + tb);
4833 gnm_float l = (yp * yp - 3) / 6;
4834 gnm_float w = yp * gnm_sqrt (h + l) / h -
4835 (tb - ta) * (l + (5 - 4 / h) / 6);
4836 return a / (a + b * gnm_exp (2 * w));
4837 }
4838
4839
4840 gnm_float
qbeta(gnm_float p,gnm_float pin,gnm_float qin,gboolean lower_tail,gboolean log_p)4841 qbeta (gnm_float p, gnm_float pin, gnm_float qin, gboolean lower_tail, gboolean log_p)
4842 {
4843 gnm_float x0, q, shape[2];
4844 gnm_float S = pin + qin;
4845
4846 #ifdef IEEE_754
4847 if (gnm_isnan (S) || gnm_isnan (p))
4848 return S + p;
4849 #endif
4850 R_Q_P01_check (p);
4851
4852 if (pin < 0. || qin < 0.) ML_ERR_return_NAN;
4853
4854 if (!log_p && p > 0.9) {
4855 /* We're far into the tail. Flip. */
4856 p = 1 - p;
4857 lower_tail = !lower_tail;
4858 }
4859
4860 if (pin >= 1 && qin >= 1)
4861 x0 = abramowitz_stegun_26_5_22 (p, pin, qin, lower_tail, log_p);
4862 else {
4863 /*
4864 * The density function is U-shaped.
4865 */
4866 gnm_float phalf = pbeta (0.5, pin, qin, lower_tail, log_p);
4867 gnm_float lb = gnm_lbeta (pin, qin);
4868
4869 if (!lower_tail == (p > phalf)) {
4870 /*
4871 * The following approximation follows from simply ignoring
4872 * the (1-t)^(qin-1) factor from the density. That works
4873 * fine as long as we stay far away from the right tail.
4874 */
4875 x0 = gnm_exp ((gnm_log (pin) + R_DT_log (p) + lb) / pin);
4876 } else {
4877 /* Mirror. */
4878 x0 = -gnm_expm1 ((gnm_log (qin) + R_DT_Clog (p) + lb) / qin);
4879 }
4880 }
4881
4882 shape[0] = pin;
4883 shape[1] = qin;
4884
4885 q = pfuncinverter (p, shape, lower_tail, log_p, 0, 1, x0,
4886 pbeta1, dbeta1);
4887 return q;
4888 }
4889
4890 gnm_float
qf(gnm_float p,gnm_float n1,gnm_float n2,gboolean lower_tail,gboolean log_p)4891 qf (gnm_float p, gnm_float n1, gnm_float n2, gboolean lower_tail, gboolean log_p)
4892 {
4893 gnm_float q, qc;
4894 #ifdef IEEE_754
4895 if (gnm_isnan(p) || gnm_isnan(n1) || gnm_isnan(n2))
4896 return p + n1 + n2;
4897 #endif
4898 if (n1 <= 0. || n2 <= 0.) ML_ERR_return_NAN;
4899
4900 R_Q_P01_check(p);
4901 if (p == R_DT_0)
4902 return 0;
4903
4904 q = qbeta (p, n2 / 2, n1 / 2, !lower_tail, log_p);
4905 if (q < 0.9)
4906 qc = 1 - q;
4907 else
4908 qc = qbeta (p, n1 / 2, n2 / 2, lower_tail, log_p);
4909
4910 return (qc / q) * (n2 / n1);
4911 }
4912
4913
4914 gnm_float
pbinom2(gnm_float x0,gnm_float x1,gnm_float n,gnm_float p)4915 pbinom2 (gnm_float x0, gnm_float x1, gnm_float n, gnm_float p)
4916 {
4917 gnm_float Pl;
4918
4919 if (x0 > n || x1 < 0 || x0 > x1)
4920 return 0;
4921
4922 if (x0 == x1)
4923 return dbinom (x0, n, p, FALSE);
4924
4925 if (x0 <= 0)
4926 return pbinom (x1, n, p, TRUE, FALSE);
4927
4928 if (x1 >= n)
4929 return pbinom (x0 - 1, n, p, FALSE, FALSE);
4930
4931 Pl = pbinom (x0 - 1, n, p, TRUE, FALSE);
4932 if (Pl > 0.75) {
4933 gnm_float PlC = pbinom (x0 - 1, n, p, FALSE, FALSE);
4934 gnm_float PrC = pbinom (x1, n, p, FALSE, FALSE);
4935 return PlC - PrC;
4936 } else {
4937 gnm_float Pr = pbinom (x1, n, p, TRUE, FALSE);
4938 return Pr - Pl;
4939 }
4940 }
4941
4942 /* ------------------------------------------------------------------------- */
4943 /**
4944 * expmx2h:
4945 * @x: a number
4946 *
4947 * Returns: The result of exp(-0.5*@x*@x) with higher accuracy than the
4948 * naive formula.
4949 */
4950 gnm_float
expmx2h(gnm_float x)4951 expmx2h (gnm_float x)
4952 {
4953 x = gnm_abs (x);
4954
4955 if (x < 5 || gnm_isnan (x))
4956 return gnm_exp (-0.5 * x * x);
4957 else if (x >= GNM_MAX_EXP * M_LN2gnum + 10)
4958 return 0;
4959 else {
4960 /*
4961 * Split x into two parts, x=x1+x2, such that |x2|<=2^-16.
4962 * Assuming that we are using IEEE doubles, that means that
4963 * x1*x1 is error free for x<1024 (above which we will underflow
4964 * anyway). If we are not using IEEE doubles then this is
4965 * still an improvement over the naive formula.
4966 */
4967 gnm_float x1 = gnm_floor (x * 65536 + 0.5) / 65536;
4968 gnm_float x2 = x - x1;
4969 return (gnm_exp (-0.5 * x1 * x1) *
4970 gnm_exp ((-0.5 * x2 - x1) * x2));
4971 }
4972 }
4973
4974 /* ------------------------------------------------------------------------- */
4975
4976 /**
4977 * gnm_agm:
4978 * @a: a number
4979 * @b: a number
4980 *
4981 * Returns: The arithmetic-geometric mean of @a and @b.
4982 */
4983 gnm_float
gnm_agm(gnm_float a,gnm_float b)4984 gnm_agm (gnm_float a, gnm_float b)
4985 {
4986 gnm_float ab = a * b;
4987 gnm_float scale = 1;
4988 int i;
4989
4990 if (a < 0 || b < 0 || gnm_isnan (ab))
4991 return gnm_nan;
4992
4993 if (a == b)
4994 return a;
4995
4996 if (ab == 0 || ab == gnm_pinf) {
4997 int ea, eb;
4998
4999 if (a == 0 || b == 0)
5000 return 0;
5001
5002 // Underflow or overflow
5003 (void)gnm_frexp (a, &ea);
5004 (void)gnm_frexp (b, &eb);
5005 scale = gnm_ldexp (1, -(ea + eb) / 2);
5006 a *= scale;
5007 b *= scale;
5008 }
5009
5010 for (i = 1; i < 20; i++) {
5011 gnm_float am = (a + b) / 2;
5012 gnm_float gm = gnm_sqrt (a * b);
5013 a = am;
5014 b = gm;
5015 if (gnm_abs (a - b) < a * GNM_EPSILON)
5016 break;
5017 }
5018 if (i == 20)
5019 g_warning ("AGM failed to converge.");
5020
5021 return a / scale;
5022 }
5023
5024 /**
5025 * gnm_lambert_w:
5026 * @x: a number
5027 * @k: branch, either 0 or -1
5028 *
5029 * Returns: The Lambert W function at @x. @k is a branch number: 0
5030 * for the primary branch and -1 for the alternate.
5031 */
5032 gnm_float
gnm_lambert_w(gnm_float x,int k)5033 gnm_lambert_w (gnm_float x, int k)
5034 {
5035 gnm_float w;
5036 static const gnm_float one_over_e = 1 / M_Egnum;
5037 const gnm_float sqrt_one_over_e = gnm_sqrt (1 / M_Egnum);
5038 static const gboolean debug = FALSE;
5039 gnm_float wmin, wmax;
5040 int i, imax = 20;
5041
5042 if (gnm_isnan (x) || x < -one_over_e)
5043 return gnm_nan;
5044 else if (x == -one_over_e) {
5045 // This is technically wrong. The mathematically right
5046 // value of 1/e is 0.367879441171442321595524...
5047 // which rounds to 0.367879441171442334024277... as a double.
5048 // Note, that the rounding went up. In other words, when we
5049 // observe "x == -one_over_e" then x is already less that the
5050 // mathematically right value of -1/e and we ought to return
5051 // NaN. That is, however, a somewhat useless behaviour so
5052 // we return -1 instead.
5053 //
5054 // NOTE 1: The analysis might be different for long double.
5055 // NOTE 2: The analysis assumes that 1/e when computed as double
5056 // division (based on rounded e) produces the right
5057 // value. It does, see goffice's test/constants
5058 return -1;
5059 }
5060
5061 if (k == 0) {
5062 if (x == gnm_pinf)
5063 return gnm_pinf;
5064 if (x < 0)
5065 w = 1.5 * (gnm_sqrt (x + one_over_e) - sqrt_one_over_e);
5066 else if (x < 10)
5067 w = gnm_sqrt (x) / 1.7;
5068 else {
5069 gnm_float l1 = gnm_log (x);
5070 gnm_float l2 = gnm_log (l1);
5071 w = l1 - l2;
5072 }
5073 wmin = -1;
5074 wmax = gnm_pinf;
5075 } else if (k == -1) {
5076 if (x >= 0)
5077 return (x == 0) ? gnm_ninf : gnm_nan;
5078 if (x < -0.1)
5079 w = -1 - 3 * gnm_sqrt (x + one_over_e);
5080 else {
5081 gnm_float l1 = gnm_log (-x);
5082 gnm_float l2 = gnm_log (-l1);
5083 w = l1 - l2;
5084 }
5085
5086 wmin = gnm_ninf;
5087 wmax = -1;
5088 } else
5089 return gnm_nan;
5090
5091 if (debug) g_printerr ("x = %.20g w=%.20g\n", x, w);
5092 for (i = 0; i < imax; i++) {
5093 gnm_float ew = gnm_exp (w);
5094 gnm_float wew = w * ew;
5095 gnm_float d1 = ew * (w + 1);
5096 gnm_float d2 = ew * (w + 2);
5097 gnm_float dw;
5098 gnm_float wold = w;
5099
5100 dw = (-2 * ((wew - x) * d1) / (2 * d1 * d1 - (wew - x) * d2));
5101 w += dw;
5102
5103 if (w <= wmin || w >= wmax) {
5104 // We overshot
5105 gnm_float l = (w < wmin ? wmin : wmax);
5106 g_printerr (" (%2d w = %.20g)\n", i, w);
5107 dw = (l - wold) * 15 / 16;
5108 w = wold + dw;
5109 }
5110
5111 if (debug) {
5112 g_printerr (" %2d w = %.20g\n", i, w);
5113 if (i == imax - 1) {
5114 g_printerr (" wew = %.20g\n", wew);
5115 g_printerr (" d1 = %.20g\n", d1);
5116 g_printerr (" d2 = %.20g\n", d2);
5117 }
5118 }
5119
5120 if (gnm_abs (dw) <= 2 * GNM_EPSILON * gnm_abs (w))
5121 break;
5122 }
5123
5124 return w;
5125 }
5126
5127
5128 /**
5129 * pow1p:
5130 * @x: a number
5131 * @y: a number
5132 *
5133 * Returns: The result of (1+@x)^@y with less rounding error than the
5134 * naive formula.
5135 */
5136 gnm_float
pow1p(gnm_float x,gnm_float y)5137 pow1p (gnm_float x, gnm_float y)
5138 {
5139 /*
5140 * We defer to the naive algorithm in two cases: (1) when 1+x
5141 * is exact (let us hope the compiler does not mess this up),
5142 * and (2) when |x|>1/2 and we have no better algorithm.
5143 */
5144
5145 if ((x + 1) - x == 1 || gnm_abs (x) > 0.5 ||
5146 gnm_isnan (x) || gnm_isnan (y))
5147 return gnm_pow (1 + x, y);
5148 else if (y < 0)
5149 return 1 / pow1p (x, -y);
5150 else {
5151 gnm_float x1 = gnm_floor (x * 65536 + 0.5) / 65536;
5152 gnm_float x2 = x - x1;
5153 gnm_float h, l;
5154 ebd0 (y, y*(x+1), &h, &l);
5155 PAIR_ADD (-y*x1, h, l);
5156 PAIR_ADD (-y*x2, h, l);
5157
5158 #if 0
5159 g_printerr ("pow1p (%.20g, %.20g)\n", x, y);
5160 #endif
5161
5162 return gnm_exp (-l) * gnm_exp (-h);
5163 }
5164 }
5165
5166 /**
5167 * pow1pm1:
5168 * @x: a number
5169 * @y: a number
5170 *
5171 * Returns: The result of (1+@x)^@y-1 with less rounding error than the
5172 * naive formula.
5173 */
5174 gnm_float
pow1pm1(gnm_float x,gnm_float y)5175 pow1pm1 (gnm_float x, gnm_float y)
5176 {
5177 if (x <= -1)
5178 return gnm_pow (1 + x, y) - 1;
5179 else
5180 return gnm_expm1 (y * gnm_log1p (x));
5181 }
5182
5183
5184 /*
5185 ---------------------------------------------------------------------
5186 Matrix functions
5187 ---------------------------------------------------------------------
5188 */
5189
5190 GType
gnm_matrix_get_type(void)5191 gnm_matrix_get_type (void)
5192 {
5193 static GType t = 0;
5194
5195 if (t == 0)
5196 t = g_boxed_type_register_static ("GnmMatrix",
5197 (GBoxedCopyFunc)gnm_matrix_ref,
5198 (GBoxedFreeFunc)gnm_matrix_unref);
5199 return t;
5200 }
5201
5202 /**
5203 * gnm_matrix_new:
5204 * @rows: Number of rows.
5205 * @cols: Number of columns.
5206 *
5207 * Returns: (transfer full): A new #GnmMatrix.
5208 */
5209 /* Note the order: y then x. */
5210 GnmMatrix *
gnm_matrix_new(int rows,int cols)5211 gnm_matrix_new (int rows, int cols)
5212 {
5213 GnmMatrix *m = g_new (GnmMatrix, 1);
5214 int r;
5215
5216 m->ref_count = 1;
5217 m->rows = rows;
5218 m->cols = cols;
5219 m->data = g_new (gnm_float *, rows);
5220 for (r = 0; r < rows; r++)
5221 m->data[r] = g_new (gnm_float, cols);
5222
5223 return m;
5224 }
5225
5226 /**
5227 * gnm_matrix_ref:
5228 * @m: (transfer none) (nullable): #GnmMatrix
5229 *
5230 * Returns: (transfer full) (nullable): a new reference to @m.
5231 */
5232 GnmMatrix *
gnm_matrix_ref(GnmMatrix * m)5233 gnm_matrix_ref (GnmMatrix *m)
5234 {
5235 if (m)
5236 m->ref_count++;
5237 return m;
5238 }
5239
5240 /**
5241 * gnm_matrix_unref:
5242 * @m: (transfer full) (nullable): #GnmMatrix
5243 */
5244 void
gnm_matrix_unref(GnmMatrix * m)5245 gnm_matrix_unref (GnmMatrix *m)
5246 {
5247 int r;
5248
5249 if (!m || m->ref_count-- > 1)
5250 return;
5251
5252 for (r = 0; r < m->rows; r++)
5253 g_free (m->data[r]);
5254 g_free (m->data);
5255 g_free (m);
5256 }
5257
5258 /**
5259 * gnm_matrix_is_empty:
5260 * @m: (nullable): A #GnmMatrix
5261 *
5262 * Returns: %TRUE if @m is empty.
5263 */
5264 gboolean
gnm_matrix_is_empty(GnmMatrix const * m)5265 gnm_matrix_is_empty (GnmMatrix const *m)
5266 {
5267 return m == NULL || m->rows <= 0 || m->cols <= 0;
5268 }
5269
5270 /**
5271 * gnm_matrix_from_value:
5272 * @v: #GnmValue
5273 * @perr: (out) (transfer full): #GnmValue with error value
5274 * @ep: Evaluation location
5275 *
5276 * Returns: (transfer full) (nullable): A new #GnmMatrix, %NULL on error.
5277 */
5278 GnmMatrix *
gnm_matrix_from_value(GnmValue const * v,GnmValue ** perr,GnmEvalPos const * ep)5279 gnm_matrix_from_value (GnmValue const *v, GnmValue **perr, GnmEvalPos const *ep)
5280 {
5281 int cols, rows;
5282 int c, r;
5283 GnmMatrix *m = NULL;
5284
5285 *perr = NULL;
5286 cols = value_area_get_width (v, ep);
5287 rows = value_area_get_height (v, ep);
5288 m = gnm_matrix_new (rows, cols);
5289 for (r = 0; r < rows; r++) {
5290 for (c = 0; c < cols; c++) {
5291 GnmValue const *v1 = value_area_fetch_x_y (v, c, r, ep);
5292 if (VALUE_IS_ERROR (v1)) {
5293 *perr = value_dup (v1);
5294 gnm_matrix_unref (m);
5295 return NULL;
5296 }
5297
5298 m->data[r][c] = value_get_as_float (v1);
5299 }
5300 }
5301 return m;
5302 }
5303
5304 /**
5305 * gnm_matrix_to_value:
5306 * @m: #GnmMatrix
5307 *
5308 * Returns: (transfer full): A #GnmValue array
5309 */
5310 GnmValue *
gnm_matrix_to_value(GnmMatrix const * m)5311 gnm_matrix_to_value (GnmMatrix const *m)
5312 {
5313 GnmValue *res = value_new_array_non_init (m->cols, m->rows);
5314 int c, r;
5315
5316 for (c = 0; c < m->cols; c++) {
5317 res->v_array.vals[c] = g_new (GnmValue *, m->rows);
5318 for (r = 0; r < m->rows; r++)
5319 res->v_array.vals[c][r] = value_new_float (m->data[r][c]);
5320 }
5321 return res;
5322 }
5323
5324 /**
5325 * gnm_matrix_multiply:
5326 * @C: Output #GnmMatrix
5327 * @A: #GnmMatrix
5328 * @B: #GnmMatrix
5329 *
5330 * Computes @A * @B and stores the result in @C. The matrices must have
5331 * suitable sizes.
5332 */
5333 void
gnm_matrix_multiply(GnmMatrix * C,const GnmMatrix * A,const GnmMatrix * B)5334 gnm_matrix_multiply (GnmMatrix *C, const GnmMatrix *A, const GnmMatrix *B)
5335 {
5336 void *state;
5337 GnmAccumulator *acc;
5338 int c, r, i;
5339
5340 g_return_if_fail (C != NULL);
5341 g_return_if_fail (A != NULL);
5342 g_return_if_fail (B != NULL);
5343 g_return_if_fail (C->rows == A->rows);
5344 g_return_if_fail (C->cols == B->cols);
5345 g_return_if_fail (A->cols == B->rows);
5346
5347 state = gnm_accumulator_start ();
5348 acc = gnm_accumulator_new ();
5349
5350 for (r = 0; r < C->rows; r++) {
5351 for (c = 0; c < C->cols; c++) {
5352 go_accumulator_clear (acc);
5353 for (i = 0; i < A->cols; ++i) {
5354 GnmQuad p;
5355 gnm_quad_mul12 (&p,
5356 A->data[r][i],
5357 B->data[i][c]);
5358 gnm_accumulator_add_quad (acc, &p);
5359 }
5360 C->data[r][c] = gnm_accumulator_value (acc);
5361 }
5362 }
5363
5364 gnm_accumulator_free (acc);
5365 gnm_accumulator_end (state);
5366 }
5367
5368 /***************************************************************************/
5369
5370 static int
gnm_matrix_eigen_max_index(gnm_float * row,guint row_n,guint size)5371 gnm_matrix_eigen_max_index (gnm_float *row, guint row_n, guint size)
5372 {
5373 guint i, res = row_n + 1;
5374 gnm_float max;
5375
5376 if (res >= size)
5377 return (size - 1);
5378
5379 max = gnm_abs (row[res]);
5380
5381 for (i = res + 1; i < size; i++)
5382 if (gnm_abs (row[i]) > max) {
5383 res = i;
5384 max = gnm_abs (row[i]);
5385 }
5386 return res;
5387 }
5388
5389 static void
gnm_matrix_eigen_rotate(gnm_float ** matrix,guint k,guint l,guint i,guint j,gnm_float c,gnm_float s)5390 gnm_matrix_eigen_rotate (gnm_float **matrix, guint k, guint l, guint i, guint j, gnm_float c, gnm_float s)
5391 {
5392 gnm_float x = c * matrix[k][l] - s * matrix[i][j];
5393 gnm_float y = s * matrix[k][l] + c * matrix[i][j];
5394
5395 matrix[k][l] = x;
5396 matrix[i][j] = y;
5397 }
5398
5399 static void
gnm_matrix_eigen_update(guint k,gnm_float t,gnm_float * eigenvalues,gboolean * changed,guint * state)5400 gnm_matrix_eigen_update (guint k, gnm_float t, gnm_float *eigenvalues, gboolean *changed, guint *state)
5401 {
5402 gnm_float y = eigenvalues[k];
5403 gboolean unchanged;
5404 eigenvalues[k] += t;
5405 unchanged = (y == eigenvalues[k]);
5406 if (changed[k] && unchanged) {
5407 changed[k] = FALSE;
5408 (*state)--;
5409 } else if (!changed[k] && !unchanged) {
5410 changed[k] = TRUE;
5411 (*state)++;
5412 }
5413 }
5414
5415 /**
5416 * gnm_matrix_eigen:
5417 * @m: Input #GnmMatrix
5418 * @EIG: Output #GnmMatrix
5419 * @eigenvalues: (out): Output location for eigen values.
5420 *
5421 * Calculates the eigenvalues and eigenvectors of a real symmetric matrix.
5422 *
5423 * This is the Jacobi iterative process in which we use a sequence of
5424 * Jacobi rotations (two-sided Givens rotations) in order to reduce the
5425 * magnitude of off-diagonal elements while preserving eigenvalues.
5426 */
5427 gboolean
gnm_matrix_eigen(GnmMatrix const * m,GnmMatrix * EIG,gnm_float * eigenvalues)5428 gnm_matrix_eigen (GnmMatrix const *m, GnmMatrix *EIG, gnm_float *eigenvalues)
5429 {
5430 guint i, state, usize, *ind;
5431 gboolean *changed;
5432 guint counter = 0;
5433 gnm_float **matrix;
5434 gnm_float **eigenvectors;
5435
5436 g_return_val_if_fail (m != NULL, FALSE);
5437 g_return_val_if_fail (m->rows == m->cols, FALSE);
5438 g_return_val_if_fail (EIG != NULL, FALSE);
5439 g_return_val_if_fail (EIG->rows == EIG->cols, FALSE);
5440 g_return_val_if_fail (EIG->rows == m->rows, FALSE);
5441
5442 matrix = m->data;
5443 eigenvectors = EIG->data;
5444 usize = m->rows;
5445
5446 state = usize;
5447
5448 ind = g_new (guint, usize);
5449 changed = g_new (gboolean, usize);
5450
5451 for (i = 0; i < usize; i++) {
5452 guint j;
5453 for (j = 0; j < usize; j++)
5454 eigenvectors[j][i] = 0.;
5455 eigenvectors[i][i] = 1.;
5456 eigenvalues[i] = matrix[i][i];
5457 ind[i] = gnm_matrix_eigen_max_index (matrix[i], i, usize);
5458 changed[i] = TRUE;
5459 }
5460
5461 while (usize > 1 && state != 0) {
5462 guint k, l, m = 0;
5463 gnm_float c, s, y, pivot, t;
5464
5465 counter++;
5466 if (counter > 400000) {
5467 g_free (ind);
5468 g_free (changed);
5469 g_print ("gnm_matrix_eigen exceeded iterations\n");
5470 return FALSE;
5471 }
5472 for (k = 1; k < usize - 1; k++)
5473 if (gnm_abs (matrix[k][ind[k]]) > gnm_abs (matrix[m][ind[m]]))
5474 m = k;
5475 l = ind[m];
5476 pivot = matrix[m][l];
5477 /* pivot is (m,l) */
5478 if (pivot == 0) {
5479 /* All remaining off-diagonal elements are zero. We're done. */
5480 break;
5481 }
5482
5483 y = (eigenvalues[l] - eigenvalues[m]) / 2;
5484 t = gnm_abs (y) + gnm_hypot (pivot, y);
5485 s = gnm_hypot (pivot, t);
5486 c = t / s;
5487 s = pivot / s;
5488 t = pivot * pivot / t;
5489 if (y < 0) {
5490 s = -s;
5491 t = -t;
5492 }
5493 matrix[m][l] = 0.;
5494 gnm_matrix_eigen_update (m, -t, eigenvalues, changed, &state);
5495 gnm_matrix_eigen_update (l, t, eigenvalues, changed, &state);
5496 for (i = 0; i < m; i++)
5497 gnm_matrix_eigen_rotate (matrix, i, m, i, l, c, s);
5498 for (i = m + 1; i < l; i++)
5499 gnm_matrix_eigen_rotate (matrix, m, i, i, l, c, s);
5500 for (i = l + 1; i < usize; i++)
5501 gnm_matrix_eigen_rotate (matrix, m, i, l, i, c, s);
5502 for (i = 0; i < usize; i++) {
5503 gnm_float x = c * eigenvectors[i][m] - s * eigenvectors[i][l];
5504 gnm_float y = s * eigenvectors[i][m] + c * eigenvectors[i][l];
5505
5506 eigenvectors[i][m] = x;
5507 eigenvectors[i][l] = y;
5508 }
5509 ind[m] = gnm_matrix_eigen_max_index (matrix[m], m, usize);
5510 ind[l] = gnm_matrix_eigen_max_index (matrix[l], l, usize);
5511 }
5512
5513 g_free (ind);
5514 g_free (changed);
5515
5516 return TRUE;
5517 }
5518
5519 /* ------------------------------------------------------------------------- */
5520
5521 static void
swap_row_and_col(GnmMatrix * M,int a,int b)5522 swap_row_and_col (GnmMatrix *M, int a, int b)
5523 {
5524 gnm_float *r;
5525 int i;
5526
5527 // Swap rows
5528 r = M->data[a];
5529 M->data[a] = M->data[b];
5530 M->data[b] = r;
5531
5532 // Swap cols
5533 for (i = 0; i < M->rows; i++) {
5534 gnm_float d = M->data[i][a];
5535 M->data[i][a] = M->data[i][b];
5536 M->data[i][b] = d;
5537 }
5538 }
5539
5540
5541 gboolean
gnm_matrix_modified_cholesky(GnmMatrix const * A,GnmMatrix * L,gnm_float * D,gnm_float * E,int * P)5542 gnm_matrix_modified_cholesky (GnmMatrix const *A,
5543 GnmMatrix *L,
5544 gnm_float *D,
5545 gnm_float *E,
5546 int *P)
5547 {
5548 int n = A->cols;
5549 gnm_float nu, xsi, gam, bsqr, delta;
5550 int i, j;
5551 GnmMatrix *G, *C;
5552
5553 g_return_val_if_fail (A->rows == A->cols, FALSE);
5554 g_return_val_if_fail (A->rows == L->rows, FALSE);
5555 g_return_val_if_fail (A->cols == L->cols, FALSE);
5556
5557 // Copy A into L; Use G and C as aliases for L.
5558 for (i = 0; i < n; i++)
5559 for (j = 0; j < n; j++)
5560 L->data[i][j] = A->data[i][j];
5561 G = L;
5562 C = G;
5563
5564 // Init permutation as identity.
5565 for (i = 0; i < n; i++)
5566 P[i] = i;
5567
5568 nu = n == 1 ? 1.0 : gnm_sqrt (n * n - 1);
5569 gam = xsi = 0;
5570 for (i = 0; i < n; i++) {
5571 gnm_float aii = gnm_abs (G->data[i][i]);
5572 gam = MAX (gam, aii);
5573 for (j = i + 1; j < n; j++) {
5574 gnm_float aij = gnm_abs (G->data[i][j]);
5575 xsi = MAX (xsi, aij);
5576 }
5577 }
5578 bsqr = MAX (MAX (gam, xsi / nu), GNM_EPSILON);
5579 delta = MAX (gam + xsi, 1.0) * GNM_EPSILON;
5580
5581 for (j = 0; j < n; j++) {
5582 int q, s;
5583 gnm_float theta_j = 0, dj;
5584
5585 q = j;
5586 for (i = j + 1; i < n; i++) {
5587 if (gnm_abs (C->data[i][i]) > gnm_abs (C->data[q][q]))
5588 q = i;
5589 }
5590
5591 if (j != q) {
5592 int a;
5593 gnm_float b;
5594
5595 swap_row_and_col (L, j, q);
5596 a = P[j]; P[j] = P[q]; P[q] = a;
5597 b = D[j]; D[j] = D[q]; D[q] = b;
5598 if (E) { b = E[j]; E[j] = E[q]; E[q] = b; }
5599 }
5600
5601 for (s = 0; s < j; s++)
5602 L->data[j][s] = C->data[j][s] / D[s];
5603
5604 for (i = j + 1; i < n; i++) {
5605 int s;
5606 gnm_float d = G->data[i][j];
5607
5608 for (s = 0; s < j; s++)
5609 d -= L->data[j][s] * C->data[i][s];
5610 C->data[i][j] = d;
5611
5612 theta_j = MAX (theta_j, gnm_abs (d));
5613 }
5614
5615 dj = MAX (theta_j * theta_j / bsqr, delta);
5616 dj = MAX (dj, gnm_abs (C->data[j][j]));
5617 D[j] = dj;
5618 if (E) E[j] = dj - C->data[j][j];
5619
5620 for (i = j + 1; i < n; i++) {
5621 gnm_float cij = C->data[i][j];
5622 C->data[i][i] -= cij * cij / D[j];
5623 }
5624 }
5625
5626 for (i = 0; i < n; i++) {
5627 for (j = i + 1; j < n; j++)
5628 L->data[i][j] = 0;
5629 L->data[i][i] = 1;
5630 }
5631
5632 return TRUE;
5633 }
5634
5635 GORegressionResult
gnm_linear_solve_posdef(GnmMatrix const * A,const gnm_float * b,gnm_float * x)5636 gnm_linear_solve_posdef (GnmMatrix const *A, const gnm_float *b,
5637 gnm_float *x)
5638 {
5639 int i, j, n;
5640 GnmMatrix *L;
5641 gnm_float *D, *E;
5642 int *P;
5643 GORegressionResult res;
5644 gboolean ok;
5645
5646 g_return_val_if_fail (A != NULL, GO_REG_invalid_dimensions);
5647 g_return_val_if_fail (A->rows == A->cols, GO_REG_invalid_dimensions);
5648 g_return_val_if_fail (b != NULL, GO_REG_invalid_dimensions);
5649 g_return_val_if_fail (x != NULL, GO_REG_invalid_dimensions);
5650
5651 n = A->cols;
5652 L = gnm_matrix_new (n, n);
5653 D = g_new (gnm_float, n);
5654 E = g_new (gnm_float, n);
5655 P = g_new (int, n);
5656
5657 ok = gnm_matrix_modified_cholesky (A, L, D, E, P);
5658 if (!ok) {
5659 res = GO_REG_invalid_data;
5660 goto done;
5661 }
5662
5663 if (gnm_debug_flag ("posdef")) {
5664 for (i = 0; i < n; i++)
5665 g_printerr ("Posdef E[i] = %g\n", E[P[i]]);
5666 }
5667
5668 // The only information from the above we use is E and P.
5669 // However, we reuse the memory for L for A+E
5670 for (i = 0; i < n; i++) {
5671 for (j = 0; j < n; j++)
5672 L->data[i][j] = A->data[i][j];
5673 L->data[i][i] += E[P[i]];
5674 }
5675
5676 res = gnm_linear_solve (L, b, x);
5677
5678 done:
5679 g_free (P);
5680 g_free (E);
5681 g_free (D);
5682 gnm_matrix_unref (L);
5683
5684 return res;
5685 }
5686
5687 /* ------------------------------------------------------------------------- */
5688
5689 GORegressionResult
gnm_linear_solve(GnmMatrix const * A,const gnm_float * b,gnm_float * x)5690 gnm_linear_solve (GnmMatrix const *A, const gnm_float *b,
5691 gnm_float *x)
5692 {
5693 g_return_val_if_fail (A != NULL, GO_REG_invalid_dimensions);
5694 g_return_val_if_fail (A->rows == A->cols, GO_REG_invalid_dimensions);
5695 g_return_val_if_fail (b != NULL, GO_REG_invalid_dimensions);
5696 g_return_val_if_fail (x != NULL, GO_REG_invalid_dimensions);
5697
5698 return
5699 #ifdef GNM_WITH_LONG_DOUBLE
5700 go_linear_solvel
5701 #else
5702 go_linear_solve
5703 #endif
5704 (A->data, b, A->rows, x);
5705 }
5706
5707 GORegressionResult
gnm_linear_solve_multiple(GnmMatrix const * A,GnmMatrix * B)5708 gnm_linear_solve_multiple (GnmMatrix const *A, GnmMatrix *B)
5709 {
5710 g_return_val_if_fail (A != NULL, GO_REG_invalid_dimensions);
5711 g_return_val_if_fail (B != NULL, GO_REG_invalid_dimensions);
5712 g_return_val_if_fail (A->rows == A->cols, GO_REG_invalid_dimensions);
5713 g_return_val_if_fail (A->rows == B->rows, GO_REG_invalid_dimensions);
5714
5715 return
5716 #ifdef GNM_WITH_LONG_DOUBLE
5717 go_linear_solve_multiplel
5718 #else
5719 go_linear_solve_multiple
5720 #endif
5721 (A->data, B->data, A->rows, B->cols);
5722 }
5723
5724 /* ------------------------------------------------------------------------- */
5725
5726 #ifdef GNM_SUPPLIES_ERFL
5727 long double
erfl(long double x)5728 erfl (long double x)
5729 {
5730 if (fabsl (x) < 0.125) {
5731 /* For small x the pnorm formula loses precision. */
5732 long double sum = 0;
5733 long double term = x * 2 / sqrtl (M_PIgnum);
5734 long double n;
5735 long double x2 = x * x;
5736
5737 for (n = 0; fabsl (term) >= fabsl (sum) * LDBL_EPSILON ; n++) {
5738 sum += term / (2 * n + 1);
5739 term *= -x2 / (n + 1);
5740 }
5741
5742 return sum;
5743 }
5744 return pnorm (x * M_SQRT2gnum, 0, 1, TRUE, FALSE) * 2 - 1;
5745 }
5746 #endif
5747
5748 /* ------------------------------------------------------------------------- */
5749
5750 #ifdef GNM_SUPPLIES_ERFCL
5751 long double
erfcl(long double x)5752 erfcl (long double x)
5753 {
5754 return 2 * pnorm (x * M_SQRT2gnum, 0, 1, FALSE, FALSE);
5755 }
5756 #endif
5757
5758 /* ------------------------------------------------------------------------- */
5759
5760 static gnm_float
gnm_owent_T1(gnm_float h,gnm_float a,int order)5761 gnm_owent_T1 (gnm_float h, gnm_float a, int order)
5762 {
5763 const gnm_float hs = -0.5 * (h * h);
5764 const gnm_float dhs = gnm_exp (hs);
5765 const gnm_float as = a * a;
5766 gnm_float aj = a / (M_PIgnum * 2);
5767 gnm_float dj = gnm_expm1 (hs);
5768 gnm_float gj = hs * dhs;
5769 gnm_float res = gnm_atan (a) / (M_PIgnum * 2);
5770 int j;
5771
5772 for (j = 1; j <= order; j++) {
5773 res += dj * aj / (j + j - 1);
5774
5775 aj *= as;
5776 dj = gj - dj;
5777 gj *= hs / (j + 1);
5778 }
5779
5780 return res;
5781 }
5782
5783 static gnm_float
gnm_owent_T2(gnm_float h,gnm_float a,int order)5784 gnm_owent_T2 (gnm_float h, gnm_float a, int order)
5785 {
5786 const gnm_float ah = a * h;
5787 const gnm_float as = -a * a;
5788 const gnm_float y = 1 / (h * h);
5789 gnm_float val = 0;
5790 gnm_float vi = a * dnorm (ah, 0, 1, FALSE);
5791 gnm_float z = gnm_erf (ah / M_SQRT2gnum) / (2 * h);
5792 int i;
5793
5794 for (i = 1; i <= 2 * order + 1; i += 2) {
5795 val += z;
5796 z = y * (vi - i * z);
5797 vi *= as;
5798 }
5799 return val * dnorm (h, 0, 1, FALSE);
5800 }
5801
5802 static gnm_float
gnm_owent_T3(gnm_float h,gnm_float a,int order)5803 gnm_owent_T3 (gnm_float h, gnm_float a, int order)
5804 {
5805 static const gnm_float c2[] = {
5806 GNM_const(+0.99999999999999987510),
5807 GNM_const(-0.99999999999988796462),
5808 GNM_const(+0.99999999998290743652),
5809 GNM_const(-0.99999999896282500134),
5810 GNM_const(+0.99999996660459362918),
5811 GNM_const(-0.99999933986272476760),
5812 GNM_const(+0.99999125611136965852),
5813 GNM_const(-0.99991777624463387686),
5814 GNM_const(+0.99942835555870132569),
5815 GNM_const(-0.99697311720723000295),
5816 GNM_const(+0.98751448037275303682),
5817 GNM_const(-0.95915857980572882813),
5818 GNM_const(+0.89246305511006708555),
5819 GNM_const(-0.76893425990463999675),
5820 GNM_const(+0.58893528468484693250),
5821 GNM_const(-0.38380345160440256652),
5822 GNM_const(+0.20317601701045299653),
5823 GNM_const(-0.82813631607004984866E-01),
5824 GNM_const(+0.24167984735759576523E-01),
5825 GNM_const(-0.44676566663971825242E-02),
5826 GNM_const(+0.39141169402373836468E-03)
5827 };
5828
5829 const gnm_float ah = a * h;
5830 const gnm_float as = a * a;
5831 const gnm_float y = 1 / (h * h);
5832 gnm_float vi = a * dnorm (ah, 0, 1, FALSE);
5833 gnm_float zi = gnm_erf (ah / M_SQRT2gnum) / (2 * h);
5834 gnm_float val = 0;
5835 int i;
5836
5837 g_return_val_if_fail (order < (int)G_N_ELEMENTS(c2), gnm_nan);
5838
5839 for (i = 0; i <= order; i++) {
5840 val += zi * c2[i];
5841 zi = y * ((i + i + 1) * zi - vi);
5842 vi *= as;
5843 }
5844 return val * dnorm (h, 0, 1, FALSE);
5845 }
5846
5847 static gnm_float
gnm_owent_T4(gnm_float h,gnm_float a,int order)5848 gnm_owent_T4 (gnm_float h, gnm_float a, int order)
5849 {
5850 const gnm_float hs = h * h;
5851 const gnm_float as = -a * a;
5852 gnm_float ai = a * gnm_exp (-0.5 * hs * (1 - as)) / (2 * M_PIgnum);
5853 gnm_float yi = 1;
5854 gnm_float val = 0;
5855 int i;
5856
5857 for (i = 1; i <= 2 * order + 1; i += 2) {
5858 val += ai * yi;
5859 yi = (1 - hs * yi) / (i + 2);
5860 ai *= as;
5861 }
5862 return val;
5863 }
5864
5865 static gnm_float
gnm_owent_T5(gnm_float h,gnm_float a,int order)5866 gnm_owent_T5 (gnm_float h, gnm_float a, int order)
5867 {
5868 static const gnm_float pts[] = {
5869 GNM_const(0.35082039676451715489E-02),
5870 GNM_const(0.31279042338030753740E-01),
5871 GNM_const(0.85266826283219451090E-01),
5872 GNM_const(0.16245071730812277011),
5873 GNM_const(0.25851196049125434828),
5874 GNM_const(0.36807553840697533536),
5875 GNM_const(0.48501092905604697475),
5876 GNM_const(0.60277514152618576821),
5877 GNM_const(0.71477884217753226516),
5878 GNM_const(0.81475510988760098605),
5879 GNM_const(0.89711029755948965867),
5880 GNM_const(0.95723808085944261843),
5881 GNM_const(0.99178832974629703586)
5882 };
5883 static const gnm_float wts[] = {
5884 0.18831438115323502887E-01,
5885 0.18567086243977649478E-01,
5886 0.18042093461223385584E-01,
5887 0.17263829606398753364E-01,
5888 0.16243219975989856730E-01,
5889 0.14994592034116704829E-01,
5890 0.13535474469662088392E-01,
5891 0.11886351605820165233E-01,
5892 0.10070377242777431897E-01,
5893 0.81130545742299586629E-02,
5894 0.60419009528470238773E-02,
5895 0.38862217010742057883E-02,
5896 0.16793031084546090448E-02
5897 };
5898 const gnm_float as = a * a;
5899 const gnm_float hs = -0.5 * h * h;
5900 gnm_float val = 0;
5901 int i;
5902
5903 g_return_val_if_fail (order <= (int)G_N_ELEMENTS(pts), gnm_nan);
5904 g_return_val_if_fail (order <= (int)G_N_ELEMENTS(wts), gnm_nan);
5905
5906 for (i = 0; i < order; i++) {
5907 gnm_float r = 1 + as * pts[i];
5908 val += wts[i] * gnm_exp (hs * r) / r;
5909 }
5910
5911 return val * a;
5912 }
5913
5914 static gnm_float
gnm_owent_T6(gnm_float h,gnm_float a,int order)5915 gnm_owent_T6 (gnm_float h, gnm_float a, int order)
5916 {
5917 const gnm_float normh = pnorm (h, 0, 1, FALSE, FALSE);
5918 const gnm_float normhC = 1 - normh;
5919 const gnm_float y = 1 - a;
5920 const gnm_float r = gnm_atan2 (y, 1 + a);
5921 gnm_float val = 0.5 * normh * normhC;
5922
5923 if (r != 0)
5924 val -= r * gnm_exp (-0.5 * y * h * h / r) / (2 * M_PIgnum);
5925
5926 return val;
5927 }
5928
5929 static gnm_float
gnm_owent_helper(gnm_float h,gnm_float a)5930 gnm_owent_helper (gnm_float h, gnm_float a)
5931 {
5932 static const double hrange[] = {
5933 0.02, 0.06, 0.09, 0.125, 0.26, 0.4, 0.6,
5934 1.6, 1.7, 2.33, 2.4, 3.36, 3.4, 4.8
5935 };
5936 static const double arange[] = {
5937 0.025, 0.09, 0.15, 0.36, 0.5, 0.9, 0.99999
5938 };
5939 static const guint8 method[] = {
5940 1, 1, 2,13,13,13,13,13,13,13,13,16,16,16, 9,
5941 1, 2, 2, 3, 3, 5, 5,14,14,15,15,16,16,16, 9,
5942 2, 2, 3, 3, 3, 5, 5,15,15,15,15,16,16,16,10,
5943 2, 2, 3, 5, 5, 5, 5, 7, 7,16,16,16,16,16,10,
5944 2, 3, 3, 5, 5, 6, 6, 8, 8,17,17,17,12,12,11,
5945 2, 3, 5, 5, 5, 6, 6, 8, 8,17,17,17,12,12,12,
5946 2, 3, 4, 4, 6, 6, 8, 8,17,17,17,17,17,12,12,
5947 2, 3, 4, 4, 6, 6,18,18,18,18,17,17,17,12,12
5948 };
5949 int ai, hi;
5950
5951 g_return_val_if_fail (h >= 0, gnm_nan);
5952 g_return_val_if_fail (a >= 0 && a <= 1, gnm_nan);
5953
5954 for (ai = 0; ai < (int)G_N_ELEMENTS(arange); ai++)
5955 if (a <= arange[ai])
5956 break;
5957 for (hi = 0; hi < (int)G_N_ELEMENTS(hrange); hi++)
5958 if (h <= hrange[hi])
5959 break;
5960
5961 switch (method[ai * (1 + G_N_ELEMENTS(hrange)) + hi]) {
5962 case 1: return gnm_owent_T1 (h, a, 2);
5963 case 2: return gnm_owent_T1 (h, a, 3);
5964 case 3: return gnm_owent_T1 (h, a, 4);
5965 case 4: return gnm_owent_T1 (h, a, 5);
5966 case 5: return gnm_owent_T1 (h, a, 7);
5967 case 6: return gnm_owent_T1 (h, a, 10);
5968 case 7: return gnm_owent_T1 (h, a, 12);
5969 case 8: return gnm_owent_T1 (h, a, 18);
5970 case 9: return gnm_owent_T2 (h, a, 10);
5971 case 10: return gnm_owent_T2 (h, a, 20);
5972 case 11: return gnm_owent_T2 (h, a, 30);
5973 case 12: return gnm_owent_T3 (h, a, 20);
5974 case 13: return gnm_owent_T4 (h, a, 4);
5975 case 14: return gnm_owent_T4 (h, a, 7);
5976 case 15: return gnm_owent_T4 (h, a, 8);
5977 case 16: return gnm_owent_T4 (h, a, 20);
5978 case 17: return gnm_owent_T5 (h, a, 13);
5979 case 18: return gnm_owent_T6 (h, a, 0);
5980 default:
5981 g_assert_not_reached ();
5982 }
5983 }
5984
5985 /*
5986 * See "Fast and Accurate Calculation of Owen’s T-Function" by
5987 * Mike Patefield and David Tandy. Journal of Statistical Software,
5988 * Volume 5, Issue 5, July 2000.
5989 *
5990 * Original code licensed under GPLv2+.
5991 */
5992 gnm_float
gnm_owent(gnm_float h,gnm_float a)5993 gnm_owent (gnm_float h, gnm_float a)
5994 {
5995 gnm_float res;
5996 gboolean neg;
5997
5998 /* Even in the "h" argument. */
5999 h = gnm_abs (h);
6000
6001 /* Odd in the "a" argument. */
6002 neg = (a < 0);
6003 a = gnm_abs (a);
6004
6005 if (a == 0)
6006 res = 0;
6007 else if (h == 0)
6008 res = gnm_atan (a) / (2 * M_PIgnum);
6009 else if (a == 1)
6010 res = 0.5 * pnorm (h, 0, 1, TRUE, FALSE) *
6011 pnorm (h, 0, 1, FALSE, FALSE);
6012 else if (a <= 1)
6013 res = gnm_owent_helper (h, a);
6014 else {
6015 gnm_float ah = h * a;
6016 /*
6017 * Use formula (2):
6018 *
6019 * T(h,a) = .5N(h) + .5N(ha) - N(h)N(ha) - T(ha,1/a)
6020 *
6021 * with care to avoid cancellation.
6022 */
6023 if (h <= 0.67) {
6024 gnm_float nh = 0.5 * gnm_erf (h / M_SQRT2gnum);
6025 gnm_float nah = 0.5 * gnm_erf (ah / M_SQRT2gnum);
6026 res = 0.25 - nh * nah -
6027 gnm_owent_helper (ah, 1 / a);
6028 } else {
6029 gnm_float nh = pnorm (h, 0, 1, FALSE, FALSE);
6030 gnm_float nah = pnorm (ah, 0, 1, FALSE, FALSE);
6031 res = 0.5 * (nh + nah) - nh * nah -
6032 gnm_owent_helper (ah, 1 / a);
6033 }
6034 }
6035
6036 /* Odd in the "a" argument. */
6037 if (neg)
6038 res = 0 - res;
6039
6040 return res;
6041 }
6042
6043 /* ------------------------------------------------------------------------- */
6044