From: Tanner Gooding Date: Mon, 17 Sep 2018 05:34:40 +0000 (-0700) Subject: Porting the Dragon4 algorithm to managed code. X-Git-Tag: submit/tizen/20210909.063632~11030^2~3866 X-Git-Url: http://review.tizen.org/git/?a=commitdiff_plain;h=830852a4787e89996bcc711d184b4882ca268bed;p=platform%2Fupstream%2Fdotnet%2Fruntime.git Porting the Dragon4 algorithm to managed code. Commit migrated from https://github.com/dotnet/coreclr/commit/b9e981c3e8fc1c718958ca3e0433ba6a6e3cf3cc --- diff --git a/src/coreclr/src/System.Private.CoreLib/src/System/Number.CoreCLR.cs b/src/coreclr/src/System.Private.CoreLib/src/System/Number.CoreCLR.cs index 30e70d0..310987a 100644 --- a/src/coreclr/src/System.Private.CoreLib/src/System/Number.CoreCLR.cs +++ b/src/coreclr/src/System.Private.CoreLib/src/System/Number.CoreCLR.cs @@ -9,9 +9,6 @@ namespace System internal static partial class Number { [MethodImpl(MethodImplOptions.InternalCall)] - public static extern void DoubleToNumber(double value, int precision, ref NumberBuffer number); - - [MethodImpl(MethodImplOptions.InternalCall)] public static extern double NumberToDouble(ref NumberBuffer number); } } diff --git a/src/libraries/System.Private.CoreLib/src/System.Private.CoreLib.Shared.projitems b/src/libraries/System.Private.CoreLib/src/System.Private.CoreLib.Shared.projitems index b1d9b4a..bd0b571 100644 --- a/src/libraries/System.Private.CoreLib/src/System.Private.CoreLib.Shared.projitems +++ b/src/libraries/System.Private.CoreLib/src/System.Private.CoreLib.Shared.projitems @@ -278,6 +278,7 @@ + diff --git a/src/libraries/System.Private.CoreLib/src/System/Double.cs b/src/libraries/System.Private.CoreLib/src/System/Double.cs index d85b4bc..a4c0119 100644 --- a/src/libraries/System.Private.CoreLib/src/System/Double.cs +++ b/src/libraries/System.Private.CoreLib/src/System/Double.cs @@ -44,6 +44,22 @@ namespace System // We use this explicit definition to avoid the confusion between 0.0 and -0.0. internal const double NegativeZero = -0.0; + [NonVersionable] + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static int GetExponent(double d) + { + var bits = BitConverter.DoubleToInt64Bits(d); + return (int)((bits >> 52) & 0x7FF); + } + + [NonVersionable] + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static long GetMantissa(double d) + { + var bits = BitConverter.DoubleToInt64Bits(d); + return (bits & 0xFFFFFFFFFFFFF); + } + /// Determines whether the specified value is finite (zero, subnormal, or normal). [NonVersionable] [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/src/libraries/System.Private.CoreLib/src/System/Number.Dragon4.cs b/src/libraries/System.Private.CoreLib/src/System/Number.Dragon4.cs new file mode 100644 index 0000000..ec2ae16 --- /dev/null +++ b/src/libraries/System.Private.CoreLib/src/System/Number.Dragon4.cs @@ -0,0 +1,257 @@ +// 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.Runtime.InteropServices; +using System.Runtime.CompilerServices; +using Internal.Runtime.CompilerServices; + +namespace System +{ + internal static partial class Number + { + private static unsafe void Dragon4(double value, int precision, ref NumberBuffer number) + { + // ======================================================================================================================================== + // This implementation is based on the paper: https://www.cs.indiana.edu/~dyb/pubs/FP-Printing-PLDI96.pdf + // Besides the paper, some of the code and ideas are modified from http://www.ryanjuckett.com/programming/printing-floating-point-numbers/ + // You must read these two materials to fully understand the code. + // + // Note: we only support fixed format input. + // ======================================================================================================================================== + // + // Overview: + // + // The input double number can be represented as: + // value = f * 2^e = r / s. + // + // f: the output mantissa. Note: f is not the 52 bits mantissa of the input double number. + // e: biased exponent. + // r: numerator. + // s: denominator. + // k: value = d0.d1d2 . . . dn * 10^k + + // Step 1: + // Extract meta data from the input double value. + // + // Refer to IEEE double precision floating point format. + ulong f = (ulong)(ExtractFractionAndBiasedExponent(value, out int e)); + int mantissaHighBitIndex = (e == -1074) ? (int)(BigInteger.LogBase2(f)) : 52; + + // Step 2: + // Estimate k. We'll verify it and fix any error later. + // + // This is an improvement of the estimation in the original paper. + // Inspired by http://www.ryanjuckett.com/programming/printing-floating-point-numbers/ + // + // LOG10V2 = 0.30102999566398119521373889472449 + // DRIFT_FACTOR = 0.69 = 1 - log10V2 - epsilon (a small number account for drift of floating point multiplication) + int k = (int)(Math.Ceiling(((mantissaHighBitIndex + e) * Log10V2) - DriftFactor)); + + // Step 3: + // Store the input double value in BigInteger format. + // + // To keep the precision, we represent the double value as r/s. + // We have several optimization based on following table in the paper. + // + // ---------------------------------------------------------------------------------------------------------- + // | e >= 0 | e < 0 | + // ---------------------------------------------------------------------------------------------------------- + // | f != b^(P - 1) | f = b^(P - 1) | e = min exp or f != b^(P - 1) | e > min exp and f = b^(P - 1) | + // -------------------------------------------------------------------------------------------------------------- + // | r | f * b^e * 2 | f * b^(e + 1) * 2 | f * 2 | f * b * 2 | + // -------------------------------------------------------------------------------------------------------------- + // | s | 2 | b * 2 | b^(-e) * 2 | b^(-e + 1) * 2 | + // -------------------------------------------------------------------------------------------------------------- + // + // Note, we do not need m+ and m- because we only support fixed format input here. + // m+ and m- are used for free format input, which need to determine the exact range of values + // that would round to value when input so that we can generate the shortest correct number.digits. + // + // In our case, we just output number.digits until reaching the expected precision. + + var r = new BigInteger(f); + var s = new BigInteger(0); + + if (e >= 0) + { + // When f != b^(P - 1): + // r = f * b^e * 2 + // s = 2 + // value = r / s = f * b^e * 2 / 2 = f * b^e / 1 + // + // When f = b^(P - 1): + // r = f * b^(e + 1) * 2 + // s = b * 2 + // value = r / s = f * b^(e + 1) * 2 / b * 2 = f * b^e / 1 + // + // Therefore, we can simply say that when e >= 0: + // r = f * b^e = f * 2^e + // s = 1 + + r.ShiftLeft((uint)(e)); + s.SetUInt64(1); + } + else + { + // When e = min exp or f != b^(P - 1): + // r = f * 2 + // s = b^(-e) * 2 + // value = r / s = f * 2 / b^(-e) * 2 = f / b^(-e) + // + // When e > min exp and f = b^(P - 1): + // r = f * b * 2 + // s = b^(-e + 1) * 2 + // value = r / s = f * b * 2 / b^(-e + 1) * 2 = f / b^(-e) + // + // Therefore, we can simply say that when e < 0: + // r = f + // s = b^(-e) = 2^(-e) + + BigInteger.ShiftLeft(1, (uint)(-e), ref s); + } + + // According to the paper, we should use k >= 0 instead of k > 0 here. + // However, if k = 0, both r and s won't be changed, we don't need to do any operation. + // + // Following are the Scheme code from the paper: + // -------------------------------------------------------------------------------- + // (if (>= est 0) + // (fixup r (* s (exptt B est)) m+ m− est B low-ok? high-ok? ) + // (let ([scale (exptt B (− est))]) + // (fixup (* r scale) s (* m+ scale) (* m− scale) est B low-ok? high-ok? )))) + // -------------------------------------------------------------------------------- + // + // If est is 0, (* s (exptt B est)) = s, (* r scale) = (* r (exptt B (− est)))) = r. + // + // So we just skip when k = 0. + + if (k > 0) + { + BigInteger poweredValue = new BigInteger(0); + BigInteger.Pow10((uint)(k), ref poweredValue); + s.Multiply(ref poweredValue); + } + else if (k < 0) + { + BigInteger poweredValue = new BigInteger(0); + BigInteger.Pow10((uint)(-k), ref poweredValue); + r.Multiply(ref poweredValue); + } + + if (BigInteger.Compare(ref r, ref s) >= 0) + { + // The estimation was incorrect. Fix the error by increasing 1. + k += 1; + } + else + { + r.Multiply10(); + } + + number.scale = (k - 1); + + // This the prerequisite of calling BigInteger.HeuristicDivide(). + BigInteger.PrepareHeuristicDivide(ref r, ref s); + + // Step 4: + // Calculate number.digits. + // + // Output number.digits until reaching the last but one precision or the numerator becomes zero. + + int digitsNum = 0; + int currentDigit = 0; + + while (true) + { + currentDigit = (int)(BigInteger.HeuristicDivide(ref r, ref s)); + + if (r.IsZero() || ((digitsNum + 1) == precision)) + { + break; + } + + number.digits[digitsNum] = (char)('0' + currentDigit); + digitsNum++; + + r.Multiply10(); + } + + // Step 5: + // Set the last digit. + // + // We round to the closest digit by comparing value with 0.5: + // compare( value, 0.5 ) + // = compare( r / s, 0.5 ) + // = compare( r, 0.5 * s) + // = compare(2 * r, s) + // = compare(r << 1, s) + + r.ShiftLeft(1); + int compareResult = BigInteger.Compare(ref r, ref s); + bool isRoundDown = compareResult < 0; + + // We are in the middle, round towards the even digit (i.e. IEEE rouding rules) + if (compareResult == 0) + { + isRoundDown = (currentDigit & 1) == 0; + } + + if (isRoundDown) + { + number.digits[digitsNum] = (char)('0' + currentDigit); + digitsNum++; + } + else + { + char* pCurrentDigit = (number.digits + digitsNum); + + // Rounding up for 9 is special. + if (currentDigit == 9) + { + // find the first non-nine prior digit + while (true) + { + // If we are at the first digit + if (pCurrentDigit == number.digits) + { + // Output 1 at the next highest exponent + *pCurrentDigit = '1'; + digitsNum++; + number.scale += 1; + break; + } + + pCurrentDigit--; + digitsNum--; + + if (*pCurrentDigit != '9') + { + // increment the digit + *pCurrentDigit += (char)(1); + digitsNum++; + break; + } + } + } + else + { + // It's simple if the digit is not 9. + *pCurrentDigit = (char)('0' + currentDigit + 1); + digitsNum++; + } + } + + while (digitsNum < precision) + { + number.digits[digitsNum] = '0'; + digitsNum++; + } + + number.digits[precision] = '\0'; + + number.scale++; + number.sign = double.IsNegative(value); + } + } +} diff --git a/src/libraries/System.Private.CoreLib/src/System/Number.Formatting.cs b/src/libraries/System.Private.CoreLib/src/System/Number.Formatting.cs index 3d3c15b..d3d9279 100644 --- a/src/libraries/System.Private.CoreLib/src/System/Number.Formatting.cs +++ b/src/libraries/System.Private.CoreLib/src/System/Number.Formatting.cs @@ -1,4 +1,4 @@ -// Licensed to the .NET Foundation under one or more agreements. +// 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. @@ -2285,5 +2285,58 @@ SkipRounding: value /= 1000000000; return rem; } + + private static unsafe void DoubleToNumber(double value, int precision, ref NumberBuffer number) + { + number.precision = precision; + + if (double.GetExponent(value) == 0x7FF) + { + number.scale = double.IsNaN(value) ? ScaleNAN : ScaleINF; + number.sign = double.IsNegative(value); + number.digits[0] = '\0'; + } + else if (value == 0.0) + { + number.scale = 0; + number.sign = double.IsNegative(value); + number.digits[0] = '\0'; + } + else + { + Dragon4(value, precision, ref number); + } + } + + private static long ExtractFractionAndBiasedExponent(double value, out int exponent) + { + long fraction = double.GetMantissa(value); + exponent = double.GetExponent(value); + + if (exponent != 0) + { + // For normalized value, according to https://en.wikipedia.org/wiki/Double-precision_floating-point_format + // value = 1.fraction * 2^(exp - 1023) + // = (1 + mantissa / 2^52) * 2^(exp - 1023) + // = (2^52 + mantissa) * 2^(exp - 1023 - 52) + // + // So f = (2^52 + mantissa), e = exp - 1075; + + fraction |= ((long)(1) << 52); + exponent -= 1075; + } + else + { + // For denormalized value, according to https://en.wikipedia.org/wiki/Double-precision_floating-point_format + // value = 0.fraction * 2^(1 - 1023) + // = (mantissa / 2^52) * 2^(-1022) + // = mantissa * 2^(-1022 - 52) + // = mantissa * 2^(-1074) + // So f = mantissa, e = -1074 + exponent = -1074; + } + + return fraction; + } } } diff --git a/src/libraries/System.Private.CoreLib/src/System/Number.NumberBuffer.cs b/src/libraries/System.Private.CoreLib/src/System/Number.NumberBuffer.cs index 1730cf1..fe9eaf6 100644 --- a/src/libraries/System.Private.CoreLib/src/System/Number.NumberBuffer.cs +++ b/src/libraries/System.Private.CoreLib/src/System/Number.NumberBuffer.cs @@ -12,6 +12,11 @@ namespace System { private const int NumberMaxDigits = 50; // needs to == NUMBER_MAXDIGITS in coreclr's src/classlibnative/bcltype/number.h. + private const double Log10V2 = 0.30102999566398119521373889472449; + + // DriftFactor = 1 - Log10V2 - epsilon (a small number account for drift of floating point multiplication) + private const double DriftFactor = 0.69; + [StructLayout(LayoutKind.Sequential, Pack = 1)] internal unsafe ref struct NumberBuffer // needs to match layout of NUMBER in coreclr's src/classlibnative/bcltype/number.h {