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

// 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);
            }
        }
    }
}