digitalmars.D - Floating point feqrel(), final submission.
- Don Clugston (86/86) Aug 17 2005 I've modified my code based on comments by Ben.
- Ben Hinkle (13/99) Aug 18 2005 Nice. Since the property mant_dig is available from variables as well as...
- Walter (1/1) Aug 21 2005 Nice work, Don. I've folded it in. Thanks!
I've modified my code based on comments by Ben.
It's now less ambitious, and doesn't try to do anything
with numbers that vary by a factor of more than 2.
This has allowed me to streamline the code somewhat.
It is generally not suitable for comparisons near zero or infinity,
unless you demand an exact match.
FWIW, this function should be almost as fast as ">"
on a PII, because it doesn't use any (slow) floating-point
compares! All branches except for the first == one are highly predictable.
Ready to be slotted into std.math or similar. Enjoy!
-Don.
========================================================
/*
int feqrel(real x, real y)
To what precision is x equal to y?
Public Domain. Author: Don Clugston, 18 Aug 2005.
Returns the number of mantissa bits which are equal in x and y.
eg, 0x1.F8p+60 and 0x1.F1p+60 are equal to 5 bits of precision.
If x == y, then feqrel(x, y) == real.mant_dig.
If x and y differ by a factor of two or more, or if one or both
is a nan, the return value is 0.
*/
int feqrel(real a, real b)
{
if (a==b) return real.mant_dig; // ensure diff!=0, cope with INF.
real diff = fabs(a-b);
ushort *pa = cast(ushort *)(&a);
ushort *pb = cast(ushort *)(&b);
ushort *pd = cast(ushort *)(&diff);
// The difference in abs(exponent) between a or b and abs(a-b)
// is equal to the number of mantissa bits of a which are
// equal to b. If negative, a and b have different exponents.
// If positive, a and b are equal to 'bitsdiff' bits.
// AND with 0x7FFF to form the absolute value.
// To avoid out-by-1 errors, we subtract 1 so it rounds down
// if the exponents were different. This means 'bitsdiff' is
// always 1 lower than we want, except that if bitsdiff==0,
// they could have 0 or 1 bits in common.
int bitsdiff = ( ((pa[4]&0x7FFF) + (pb[4]&0x7FFF)-1)>>1) - pd[4];
if (pd[4]== 0) { // Difference is denormal
// For denormals, we need to add the number of zeros that
// lie at the start of diff's mantissa.
// We do this by multiplying by 2^real.mant_dig
diff*=0x1p+63;
return bitsdiff + real.mant_dig - pd[4];
}
if (bitsdiff>0) return bitsdiff+1; // add the 1 we subtracted before
// Avoid out-by-1 errors when factor is almost 2.
return bitsdiff==0 ? pa[4]==pb[4] : 0;
}
unittest {
// Exact equality
assert(feqrel(real.max,real.max)==real.mant_dig);
assert(feqrel(0,0)==real.mant_dig);
assert(feqrel(7.1824,7.1824)==real.mant_dig);
assert(feqrel(real.infinity,real.infinity)==real.mant_dig);
// a few bits away from exact equality
real w=1;
for (int i=1; i<real.mant_dig-1; ++i) {
assert(feqrel(1+w*real.epsilon,1)==real.mant_dig-i);
assert(feqrel(1-w*real.epsilon,1)==real.mant_dig-i);
assert(feqrel(1,1+(w-1)*real.epsilon)==real.mant_dig-i+1);
w*=2;
}
assert(feqrel(1.5+real.epsilon,1.5)==real.mant_dig-1);
assert(feqrel(1.5-real.epsilon,1.5)==real.mant_dig-1);
assert(feqrel(1.5-real.epsilon,1.5+real.epsilon)==real.mant_dig-2);
// Numbers that are close
assert(feqrel(0x1.Bp+84, 0x1.B8p+84)==5);
assert(feqrel(0x1.8p+10, 0x1.Cp+10)==2);
assert(feqrel(1.5*(1-real.epsilon), 1)==2);
assert(feqrel(1.5, 1)==1);
assert(feqrel(2*(1-real.epsilon), 1)==1);
// Factors of 2
assert(feqrel(real.max,real.infinity)==0);
assert(feqrel(2*(1-real.epsilon), 1)==1);
assert(feqrel(1, 2)==0);
assert(feqrel(4, 1)==0);
// Extreme inequality
assert(feqrel(real.nan,real.nan)==0);
assert(feqrel(0,-real.nan)==0);
assert(feqrel(real.nan,real.infinity)==0);
assert(feqrel(real.infinity,-real.infinity)==0);
assert(feqrel(-real.max,real.infinity)==0);
assert(feqrel(real.max,-real.max)==0);
}
Aug 17 2005
"Ben Hinkle" <bhinkle mathworks.com>
Nice. Since the property mant_dig is available from variables as well as
types I see typical user code as looking like
double x = 1.0;
double y = ...;
... perform some computation ...
if (feqrel(x,y) > x.mant_dig/2) {
... x and y agree to at least 1/2 of x's precision
}
It would be great if it went into std.math. I'd add it to the phobos doc for
std.math and send Walter the code+doc. I'm guessing removing feq from
std.math and replacing with feqrel would be nice, too.
"Don Clugston" <dac nospam.com.au> wrote in message
news:de0o62$pbl$1 digitaldaemon.com...
I've modified my code based on comments by Ben.
It's now less ambitious, and doesn't try to do anything
with numbers that vary by a factor of more than 2.
This has allowed me to streamline the code somewhat.
It is generally not suitable for comparisons near zero or infinity,
unless you demand an exact match.
FWIW, this function should be almost as fast as ">"
on a PII, because it doesn't use any (slow) floating-point
compares! All branches except for the first == one are highly predictable.
Ready to be slotted into std.math or similar. Enjoy!
-Don.
========================================================
/*
int feqrel(real x, real y)
To what precision is x equal to y?
Public Domain. Author: Don Clugston, 18 Aug 2005.
Returns the number of mantissa bits which are equal in x and y.
eg, 0x1.F8p+60 and 0x1.F1p+60 are equal to 5 bits of precision.
If x == y, then feqrel(x, y) == real.mant_dig.
If x and y differ by a factor of two or more, or if one or both
is a nan, the return value is 0.
*/
int feqrel(real a, real b)
{
if (a==b) return real.mant_dig; // ensure diff!=0, cope with INF.
real diff = fabs(a-b);
ushort *pa = cast(ushort *)(&a);
ushort *pb = cast(ushort *)(&b);
ushort *pd = cast(ushort *)(&diff);
// The difference in abs(exponent) between a or b and abs(a-b)
// is equal to the number of mantissa bits of a which are
// equal to b. If negative, a and b have different exponents.
// If positive, a and b are equal to 'bitsdiff' bits.
// AND with 0x7FFF to form the absolute value.
// To avoid out-by-1 errors, we subtract 1 so it rounds down
// if the exponents were different. This means 'bitsdiff' is
// always 1 lower than we want, except that if bitsdiff==0,
// they could have 0 or 1 bits in common.
int bitsdiff = ( ((pa[4]&0x7FFF) + (pb[4]&0x7FFF)-1)>>1) - pd[4];
if (pd[4]== 0) { // Difference is denormal
// For denormals, we need to add the number of zeros that
// lie at the start of diff's mantissa.
// We do this by multiplying by 2^real.mant_dig
diff*=0x1p+63;
return bitsdiff + real.mant_dig - pd[4];
}
if (bitsdiff>0) return bitsdiff+1; // add the 1 we subtracted before
// Avoid out-by-1 errors when factor is almost 2.
return bitsdiff==0 ? pa[4]==pb[4] : 0;
}
unittest {
// Exact equality
assert(feqrel(real.max,real.max)==real.mant_dig);
assert(feqrel(0,0)==real.mant_dig);
assert(feqrel(7.1824,7.1824)==real.mant_dig);
assert(feqrel(real.infinity,real.infinity)==real.mant_dig);
// a few bits away from exact equality
real w=1;
for (int i=1; i<real.mant_dig-1; ++i) {
assert(feqrel(1+w*real.epsilon,1)==real.mant_dig-i);
assert(feqrel(1-w*real.epsilon,1)==real.mant_dig-i);
assert(feqrel(1,1+(w-1)*real.epsilon)==real.mant_dig-i+1);
w*=2;
}
assert(feqrel(1.5+real.epsilon,1.5)==real.mant_dig-1);
assert(feqrel(1.5-real.epsilon,1.5)==real.mant_dig-1);
assert(feqrel(1.5-real.epsilon,1.5+real.epsilon)==real.mant_dig-2);
// Numbers that are close
assert(feqrel(0x1.Bp+84, 0x1.B8p+84)==5);
assert(feqrel(0x1.8p+10, 0x1.Cp+10)==2);
assert(feqrel(1.5*(1-real.epsilon), 1)==2);
assert(feqrel(1.5, 1)==1);
assert(feqrel(2*(1-real.epsilon), 1)==1);
// Factors of 2
assert(feqrel(real.max,real.infinity)==0);
assert(feqrel(2*(1-real.epsilon), 1)==1);
assert(feqrel(1, 2)==0);
assert(feqrel(4, 1)==0);
// Extreme inequality
assert(feqrel(real.nan,real.nan)==0);
assert(feqrel(0,-real.nan)==0);
assert(feqrel(real.nan,real.infinity)==0);
assert(feqrel(real.infinity,-real.infinity)==0);
assert(feqrel(-real.max,real.infinity)==0);
assert(feqrel(real.max,-real.max)==0);
}
Aug 18 2005
