sifakis@UUNET.UU.NET (George Sifakis) (09/14/89)
The source that follows consists of a main program and a function which
computes square roots. The main program feeds the function random real
numbers and compares the results with the local C library "sqrt" function.
If compiled without optimization, the results agree to 10**-12. If compiled
with optimization ("-O") the results are grossly different. This indicates
to me a problem with optimization.
Here is the environment:
gcc version 1.35
config.h -> config/xm-m68k.h
tm.h -> config/tm-sun3-nfp.h
Compilation command: gcc -msoft-float [-O] -Blib/ file.c -lm
To run it: a.out 5 5 1 (5 rounds, lots of output)
I would be interested to hear any comments/results. Thanks for a great
compiler.
George Sifakis
ACUSON
1365 Shorebird
Mtn. View, CA 94043
(415)969-9112 x2635
----------------------- CUT HERE -----------------------------------
/*
** Tests delta between versions of a math function
*/
#include <stdio.h>
double delta = .00000000001;
double sqrt(),
sqrt1();
double (*funct[])() = {
sqrt, /* First entry is the comparison base */
sqrt1
};
main(argc, argv)
int argc;
char **argv;
{
int cnt,
ckp,
ckp_int,
max,
f,
report = 0;
long random();
double dval,
res0,
resc;
char s[128];
if (argc < 3) {
printf("Usage is: %s number_of_random_tests ckp_interval [report]\n", argv[0]);
return(-1);
}
sscanf(argv[1], "%d", &max);
sscanf(argv[2], "%d", &ckp_int);
if (argc == 4)
report = 1;
srandom(time((char *)0));
for (cnt = 0, ckp = ckp_int; cnt < max; cnt++, ckp--) {
if (ckp <= 0) {
printf("CHECKPOINT: round %d/%d\n", cnt, max);
ckp = ckp_int;
}
sprintf(s, "%d.%d", random(), random());
sscanf(s, "%lf", &dval);
res0 = (*funct[0])(dval);
if (report)
printf("\nval=%10g\tf0=%10g\n", dval, res0);
for (f=1; f < sizeof(funct) / sizeof(double *); f++) {
resc = (*funct[f])(dval);
if (report)
printf("f%d=%10g\tdiff=%10g\n", f, resc, res0 - resc);
if (res0 - resc > delta || resc - res0 > delta)
printf("ERROR: f0(%f)=%f, f%d(%f)=%f\n",
dval, res0, f, dval, resc);
}
}
}
/* SQRT
* RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT
* FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE
* CODED IN C BY K.C. NG, 3/22/85.
*
* Warning: this code should not get compiled in unless ALL of
* the following machine-dependent routines are supplied.
*
* Required machine dependent functions:
* swapINX(i) ...return the status of INEXACT flag and reset it to "i"
* swapRM(r) ...return the current Rounding Mode and reset it to "r"
* swapENI(e) ...return the status of inexact enable and reset it to "e"
* addc(t) ...perform t=t+1 regarding t as a 64 bit unsigned integer
* subc(t) ...perform t=t-1 regarding t as a 64 bit unsigned integer
*/
#define swapINX(i) 0
#define swapRM(r) 0
#define swapENI(e) 0
#define addc(t) (t=t+1)
#define subc(t) (t=t-1)
static unsigned long table[] = {
0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740,
58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478,
21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, };
double sqrt1(x)
double x;
{
/**double y,z,t,addc(),subc(),b54=134217728.*134217728.;**/ /* b54=2**54 */
double y,z,t,b54=134217728.*134217728.; /* b54=2**54 */
long mx,scalx,mexp=0x7ff00000;
/**int i,j,r,e,swapINX(),swapRM(),swapENI(); **/
int i,j,r,e;
unsigned long *py=(unsigned long *) &y ,
*pt=(unsigned long *) &t ,
*px=(unsigned long *) &x ;
#ifdef national /* ordering of word in a floating point number */
int n0=1, n1=0;
#else
int n0=0, n1=1;
#endif
/* Rounding Mode: RN ...round-to-nearest
* RZ ...round-towards 0
* RP ...round-towards +INF
* RM ...round-towards -INF
*/
int RN=0,RZ=1,RP=2,RM=3;/* machine dependent: work on a Zilog Z8070
* and a National 32081 & 16081
*/
/* exceptions */
if(x!=x||x==0.0) return(x); /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */
if((mx=px[n0]&mexp)==mexp) return(x); /* sqrt(+INF) is +INF */
/* save, reset, initialize */
e=swapENI(0); /* ...save and reset the inexact enable */
i=swapINX(0); /* ...save INEXACT flag */
r=swapRM(RN); /* ...save and reset the Rounding Mode to RN */
scalx=0;
/* subnormal number, scale up x to x*2**54 */
if(mx==0) {x *= b54 ; scalx-=0x01b00000;}
/* scale x to avoid intermediate over/underflow:
* if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */
if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;}
if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;}
/* magic initial approximation to almost 8 sig. bits */
py[n0]=(px[n0]>>1)+0x1ff80000;
py[n0]=py[n0]-table[(py[n0]>>15)&31];
/* Heron's rule once with correction to improve y to almost 18 sig. bits */
t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0;
/* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */
t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y;
t=z/(t+x) ; pt[n0]+=0x00100000; y+=t;
/* twiddle last bit to force y correctly rounded */
swapRM(RZ); /* ...set Rounding Mode to round-toward-zero */
swapINX(0); /* ...clear INEXACT flag */
swapENI(e); /* ...restore inexact enable status */
t=x/y; /* ...chopped quotient, possibly inexact */
j=swapINX(i); /* ...read and restore inexact flag */
if(j==0) { if(t==y) goto end; else t=subc(t); } /* ...t=t-ulp */
b54+0.1; /* ..trigger inexact flag, sqrt(x) is inexact */
if(r==RN) t=addc(t); /* ...t=t+ulp */
else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */
y=y+t; /* ...chopped sum */
py[n0]=py[n0]-0x00100000; /* ...correctly rounded sqrt(x) */
end: py[n0]=py[n0]+scalx; /* ...scale back y */
swapRM(r); /* ...restore Rounding Mode */
return(y);
}