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