www.digitalmars.com         C & C++   DMDScript  

digitalmars.D - Floating point feqrel(), final submission.

reply Don Clugston <dac nospam.com.au> writes:
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
next sibling parent "Ben Hinkle" <bhinkle mathworks.com> writes:
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
prev sibling parent "Walter" <newshound digitalmars.com> writes:
Nice work, Don. I've folded it in. Thanks!
Aug 21 2005