xref: /openbsd-src/regress/lib/libm/fpaccuracy/trailer.h (revision 2b0358df1d88d06ef4139321dd05bd5e05d91eaf)
1*2b0358dfSmartynas /*	$OpenBSD: trailer.h,v 1.1 2009/04/09 01:24:43 martynas Exp $	*/
2*2b0358dfSmartynas 
3*2b0358dfSmartynas /*
4*2b0358dfSmartynas  * Copyright (c) 2009 Gaston H. Gonnet <gonnet@inf.ethz.ch>
5*2b0358dfSmartynas  *
6*2b0358dfSmartynas  * Permission to use, copy, modify, and distribute this software for any
7*2b0358dfSmartynas  * purpose with or without fee is hereby granted, provided that the above
8*2b0358dfSmartynas  * copyright notice and this permission notice appear in all copies.
9*2b0358dfSmartynas  *
10*2b0358dfSmartynas  * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11*2b0358dfSmartynas  * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12*2b0358dfSmartynas  * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13*2b0358dfSmartynas  * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14*2b0358dfSmartynas  * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15*2b0358dfSmartynas  * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16*2b0358dfSmartynas  * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17*2b0358dfSmartynas  */
18*2b0358dfSmartynas 
19*2b0358dfSmartynas static struct point { double arg, res, err1, err2; } points[N+1];
20*2b0358dfSmartynas 
21*2b0358dfSmartynas #include <stdio.h>
22*2b0358dfSmartynas #include <math.h>
23*2b0358dfSmartynas #ifdef FORCE_FPU_DOUBLE
24*2b0358dfSmartynas #include <fpu_control.h>
25*2b0358dfSmartynas #endif
26*2b0358dfSmartynas 
27*2b0358dfSmartynas #define SHOW 20
28*2b0358dfSmartynas 
29*2b0358dfSmartynas int
Fn(FILE * out)30*2b0358dfSmartynas Fn(FILE *out) {
31*2b0358dfSmartynas   struct point *by_err1[N], *by_err2[N], *q;
32*2b0358dfSmartynas   struct input_point *p;
33*2b0358dfSmartynas   double t, t1, t2;
34*2b0358dfSmartynas   int i, tot1={0}, tot2={0}, tot3={0}, tot4={0}, tot5={0};
35*2b0358dfSmartynas 
36*2b0358dfSmartynas #ifdef FORCE_FPU_DOUBLE
37*2b0358dfSmartynas #define FPU_CW  (_FPU_MASK_IM+_FPU_MASK_ZM+_FPU_MASK_OM+_FPU_MASK_UM+ \
38*2b0358dfSmartynas 	_FPU_MASK_PM+_FPU_MASK_DM+_FPU_DOUBLE+_FPU_RC_NEAREST)
39*2b0358dfSmartynas 
40*2b0358dfSmartynas   fpu_control_t cw = FPU_CW;
41*2b0358dfSmartynas   _FPU_SETCW (cw);
42*2b0358dfSmartynas #endif
43*2b0358dfSmartynas 
44*2b0358dfSmartynas   /* compute the results and store them in points structure */
45*2b0358dfSmartynas   for( p=input_points, q=points; q-points<N; p++,q++ ) {
46*2b0358dfSmartynas 	q->arg = scalb(p->arg_m,p->arg_e);
47*2b0358dfSmartynas 	q->res = F( q->arg );
48*2b0358dfSmartynas 	}
49*2b0358dfSmartynas 
50*2b0358dfSmartynas   for( i=0; i<N; i++ ) by_err1[i] = by_err2[i] = points+i;
51*2b0358dfSmartynas 
52*2b0358dfSmartynas   /* do something silly to avoid loop merging optimization */
53*2b0358dfSmartynas   q->res = points[N/2].res+points[N/3].res+points[N/4].res;
54*2b0358dfSmartynas 
55*2b0358dfSmartynas   /* Now compute the error directly, for the benefit of more precise hw */
56*2b0358dfSmartynas   for( p=input_points, q=points; q-points<N; p++,q++ ) {
57*2b0358dfSmartynas 	if( p->val_e >= DBL_MAX_EXP-2 ) {
58*2b0358dfSmartynas 	    t1 = scalb(1.0,p->val_e/2);
59*2b0358dfSmartynas 	    t2 = scalb(1.0,p->val_e-p->val_e/2);
60*2b0358dfSmartynas 	    t = F(q->arg)*t1*t2 - p->val;
61*2b0358dfSmartynas 	    }
62*2b0358dfSmartynas 	else {
63*2b0358dfSmartynas 	    t1 = scalb(1.0,p->val_e);
64*2b0358dfSmartynas 	    t = F(q->arg)*t1 - p->val;
65*2b0358dfSmartynas 	    }
66*2b0358dfSmartynas 	t -= p->eps;
67*2b0358dfSmartynas 	q->err2 = fabs(t);
68*2b0358dfSmartynas 	t = scalb(q->res,p->val_e) - p->val;
69*2b0358dfSmartynas 	t -= p->eps;
70*2b0358dfSmartynas 	q->err1 = fabs(t);
71*2b0358dfSmartynas 	}
72*2b0358dfSmartynas 
73*2b0358dfSmartynas   /* Sort by errors in decreasing order */
74*2b0358dfSmartynas #define KEY(x) (-(x)->err1)
75*2b0358dfSmartynas   SORT(by_err1,0,N-1,q);
76*2b0358dfSmartynas #undef KEY
77*2b0358dfSmartynas #define KEY(x) (-(x)->err2)
78*2b0358dfSmartynas   SORT(by_err2,0,N-1,q);
79*2b0358dfSmartynas 
80*2b0358dfSmartynas   /* count the number of differences in errors */
81*2b0358dfSmartynas   for( q=points; q-points < N; q++ ) {
82*2b0358dfSmartynas 	if( q->err1 > q->err2 ) tot1++;
83*2b0358dfSmartynas 	else if( 2*q->err1 < q->err2 ) tot5++;
84*2b0358dfSmartynas 	else if( 1.1*q->err1 < q->err2 ) tot4++;
85*2b0358dfSmartynas 	else if( 1.01*q->err1 < q->err2 ) tot3++;
86*2b0358dfSmartynas 	else if( q->err1 < q->err2 ) tot2++;
87*2b0358dfSmartynas 	}
88*2b0358dfSmartynas 
89*2b0358dfSmartynas   printf( "result of %s is ", Fs );
90*2b0358dfSmartynas   if( tot1==0 ) printf( "never more precise than double\n\n" );
91*2b0358dfSmartynas   else printf( "more precise than double %d out of %d times\n\n", tot1, N );
92*2b0358dfSmartynas 
93*2b0358dfSmartynas   if( tot2+tot3+tot4+tot5 )
94*2b0358dfSmartynas 	printf( "%d errors <= 1%% worse in accum, %d <= 10%%,"
95*2b0358dfSmartynas 	    "%d <= 100%%, %d > 100%%\n\n", tot2, tot3, tot4, tot5 );
96*2b0358dfSmartynas 
97*2b0358dfSmartynas   for( i=N-1; i>=0 && by_err2[i]->err2 == 0; i-- );
98*2b0358dfSmartynas   if( N-1-i > 0 )
99*2b0358dfSmartynas 	printf( "%d results were exact to double the precision\n\n", N-1-i );
100*2b0358dfSmartynas 
101*2b0358dfSmartynas   if( tot1 > 10 ) {
102*2b0358dfSmartynas 	printf( "%d largest ulp errors (from result in accumulator)\n", SHOW );
103*2b0358dfSmartynas 	for( i=0; i<SHOW; i++ )
104*2b0358dfSmartynas 	    printf( "   %.5f ulp for %s(%.18g) = %.18g)\n",
105*2b0358dfSmartynas 		by_err2[i]->err2, Fs, by_err2[i]->arg, by_err2[i]->res );
106*2b0358dfSmartynas 	printf( "\n" );
107*2b0358dfSmartynas 	}
108*2b0358dfSmartynas 
109*2b0358dfSmartynas   printf( "%d largest ulp errors (stored in a double)\n", SHOW );
110*2b0358dfSmartynas   for( i=0; i<SHOW; i++ )
111*2b0358dfSmartynas 	printf( "   %.5f ulp for %s(%.18g) = %.18g)\n",
112*2b0358dfSmartynas 	    by_err1[i]->err1, Fs, by_err1[i]->arg, by_err1[i]->res );
113*2b0358dfSmartynas   printf( "\n" );
114*2b0358dfSmartynas 
115*2b0358dfSmartynas   fprintf( out, "%8s %5d %27.5f %26.18g %25.18g\n", Fs, N,
116*2b0358dfSmartynas 	by_err1[0]->err1, by_err1[0]->arg, by_err1[0]->res );
117*2b0358dfSmartynas 
118*2b0358dfSmartynas   return (0);
119*2b0358dfSmartynas }
120