Porting bcltype/diyfp.cpp to managed code as shared/System/Number.DiyFp.cs

1 files changed, 129 insertions, 0 deletions

diff --git a/src/System.Private.CoreLib/shared/System/Number.DiyFp.cs b/src/System.Private.CoreLib/shared/System/Number.DiyFp.cs
new file mode 100644
index 0000000000..c6cab0d0b4
--- /dev/null
+++ b/src/System.Private.CoreLib/shared/System/Number.DiyFp.cs
@@ -0,0 +1,129 @@
+// Licensed to the .NET Foundation under one or more agreements.
+// The .NET Foundation licenses this file to you under the MIT license.
+// See the LICENSE file in the project root for more information.
+
+using System.Diagnostics;
+
+namespace System
+{
+    internal static partial class Number
+    {
+        // An exteneded floating-point data structure which is required by Grisu3 algorithm.
+        // It defines a 64-bit significand and a 32-bit exponent,
+        // which is EXTENDED compare to IEEE double precision floating-point number (53 bits significand and 11 bits exponent).
+        //
+        // Original Grisu algorithm produces suboptimal results. To shorten the output (which is part of Grisu2/Grisu3's job),
+        // we need additional 11 bits of the significand f and larger exponent e (A larger exponent range is used to avoid overflow. A 32-bit exponent is far big enough).
+        // To fully understand how Grisu3 uses this data structure to get better result, please read http://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf
+        internal struct DiyFp
+        {
+            public const int SignificandLength = 64;
+
+            // Extended significand.
+            // IEEE 754 double-precision numbers only require 53 bits significand.
+            // However, in Grisu3 we need additional 11 bits so we declare _f as ulong.
+            // Please note _f does not include sign bit.
+            private ulong _f;
+
+            // Extended exponent.
+            // IEEE 754 double-precision numbers only require 11 bits exponent.
+            // However, in Grisu3 we need extra space to avoid overflow so we declare _e as int.
+            // Please note _e is a biased exponent if the DiyFp instance is generated by GenerateNormalizedDiyFp().
+            private int _e;
+
+            public DiyFp(ulong f, int e)
+            {
+                _f = f;
+                _e = e;
+            }
+
+            public ulong f
+            {
+                get
+                {
+                    return _f;
+                }
+
+                set
+                {
+                    _f = value;
+                }
+            }
+
+            public int e
+            {
+                get
+                {
+                    return _e;
+                }
+
+                set
+                {
+                    _e = value;
+                }
+            }
+
+            public static void GenerateNormalizedDiyFp(double value, out DiyFp result)
+            {
+                Debug.Assert(value > 0.0);
+
+                long f = ExtractFractionAndBiasedExponent(value, out int e);
+
+                while ((f & (1L << 52)) == 0)
+                {
+                    f <<= 1;
+                    e--;
+                }
+
+                int lengthDelta = (SignificandLength - 53);
+                f <<= lengthDelta;
+                e -= lengthDelta;
+
+                result = new DiyFp((ulong)(f), e);
+            }
+
+            public static void Minus(ref DiyFp lhs, ref DiyFp rhs, out DiyFp result)
+            {
+                result = lhs;
+                result.Minus(ref rhs);
+            }
+
+            public static void Multiply(ref DiyFp lhs, ref DiyFp rhs, out DiyFp result)
+            {
+                result = lhs;
+                result.Multiply(ref rhs);
+            }
+
+            public void Minus(ref DiyFp rhs)
+            {
+                Debug.Assert(_e == rhs._e);
+                Debug.Assert(_f >= rhs._f);
+
+                _f -= rhs._f;
+            }
+
+            public void Multiply(ref DiyFp rhs)
+            {
+                ulong lf = _f;
+                ulong rf = rhs._f;
+
+                uint a = (uint)(lf >> 32);
+                uint b = (uint)(lf);
+
+                uint c = (uint)(rf >> 32);
+                uint d = (uint)(rf);
+
+                ulong ac = ((ulong)(a) * c);
+                ulong bc = ((ulong)(b) * c);
+                ulong ad = ((ulong)(a) * d);
+                ulong bd = ((ulong)(b) * d);
+
+                ulong tmp = (bd >> 32) + (uint)(ad) + (uint)(bc);
+                tmp += (1UL << 31);
+
+                _f = ac + (ad >> 32) + (bc >> 32) + (tmp >> 32);
+                _e += (rhs._e + SignificandLength);
+            }
+        }
+    }
+}