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 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 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