netcore/System.Private.CoreLib/shared/System/Number.DiyFp.cs

   1 // Licensed to the .NET Foundation under one or more agreements.
   2 // The .NET Foundation licenses this file to you under the MIT license.
   3 // See the LICENSE file in the project root for more information.
   4
   5 using System.Diagnostics;
   6 using System.Numerics;
   7
   8 namespace System
   9 {
  10     internal static partial class Number
  11     {
  12         // This is a port of the `DiyFp` implementation here: https://github.com/google/double-conversion/blob/a711666ddd063eb1e4b181a6cb981d39a1fc8bac/double-conversion/diy-fp.h
  13         // The backing structure and how it is used is described in more detail here: http://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf
  14
  15         // This "Do It Yourself Floating Point" class implements a floating-point number with a ulong significand and an int exponent.
  16         // Normalized DiyFp numbers will have the most significant bit of the significand set.
  17         // Multiplication and Subtraction do not normalize their results.
  18         // DiyFp are not designed to contain special doubles (NaN and Infinity).
  19         internal readonly ref struct DiyFp
  20         {
  21             public const int DoubleImplicitBitIndex = 52;
  22             public const int SingleImplicitBitIndex = 23;
  23
  24             public const int SignificandSize = 64;
  25
  26             public readonly ulong f;
  27             public readonly int e;
  28
  29             // Computes the two boundaries of value.
  30             //
  31             // The bigger boundary (mPlus) is normalized.
  32             // The lower boundary has the same exponent as mPlus.
  33             //
  34             // Precondition:
  35             //  The value encoded by value must be greater than 0.
  36             public static DiyFp CreateAndGetBoundaries(double value, out DiyFp mMinus, out DiyFp mPlus)
  37             {
  38                 var result = new DiyFp(value);
  39                 result.GetBoundaries(DoubleImplicitBitIndex, out mMinus, out mPlus);
  40                 return result;
  41             }
  42
  43             // Computes the two boundaries of value.
  44             //
  45             // The bigger boundary (mPlus) is normalized.
  46             // The lower boundary has the same exponent as mPlus.
  47             //
  48             // Precondition:
  49             //  The value encoded by value must be greater than 0.
  50             public static DiyFp CreateAndGetBoundaries(float value, out DiyFp mMinus, out DiyFp mPlus)
  51             {
  52                 var result = new DiyFp(value);
  53                 result.GetBoundaries(SingleImplicitBitIndex, out mMinus, out mPlus);
  54                 return result;
  55             }
  56
  57             public DiyFp(double value)
  58             {
  59                 Debug.Assert(double.IsFinite(value));
  60                 Debug.Assert(value > 0.0);
  61                 f = ExtractFractionAndBiasedExponent(value, out e);
  62             }
  63
  64             public DiyFp(float value)
  65             {
  66                 Debug.Assert(float.IsFinite(value));
  67                 Debug.Assert(value > 0.0f);
  68                 f = ExtractFractionAndBiasedExponent(value, out e);
  69             }
  70
  71             public DiyFp(ulong f, int e)
  72             {
  73                 this.f = f;
  74                 this.e = e;
  75             }
  76
  77             public DiyFp Multiply(in DiyFp other)
  78             {
  79                 // Simply "emulates" a 128-bit multiplication
  80                 //
  81                 // However: the resulting number only contains 64-bits. The least
  82                 // signficant 64-bits are only used for rounding the most significant
  83                 // 64-bits.
  84
  85                 uint a = (uint)(f >> 32);
  86                 uint b = (uint)(f);
  87
  88                 uint c = (uint)(other.f >> 32);
  89                 uint d = (uint)(other.f);
  90
  91                 ulong ac = ((ulong)(a) * c);
  92                 ulong bc = ((ulong)(b) * c);
  93                 ulong ad = ((ulong)(a) * d);
  94                 ulong bd = ((ulong)(b) * d);
  95
  96                 ulong tmp = (bd >> 32) + (uint)(ad) + (uint)(bc);
  97
  98                 // By adding (1UL << 31) to tmp, we round the final result.
  99                 // Halfway cases will be rounded up.
 100
 101                 tmp += (1U << 31);
 102
 103                 return new DiyFp(ac + (ad >> 32) + (bc >> 32) + (tmp >> 32), e + other.e + SignificandSize);
 104             }
 105
 106             public DiyFp Normalize()
 107             {
 108                 // This method is mainly called for normalizing boundaries.
 109                 //
 110                 // We deviate from the reference implementation by just using
 111                 // our LeadingZeroCount function so that we only need to shift
 112                 // and subtract once.
 113
 114                 Debug.Assert(f != 0);
 115                 int lzcnt = BitOperations.LeadingZeroCount(f);
 116                 return new DiyFp(f << lzcnt, e - lzcnt);
 117             }
 118
 119             // The exponents of both numbers must be the same.
 120             // The significand of 'this' must be bigger than the significand of 'other'.
 121             // The result will not be normalized.
 122             public DiyFp Subtract(in DiyFp other)
 123             {
 124                 Debug.Assert(e == other.e);
 125                 Debug.Assert(f >= other.f);
 126                 return new DiyFp(f - other.f, e);
 127             }
 128
 129             private void GetBoundaries(int implicitBitIndex, out DiyFp mMinus, out DiyFp mPlus)
 130             {
 131                 mPlus = new DiyFp((f << 1) + 1, e - 1).Normalize();
 132
 133                 // The boundary is closer if the sigificand is of the form:
 134                 //      f == 2^p-1
 135                 //
 136                 // Think of v = 1000e10 and v- = 9999e9
 137                 // Then the boundary == (v - v-) / 2 is not just at a distance of 1e9 but at a distance of 1e8.
 138                 // The only exception is for the smallest normal, where the largest denormal is at the same distance as its successor.
 139                 //
 140                 // Note: denormals have the same exponent as the smallest normals.
 141
 142                 // We deviate from the reference implementation by just checking if the significand has only the implicit bit set.
 143                 // In this scenario, we know that all the explicit bits are 0 and that the unbiased exponent is non-zero.
 144                 if (f == (1UL << implicitBitIndex))
 145                 {
 146                     mMinus = new DiyFp((f << 2) - 1, e - 2);
 147                 }
 148                 else
 149                 {
 150                     mMinus = new DiyFp((f << 1) - 1, e - 1);
 151                 }
 152
 153                 mMinus = new DiyFp(mMinus.f << (mMinus.e - mPlus.e), mPlus.e);
 154             }
 155         }
 156     }
 157 }