1 /* $OpenBSD: trailer.h,v 1.1 2009/04/09 01:24:43 martynas Exp $ */
2
3 /*
4 * Copyright (c) 2009 Gaston H. Gonnet <gonnet@inf.ethz.ch>
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 static struct point { double arg, res, err1, err2; } points[N+1];
20
21 #include <stdio.h>
22 #include <math.h>
23 #ifdef FORCE_FPU_DOUBLE
24 #include <fpu_control.h>
25 #endif
26
27 #define SHOW 20
28
29 int
Fn(FILE * out)30 Fn(FILE *out) {
31 struct point *by_err1[N], *by_err2[N], *q;
32 struct input_point *p;
33 double t, t1, t2;
34 int i, tot1={0}, tot2={0}, tot3={0}, tot4={0}, tot5={0};
35
36 #ifdef FORCE_FPU_DOUBLE
37 #define FPU_CW (_FPU_MASK_IM+_FPU_MASK_ZM+_FPU_MASK_OM+_FPU_MASK_UM+ \
38 _FPU_MASK_PM+_FPU_MASK_DM+_FPU_DOUBLE+_FPU_RC_NEAREST)
39
40 fpu_control_t cw = FPU_CW;
41 _FPU_SETCW (cw);
42 #endif
43
44 /* compute the results and store them in points structure */
45 for( p=input_points, q=points; q-points<N; p++,q++ ) {
46 q->arg = scalb(p->arg_m,p->arg_e);
47 q->res = F( q->arg );
48 }
49
50 for( i=0; i<N; i++ ) by_err1[i] = by_err2[i] = points+i;
51
52 /* do something silly to avoid loop merging optimization */
53 q->res = points[N/2].res+points[N/3].res+points[N/4].res;
54
55 /* Now compute the error directly, for the benefit of more precise hw */
56 for( p=input_points, q=points; q-points<N; p++,q++ ) {
57 if( p->val_e >= DBL_MAX_EXP-2 ) {
58 t1 = scalb(1.0,p->val_e/2);
59 t2 = scalb(1.0,p->val_e-p->val_e/2);
60 t = F(q->arg)*t1*t2 - p->val;
61 }
62 else {
63 t1 = scalb(1.0,p->val_e);
64 t = F(q->arg)*t1 - p->val;
65 }
66 t -= p->eps;
67 q->err2 = fabs(t);
68 t = scalb(q->res,p->val_e) - p->val;
69 t -= p->eps;
70 q->err1 = fabs(t);
71 }
72
73 /* Sort by errors in decreasing order */
74 #define KEY(x) (-(x)->err1)
75 SORT(by_err1,0,N-1,q);
76 #undef KEY
77 #define KEY(x) (-(x)->err2)
78 SORT(by_err2,0,N-1,q);
79
80 /* count the number of differences in errors */
81 for( q=points; q-points < N; q++ ) {
82 if( q->err1 > q->err2 ) tot1++;
83 else if( 2*q->err1 < q->err2 ) tot5++;
84 else if( 1.1*q->err1 < q->err2 ) tot4++;
85 else if( 1.01*q->err1 < q->err2 ) tot3++;
86 else if( q->err1 < q->err2 ) tot2++;
87 }
88
89 printf( "result of %s is ", Fs );
90 if( tot1==0 ) printf( "never more precise than double\n\n" );
91 else printf( "more precise than double %d out of %d times\n\n", tot1, N );
92
93 if( tot2+tot3+tot4+tot5 )
94 printf( "%d errors <= 1%% worse in accum, %d <= 10%%,"
95 "%d <= 100%%, %d > 100%%\n\n", tot2, tot3, tot4, tot5 );
96
97 for( i=N-1; i>=0 && by_err2[i]->err2 == 0; i-- );
98 if( N-1-i > 0 )
99 printf( "%d results were exact to double the precision\n\n", N-1-i );
100
101 if( tot1 > 10 ) {
102 printf( "%d largest ulp errors (from result in accumulator)\n", SHOW );
103 for( i=0; i<SHOW; i++ )
104 printf( " %.5f ulp for %s(%.18g) = %.18g)\n",
105 by_err2[i]->err2, Fs, by_err2[i]->arg, by_err2[i]->res );
106 printf( "\n" );
107 }
108
109 printf( "%d largest ulp errors (stored in a double)\n", SHOW );
110 for( i=0; i<SHOW; i++ )
111 printf( " %.5f ulp for %s(%.18g) = %.18g)\n",
112 by_err1[i]->err1, Fs, by_err1[i]->arg, by_err1[i]->res );
113 printf( "\n" );
114
115 fprintf( out, "%8s %5d %27.5f %26.18g %25.18g\n", Fs, N,
116 by_err1[0]->err1, by_err1[0]->arg, by_err1[0]->res );
117
118 return (0);
119 }
120