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);
> }
| 
Permalink
Reply
