Posted by: leppie | March 5, 2009

System.Double deconstruction extension methods

Recently I discovered an issue with the = procedure in IronScheme, anyways, I was interpreting the semantics incorrectly, and realised I had to deconstruct my floating point representation (System.Double) to get the exact representation of the number, then present it as a rational to check for mathematically correct equality.

NOTE: this code needs some form of big number (the one used here is from the DLR) that can do bitshifting or big powers (up to 1023) of 2.

static class DoubleExtensions
{
  const long MANTISSA_MASK = (long)(ulong.MaxValue >> (64 - 52));
  const long EXPONENT_MASK = (long)(ulong.MaxValue >> (64 - 11));
  const long SIGN_MASK = 0x1;

  public static long GetMantissa(this double d)
  {
    var r = BitConverter.DoubleToInt64Bits(d);
    return GetMantissa(r);
  }

  static long GetMantissa(long r)
  {
    return r & MANTISSA_MASK;
  }

  public static int GetExponent(this double d)
  {
    var r = BitConverter.DoubleToInt64Bits(d);
    return GetExponent(r);
  }

  static int GetExponent(long r)
  {
    r >>= 52;
    return (int)(r & EXPONENT_MASK);
  }

  public static int GetSign(this double d)
  {
    var r = BitConverter.DoubleToInt64Bits(d);
    return GetSign(r);
  }

  static int GetSign(long r)
  {
    r >>= 63;
    return (int)(r & SIGN_MASK);
  }

  readonly static BigInteger MANTISSA = (ulong.MaxValue >> (64 - 52)) + 1;

  public static bool GetComponents(this double d, out BigInteger numerator, out BigInteger denominator)
  {
    if (double.IsNaN(d))
    {
      numerator = denominator = 0;
      return false;
    }

    if (double.IsNegativeInfinity(d))
    {
      numerator = -1;
      denominator = 0;
      return false;
    }

    if (double.IsPositiveInfinity(d))
    {
      numerator = 1;
      denominator = 0;
      return false;
    }

    if (d == 0.0)
    {
      numerator = 0;
      denominator = 1;
      return true;
    }

    var r = BitConverter.DoubleToInt64Bits(d);

    var m = GetMantissa(r);
    var e = GetExponent(r);
    var s = GetSign(r) == 0 ? 1 : -1;

    const int BIAS = 1023;
    var re = e - BIAS;

    var exp = (((BigInteger)1) << Math.Abs(re));

    if (re < 0)
    {
      denominator = MANTISSA * s * exp;
      numerator = (MANTISSA + m);
    }
    else
    {
      denominator = MANTISSA * s;
      numerator = exp * (MANTISSA + m);
    }

    return true;
  }
}
Advertisements

Leave a Reply

Please log in using one of these methods to post your comment:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Categories

%d bloggers like this: