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 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