1 /*
2 * Copyright (c) 1994 Cygnus Support.
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms are permitted
6 * provided that the above copyright notice and this paragraph are
7 * duplicated in all such forms and that any documentation,
8 * and/or other materials related to such
9 * distribution and use acknowledge that the software was developed
10 * at Cygnus Support, Inc. Cygnus Support, Inc. may not be used to
11 * endorse or promote products derived from this software without
12 * specific prior written permission.
13 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
14 * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
15 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
16 */
17 /*
18 Test the library maths functions using trusted precomputed test
19 vectors.
20
21 These vectors were originally generated on a sun3 with a 68881 using
22 80 bit precision, but ...
23
24 Each function is called with a variety of interesting arguments.
25 Note that many of the polynomials we use behave badly when the
26 domain is stressed, so the numbers in the vectors depend on what is
27 useful to test - eg sin(1e30) is pointless - the arg has to be
28 reduced modulo pi, and after that there's no bits of significance
29 left to evaluate with - any number would be just as precise as any
30 other.
31
32
33 */
34
35 #include "test.h"
36 #include <math.h>
37 #include <float.h>
38 #include <errno.h>
39 #include <stdio.h>
40 #include <string.h>
41 #include <stdbool.h>
42
43 extern int inacc;
44
45 int merror;
46 double mretval = 64;
47 int traperror = 1;
48 char *mname;
49
50 void
translate_to(FILE * file,double r)51 translate_to (FILE *file,
52 double r)
53 {
54 __ieee_double_shape_type bits;
55 bits.value = r;
56 fprintf(file, "0x%08lx, 0x%08lx", (unsigned long) bits.parts.msw, (unsigned long) bits.parts.lsw);
57 }
58
59 /* Convert double to float, preserving issignaling status and NaN sign */
60
61 #define to_float(f,d) do { \
62 if (__issignaling(d)) { \
63 if (signbit(d)) \
64 (f) = -__builtin_nansf(""); \
65 else \
66 (f) = __builtin_nansf(""); \
67 } else if (isnan(d)) { \
68 if (signbit(d)) \
69 (f) = -__builtin_nanf(""); \
70 else \
71 (f) = __builtin_nanf(""); \
72 } else \
73 (f) = (float) (d); \
74 } while(0)
75
76 #define to_double(d,f) do { \
77 if (__issignalingf(f)) { \
78 if (signbit(f)) \
79 (d) = -__builtin_nans(""); \
80 else \
81 (d) = __builtin_nans(""); \
82 } else if (isnanf(f)) { \
83 if (signbit(f)) \
84 (d) = -__builtin_nan(""); \
85 else \
86 (d) = __builtin_nan(""); \
87 } else \
88 (d) = (double) (f); \
89 } while(0)
90
91 int
ffcheck_id(double is,one_line_type * p,char * name,int serrno,int merror,int id)92 ffcheck_id(double is,
93 one_line_type *p,
94 char *name,
95 int serrno,
96 int merror,
97 int id)
98 {
99 /* Make sure the answer isn't to far wrong from the correct value */
100 __ieee_double_shape_type correct, isbits;
101 int mag;
102 isbits.value = is;
103
104 (void) serrno;
105 (void) merror;
106 correct.parts.msw = p->qs[id].msw;
107 correct.parts.lsw = p->qs[id].lsw;
108
109 int error_bit = p->error_bit;
110
111 mag = mag_of_error(correct.value, is);
112
113 if (error_bit > 63)
114 error_bit = 63;
115
116 /* All signaling nans are the "same", as are all quiet nans.
117 * On i386, just returning a value converts signaling nans to quiet
118 * nans, let that pass
119 */
120 if (isnan(correct.value) && isnan(is)
121 #ifndef __i386__
122 && (issignaling(correct.value) == issignaling(is))
123 #endif
124 )
125 {
126 mag = 64;
127 }
128
129 if (mag < error_bit)
130 {
131 inacc ++;
132
133 printf("%s:%d, inaccurate answer: bit %d (%08lx%08lx %08lx%08lx) (%g %g)\n",
134 name, p->line, mag,
135 (unsigned long) correct.parts.msw,
136 (unsigned long) correct.parts.lsw,
137 (unsigned long) isbits.parts.msw,
138 (unsigned long) isbits.parts.lsw,
139 correct.value, is);
140 }
141
142 #if 0
143 if (p->qs[id].merror != merror)
144 {
145 /* Beware, matherr doesn't exist anymore. */
146 printf("testing %s_vec.c:%d, matherr wrong: %d %d\n",
147 name, p->line, merror, p->qs[id].merror);
148 }
149
150 if (p->qs[id].errno_val != errno)
151 {
152 printf("testing %s_vec.c:%d, errno wrong: %d %d\n",
153 name, p->line, errno, p->qs[id].errno_val);
154
155 }
156 #endif
157 return mag;
158 }
159
160 int
ffcheck(double is,one_line_type * p,char * name,int serrno,int merror)161 ffcheck(double is,
162 one_line_type *p,
163 char *name,
164 int serrno,
165 int merror)
166 {
167 return ffcheck_id(is, p, name, serrno, merror, 0);
168 }
169
170 int
fffcheck_id(float is,one_line_type * p,char * name,int serrno,int merror,int id)171 fffcheck_id (float is,
172 one_line_type *p,
173 char *name,
174 int serrno,
175 int merror,
176 int id)
177 {
178 /* Make sure the answer isn't to far wrong from the correct value */
179 __ieee_float_shape_type correct, isbits;
180 __ieee_double_shape_type correct_double;
181 __ieee_double_shape_type is_double;
182 int mag;
183
184 (void) serrno;
185 (void) merror;
186
187 isbits.value = is;
188 to_double(is_double.value, is);
189
190 correct_double.parts.msw = p->qs[id].msw;
191 correct_double.parts.lsw = p->qs[id].lsw;
192 to_float(correct.value, correct_double.value);
193 // printf("%s got 0x%08x want 0x%08x\n", name, isbits.p1, correct.p1);
194
195 int error_bit = p->error_bit;
196
197 if (error_bit > 31)
198 error_bit = 31;
199
200 mag = fmag_of_error(correct.value, is);
201
202 /* All signaling nans are the "same", as are all quiet nans */
203 if (isnan(correct.value) && isnan(is)
204 #ifndef __i386__
205 /* i386 calling convention ends up always converting snan into qnan */
206 && (__issignalingf(correct.value) == __issignalingf(is))
207 #endif
208 )
209 {
210 mag = 32;
211 }
212
213 if (mag < error_bit)
214 {
215 inacc ++;
216
217 printf("%s:%d, inaccurate answer: bit %d (%08lx %08lx) (%g %g) (is double 0x%08lx, 0x%08lx)\n",
218 name, p->line, mag,
219 (unsigned long) (uint32_t) correct.p1,
220 (unsigned long) (uint32_t) isbits.p1,
221 (double) correct.value,
222 (double) is,
223 (unsigned long) (uint32_t) is_double.parts.msw,
224 (unsigned long) (uint32_t) is_double.parts.lsw);
225 }
226
227 #if 0
228 if (p->qs[id].merror != merror)
229 {
230 /* Beware, matherr doesn't exist anymore. */
231 printf("testing %s_vec.c:%d, matherr wrong: %d %d\n",
232 name, p->line, merror, p->qs[id].merror);
233 }
234
235 if (p->qs[id].errno_val != errno)
236 {
237 printf("testing %s_vec.c:%d, errno wrong: %d %d\n",
238 name, p->line, errno, p->qs[id].errno_val);
239
240 }
241 #endif
242 return mag;
243 }
244
245 int
fffcheck(float is,one_line_type * p,char * name,int serrno,int merror)246 fffcheck (float is,
247 one_line_type *p,
248 char *name,
249 int serrno,
250 int merror)
251 {
252 return fffcheck_id(is, p, name, serrno, merror, 0);
253 }
254
255 double
thedouble(uint32_t msw,uint32_t lsw,double * r)256 thedouble (uint32_t msw,
257 uint32_t lsw, double *r)
258 {
259 __ieee_double_shape_type x;
260
261 x.parts.msw = msw;
262 x.parts.lsw = lsw;
263 if (r)
264 memcpy(r, &x.value, sizeof(double));
265 return x.value;
266 }
267
268 int calc;
269 extern int reduce;
270
271
272 void
frontline(FILE * f,int mag,one_line_type * p,double result,int merror,int my_errno,char * args,char * name)273 frontline (FILE *f,
274 int mag,
275 one_line_type *p,
276 double result,
277 int merror,
278 int my_errno,
279 char *args,
280 char *name)
281 {
282 (void) name;
283 /* float returns can never have more than 32 bits of accuracy */
284 if (*args == 'f' && mag > 32)
285 mag = 32;
286
287 if (reduce && p->error_bit < mag)
288 {
289 fprintf(f, "{%2d,", p->error_bit);
290 }
291 else
292 {
293 fprintf(f, "{%2d,",mag);
294 }
295
296
297 fprintf(f,"%2d,%3d,", merror,my_errno);
298 fprintf(f, "__LINE__, ");
299
300 if (calc)
301 {
302 translate_to(f, result);
303 }
304 else
305 {
306 translate_to(f, thedouble(p->qs[0].msw, p->qs[0].lsw, NULL));
307 }
308
309 fprintf(f, ", ");
310
311 fprintf(f,"0x%08lx, 0x%08lx", (unsigned long) p->qs[1].msw, (unsigned long) p->qs[1].lsw);
312
313
314 if (args[2])
315 {
316 fprintf(f, ", ");
317 fprintf(f,"0x%08lx, 0x%08lx", (unsigned long) p->qs[2].msw, (unsigned long) p->qs[2].lsw);
318 }
319
320 fprintf(f,"}, /* %g=f(%g",result,
321 thedouble(p->qs[1].msw, p->qs[1].lsw, NULL));
322
323 if (args[2])
324 {
325 fprintf(f,", %g", thedouble(p->qs[2].msw,p->qs[2].lsw, NULL));
326 }
327 fprintf(f, ")*/\n");
328 }
329
330 void
finish(FILE * f,int vector,double result,one_line_type * p,char * args,char * name)331 finish (FILE *f,
332 int vector,
333 double result,
334 one_line_type *p,
335 char *args,
336 char *name)
337 {
338 int mag;
339
340 mag = ffcheck(result, p,name, merror, errno);
341 if (vector)
342 {
343 frontline(f, mag, p, result, merror, errno, args , name);
344 }
345 }
346
347 void
finish2(FILE * f,int vector,double result,double result2,one_line_type * p,char * args,char * name)348 finish2 (FILE *f,
349 int vector,
350 double result,
351 double result2,
352 one_line_type *p,
353 char *args,
354 char *name)
355 {
356 int mag, mag2;
357
358 mag = ffcheck(result, p,name, merror, errno);
359 mag2 = ffcheck_id(result2, p,name, merror, errno, 2);
360 if (mag2 < mag)
361 mag = mag2;
362 if (vector)
363 {
364 __ieee_double_shape_type result2_double;
365 result2_double.value = result2;
366 p->qs[2].msw = result2_double.parts.msw;
367 p->qs[2].lsw = result2_double.parts.lsw;
368 frontline(f, mag, p, result, merror, errno, args , name);
369 }
370 }
371
372 void
ffinish(FILE * f,int vector,float fresult,one_line_type * p,char * args,char * name)373 ffinish (FILE *f,
374 int vector,
375 float fresult,
376 one_line_type *p,
377 char *args,
378 char *name)
379 {
380 int mag;
381
382 mag = fffcheck(fresult, p,name, merror, errno);
383 if (vector)
384 {
385 frontline(f, mag, p, (double) fresult, merror, errno, args , name);
386 }
387 }
388
389 void
ffinish2(FILE * f,int vector,float fresult,float fresult2,one_line_type * p,char * args,char * name)390 ffinish2 (FILE *f,
391 int vector,
392 float fresult,
393 float fresult2,
394 one_line_type *p,
395 char *args,
396 char *name)
397 {
398 int mag, mag2;
399
400 mag = fffcheck(fresult, p,name, merror, errno);
401 mag2 = fffcheck_id(fresult2, p, name, merror, errno, 2);
402 if (mag2 < mag)
403 mag = mag2;
404 if (vector)
405 {
406 __ieee_double_shape_type result2_double;
407 to_double(result2_double.value, (double) fresult2);
408 p->qs[2].msw = result2_double.parts.msw;
409 p->qs[2].lsw = result2_double.parts.lsw;
410 frontline(f, mag, p, (double) fresult, merror, errno, args , name);
411 }
412 }
413
414 extern int redo;
415
416 bool
in_float_range(double arg)417 in_float_range(double arg)
418 {
419 if (isinf(arg) || isnan(arg))
420 return true;
421 return -(double) FLT_MAX <= arg && arg < (double) FLT_MAX;
422 }
423
424 void
run_vector_1(int vector,one_line_type * p,char * func,char * name,char * args)425 run_vector_1 (int vector,
426 one_line_type *p,
427 char *func,
428 char *name,
429 char *args)
430 {
431 FILE *f = NULL;
432 double result;
433 float fresult;
434
435 #if 0
436 if (vector)
437 {
438
439 VECOPEN(name, f);
440
441 if (redo)
442 {
443 double k;
444 union {
445 struct { uint32_t a, b; };
446 double d;
447 } d, d4;
448
449 for (k = -.2; k < .2; k+= 0.00132)
450 {
451 d.d = k;
452 d4.d = k + 4;
453 fprintf(f,"{1,1, 1,1, 0,0,0x%08lx,0x%08lx, 0x%08lx, 0x%08lx},\n",
454 (unsigned long) d.a, (unsigned long) d.b, (unsigned long) d4.a, (unsigned long) d4.b);
455
456 }
457
458 for (k = -1.2; k < 1.2; k+= 0.01)
459 {
460 d.d = k;
461 d4.d = k + 4;
462 fprintf(f,"{1,1, 1,1, 0,0,0x%08lx,0x%08lx, 0x%08lx, 0x%08lx},\n",
463 (unsigned long) d.a, (unsigned long) d.b, (unsigned long) d4.a, (unsigned long) d4.b);
464
465 }
466 for (k = -M_PI *2; k < M_PI *2; k+= M_PI/2)
467 {
468 d.d = k;
469 d4.d = k + 4;
470 fprintf(f,"{1,1, 1,1, 0,0,0x%08lx,0x%08lx, 0x%08lx, 0x%08lx},\n",
471 (unsigned long) d.a, (unsigned long) d.b, (unsigned long) d4.a, (unsigned long) d4.b);
472
473 }
474
475 for (k = -30; k < 30; k+= 1.7)
476 {
477 d.d = k;
478 d4.d = k + 4;
479 fprintf(f,"{2,2, 1,1, 0,0, 0x%08lx,0x%08lx, 0x%08lx, 0x%08lx},\n",
480 (unsigned long) d.a, (unsigned long) d.b, (unsigned long) d4.a, (unsigned long) d4.b);
481
482 }
483 VECCLOSE(f, name, args);
484 return;
485 }
486 }
487 #endif
488
489 newfunc(name);
490 while (p->line)
491 {
492 double arg1;
493 thedouble(p->qs[1].msw, p->qs[1].lsw, &arg1);
494 double arg2;
495 thedouble(p->qs[2].msw, p->qs[2].lsw, &arg2);
496
497 errno = 0;
498 merror = 0;
499 mname = 0;
500
501
502 line(p->line);
503
504 merror = 0;
505 errno = 123;
506
507 if (strcmp(args,"dd")==0)
508 {
509 typedef double (*pdblfunc) (double);
510
511 /* Double function returning a double */
512
513 result = ((pdblfunc)(func))(arg1);
514
515 finish(f,vector, result, p, args, name);
516 }
517 else if (strcmp(args,"ff")==0)
518 {
519 float arga;
520
521 typedef float (*pdblfunc) (float);
522
523 /* Double function returning a double */
524
525 if (in_float_range(arg1))
526 {
527 to_float(arga, arg1);
528 fresult = ((pdblfunc)(func))(arga);
529 ffinish(f, vector, fresult, p,args, name);
530 }
531 }
532 else if (strcmp(args,"ddd")==0)
533 {
534 typedef double (*pdblfunc) (double,double);
535
536 result = ((pdblfunc)(func))(arg1,arg2);
537 finish(f, vector, result, p,args, name);
538 }
539 else if (strcmp(args,"fff")==0)
540 {
541 float arga;
542 float argb;
543
544 typedef float (*pdblfunc) (float,float);
545
546 if (in_float_range(arg1) && in_float_range(arg2))
547 {
548 to_float(arga, arg1);
549 to_float(argb, arg2);
550 fresult = ((pdblfunc)(func))(arga, argb);
551 // printf("0x%08x = %s(0x%08x, 0x%08x)\n",
552 // float_bits(fresult), name, float_bits(arga), float_bits(argb));
553 ffinish(f, vector, fresult, p,args, name);
554 }
555 }
556 else if (strcmp(args,"did")==0)
557 {
558 typedef double (*pdblfunc) (int,double);
559
560 result = ((pdblfunc)(func))((int)arg1,arg2);
561 finish(f, vector, result, p,args, name);
562 }
563 else if (strcmp(args,"fif")==0)
564 {
565 float argb;
566
567 typedef float (*pdblfunc) (int,float);
568
569
570 if (in_float_range(arg2))
571 {
572 to_float(argb, arg2);
573 fresult = ((pdblfunc)(func))((int)arg1, argb);
574 ffinish(f, vector, fresult, p,args, name);
575 }
576 }
577 else if (strcmp(args,"id")==0)
578 {
579 typedef int (*pdblfunc)(double);
580
581 result = (double) ((pdblfunc)(func))(arg1);
582 finish(f, vector, result, p, args, name);
583 }
584 else if (strcmp(args,"ddp") == 0)
585 {
586 typedef double (*pdblfunc)(double, double *);
587 double result2;
588
589 result = (double) ((pdblfunc)(func))(arg1, &result2);
590 finish2(f, vector, result, result2, p, args, name);
591 }
592 else if (strcmp(args,"ffp") == 0)
593 {
594 typedef float (*pdblfunc)(float, float *);
595 float result2;
596
597 fresult = (float) ((pdblfunc)(func))(arg1, &result2);
598 ffinish2(f, vector, fresult, result2, p, args, name);
599 }
600 else if (strcmp(args,"ddi") == 0)
601 {
602 typedef double (*pdblfunc)(double, int);
603
604 result = ((pdblfunc)func) (arg1, (int) arg2);
605 finish(f, vector, result, p, args, name);
606 }
607 else if (strcmp(args,"ffi") == 0)
608 {
609 typedef float (*pdblfunc)(float, int);
610
611 fresult = ((pdblfunc)func) (arg1, (int) arg2);
612 ffinish(f, vector, fresult, p, args, name);
613 }
614 p++;
615 }
616 if (vector && f)
617 {
618 VECCLOSE(f, name, args);
619 }
620 }
621
622 void
test_math(int vector)623 test_math (int vector)
624 {
625 test_acos(vector);
626 test_acosf(vector);
627 test_acosh(vector);
628 test_acoshf(vector);
629 test_asin(vector);
630 test_asinf(vector);
631 test_asinh(vector);
632 test_asinhf(vector);
633 test_atan(vector);
634 test_atan2(vector);
635 test_atan2f(vector);
636 test_atanf(vector);
637 test_atanh(vector);
638 test_atanhf(vector);
639 test_ceil(vector);
640 test_ceilf(vector);
641 test_copysign(vector);
642 test_copysignf(vector);
643 test_cos(vector);
644 test_cosf(vector);
645 test_cosh(vector);
646 test_coshf(vector);
647 test_erf(vector);
648 test_erfc(vector);
649 test_erfcf(vector);
650 test_erff(vector);
651 test_exp(vector);
652 test_expf(vector);
653 test_fabs(vector);
654 test_fabsf(vector);
655 test_floor(vector);
656 test_floorf(vector);
657 test_fmod(vector);
658 test_fmodf(vector);
659 test_gamma(vector);
660 test_gammaf(vector);
661 test_hypot(vector);
662 test_hypotf(vector);
663 test_issignaling(vector);
664 test_j0(vector);
665 test_j0f(vector);
666 test_j1(vector);
667 test_j1f(vector);
668 test_jn(vector);
669 test_jnf(vector);
670 test_log(vector);
671 test_log10(vector);
672 test_log10f(vector);
673 test_log1p(vector);
674 test_log1pf(vector);
675 test_log2(vector);
676 test_log2f(vector);
677 test_logf(vector);
678 test_modf(vector);
679 test_modff(vector);
680 test_pow_vec(vector);
681 test_powf_vec(vector);
682 test_scalb(vector);
683 test_scalbf(vector);
684 test_scalbn(vector);
685 test_scalbnf(vector);
686 test_sin(vector);
687 test_sinf(vector);
688 test_sinh(vector);
689 test_sinhf(vector);
690 test_sqrt(vector);
691 test_sqrtf(vector);
692 test_tan(vector);
693 test_tanf(vector);
694 test_tanh(vector);
695 test_tanhf(vector);
696 test_trunc(vector);
697 test_truncf(vector);
698 test_y0(vector);
699 test_y0f(vector);
700 test_y1(vector);
701 test_y1f(vector);
702 test_y1f(vector);
703 test_ynf(vector);
704 }
705
706 /* These have to be played with to get to compile on machines which
707 don't have the fancy <foo>f entry points
708 */
709
710 #if 0
711 float cosf (float a) { return cos((double)a); }
712 float sinf (float a) { return sin((double)a); }
713 float log1pf (float a) { return log1p((double)a); }
714 float tanf (float a) { return tan((double)a); }
715 float ceilf (float a) { return ceil(a); }
716 float floorf (float a) { return floor(a); }
717 #endif
718
719 /*ndef HAVE_FLOAT*/
720 #if 0
721
722 float fmodf(a,b) float a,b; { return fmod(a,b); }
723 float hypotf(a,b) float a,b; { return hypot(a,b); }
724
725 float acosf(a) float a; { return acos(a); }
726 float acoshf(a) float a; { return acosh(a); }
727 float asinf(a) float a; { return asin(a); }
728 float asinhf(a) float a; { return asinh(a); }
729 float atanf(a) float a; { return atan(a); }
730 float atanhf(a) float a; { return atanh(a); }
731
732 float coshf(a) float a; { return cosh(a); }
733 float erff(a) float a; { return erf(a); }
734 float erfcf(a) float a; { return erfc(a); }
735 float expf(a) float a; { return exp(a); }
736 float fabsf(a) float a; { return fabs(a); }
737
738 float gammaf(a) float a; { return gamma(a); }
739 float j0f(a) float a; { return j0(a); }
740 float j1f(a) float a; { return j1(a); }
741 float log10f(a) float a; { return log10(a); }
742
743 float logf(a) float a; { return log(a); }
744
745 float sinhf(a) float a; { return sinh(a); }
746 float sqrtf(a) float a; { return sqrt(a); }
747
748 float tanhf(a) float a; { return tanh(a); }
749 float y0f(a) float a; { return y0(a); }
750 float y1f(a) float a; { return y1(a); }
751 #endif
752
