[444] | 1 | /* |
---|
| 2 | Test the library maths functions using trusted precomputed test |
---|
| 3 | vectors. |
---|
| 4 | |
---|
| 5 | These vectors were originally generated on a sun3 with a 68881 using |
---|
| 6 | 80 bit precision, but ... |
---|
| 7 | |
---|
| 8 | Each function is called with a variety of interesting arguments. |
---|
| 9 | Note that many of the polynomials we use behave badly when the |
---|
| 10 | domain is stressed, so the numbers in the vectors depend on what is |
---|
| 11 | useful to test - eg sin(1e30) is pointless - the arg has to be |
---|
| 12 | reduced modulo pi, and after that there's no bits of significance |
---|
| 13 | left to evaluate with - any number would be just as precise as any |
---|
| 14 | other. |
---|
| 15 | |
---|
| 16 | |
---|
| 17 | */ |
---|
| 18 | |
---|
| 19 | #include "test.h" |
---|
| 20 | #include <math.h> |
---|
| 21 | #include <ieeefp.h> |
---|
| 22 | #include <float.h> |
---|
| 23 | #include <math.h> |
---|
| 24 | #include <errno.h> |
---|
| 25 | #include <stdio.h> |
---|
| 26 | |
---|
| 27 | int inacc; |
---|
| 28 | |
---|
| 29 | int merror; |
---|
| 30 | double mretval = 64; |
---|
| 31 | int traperror = 1; |
---|
| 32 | char *mname; |
---|
| 33 | |
---|
| 34 | int verbose; |
---|
| 35 | |
---|
| 36 | /* To test exceptions - we trap them all and return a known value */ |
---|
| 37 | int |
---|
| 38 | matherr (struct exception *e) |
---|
| 39 | { |
---|
| 40 | if (traperror) |
---|
| 41 | { |
---|
| 42 | merror = e->type + 12; |
---|
| 43 | mname = e->name; |
---|
| 44 | e->retval = mretval; |
---|
| 45 | errno = merror + 24; |
---|
| 46 | return 1; |
---|
| 47 | } |
---|
| 48 | return 0; |
---|
| 49 | } |
---|
| 50 | |
---|
| 51 | |
---|
| 52 | void translate_to (FILE *file, |
---|
| 53 | double r) |
---|
| 54 | { |
---|
| 55 | __ieee_double_shape_type bits; |
---|
| 56 | bits.value = r; |
---|
| 57 | fprintf(file, "0x%08x, 0x%08x", bits.parts.msw, bits.parts.lsw); |
---|
| 58 | } |
---|
| 59 | |
---|
| 60 | int |
---|
| 61 | ffcheck (double is, |
---|
| 62 | one_line_type *p, |
---|
| 63 | char *name, |
---|
| 64 | int serrno, |
---|
| 65 | int merror) |
---|
| 66 | { |
---|
| 67 | /* Make sure the answer isn't to far wrong from the correct value */ |
---|
| 68 | __ieee_double_shape_type correct, isbits; |
---|
| 69 | int mag; |
---|
| 70 | isbits.value = is; |
---|
| 71 | |
---|
| 72 | correct.parts.msw = p->qs[0].msw; |
---|
| 73 | correct.parts.lsw = p->qs[0].lsw; |
---|
| 74 | |
---|
| 75 | mag = mag_of_error(correct.value, is); |
---|
| 76 | |
---|
| 77 | if (mag < p->error_bit) |
---|
| 78 | { |
---|
| 79 | inacc ++; |
---|
| 80 | |
---|
| 81 | printf("%s:%d, inaccurate answer: bit %d (%08x%08x %08x%08x) (%g %g)\n", |
---|
| 82 | name, p->line, mag, |
---|
| 83 | correct.parts.msw, |
---|
| 84 | correct.parts.lsw, |
---|
| 85 | isbits.parts.msw, |
---|
| 86 | isbits.parts.lsw, |
---|
| 87 | correct.value, is); |
---|
| 88 | } |
---|
| 89 | |
---|
| 90 | #if 0 |
---|
| 91 | if (p->qs[0].merror != merror) |
---|
| 92 | { |
---|
| 93 | printf("testing %s_vec.c:%d, matherr wrong: %d %d\n", |
---|
| 94 | name, p->line, merror, p->qs[0].merror); |
---|
| 95 | } |
---|
| 96 | |
---|
| 97 | if (p->qs[0].errno_val != errno) |
---|
| 98 | { |
---|
| 99 | printf("testing %s_vec.c:%d, errno wrong: %d %d\n", |
---|
| 100 | name, p->line, errno, p->qs[0].errno_val); |
---|
| 101 | |
---|
| 102 | } |
---|
| 103 | #endif |
---|
| 104 | return mag; |
---|
| 105 | } |
---|
| 106 | |
---|
| 107 | double |
---|
| 108 | thedouble (long msw, |
---|
| 109 | long lsw) |
---|
| 110 | { |
---|
| 111 | __ieee_double_shape_type x; |
---|
| 112 | |
---|
| 113 | x.parts.msw = msw; |
---|
| 114 | x.parts.lsw = lsw; |
---|
| 115 | return x.value; |
---|
| 116 | } |
---|
| 117 | |
---|
| 118 | int calc; |
---|
| 119 | int reduce; |
---|
| 120 | |
---|
| 121 | |
---|
| 122 | frontline (FILE *f, |
---|
| 123 | int mag, |
---|
| 124 | one_line_type *p, |
---|
| 125 | double result, |
---|
| 126 | int merror, |
---|
| 127 | int errno, |
---|
| 128 | char *args, |
---|
| 129 | char *name) |
---|
| 130 | { |
---|
| 131 | if (reduce && p->error_bit < mag) |
---|
| 132 | { |
---|
| 133 | fprintf(f, "{%2d,", p->error_bit); |
---|
| 134 | } |
---|
| 135 | else |
---|
| 136 | { |
---|
| 137 | fprintf(f, "{%2d,",mag); |
---|
| 138 | } |
---|
| 139 | |
---|
| 140 | |
---|
| 141 | fprintf(f,"%2d,%3d,", merror,errno); |
---|
| 142 | fprintf(f, "__LINE__, "); |
---|
| 143 | |
---|
| 144 | if (calc) |
---|
| 145 | { |
---|
| 146 | translate_to(f, result); |
---|
| 147 | } |
---|
| 148 | else |
---|
| 149 | { |
---|
| 150 | translate_to(f, thedouble(p->qs[0].msw, p->qs[0].lsw)); |
---|
| 151 | } |
---|
| 152 | |
---|
| 153 | fprintf(f, ", "); |
---|
| 154 | |
---|
| 155 | fprintf(f,"0x%08x, 0x%08x", p->qs[1].msw, p->qs[1].lsw); |
---|
| 156 | |
---|
| 157 | |
---|
| 158 | if (args[2]) |
---|
| 159 | { |
---|
| 160 | fprintf(f, ", "); |
---|
| 161 | fprintf(f,"0x%08x, 0x%08x", p->qs[2].msw, p->qs[2].lsw); |
---|
| 162 | } |
---|
| 163 | |
---|
| 164 | fprintf(f,"}, /* %g=f(%g",result, |
---|
| 165 | thedouble(p->qs[1].msw, p->qs[1].lsw)); |
---|
| 166 | |
---|
| 167 | if (args[2]) |
---|
| 168 | { |
---|
| 169 | fprintf(f,", %g", thedouble(p->qs[2].msw,p->qs[2].lsw)); |
---|
| 170 | } |
---|
| 171 | fprintf(f, ")*/\n"); |
---|
| 172 | } |
---|
| 173 | |
---|
| 174 | finish (FILE *f, |
---|
| 175 | int vector, |
---|
| 176 | double result, |
---|
| 177 | one_line_type *p, |
---|
| 178 | char *args, |
---|
| 179 | char *name) |
---|
| 180 | { |
---|
| 181 | int mag; |
---|
| 182 | |
---|
| 183 | mag = ffcheck(result, p,name, merror, errno); |
---|
| 184 | if (vector) |
---|
| 185 | { |
---|
| 186 | frontline(f, mag, p, result, merror, errno, args , name); |
---|
| 187 | } |
---|
| 188 | } |
---|
| 189 | int redo; |
---|
| 190 | |
---|
| 191 | run_vector_1 (int vector, |
---|
| 192 | one_line_type *p, |
---|
| 193 | char *func, |
---|
| 194 | char *name, |
---|
| 195 | char *args) |
---|
| 196 | { |
---|
| 197 | FILE *f; |
---|
| 198 | int mag; |
---|
| 199 | double result; |
---|
| 200 | |
---|
| 201 | if (vector) |
---|
| 202 | { |
---|
| 203 | |
---|
| 204 | VECOPEN(name, f); |
---|
| 205 | |
---|
| 206 | if (redo) |
---|
| 207 | { |
---|
| 208 | double k; |
---|
| 209 | |
---|
| 210 | for (k = -.2; k < .2; k+= 0.00132) |
---|
| 211 | { |
---|
| 212 | |
---|
| 213 | fprintf(f,"{1,1, 1,1, 0,0,0x%08x,0x%08x, 0x%08x, 0x%08x},\n", |
---|
| 214 | k,k+4); |
---|
| 215 | |
---|
| 216 | } |
---|
| 217 | |
---|
| 218 | for (k = -1.2; k < 1.2; k+= 0.01) |
---|
| 219 | { |
---|
| 220 | |
---|
| 221 | fprintf(f,"{1,1, 1,1, 0,0,0x%08x,0x%08x, 0x%08x, 0x%08x},\n", |
---|
| 222 | k,k+4); |
---|
| 223 | |
---|
| 224 | } |
---|
| 225 | for (k = -M_PI *2; k < M_PI *2; k+= M_PI/2) |
---|
| 226 | { |
---|
| 227 | |
---|
| 228 | fprintf(f,"{1,1, 1,1, 0,0,0x%08x,0x%08x, 0x%08x, 0x%08x},\n", |
---|
| 229 | k,k+4); |
---|
| 230 | |
---|
| 231 | } |
---|
| 232 | |
---|
| 233 | for (k = -30; k < 30; k+= 1.7) |
---|
| 234 | { |
---|
| 235 | |
---|
| 236 | fprintf(f,"{2,2, 1,1, 0,0, 0x%08x,0x%08x, 0x%08x, 0x%08x},\n", |
---|
| 237 | k,k+4); |
---|
| 238 | |
---|
| 239 | } |
---|
| 240 | VECCLOSE(f, name, args); |
---|
| 241 | return; |
---|
| 242 | } |
---|
| 243 | } |
---|
| 244 | |
---|
| 245 | newfunc(name); |
---|
| 246 | while (p->line) |
---|
| 247 | { |
---|
| 248 | double arg1 = thedouble(p->qs[1].msw, p->qs[1].lsw); |
---|
| 249 | double arg2 = thedouble(p->qs[2].msw, p->qs[2].lsw); |
---|
| 250 | |
---|
| 251 | double r; |
---|
| 252 | double rf; |
---|
| 253 | |
---|
| 254 | errno = 0; |
---|
| 255 | merror = 0; |
---|
| 256 | mname = 0; |
---|
| 257 | |
---|
| 258 | |
---|
| 259 | line(p->line); |
---|
| 260 | |
---|
| 261 | merror = 0; |
---|
| 262 | errno = 123; |
---|
| 263 | |
---|
| 264 | if (strcmp(args,"dd")==0) |
---|
| 265 | { |
---|
| 266 | typedef double (*pdblfunc) (double); |
---|
| 267 | |
---|
| 268 | /* Double function returning a double */ |
---|
| 269 | |
---|
| 270 | result = ((pdblfunc)(func))(arg1); |
---|
| 271 | |
---|
| 272 | finish(f,vector, result, p, args, name); |
---|
| 273 | } |
---|
| 274 | else if (strcmp(args,"ff")==0) |
---|
| 275 | { |
---|
| 276 | float arga; |
---|
| 277 | double a; |
---|
| 278 | |
---|
| 279 | typedef float (*pdblfunc) (float); |
---|
| 280 | |
---|
| 281 | /* Double function returning a double */ |
---|
| 282 | |
---|
| 283 | if (arg1 < FLT_MAX ) |
---|
| 284 | { |
---|
| 285 | arga = arg1; |
---|
| 286 | result = ((pdblfunc)(func))(arga); |
---|
| 287 | finish(f, vector, result, p,args, name); |
---|
| 288 | } |
---|
| 289 | } |
---|
| 290 | else if (strcmp(args,"ddd")==0) |
---|
| 291 | { |
---|
| 292 | typedef double (*pdblfunc) (double,double); |
---|
| 293 | |
---|
| 294 | result = ((pdblfunc)(func))(arg1,arg2); |
---|
| 295 | finish(f, vector, result, p,args, name); |
---|
| 296 | } |
---|
| 297 | else if (strcmp(args,"fff")==0) |
---|
| 298 | { |
---|
| 299 | double a,b; |
---|
| 300 | |
---|
| 301 | float arga; |
---|
| 302 | float argb; |
---|
| 303 | |
---|
| 304 | typedef float (*pdblfunc) (float,float); |
---|
| 305 | |
---|
| 306 | |
---|
| 307 | if (arg1 < FLT_MAX && arg2 < FLT_MAX) |
---|
| 308 | { |
---|
| 309 | arga = arg1; |
---|
| 310 | argb = arg2; |
---|
| 311 | result = ((pdblfunc)(func))(arga, argb); |
---|
| 312 | finish(f, vector, result, p,args, name); |
---|
| 313 | } |
---|
| 314 | } |
---|
| 315 | else if (strcmp(args,"did")==0) |
---|
| 316 | { |
---|
| 317 | typedef double (*pdblfunc) (int,double); |
---|
| 318 | |
---|
| 319 | result = ((pdblfunc)(func))((int)arg1,arg2); |
---|
| 320 | finish(f, vector, result, p,args, name); |
---|
| 321 | } |
---|
| 322 | else if (strcmp(args,"fif")==0) |
---|
| 323 | { |
---|
| 324 | double a,b; |
---|
| 325 | |
---|
| 326 | float arga; |
---|
| 327 | float argb; |
---|
| 328 | |
---|
| 329 | typedef float (*pdblfunc) (int,float); |
---|
| 330 | |
---|
| 331 | |
---|
| 332 | if (arg1 < FLT_MAX && arg2 < FLT_MAX) |
---|
| 333 | { |
---|
| 334 | arga = arg1; |
---|
| 335 | argb = arg2; |
---|
| 336 | result = ((pdblfunc)(func))((int)arga, argb); |
---|
| 337 | finish(f, vector, result, p,args, name); |
---|
| 338 | } |
---|
| 339 | } |
---|
| 340 | |
---|
| 341 | p++; |
---|
| 342 | } |
---|
| 343 | if (vector) |
---|
| 344 | { |
---|
| 345 | VECCLOSE(f, name, args); |
---|
| 346 | } |
---|
| 347 | } |
---|
| 348 | |
---|
| 349 | void |
---|
| 350 | test_math (void) |
---|
| 351 | { |
---|
| 352 | test_acos(0); |
---|
| 353 | test_acosf(0); |
---|
| 354 | test_acosh(0); |
---|
| 355 | test_acoshf(0); |
---|
| 356 | test_asin(0); |
---|
| 357 | test_asinf(0); |
---|
| 358 | test_asinh(0); |
---|
| 359 | test_asinhf(0); |
---|
| 360 | test_atan(0); |
---|
| 361 | test_atan2(0); |
---|
| 362 | test_atan2f(0); |
---|
| 363 | test_atanf(0); |
---|
| 364 | test_atanh(0); |
---|
| 365 | test_atanhf(0); |
---|
| 366 | test_ceil(0); |
---|
| 367 | test_ceilf(0); |
---|
| 368 | test_cos(0); |
---|
| 369 | test_cosf(0); |
---|
| 370 | test_cosh(0); |
---|
| 371 | test_coshf(0); |
---|
| 372 | test_erf(0); |
---|
| 373 | test_erfc(0); |
---|
| 374 | test_erfcf(0); |
---|
| 375 | test_erff(0); |
---|
| 376 | test_exp(0); |
---|
| 377 | test_expf(0); |
---|
| 378 | test_fabs(0); |
---|
| 379 | test_fabsf(0); |
---|
| 380 | test_floor(0); |
---|
| 381 | test_floorf(0); |
---|
| 382 | test_fmod(0); |
---|
| 383 | test_fmodf(0); |
---|
| 384 | test_gamma(0); |
---|
| 385 | test_gammaf(0); |
---|
| 386 | test_hypot(0); |
---|
| 387 | test_hypotf(0); |
---|
| 388 | test_j0(0); |
---|
| 389 | test_j0f(0); |
---|
| 390 | test_j1(0); |
---|
| 391 | test_j1f(0); |
---|
| 392 | test_jn(0); |
---|
| 393 | test_jnf(0); |
---|
| 394 | test_log(0); |
---|
| 395 | test_log10(0); |
---|
| 396 | test_log10f(0); |
---|
| 397 | test_log1p(0); |
---|
| 398 | test_log1pf(0); |
---|
| 399 | test_log2(0); |
---|
| 400 | test_log2f(0); |
---|
| 401 | test_logf(0); |
---|
| 402 | test_sin(0); |
---|
| 403 | test_sinf(0); |
---|
| 404 | test_sinh(0); |
---|
| 405 | test_sinhf(0); |
---|
| 406 | test_sqrt(0); |
---|
| 407 | test_sqrtf(0); |
---|
| 408 | test_tan(0); |
---|
| 409 | test_tanf(0); |
---|
| 410 | test_tanh(0); |
---|
| 411 | test_tanhf(0); |
---|
| 412 | test_y0(0); |
---|
| 413 | test_y0f(0); |
---|
| 414 | test_y1(0); |
---|
| 415 | test_y1f(0); |
---|
| 416 | test_y1f(0); |
---|
| 417 | test_ynf(0); |
---|
| 418 | } |
---|
| 419 | |
---|
| 420 | /* These have to be played with to get to compile on machines which |
---|
| 421 | don't have the fancy <foo>f entry points |
---|
| 422 | */ |
---|
| 423 | |
---|
| 424 | #if 0 |
---|
| 425 | float cosf (float a) { return cos((double)a); } |
---|
| 426 | float sinf (float a) { return sin((double)a); } |
---|
| 427 | float log1pf (float a) { return log1p((double)a); } |
---|
| 428 | float tanf (float a) { return tan((double)a); } |
---|
| 429 | float ceilf (float a) { return ceil(a); } |
---|
| 430 | float floorf (float a) { return floor(a); } |
---|
| 431 | #endif |
---|
| 432 | |
---|
| 433 | /*ndef HAVE_FLOAT*/ |
---|
| 434 | #if 0 |
---|
| 435 | |
---|
| 436 | float fmodf(a,b) float a,b; { return fmod(a,b); } |
---|
| 437 | float hypotf(a,b) float a,b; { return hypot(a,b); } |
---|
| 438 | |
---|
| 439 | float acosf(a) float a; { return acos(a); } |
---|
| 440 | float acoshf(a) float a; { return acosh(a); } |
---|
| 441 | float asinf(a) float a; { return asin(a); } |
---|
| 442 | float asinhf(a) float a; { return asinh(a); } |
---|
| 443 | float atanf(a) float a; { return atan(a); } |
---|
| 444 | float atanhf(a) float a; { return atanh(a); } |
---|
| 445 | |
---|
| 446 | float coshf(a) float a; { return cosh(a); } |
---|
| 447 | float erff(a) float a; { return erf(a); } |
---|
| 448 | float erfcf(a) float a; { return erfc(a); } |
---|
| 449 | float expf(a) float a; { return exp(a); } |
---|
| 450 | float fabsf(a) float a; { return fabs(a); } |
---|
| 451 | |
---|
| 452 | float gammaf(a) float a; { return gamma(a); } |
---|
| 453 | float j0f(a) float a; { return j0(a); } |
---|
| 454 | float j1f(a) float a; { return j1(a); } |
---|
| 455 | float log10f(a) float a; { return log10(a); } |
---|
| 456 | |
---|
| 457 | float logf(a) float a; { return log(a); } |
---|
| 458 | |
---|
| 459 | float sinhf(a) float a; { return sinh(a); } |
---|
| 460 | float sqrtf(a) float a; { return sqrt(a); } |
---|
| 461 | |
---|
| 462 | float tanhf(a) float a; { return tanh(a); } |
---|
| 463 | float y0f(a) float a; { return y0(a); } |
---|
| 464 | float y1f(a) float a; { return y1(a); } |
---|
| 465 | #endif |
---|