be careful about optimizing a+b-a-b !
David Hough
uunet!Eng.Sun.COM!David.Hough
Fri May 8 14:59:34 PDT 1992
Many of you will recognize loops like this one from Cody's cmachar:
a = two;
b = a;
zero = 0.0e0;
tmp = ((a+one)-a)-one;
while (tmp == zero) {
a = a+a;
tmp = a+one;
tmp1 = tmp-a;
tmp = tmp1-one;
}
as typical in numerical programs that dynamically determine precision.
While there are other methods (float.h) that may be helpful toward that
end in ANSI-C, and environmental inquiries in Fortran-90,
there are a lot of programs written in traditional C,
or translated from Fortran-77 via F2C, that do it dynamically. And it may
not be easy to recognize them if you're not used to thinking this way.
Compiler writers would do us a favor by NOT doing us a favor by optimizing
away a+1-a-1 in the foregoing to zero. The whole point of the computation
is to detect the roundoff error. Unfortunately in the rush to better
SPECmarks, some people may succumb to temptation. Just don't make it
the default, please. For instance, on SPARC running SunOS 4.1.2,
Lucid C 2.1 and its predecessor 2.0
both goof on the foregoing loop in the default compilation mode (also known
as -O0). Correct compilation can be obtained by compiling -g or -O-,
or by declaring "tmp" volatile.
It's not just roundoff inquiries, either; the following part of a GCD
computation in floating point also needed volatile declarations for
some of the variables:
do {
g_0 = g; g = g_1; j_0 = j; j = j_1; k_0 = k; k = k_1;
m = rint(g_0/g);
g_1 = g_0 - m*g; j_1 = j_0 - m*j; k_1 = k_0 - m*k;
} while (g_1 != 0);
as did the following random number computation from another
program (for Lo):
real ranf()
{
short *p,*q,k; long Hi;
#ifdef __lucid
volatile long Lo;
#else
long Lo;
#endif
/* generate product using double precision simulation (comments */
/* refer to In's lower 16 bits as "L", its upper 16 bits as "H") */
p=(short *)&In[strm]; Hi= *(p)*A; /* 16807*H->Hi */
*(p)=0; Lo=In[strm]*A; /* 16807*L->Lo */
p=(short *)&Lo; Hi+= *(p); /* add high-order bits of Lo to Hi */
q=(short *)&Hi; /* low-order bits of Hi->LO */
*(p)= *(q+1)&0X7FFF; /* clear sign bit */
k= *(q)<<1; if (*(q+1)&0X8000) then k++; /* Hi bits 31-45->K */
/* form Z + K [- M] (where Z=Lo): presubtract M to avoid overflow */
Lo-=M; Lo+=k; if (Lo<0) then Lo+=M;
In[strm]=Lo;
return((real)Lo*4.656612875E-10); /* Lo x 1/(2**31-1) */
}
Speaking of F2C, another problem with Lucid C 2.1
requires that certain static local variables be moved outside the function
in the F2C source file putpcc.c:
#ifdef __lucid
static int initargs[4] = {TYCOMPLEX, TYDCOMPLEX, TYCHAR, TYFTNLEN+100};
static int *init_ap[TYSUBR+1] = {0,0,0,0,0,0,
initargs, initargs+1,0,initargs+2};
#endif
void
save_argtypes(arglist, at0, at1, ccall, fname, stg, nchargs, type, zap)
chainp arglist;
Argtypes **at0, **at1;
int ccall, stg, nchargs, type, zap;
char *fname;
{
Argtypes *at;
chainp cp;
int i, i0, j, k, nargs, *t, *te;
Atype *atypes;
expptr q;
char buf[208];
#ifndef __lucid
static int initargs[4] = {TYCOMPLEX, TYDCOMPLEX, TYCHAR, TYFTNLEN+100};
static int *init_ap[TYSUBR+1] = {0,0,0,0,0,0,
initargs, initargs+1,0,initargs+2};
#endif
No doubt Lucid customers could obtain fixes from Lucid once they had isolated
the problems. The version of Lucid C 2.1 was the one released last month.
I'm only pointing out these trouble areas as a warning to the many compiler
writers in NCEG to be careful how they optimize floating-point codes.
For everybody's benefit, I'll include the source for a test program by
Bonnie Toy for testing floating-point multiplication, based on ideas
of W. Kahan. It's in a form for testing on SunOS 4.x but might be portable
to other systems with some effort.
Compile with -DSP for single precision and -DDP for double precision.
It turned up a couple of the problems mentioned above.
The header comments are somewhat out of date.
If it doesn't work for you, please don't complain to us about it!
-------------------------------------------------------------------------------
/* CMULT, a C program to test whether multiplication
* is correctly rounded according to the IEEE standard 754.
*
* Based on W. Kahan's "To Test Whether Binary Floating-Point Multiplication is
* Correctly Rounded", July 13, 1988.
*
* Implemented by Bonnie Toy.
* Modified by Shing Ma, K.C. Ng, M. Mueller.
* Last change: 1/18/89.
*
* Required system supported functions:
* ieee_flags, nextafter, floor, scalbn, rint, ceil,
* copysign, infinity.
*
* Compile switches:
* -DN_RAND=n Test n random arguments
* -DSINGLE Test single precision multiplication
* -DDOUBLE Test double precision multiplication
* -Ddebug Print X, Y, computed X*Y, expected X*Y on failure
*
* Assumptions:
* Addition and subtraction are correctly rounded.
* Multiplication and division are in error by less than one ulp.
*
* Idea:
* Generate arguments X and Y whose product lies halfway between
* two representable numbers, or nearly halfway. Check if the
* rounded product [X*Y] is CORRECTLY rounded, and if rounding
* is CONSISTENT.
*
* A floating point arithmetic operation is correctly rounded
* when its rounding error does not exceed half an ulp (unit
* in the last place.) When the error is exactly half an ulp,
* the IEEE standard requires rounding to nearest EVEN (the
* nearest representable number whose last bit is zero.)
*
* Functions/procedures contained in this file:
* rand Random number generator
* cmultrand Seeds random number generator
* getit Determines number of significant bits being used
* round2 Series of clever *, / that returns 2 if *, / are
* correctly rounded, 1 otherwise
* halfway Drives halfway tests
* test_half Tests halfway cases
* nearhalf Drives hear-halfway tests
* test_nearhalf Tests near-halfway cases
* main Drive halfway and nearhalf
* gcd Find j and k so that 1= gcd(i,4t) = i*j - 4*t*k
* prhex Print hex representation of argument (for debug)
*/
/* Global definitions */
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
/*
* Need to define ASSIGNFLOAT?
*/
#ifndef ASSIGNFLOAT
#if defined(FLOATFUNCTIONTYPE)
#define ASSIGNFLOAT(x,y) { union {FLOATFUNCTIONTYPE _d; float _f; } _kluge; _kluge._d = (y); x = _kluge._f; }
#else
#define ASSIGNFLOAT(x,y) x=(y)
#endif
#endif
#ifdef SP
#define GENERIC float
#else
#ifdef DP
#define GENERIC double
#else
#define GENERIC double /* default is IEEE double precision */
#endif
#endif
#ifdef f68881
#define EXTENDED
#endif
#ifdef f80387
#define EXTENDED
#endif
#ifndef N_RAND
#define N_RAND 10000
#endif
#ifndef ROUND2_A
#define ROUND2_A (GENERIC)1717.0
#endif
#define seed 0 /* seed for random number generator */
#define ZERO (GENERIC)0
#define ONE (GENERIC)1
#define TWO (GENERIC)2
#define FOUR (GENERIC)4
struct detailed_info{
GENERIC x, y; /*numbers x & y whose product is incorrectly rounded*/
GENERIC expected; /* expected result of x*y */
GENERIC computed; /* computed result of x*y */
};
struct info {
int nf; /* number of failures */
struct detailed_info first, last, worst;
};
static int firstFailure_consistent, firstFailure_correct;
GENERIC worstSoFar_consistent, worstSoFar_correct;
int count; /* count how many times each test is performed */
static int bits, /* number of bits in GENERIC type */
getit(); /* get number of significant bits in a GENERIC type */
static GENERIC two_to_bits, twobp1, twobm2;
void prhex(); /* print hex value of argument */
void print_results(); /* print x,y, computed x*y, expected x*y, in decimal*/
/*
* Random number algorithm is based on the additive number generator described
* in Knuth's The Art of Computing, Vol II (second edition) section 3.2.2.
* It returns a positive random integer.
*/
static unsigned long Y[] = {
0x8ca0df45, 0x37334f23, 0x4a5901d2, 0xaeede075, 0xd84bd3cf,
0xa1ce3350, 0x35074a8f, 0xfd4e6da0, 0xe2c22e6f, 0x045de97e,
0x0e6d45b9, 0x201624a2, 0x01e10dca, 0x2810aef2, 0xea0be721,
0x3a3781e4, 0xa3602009, 0xd2ffcf69, 0xff7102e9, 0x36fab972,
0x5c3650ff, 0x8cd44c9c, 0x25a4a676, 0xbd6385ce, 0xcd55c306,
0xec8a31f5, 0xa87b24ce, 0x1e025786, 0x53d713c9, 0xb29d308f,
0x0dc6cf3f, 0xf11139c9, 0x3afb3780, 0x0ed6b24c, 0xef04c8fe,
0xab53d825, 0x3ca69893, 0x35460fb1, 0x058ead73, 0x0b567c59,
0xfdddca3f, 0x6317e77d, 0xaa5febe5, 0x655f73e2, 0xd42455bb,
0xe845a8bb, 0x351e4a67, 0xa36a9dfb, 0x3e0ac91d, 0xbaa0de01,
0xec60dc66, 0xdb29309e, 0xcfa52971, 0x1f3eddaf, 0xe14aae61,
};
static long j=23, k=54;
rand()
{
unsigned long m;
m = Y[j]+Y[k];
Y[k] = m;
j=j-1;
k=k-1;
if(j<0) j=54;
if(k<0) k=54;
return m&0x7fffffff;
}
/*
* cmultrand() uses ax+c mod 2**32 to generate seeds for rand(). Here
* a=8*(10**8-29)+5, c=10**9-63.
*/
static unsigned long Z[] = {
0x8ca0df45, 0x37334f23, 0x4a5901d2, 0xaeede075, 0xd84bd3cf,
0xa1ce3350, 0x35074a8f, 0xfd4e6da0, 0xe2c22e6f, 0x045de97e,
0x0e6d45b9, 0x201624a2, 0x01e10dca, 0x2810aef2, 0xea0be721,
0x3a3781e4, 0xa3602009, 0xd2ffcf69, 0xff7102e9, 0x36fab972,
0x5c3650ff, 0x8cd44c9c, 0x25a4a676, 0xbd6385ce, 0xcd55c306,
0xec8a31f5, 0xa87b24ce, 0x1e025786, 0x53d713c9, 0xb29d308f,
0x0dc6cf3f, 0xf11139c9, 0x3afb3780, 0x0ed6b24c, 0xef04c8fe,
0xab53d825, 0x3ca69893, 0x35460fb1, 0x058ead73, 0x0b567c59,
0xfdddca3f, 0x6317e77d, 0xaa5febe5, 0x655f73e2, 0xd42455bb,
0xe845a8bb, 0x351e4a67, 0xa36a9dfb, 0x3e0ac91d, 0xbaa0de01,
0xec60dc66, 0xdb29309e, 0xcfa52971, 0x1f3eddaf, 0xe14aae61,
};
static unsigned long a= 0x2faf071d, /* a = 8*(10**8-29)+5 */
c= 0x3b9ac9c1; /* c = 10**9-63 */
cmultrand(iseed)
unsigned long iseed;
{
long i;
j = 23; k=54;
if(iseed==0) for (i=0;i<55;i++) Y[i] = Z[i];
else {
Y[0] = (a*iseed+c)>>1;
for (i=1;i<55;i++) Y[i] = (a*Y[i-1]+c)>>1;
}
}
static int
getit() /* get number of significant bits */
{
int it;
#ifdef __lucid
volatile GENERIC b,t;
#else
GENERIC b,t;
#endif
it = 0;
b = ONE;
do {
it++;
b += b;
t = b+ONE;
t -= b;
t -= ONE;
} while (t == ZERO);
return it;
}
main()
{
void halfway(), nearhalf(), round2();
printf("\n CMULT tests whether multiplication is correctly rounded");
printf("\n according to ANSI/IEEE Std 754-1985, 'The IEEE Standard");
printf("\n for Binary Floating-Point Arithmetic.' ");
printf("\n\n There are five tests." );
printf("\n The first is a one-shot test of multiplication and division ");
printf("\n that results in the answer 2.0 if IEEE rounding is in effect, ");
printf("\n and the answer 1.0 otherwise. ");
printf("\n\n The next two tests check whether the product x*y is rounded ");
printf("\n consistently and correctly, when x and y are chosen so that ");
printf("\n their product is a number that falls halfway between two ");
printf("\n representable numbers. IEEE specifies rounding towards even.");
printf("\n\n The last two tests check whether the product x*y is rounded ");
printf("\n consistently and correctly, when x and y are chosen so that ");
printf("\n their product is very nearly a halfway case. IEEE specifies ");
printf("\n rounding towards nearest.");
printf("\n\n By default, double precision multiplication is tested ");
printf("\n 10,000 times. ");
printf("\n\n To change the defaults, compile CMULT with the switches:");
printf("\n -DROUND2_A=n Use n as the input to the first test.");
printf("\n Choose n so that 1000 < n < 8,000,000.");
printf("\n -DN_RAND=n Test n random arguments.");
printf("\n -DSP IEEE 754 Single Precision or VAX F; 24 bits ");
printf("\n -DDP IEEE 754 Double Precision or VAX G; 53 bits ");
printf("\n -------------------------------------------------------------\n");
/* set up some useful global variables */
/*get number of bits in GENERIC type*/
bits = getit();
#ifdef EXTENDED
#ifdef SP
bits = 24;
#else
bits = 53;
#endif
#endif
two_to_bits = scalbn(1.0, bits);
twobp1 = scalbn(1.0, bits + 1 );
twobm2 = scalbn(1.0, bits - 2 );
cmultrand(seed); /*seed random number generator*/
printf("\n\n1. One-shot test of multiplication & division: ");
round2();
halfway();
nearhalf();
ieee_retrospective_(); /* Bugs cause IEEE exceptions. */
exit(0); /* Must explicitly return zero status. */
}
void round2()
{
GENERIC A, One, Two, Half, Three, TwoThirds, Roundoff, C, S, I, D, Q, X, E, Z, first, last;
long i, j;
char *out;
/* This is a cute little test that will pick up discrepancies in
rounding. If IEEE rounding is used, the series of multiplications
and divisions in the main "for" loop, below, will always result
in Z = 2. If IEEE rounding is not used, those same operations
result in Z = 1.
/* A, the input to the series of multiplications and division, may take
any value between 1,000 and 8,000,000. Test for different input
values by defining -DROUND2_A=n during compilation. */
/* Translated into C from the program written by W. Kahan in his
paper "A Computer Program with Almost No Significance",
November 16, 1988. */
ieee_flags("set", "direction", "nearest", &out);
first = ZERO; last = ZERO;
One = ONE; /* (GENERIC)1 */
Two = One + One;
Half = One/Two;
Two = One + One;
Three = One + Two;
TwoThirds = Two/Three;
/* Roundoff = ( 3 * roundoff in TwoThirds) */
Roundoff = (((TwoThirds - Half) - Half) + (TwoThirds - Half)) +
(TwoThirds - Half);
/* C is a huge number, normally like (1/Roundoff)^2 */
C = (Roundoff == ZERO) ? 1.0e36 : One/(Roundoff*Roundoff) ;
S = One;
I = One; /* I takes on values 1,3,5,7,... */
while (I < ROUND2_A) {
D = Three; /* D takes on values D = 1 + 2^j */
for (j=1; j <= 15; j++) {
Q = I / D; /* Q = I/D rounded */
X = Q*D; /* X = I + roundoff */
E = (X-I)*C; /* E = roundoff*C */
S = E*E + S; /* S = 1 + summation of E^2 */
D = D - One + D;
}
I += Two; /* Now S=1 + summation of (roundoff * C)^2 */
}
Z = One + One/S; /* If all roundoffs = 0, then Z = 2 */
if (Z == TWO)
printf("passed\n");
else {
printf("failed for A = %f.\n", ROUND2_A);
#ifdef debug
printf("\n failed for A = %f", ROUND2_A);
printf("A has hex representation:");
prhex(ROUND2_A);
#endif
} /* end of if(Z==TWO)/else */
} /* end of round2 routine */
void halfway()
/* This method generates products X*Y that are all half-odd-integers
* in the binade 2^(bits-1) < X*Y < 2^bits, for which [X*Y] should
* round to the nearest even integer for IEEE 754, to the next larger
* integer for the VAX.
*/
{
GENERIC jl, ju, j, x, temp;
char *out, *in;
long i, m;
struct info consistent, correct;
void test_half();
/* initializations */
consistent.nf = 0;
correct.nf = 0;
firstFailure_consistent=1;
firstFailure_correct = 1;
worstSoFar_consistent=ZERO;
worstSoFar_correct=ZERO;
count = 0;
#ifdef EXTENDED
ieee_flags("set", "precision", "double", &out);
#endif
/* Generate a random odd integer in the interval (2, 2^bits) */
jl = ju = ZERO;
#ifdef SP
x = (float)3;
#endif
#ifdef DP
x = ZERO;
while (( (x < TWO) || (x > twobm2)) || (fmod(x,2.0)==0.0))
{ x = rand();
}
#endif
/* Compute in floating-point two integers
* Jl := ceil((2^pi - (x-1))/(2x))
* Ju := floor((2^(pi+1) - (x+1))/(2x)).
* Each quotient can be computed with only one rounding error which,
* ideally, should be directed upward for Jl, downward for Ju.
*/
ieee_flags("set", "direction", "positive", &out);
temp = x - ONE;
jl = ceil((two_to_bits - temp)/(TWO*x));
ieee_flags("set", "direction", "tozero", &out);
temp = x + ONE;
ju = floor((twobp1 - temp)/(TWO*x));
ieee_flags("set", "direction", "nearest", &out);
/* Then choose at random any integers J between Jl and Ju, as well
* as Jl and Ju, from which to construct test arguments Y := J + 1/2
* representable exactly in floating-point.
*/
consistent.nf = 0; correct.nf = 0;
test_half(jl, x, &consistent, &correct);
test_half(ju, x, &consistent, &correct);
if ((ju -jl) < N_RAND)
for (j=jl+1; j<ju; j++)
test_half(j,x, &consistent, &correct);
else
for(i = 2; i < N_RAND; i++){
j = (GENERIC) rand();
while ((j <= jl) || (j >= ju)) j = (GENERIC) rand();
test_half(j, x, &consistent, &correct);
}
/* print results of testing halfway cases */
printf("\n2. Halfway cases, %d cases tested:\n\n",
count);
if (consistent.nf == 0)
printf(" a. Passed consistency tests for halfway cases.\n\n");
else{
printf(" a. Failed consistency tests for halfway cases %d times. \n\n",
consistent.nf);
print_results(&consistent);
}
if (correct.nf == 0)
printf(" b. Passed correctness tests for halfway cases.\n\n");
else{
printf(" b. Failed correctness tests for halfway cases %d times. \n\n",
correct.nf);
print_results(&correct);
}
}
void test_half(j,x,consistent,correct)
GENERIC j, x;
struct info *consistent, *correct;
{
GENERIC y, u, uu, xy, xj, badNews;
y = j + 0.5;
/* Now every product X*Y turns out to be half an odd integer in a
* binade where it must round to the nearest even integer for
* IEEE 754, the next larger integer for a VAX. The rounded
* product [X*Y] should match the rounded sum [X*J + X/2] of the
* two terms each of which is computable exactly; otherwise
* multiplication and addition do not round consistently. To
* test both operations for correct rounding, calculate
* U := (X-1)*J + (X-1)/2 + J
* exactly and expect to find [X*Y] = U if U is even and IEEE
* 754 is in force, otherwise [X*Y] = U+1, or else conclude that
* rounding has failed the test.
*/
count++; /* record how many times this test is REALLY performed */
xj = x*j + x/2;
u = (x-1.0)*j+(x-1.0)/2.0+j;
uu = (fmod(u,2.0) == 0.0)?u:u+1.0;
xy = x*y;
badNews = fabs(xy-xj);
if (xy != xj){
consistent->nf++;
consistent->last.x = x;
consistent->last.y = y;
consistent->last.expected = xj;
consistent->last.computed = xy;
if (firstFailure_consistent){
consistent->first.x = x;
consistent->first.y = y;
consistent->first.expected = xj;
consistent->first.computed = xy;
firstFailure_consistent = 0;
worstSoFar_consistent = badNews;
}
if (badNews >= worstSoFar_consistent) {
consistent->worst.x = x;
consistent->worst.y = y;
consistent->worst.expected = xj;
consistent->worst.computed = xy;
worstSoFar_consistent = badNews;
}
#ifdef debug
printf("Consistent Rounding Failed\n");
printf("X = "); prhex(x); printf(", Y = "); prhex(y);
printf("\n");
printf("U = "); prhex(u); printf(", J = "); prhex(j);
printf("\n");
printf("[X*Y] = "); prhex(xy);
printf(", [X*J+X/2] = "); prhex(xj);
printf("\n");
#endif
}
badNews = fabs(xy-uu);
if (xy != uu){
correct->nf++;
correct->last.x = x;
correct->last.y = y;
correct->last.expected = uu;
correct->last.computed = xy;
if (firstFailure_correct){
correct->first.x = x;
correct->first.y = y;
correct->first.expected = uu;
correct->first.computed = xy;
firstFailure_correct = 0;
worstSoFar_correct = badNews;
}
if (badNews >= worstSoFar_correct) {
correct->worst.x = x;
correct->worst.y = y;
correct->worst.expected = xj;
correct->worst.computed = xy;
worstSoFar_correct = badNews;
}
#ifdef debug
printf("Correct Rounding Failed\n");
printf("X = "); prhex(x); printf(", Y = "); prhex(y);
printf("\n");
printf("U = "); prhex(u); printf(", J = "); prhex(j);
printf("\n");
printf("[X*Y] = "); prhex(xy);
printf(", [X*J+X/2] = "); prhex(xj);
printf("\n");
#endif
}
}
void nearhalf()
/* Generate odd integers X and Y at random, in the binade between
* 2^(pi-1) and 2^pi, of which many will satisfy either
* 2^(2*pi-1) < X*Y < 2^(2*pi) and X*Y = [X*Y] (+|-) (2^(pi-1) - 1), or
* 2^(2*pi-2) < X*Y < 2^(2*pi-1) and X*Y = [X*Y] (+|-) (2^(pi-2) - 1).
* These products come as close as possible to half-way cases without
* hitting one. IEEE 754 and VAX round them in the same way.
*/
{
GENERIC t, i;
int n;
struct info consistent, correct;
void test_nearhalf();
/* initializations */
consistent.nf = 0; correct.nf = 0;
firstFailure_consistent = 1; firstFailure_correct = 1;
worstSoFar_consistent = ZERO; worstSoFar_correct = ZERO;
count = 0;
/* Abbreviate t := 2^(pi-3), so that all integers between (+|-)8t
are representable exactly in floating-point. The first step is
to choose at random an odd integer i in the interval 0 < i < t;
i := 1, i := 3, i := t-1 and i := t-3 are good choices too.
*/
t = scalbn(1.0, bits - 3);
test_nearhalf(1.0, t, &correct, &consistent);
test_nearhalf(3.0, t, &correct, &consistent);
test_nearhalf(t-1.0, t, &correct, &consistent);
test_nearhalf(t-3.0, t, &correct, &consistent);
for (n = 4; n < N_RAND;n++) {
i = fmod((GENERIC) rand(), t); /* i < t */
i = (fmod(i, 2.0) == ONE)? i : i + ONE; /* odd(i) */
test_nearhalf(i, t, &correct, &consistent);
}
/* print results from testing nearly halfway cases */
printf("\n3. Nearly halfway cases, %d cases tested::\n\n",
count);
if (consistent.nf==0)
printf(" a. Passed consistency tests for nearly halfway cases.\n\n");
else{
printf(" a. Failed consistency tests for nearly halfway cases %d times. \n\n", consistent.nf);
print_results(&consistent);
}
if (correct.nf==0)
printf(" b. Passed correctness tests for nearly halfway cases. \n\n");
else{
printf(" b. Failed correctness tests for nearly halfway cases %d times. \n\n", correct.nf);
print_results(&correct);
}
} /* end of nearhalf */
void test_nearhalf(i, t, correct, consistent)
GENERIC i,t;
struct info *consistent, *correct;
{
GENERIC j, k, ii[4],jj[4],kk[4][4],X[4],Y[4];
GENERIC ell[4][4], lambda_1[4][4], lambda_2[4][4];
GENERIC S[4][4], D[4][4], H[4][4], T[4][4];
GENERIC sum, prod, sign, ulp_of_Hlm, two_Tlm;
GENERIC k_temp, ij_temp, temp, tempo, tempSign, badNews;
int d, L, M;
char *out;
void gcd(), prhex();
count++; /* record how many times this test is REALLY performed */
/* Compute the greatest common divisor of i and 4t using the
* Euclidean algorithm to get j and k too:
* 1 = gcd(i, 4t) = i*j - 4t*k with 0 < j < 4t and 0 <= k < i.
*/
gcd(i,t,&j,&k);
/* From this trio (i, j, k) derive a collection of quantities all
* computed exactly in floating-point thus:
*
* First let delta := sign(2t - j) = (+|-)1 according as 2t (>|<) j.
*/
d = (2*t > j)?1:-1;
/* Then compute in turn
*/
ii[0] = i;
X[0] = 4*t+ii[0]; jj[0] = j; Y[0] = 4*t+jj[0];
ii[1] = 2*t+i;
X[1] = 4*t+ii[1]; jj[1] = 2*d*t+j; Y[1] = 4*t+jj[1];
ii[2] = 4*t-ii[0];
X[2] = 4*t+ii[2]; jj[2] = 4*t-jj[0]; Y[2] = 4*t+jj[2];
ii[3] = 4*t-ii[1];
X[3] = 4*t+ii[3]; jj[3] = 4*t-jj[1]; Y[3] = 4*t+jj[3];
kk[0][0] = k;
kk[0][1] = k+i*d/2.0;
kk[1][0] = k+j/2.0;
k_temp = kk[0][1] + kk[1][0];
kk[1][1] = k_temp - k + d*t;
for (L=0;L<2;L++){
for (M=0;M<2;M++){
kk[L][M+2] = ii[L] - kk[L][M];
kk[L+2][M] = jj[M] - kk[L][M];
ij_temp = ii[L] + jj[M];
kk[L+2][M+2] = 4*t - ij_temp + kk[L][M];
}
}
for (L=0;L<4;L++){
for (M=0;M<4;M++){
ell[L][M] = copysign(1.0,(L-1.5)*(M-1.5));
}
}
/* At this point each of the intervals (0, t), (t, 2t), (2t, 3t) and
* (3t, 4t) contains just one i[l] and just one j[m], whereupon
* each interval (4t, 5t), (5t, 6t), (6t, 7t) and (7t, 8t) must
* contain just one X[L] and just one Y[M], all of them odd integers
* that can be represented exactly in floating-point. Moreover, all
* 32 products i[L]*j[M] and X[L]*Y[M] differ from multiples of 2t
* by l[L][M] = (+|-)1. In fact,
* i[L]*j[M] = 4t*k[L][M] + l[L][M] and
* x[L]*y[M] = 4t*(4t + (i[L]+j[M]) + k[L][M]) + l[L][M];
* but note that k[L][M] is a half-integer if L-M is odd, an integer
* if L-M is even.
*
* Therefore the 16 products X[L]*Y[M] are all 2*pi or 2*(pi-1) bits
* wide with pi-2 trailing bits ...00001 or ...11111; about half of
* the products will have the property needed to test multiplication
* for correct rounding, namely that they come as close as possible
* to half-way cases without hitting one. All that remains is to
* show how to compute [X[L]*Y[M]] in another way that tests whether
* multiplication is rounded correctly, consistently with addition.
*
* We need formulas that express each product as a sum of two terms,
* X[L]*Y[M] = H[L][M] + T[L][M],
* each of which is computable exactly. Here they are:
*/
#ifdef EXTENDED
ieee_flags("set", "precision", "double", &out);
#endif
for (L=0;L<4;L++){
for (M=0;M<4;M++){
temp = floor(kk[L][M]/2.0);
lambda_1[L][M] = kk[L][M] - 2.0*temp;
lambda_2[L][M] = kk[L][M] - lambda_1[L][M];
sum = ii[L] + jj[M];
sum = 4*t + sum + lambda_2[L][M];
S[L][M] = 4*t*sum;
D[L][M] = 4*t*lambda_1[L][M] + ell[L][M];
H[L][M] = S[L][M] + D[L][M];
/* sum = S[L][M] - H[L][M];*/
T[L][M] = (S[L][M] - H[L][M]) + D[L][M];
#ifdef DEBUG
printf( "T[%1d][%1d] = ", L, M );
prhex(T[L][M]);
printf("\n");
#endif
}
}
for (L=0;L<4;L++){
for (M=0;M<4;M++){
prod = X[L]*Y[M];
badNews = fabs(prod-H[L][M]);
if (prod != H[L][M]) {
consistent->nf++;
consistent->last.x = X[L];
consistent->last.y = Y[M];
consistent->last.computed = prod;
consistent->last.expected = H[L][M] + T[L][M];
if (firstFailure_consistent){
consistent->first.x = X[L];
consistent->first.y = Y[M];
consistent->first.computed = prod;
consistent->first.expected = H[L][M] + T[L][M];
firstFailure_consistent = 0;
worstSoFar_consistent=badNews;
}
if (badNews >= worstSoFar_consistent){
consistent->worst.x = X[L];
consistent->worst.y = Y[M];
consistent->worst.computed = prod;
consistent->worst.expected = H[L][M] + T[L][M];
worstSoFar_consistent = badNews;
}
#ifdef debug
printf("Inconsistent Rounding\n");
printf("L=%d, M=%d\n", L, M);
printf("H[L][M] = "); prhex(H[L][M]);
printf(", X[L]*Y[M] = "); prhex(prod);
printf("\n");
printf("X[L] = "); prhex(X[L]);
printf(", Y[M] = "); prhex(Y[M]);
printf(", T[L][M] = "); prhex(T[L][M]);
printf("\n");
#endif
}
sign = copysign(1.0,T[L][M]);
#ifdef SP
tempo = H[L][M];
tempSign = sign*infinity();
ASSIGNFLOAT(ulp_of_Hlm, r_nextafter_(&tempo,&tempSign));
#endif
#ifdef DP
ulp_of_Hlm = nextafter(H[L][M], sign*infinity());
#endif
ulp_of_Hlm = fabs(ulp_of_Hlm - H[L][M]);
two_Tlm = 2.0 * fabs(T[L][M]);
if (two_Tlm >= ulp_of_Hlm) {
correct->nf++;
correct->last.x = X[L];
correct->last.y = Y[M];
correct->last.computed = prod;
correct->last.expected = H[L][M] + T[L][M];
if (firstFailure_correct){
correct->first.x = X[L];
correct->first.y = Y[M];
correct->first.computed = prod;
correct->first.expected = H[L][M] + T[L][M];
firstFailure_correct = 0;
worstSoFar_correct=badNews;
}
if (badNews >= worstSoFar_correct){
correct->worst.x = X[L];
correct->worst.y = Y[M];
correct->worst.computed = prod;
correct->worst.expected = H[L][M] + T[L][M];
worstSoFar_correct = badNews;
}
#ifdef debug
printf("Incorrect Rounding\n");
printf("2|T[L][M]| = "); prhex(two_Tlm);
printf(", 1 ulp of H[L][M] = "); prhex(ulp_of_Hlm);
printf("\n");
#endif
} /* end of if(two_Tlm <= ulp_of_Hlm) */
}/* end of for(M=0;M<4;M++) */
} /*end of for(L=0,L<4,L++) */
} /* end of test_nearhalf */
void gcd(i, t, j_ret, k_ret)
GENERIC i, t, *j_ret, *k_ret;
/* Based upon the Euclidean algorithm for a Greatest Common Divisor,
* this program starts from a given t = 2^(bits-3) and a given randomly
* chosen positive integer i < t, and yields a sequence of triples
* <g_n, j_n, k_n> that all satisfy g_n = i*j_n - 4t*k_n while each
* new |g_n| <= |g_n|/2 until at last |g_n| = GCD(i, 4t) = 1.
*/
{
#ifdef __lucid
volatile GENERIC j,k;
volatile GENERIC g, g_0, g_1, j_0, j_1, k_0, k_1, m;
#else
GENERIC j,k;
GENERIC g, g_0, g_1, j_0, j_1, k_0, k_1, m;
#endif
/* initialize */
g = FOUR*t; g_1 = i; j = ZERO; k_1 = ZERO; j_1 = ONE; k = -ONE;
/* find factors */
do {
g_0 = g; g = g_1; j_0 = j; j = j_1; k_0 = k; k = k_1;
m = rint(g_0/g);
g_1 = g_0 - m*g; j_1 = j_0 - m*j; k_1 = k_0 - m*k;
} while (g_1 != 0);
/* adjust factors if nec. */
if (g < ZERO) {g = -g; j = -j; k = -k;}
if (j < ZERO) { j = FOUR*t + j; k = k+i;}
#ifdef debug
/* check factors */
if ((j >= FOUR*t) || (j <= ZERO)){
printf("bad value for j = %f\n", j);
}
if ((k >= i) || (k < ZERO)) {
printf(" bad value for k = %f\n", k);
}
if ((i*j - FOUR*t*k) != ONE) {
printf(" bad factors i, j, k: \n");
printf(" %f %f %f\n", i, j, k);
badFactors++;
if ((i*j - FOUR*t*k) == ZERO) zeroed++;
}
#endif
*j_ret = j; *k_ret = k;
}
void prhex(x)
double x;
{
#ifdef SP
union print_hex{
unsigned i;
float flt;
}px;
px.flt = x;
printf("0x%08X", px.i);
#endif
#ifdef DP
union print_hex{
unsigned i[2];
double dbl;
} px;
px.dbl = x;
printf("0x%08x 0x%08x", px.i[0], px.i[1]);
#endif
}
void print_results(ans)
struct info *ans;
{
#ifdef SP
printf(" First: ");
printf(" x = "); prhex(ans->first.x); printf(" = %8.7e\n", ans->first.x);
printf(" y = "); prhex(ans->first.y); printf(" = %8.7e\n", ans->first.y);
printf(" Computed x*y = "); prhex(ans->first.computed); printf(" = %8.7e\n", ans->first.computed);
printf(" Expected x*y = "); prhex(ans->first.expected); printf(" = %8.7e\n\n", ans->first.expected);
printf(" Worst: ");
printf(" x = "); prhex(ans->worst.x); printf(" = %8.7e\n", ans->worst.x);
printf(" y = "); prhex(ans->worst.y); printf(" = %8.7e\n", ans->worst.y);
printf(" Computed x*y = "); prhex(ans->worst.computed); printf(" = %8.7e\n", ans->worst.computed);
printf(" Expected x*y = "); prhex(ans->worst.expected); printf(" = %8.7e\n\n", ans->worst.expected);
printf(" Last : ");
printf(" x = "); prhex(ans->last.x); printf(" = %8.7e\n", ans->last.x);
printf(" y = "); prhex(ans->last.y); printf(" = %8.7e\n", ans->last.y);
printf(" Computed x*y = "); prhex(ans->last.computed); printf(" = %8.7e\n", ans->last.computed);
printf(" Expected x*y = "); prhex(ans->last.expected); printf(" = %8.7e\n\n", ans->last.expected);
#endif
#ifdef DP
printf(" First: ");
printf(" x = "); prhex(ans->first.x); printf(" = %16.15e\n", ans->first.x);
printf(" y = "); prhex(ans->first.y); printf(" = %16.15e\n", ans->first.y);
printf(" Computed x*y = "); prhex(ans->first.computed); printf(" = %16.15e\n", ans->first.computed);
printf(" Expected x*y = "); prhex(ans->first.expected); printf(" = %16.15e\n\n", ans->first.expected);
printf(" Worst: ");
printf(" x = "); prhex(ans->worst.x); printf(" = %16.15e\n", ans->worst.x);
printf(" y = "); prhex(ans->worst.y); printf(" = %16.15e\n", ans->worst.y);
printf(" Computed x*y = "); prhex(ans->worst.computed); printf(" = %16.15e\n", ans->worst.computed);
printf(" Expected x*y = "); prhex(ans->worst.expected); printf(" = %16.15e\n\n", ans->worst.expected);
printf(" Last : ");
printf(" x = "); prhex(ans->last.x); printf(" = %16.15e\n", ans->last.x);
printf(" y = "); prhex(ans->last.y); printf(" = %16.15e\n", ans->last.y);
printf(" Computed x*y = "); prhex(ans->last.computed); printf(" = %16.15e\n", ans->last.computed);
printf(" Expected x*y = "); prhex(ans->last.expected); printf(" = %16.15e\n\n", ans->last.expected);
#endif
}
More information about the Numeric-interest
mailing list