1 /* $OpenBSD: monot.c,v 1.2 2021/05/29 10:35:56 bluhm Exp $ */
2
3 /*
4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19 /* monot.c
20 Floating point function test vectors.
21
22 Arguments and function values are synthesized for NPTS points in
23 the vicinity of each given tabulated test point. The points are
24 chosen to be near and on either side of the likely function algorithm
25 domain boundaries. Since the function programs change their methods
26 at these points, major coding errors or monotonicity failures might be
27 detected.
28
29 August, 1998
30 S. L. Moshier */
31
32
33 #include <stdio.h>
34
35 /* Avoid including math.h. */
36 double frexp (double, int *);
37 double ldexp (double, int);
38
39 /* Number of test points to generate on each side of tabulated point. */
40 #define NPTS 100
41
42 /* Functions of one variable. */
43 double exp (double);
44 double log (double);
45 double sin (double);
46 double cos (double);
47 double tan (double);
48 double atan (double);
49 double asin (double);
50 double acos (double);
51 double sinh (double);
52 double cosh (double);
53 double tanh (double);
54 double asinh (double);
55 double acosh (double);
56 double atanh (double);
57 double gamma (double);
58 double fabs (double);
59 double floor (double);
60 #if 0
61 double polylog (int, double);
62 #endif
63
64 struct oneargument
65 {
66 char *name; /* Name of the function. */
67 double (*func) (double);
68 double arg1; /* Function argument, assumed exact. */
69 double answer1; /* Exact, close to function value. */
70 double answer2; /* answer1 + answer2 has extended precision. */
71 double derivative; /* dy/dx evaluated at x = arg1. */
72 int thresh; /* Error report threshold. 2 = 1 ULP approx. */
73 };
74
75 /* Add this to error threshold test[i].thresh. */
76 #define OKERROR 0
77
78 /* Unit of relative error in test[i].thresh. */
79 static double MACHEP = 1.1102230246251565404236316680908203125e-16;
80
81
82 /* extern double MACHEP; */
83
84 #if 0
85 double w_polylog_3 (double x) {return polylog (3, x);}
86 #endif
87
88 static struct oneargument test1[] =
89 {
90 {"exp", exp, 1.0, 2.7182769775390625,
91 4.85091998273536028747e-6, 2.71828182845904523536, 2},
92 {"exp", exp, -1.0, 3.678741455078125e-1,
93 5.29566362982159552377e-6, 3.678794411714423215955e-1, 2},
94 {"exp", exp, 0.5, 1.648712158203125,
95 9.1124970031468486507878e-6, 1.64872127070012814684865, 2},
96 {"exp", exp, -0.5, 6.065216064453125e-1,
97 9.0532673209236037995e-6, 6.0653065971263342360e-1, 2},
98 {"exp", exp, 2.0, 7.3890533447265625,
99 2.75420408772723042746e-6, 7.38905609893065022723, 2},
100 {"exp", exp, -2.0, 1.353302001953125e-1,
101 5.08304130019189399949e-6, 1.3533528323661269189e-1, 2},
102 {"log", log, 1.41421356237309492343, 3.465728759765625e-1,
103 7.1430341006605745676897e-7, 7.0710678118654758708668e-1, 2},
104 {"log", log, 7.07106781186547461715e-1, -3.46588134765625e-1,
105 1.45444856522566402246e-5, 1.41421356237309517417, 2},
106 {"sin", sin, 7.85398163397448278999e-1, 7.0709228515625e-1,
107 1.4496030297502751942956e-5, 7.071067811865475460497e-1, 2},
108 {"sin", sin, -7.85398163397448501044e-1, -7.071075439453125e-1,
109 7.62758764840238811175e-7, 7.07106781186547389040e-1, 2},
110 {"sin", sin, 1.570796326794896558, 9.999847412109375e-1,
111 1.52587890625e-5, 6.12323399573676588613e-17, 2},
112 {"sin", sin, -1.57079632679489678004, -1.0,
113 1.29302922820150306903e-32, -1.60812264967663649223e-16, 2},
114 {"sin", sin, 4.712388980384689674, -1.0,
115 1.68722975549458979398e-32, -1.83697019872102976584e-16, 2},
116 {"sin", sin, -4.71238898038468989604, 9.999847412109375e-1,
117 1.52587890625e-5, 3.83475850529283315008e-17, 2},
118 {"cos", cos, 3.92699081698724139500E-1, 9.23873901367187500000E-1,
119 5.63114409926198633370E-6, -3.82683432365089757586E-1, 2},
120 {"cos", cos, 7.85398163397448278999E-1, 7.07092285156250000000E-1,
121 1.44960302975460497458E-5, -7.07106781186547502752E-1, 2},
122 {"cos", cos, 1.17809724509617241850E0, 3.82675170898437500000E-1,
123 8.26146665231415693919E-6, -9.23879532511286738554E-1, 2},
124 {"cos", cos, 1.96349540849362069750E0, -3.82690429687500000000E-1,
125 6.99732241029898567203E-6, -9.23879532511286785419E-1, 2},
126 {"cos", cos, 2.35619449019234483700E0, -7.07107543945312500000E-1,
127 7.62758765040545859856E-7, -7.07106781186547589348E-1, 2},
128 {"cos", cos, 2.74889357189106897650E0, -9.23889160156250000000E-1,
129 9.62764496328487887036E-6, -3.82683432365089870728E-1, 2},
130 {"cos", cos, 3.14159265358979311600E0, -1.00000000000000000000E0,
131 7.49879891330928797323E-33, -1.22464679914735317723E-16, 2},
132 {"tan", tan, 7.85398163397448278999E-1, 9.999847412109375e-1,
133 1.52587890624387676600E-5, 1.99999999999999987754E0, 2},
134 {"tan", tan, 1.17809724509617241850E0, 2.41419982910156250000E0,
135 1.37332715322352112604E-5, 6.82842712474618858345E0, 2},
136 {"tan", tan, 1.96349540849362069750E0, -2.41421508789062500000E0,
137 1.52551752942854759743E-6, 6.82842712474619262118E0, 2},
138 {"tan", tan, 2.35619449019234483700E0, -1.00001525878906250000E0,
139 1.52587890623163029801E-5, 2.00000000000000036739E0, 2},
140 {"tan", tan, 2.74889357189106897650E0, -4.14215087890625000000E-1,
141 1.52551752982565655126E-6, 1.17157287525381000640E0, 2},
142 {"atan", atan, 4.14213562373094923430E-1, 3.92684936523437500000E-1,
143 1.41451752865477964149E-5, 8.53553390593273837869E-1, 2},
144 {"atan", atan, 1.0, 7.85385131835937500000E-1,
145 1.30315615108096156608E-5, 0.5, 2},
146 {"atan", atan, 2.41421356237309492343E0, 1.17808532714843750000E0,
147 1.19179477349460632350E-5, 1.46446609406726250782E-1, 2},
148 {"atan", atan, -2.41421356237309514547E0, -1.17810058593750000000E0,
149 3.34084132752141908545E-6, 1.46446609406726227789E-1, 2},
150 {"atan", atan, -1.0, -7.85400390625000000000E-1,
151 2.22722755169038433915E-6, 0.5, 2},
152 {"atan", atan, -4.14213562373095145475E-1, -3.92700195312500000000E-1,
153 1.11361377576267665972E-6, 8.53553390593273703853E-1, 2},
154 {"asin", asin, 3.82683432365089615246E-1, 3.92684936523437500000E-1,
155 1.41451752864854321970E-5, 1.08239220029239389286E0, 2},
156 {"asin", asin, 0.5, 5.23590087890625000000E-1,
157 8.68770767387307710723E-6, 1.15470053837925152902E0, 2},
158 {"asin", asin, 7.07106781186547461715E-1, 7.85385131835937500000E-1,
159 1.30315615107209645016E-5, 1.41421356237309492343E0, 2},
160 {"asin", asin, 9.23879532511286738483E-1, 1.17808532714843750000E0,
161 1.19179477349183147612E-5, 2.61312592975275276483E0, 2},
162 {"asin", asin, -0.5, -5.23605346679687500000E-1,
163 6.57108138862692289277E-6, 1.15470053837925152902E0, 2},
164 {"acos", acos, 1.95090322016128192573E-1, 1.37443542480468750000E0,
165 1.13611408471185777914E-5, -1.01959115820831832232E0, 2},
166 {"acos", acos, 3.82683432365089615246E-1, 1.17808532714843750000E0,
167 1.19179477351337991247E-5, -1.08239220029239389286E0, 2},
168 {"acos", acos, 0.5, 1.04719543457031250000E0,
169 2.11662628524615421446E-6, -1.15470053837925152902E0, 2},
170 {"acos", acos, 7.07106781186547461715E-1, 7.85385131835937500000E-1,
171 1.30315615108982668201E-5, -1.41421356237309492343E0, 2},
172 {"acos", acos, 9.23879532511286738483E-1, 3.92684936523437500000E-1,
173 1.41451752867009165605E-5, -2.61312592975275276483E0, 2},
174 {"acos", acos, 9.80785280403230430579E-1, 1.96334838867187500000E-1,
175 1.47019821746724723933E-5, -5.12583089548300990774E0, 2},
176 {"acos", acos, -0.5, 2.09439086914062500000E0,
177 4.23325257049230842892E-6, -1.15470053837925152902E0, 2},
178 {"sinh", sinh, 1.0, 1.17518615722656250000E0,
179 1.50364172389568823819E-5, 1.54308063481524377848E0, 2},
180 {"sinh", sinh, 7.09089565712818057364E2, 4.49423283712885057274E307,
181 4.25947714184369757620E208, 4.49423283712885057274E307, 2},
182 {"sinh", sinh, 2.22044604925031308085E-16, 0.00000000000000000000E0,
183 2.22044604925031308085E-16, 1.00000000000000000000E0, 2},
184 {"cosh", cosh, 7.09089565712818057364E2, 4.49423283712885057274E307,
185 4.25947714184369757620E208, 4.49423283712885057274E307, 2},
186 {"cosh", cosh, 1.0, 1.54307556152343750000E0,
187 5.07329180627847790562E-6, 1.17520119364380145688E0, 2},
188 {"cosh", cosh, 0.5, 1.12762451171875000000E0,
189 1.45348763078522622516E-6, 5.21095305493747361622E-1, 2},
190 /* tanh in OpenBSD libm has less precession, adapt expectation from 2 to 3 */
191 {"tanh", tanh, 0.5, 4.62112426757812500000E-1,
192 4.73050219725850231848E-6, 7.86447732965927410150E-1, 3},
193 {"tanh", tanh, 5.49306144334054780032E-1, 4.99984741210937500000E-1,
194 1.52587890624507506378E-5, 7.50000000000000049249E-1, 2},
195 {"tanh", tanh, 0.625, 5.54595947265625000000E-1,
196 3.77508375729399903910E-6, 6.92419147969988069631E-1, 3},
197 {"asinh", asinh, 0.5, 4.81201171875000000000E-1,
198 1.06531846034474977589E-5, 8.94427190999915878564E-1, 2},
199 {"asinh", asinh, 1.0, 8.81362915039062500000E-1,
200 1.06719804805252326093E-5, 7.07106781186547524401E-1, 2},
201 {"asinh", asinh, 2.0, 1.44363403320312500000E0,
202 1.44197568534249327674E-6, 4.47213595499957939282E-1, 2},
203 {"acosh", acosh, 2.0, 1.31695556640625000000E0,
204 2.33051856670862504635E-6, 5.77350269189625764509E-1, 2},
205 {"acosh", acosh, 1.5, 9.62417602539062500000E-1,
206 6.04758014439499551783E-6, 8.94427190999915878564E-1, 2},
207 {"acosh", acosh, 1.03125, 2.49343872070312500000E-1,
208 9.62177257298785143908E-6, 3.96911150685467059809E0, 2},
209 {"atanh", atanh, 0.5, 5.49301147460937500000E-1,
210 4.99687311734569762262E-6, 1.33333333333333333333E0, 2},
211 #if 0
212 {"polylog_3", w_polylog_3, 0.875, 1.01392710208892822265625,
213 1.21106784918910854520955967e-8, 1.4149982649102631253520580, 2},
214 {"polylog_3", w_polylog_3, 8.0000000000000004440892e-1,
215 9.10605847835540771484375e-1, 7.570301035551561731571421485488e-9,
216 1.343493250010310486342647, 2},
217 #endif
218 #if 0
219 {"gamma", gamma, 1.0, 1.0,
220 0.0, -5.772156649015328606e-1, 2},
221 {"gamma", gamma, 2.0, 1.0,
222 0.0, 4.2278433509846713939e-1, 2},
223 {"gamma", gamma, 3.0, 2.0,
224 0.0, 1.845568670196934279, 2},
225 {"gamma", gamma, 4.0, 6.0,
226 0.0, 7.536706010590802836, 2},
227 #endif
228 {"null", NULL, 0.0, 0.0, 0.0, 2},
229 };
230
231 /* These take care of extra-precise floating point register problems. */
232 static volatile double volat1;
233 static volatile double volat2;
234
235
236 /* Return the next nearest floating point value to X
237 in the direction of UPDOWN (+1 or -1).
238 (Fails if X is denormalized.) */
239
240 static double
nextval(x,updown)241 nextval (x, updown)
242 double x;
243 int updown;
244 {
245 double m;
246 int i;
247
248 volat1 = x;
249 m = 0.25 * MACHEP * volat1 * updown;
250 volat2 = volat1 + m;
251 if (volat2 != volat1)
252 printf ("successor failed\n");
253
254 for (i = 2; i < 10; i++)
255 {
256 volat2 = volat1 + i * m;
257 if (volat1 != volat2)
258 return volat2;
259 }
260
261 printf ("nextval failed\n");
262 return volat1;
263 }
264
265
266
267
268 int
monot()269 monot ()
270 {
271 double (*fun1) (double);
272 int i, j, errs, tests;
273 double x, x0, y, dy, err;
274
275 /* Set math coprocessor to double precision. */
276 /* dprec (); */
277 errs = 0;
278 tests = 0;
279 i = 0;
280
281 for (;;)
282 {
283 fun1 = test1[i].func;
284 if (fun1 == NULL)
285 break;
286 volat1 = test1[i].arg1;
287 x0 = volat1;
288 x = volat1;
289 for (j = 0; j <= NPTS; j++)
290 {
291 volat1 = x - x0;
292 dy = volat1 * test1[i].derivative;
293 dy = test1[i].answer2 + dy;
294 volat1 = test1[i].answer1 + dy;
295 volat2 = (*(fun1)) (x);
296 if (volat2 != volat1)
297 {
298 /* Report difference between program result
299 and extended precision function value. */
300 err = volat2 - test1[i].answer1;
301 err = err - dy;
302 err = err / volat1;
303 if (fabs (err) > ((OKERROR + test1[i].thresh) * MACHEP))
304 {
305 printf ("%d %s(%.16e) = %.16e, rel err = %.3e\n",
306 j, test1[i].name, x, volat2, err);
307 errs += 1;
308 }
309 }
310 x = nextval (x, 1);
311 tests += 1;
312 }
313
314 x = x0;
315 x = nextval (x, -1);
316 for (j = 1; j < NPTS; j++)
317 {
318 volat1 = x - x0;
319 dy = volat1 * test1[i].derivative;
320 dy = test1[i].answer2 + dy;
321 volat1 = test1[i].answer1 + dy;
322 volat2 = (*(fun1)) (x);
323 if (volat2 != volat1)
324 {
325 err = volat2 - test1[i].answer1;
326 err = err - dy;
327 err = err / volat1;
328 if (fabs (err) > ((OKERROR + test1[i].thresh) * MACHEP))
329 {
330 printf ("%d %s(%.16e) = %.16e, rel err = %.3e\n",
331 j, test1[i].name, x, volat2, err);
332 errs += 1;
333 }
334 }
335 x = nextval (x, -1);
336 tests += 1;
337 }
338 i += 1;
339 }
340 printf ("%d errors in %d tests\n", errs, tests);
341 return (errs);
342 }
343