src/System.Private.CoreLib/shared/System/Number.NumberToDouble.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
   7 namespace System
   8 {
   9     internal unsafe partial class Number
  10     {
  11         // precomputed tables with powers of 10. These allows us to do at most
  12         // two Mul64 during the conversion. This is important not only
  13         // for speed, but also for precision because of Mul64 computes with 1 bit error.
  14
  15         private static readonly ulong[] s_Pow10MantissaTable = new ulong[]
  16         {
  17             // powers of 10
  18             0XA0000000_00000000,     // 1
  19             0XC8000000_00000000,     // 2
  20             0XFA000000_00000000,     // 3
  21             0X9C400000_00000000,     // 4
  22             0XC3500000_00000000,     // 5
  23             0XF4240000_00000000,     // 6
  24             0X98968000_00000000,     // 7
  25             0XBEBC2000_00000000,     // 8
  26             0XEE6B2800_00000000,     // 9
  27             0X9502F900_00000000,     // 10
  28             0XBA43B740_00000000,     // 11
  29             0XE8D4A510_00000000,     // 12
  30             0X9184E72A_00000000,     // 13
  31             0XB5E620F4_80000000,     // 14
  32             0XE35FA931_A0000000,     // 15
  33
  34             // powers of 0.1
  35             0xCCCCCCCC_CCCCCCCD,     // 1
  36             0xA3D70A3D_70A3D70B,     // 2
  37             0x83126E97_8D4FDF3C,     // 3
  38             0xD1B71758_E219652E,     // 4
  39             0xA7C5AC47_1B478425,     // 5
  40             0x8637BD05_AF6C69B7,     // 6
  41             0xD6BF94D5_E57A42BE,     // 7
  42             0xABCC7711_8461CEFF,     // 8
  43             0x89705F41_36B4A599,     // 9
  44             0xDBE6FECE_BDEDD5C2,     // 10
  45             0xAFEBFF0B_CB24AB02,     // 11
  46             0x8CBCCC09_6F5088CF,     // 12
  47             0xE12E1342_4BB40E18,     // 13
  48             0xB424DC35_095CD813,     // 14
  49             0x901D7CF7_3AB0ACDC,     // 15
  50         };
  51
  52         private static readonly short[] s_Pow10ExponentTable = new short[]
  53         {
  54             // exponents for both powers of 10 and 0.1
  55             4,      // 1
  56             7,      // 2
  57             10,     // 3
  58             14,     // 4
  59             17,     // 5
  60             20,     // 6
  61             24,     // 7
  62             27,     // 8
  63             30,     // 9
  64             34,     // 10
  65             37,     // 11
  66             40,     // 12
  67             44,     // 13
  68             47,     // 14
  69             50,     // 15
  70         };
  71
  72         private static readonly ulong[] s_Pow10By16MantissaTable = new ulong[]
  73         {
  74             // powers of 10^16
  75             0x8E1BC9BF_04000000,     // 1
  76             0x9DC5ADA8_2B70B59E,     // 2
  77             0xAF298D05_0E4395D6,     // 3
  78             0xC2781F49_FFCFA6D4,     // 4
  79             0xD7E77A8F_87DAF7FA,     // 5
  80             0xEFB3AB16_C59B14A0,     // 6
  81             0x850FADC0_9923329C,     // 7
  82             0x93BA47C9_80E98CDE,     // 8
  83             0xA402B9C5_A8D3A6E6,     // 9
  84             0xB616A12B_7FE617A8,     // 10
  85             0xCA28A291_859BBF90,     // 11
  86             0xE070F78D_39275566,     // 12
  87             0xF92E0C35_37826140,     // 13
  88             0x8A5296FF_E33CC92C,     // 14
  89             0x9991A6F3_D6BF1762,     // 15
  90             0xAA7EEBFB_9DF9DE8A,     // 16
  91             0xBD49D14A_A79DBC7E,     // 17
  92             0xD226FC19_5C6A2F88,     // 18
  93             0xE950DF20_247C83F8,     // 19
  94             0x81842F29_F2CCE373,     // 20
  95             0x8FCAC257_558EE4E2,     // 21
  96
  97             // powers of 0.1^16
  98             0xE69594BE_C44DE160,     // 1
  99             0xCFB11EAD_453994C3,     // 2
 100             0xBB127C53_B17EC165,     // 3
 101             0xA87FEA27_A539E9B3,     // 4
 102             0x97C560BA_6B0919B5,     // 5
 103             0x88B402F7_FD7553AB,     // 6
 104             0xF64335BC_F065D3A0,     // 7
 105             0xDDD0467C_64BCE4C4,     // 8
 106             0xC7CABA6E_7C5382ED,     // 9
 107             0xB3F4E093_DB73A0B7,     // 10
 108             0xA21727DB_38CB0053,     // 11
 109             0x91FF8377_5423CC29,     // 12
 110             0x8380DEA9_3DA4BC82,     // 13
 111             0xECE53CEC_4A314F00,     // 14
 112             0xD5605FCD_CF32E217,     // 15
 113             0xC0314325_637A1978,     // 16
 114             0xAD1C8EAB_5EE43BA2,     // 17
 115             0x9BECCE62_836AC5B0,     // 18
 116             0x8C71DCD9_BA0B495C,     // 19
 117             0xFD00B897_47823938,     // 20
 118             0xE3E27A44_4D8D991A,     // 21
 119         };
 120
 121         private static readonly short[] s_Pow10By16ExponentTable = new short[]
 122         {
 123             // exponents for both powers of 10^16 and 0.1^16
 124             54,     // 1
 125             107,    // 2
 126             160,    // 3
 127             213,    // 4
 128             266,    // 5
 129             319,    // 6
 130             373,    // 7
 131             426,    // 8
 132             479,    // 9
 133             532,    // 10
 134             585,    // 11
 135             638,    // 12
 136             691,    // 13
 137             745,    // 14
 138             798,    // 15
 139             851,    // 16
 140             904,    // 17
 141             957,    // 18
 142             1010,   // 19
 143             1064,   // 20
 144             1117,   // 21
 145         };
 146
 147 #if DEBUG
 148         private static bool s_CheckedTables = false;
 149
 150         private static void CheckTables()
 151         {
 152             ulong val;
 153             int exp;
 154
 155             val = 0xA0000000_00000000;
 156             exp = 4; // 10
 157
 158             CheckPow10MantissaTable(val, exp, new Span<ulong>(s_Pow10MantissaTable, 0, 15), nameof(s_Pow10MantissaTable));
 159             CheckPow10ExponentTable(val, exp, new Span<short>(s_Pow10ExponentTable, 0, 15), nameof(s_Pow10ExponentTable));
 160
 161             val = 0x8E1BC9BF_04000000;
 162             exp = 54; //10^16
 163
 164             CheckPow10MantissaTable(val, exp, new Span<ulong>(s_Pow10By16MantissaTable, 0, 21), nameof(s_Pow10By16MantissaTable));
 165             CheckPow10ExponentTable(val, exp, new Span<short>(s_Pow10By16ExponentTable, 0, 21), nameof(s_Pow10By16ExponentTable));
 166
 167             val = 0xCCCCCCCC_CCCCCCCD;
 168             exp = -3; // 0.1
 169
 170             CheckPow10MantissaTable(val, exp, new Span<ulong>(s_Pow10MantissaTable, 15, 15), nameof(s_Pow10MantissaTable) + " - inv");
 171
 172             val = 0xE69594BE_C44DE160;
 173             exp = -53; // 0.1^16
 174
 175             CheckPow10MantissaTable(val, exp, new Span<ulong>(s_Pow10By16MantissaTable, 21, 21), nameof(s_Pow10By16MantissaTable) + " - inv");
 176         }
 177
 178         // debug-only verification of the precomputed tables
 179         private static void CheckPow10MantissaTable(ulong val, int exp, Span<ulong> table, string name)
 180         {
 181             ulong multval = val;
 182             int mulexp = exp;
 183             bool isBad = false;
 184
 185             for (int i = 0; i < table.Length; i++)
 186             {
 187                 if (table[i] != val)
 188                 {
 189                     if (!isBad)
 190                     {
 191                         Debug.WriteLine(name);
 192                         isBad = true;
 193                     }
 194                     Debug.WriteLine($"0x{val:X16},     // {i + 1}");
 195                 }
 196
 197                 exp += mulexp;
 198                 val = Mul64Precise(val, multval, ref exp);
 199             }
 200
 201             Debug.Assert(!isBad, "NumberToDouble table not correct. Correct version dumped to debug output.");
 202         }
 203
 204         // debug-only verification of the precomputed tables
 205         private static void CheckPow10ExponentTable(ulong val, int exp, Span<short> table, string name)
 206         {
 207             ulong multval = val;
 208             int mulexp = exp;
 209             bool isBad = false;
 210
 211             for (int i = 0; i < table.Length; i++)
 212             {
 213                 if (table[i] != exp)
 214                 {
 215                     if (!isBad)
 216                     {
 217                         Debug.WriteLine(name);
 218                         isBad = true;
 219                     }
 220                     Debug.WriteLine($"{val}, // {i + 1}");
 221                 }
 222
 223                 exp += mulexp;
 224                 val = Mul64Precise(val, multval, ref exp);
 225             }
 226
 227             Debug.Assert(!isBad, "NumberToDouble table not correct. Correct version dumped to debug output.");
 228         }
 229
 230         // slower high precision version of Mul64 for computation of the tables
 231         private static ulong Mul64Precise(ulong a, ulong b, ref int exp)
 232         {
 233             ulong hilo = ((Mul32x32To64((uint)(a >> 32), (uint)(b)) >> 1)
 234                        + (Mul32x32To64((uint)(a), (uint)(b >> 32)) >> 1)
 235                        + (Mul32x32To64((uint)(a), (uint)(b)) >> 33)) >> 30;
 236
 237             ulong val = Mul32x32To64((uint)(a >> 32), (uint)(b >> 32))
 238                       + (hilo >> 1)
 239                       + (hilo & 1);
 240
 241             // normalize
 242             if ((val & 0x80000000_00000000) == 0)
 243             {
 244                 val <<= 1;
 245                 exp--;
 246             }
 247
 248             return val;
 249         }
 250 #endif
 251
 252         // get 32-bit integer from at most 9 digits
 253         private static uint DigitsToInt(char* p, int count)
 254         {
 255             Debug.Assert((1 <= count) && (count <= 9));
 256
 257             char* end = (p + count);
 258             uint res = (uint)(p[0] - '0');
 259
 260             for (p++; p < end; p++)
 261             {
 262                 res = (10 * res) + p[0] - '0';
 263             }
 264
 265             return res;
 266         }
 267
 268         private static int GetLength(char* src)
 269         {
 270             int length = 0;
 271
 272             while (src[length] != '\0')
 273             {
 274                 length++;
 275             }
 276
 277             return length;
 278         }
 279
 280         // helper function to multiply two 32-bit uints
 281         private static ulong Mul32x32To64(uint a, uint b)
 282         {
 283             return (ulong)(a) * b;
 284         }
 285
 286         // multiply two numbers in the internal integer representation
 287         private static ulong Mul64Lossy(ulong a, ulong b, ref int exp)
 288         {
 289             // it's ok to lose some precision here - Mul64 will be called
 290             // at most twice during the conversion, so the error won't propagate
 291             // to any of the 53 significant bits of the result
 292             ulong val = Mul32x32To64((uint)(a >> 32), (uint)(b >> 32))
 293                       + (Mul32x32To64((uint)(a >> 32), (uint)(b)) >> 32)
 294                       + (Mul32x32To64((uint)(a), (uint)(b >> 32)) >> 32);
 295
 296             // normalize
 297             if ((val & 0x80000000_00000000) == 0)
 298             {
 299                 val <<= 1;
 300                 exp--;
 301             }
 302
 303             return val;
 304         }
 305
 306         private static double NumberToDouble(ref NumberBuffer number)
 307         {
 308 #if DEBUG
 309
 310             if (!s_CheckedTables)
 311             {
 312                 CheckTables();
 313                 s_CheckedTables = true;
 314             }
 315 #endif
 316
 317             char* src = number.GetDigitsPointer();
 318             int total = GetLength(src);
 319             int remaining = total;
 320
 321             // skip the leading zeros
 322             while (src[0] == '0')
 323             {
 324                 remaining--;
 325                 src++;
 326             }
 327
 328             if (remaining == 0)
 329             {
 330                 return number.sign ? -0.0 : 0.0;
 331             }
 332
 333             int count = Math.Min(remaining, 9);
 334             remaining -= count;
 335             ulong val = DigitsToInt(src, count);
 336
 337             if (remaining > 0)
 338             {
 339                 count = Math.Min(remaining, 9);
 340                 remaining -= count;
 341
 342                 // get the denormalized power of 10
 343                 uint mult = (uint)(s_Pow10MantissaTable[count - 1] >> (64 - s_Pow10ExponentTable[count - 1]));
 344                 val = Mul32x32To64((uint)(val), mult) + DigitsToInt(src + 9, count);
 345             }
 346
 347             int scale = number.scale - (total - remaining);
 348             int absscale = Math.Abs(scale);
 349             if (absscale >= 22 * 16)
 350             {
 351                 // overflow / underflow
 352                 if (scale > 0)
 353                 {
 354                     return number.sign ? double.NegativeInfinity : double.PositiveInfinity;
 355                 }
 356                 else
 357                 {
 358                     return number.sign ? -0.0 : 0.0;
 359                 }
 360             }
 361
 362             int exp = 64;
 363
 364             // normalize the mantissa
 365             if ((val & 0xFFFFFFFF_00000000) == 0) { val <<= 32; exp -= 32; }
 366             if ((val & 0xFFFF0000_00000000) == 0) { val <<= 16; exp -= 16; }
 367             if ((val & 0xFF000000_00000000) == 0) { val <<= 8; exp -= 8; }
 368             if ((val & 0xF0000000_00000000) == 0) { val <<= 4; exp -= 4; }
 369             if ((val & 0xC0000000_00000000) == 0) { val <<= 2; exp -= 2; }
 370             if ((val & 0x80000000_00000000) == 0) { val <<= 1; exp -= 1; }
 371
 372             int index = absscale & 15;
 373             if (index != 0)
 374             {
 375                 int multexp = s_Pow10ExponentTable[index - 1];
 376                 // the exponents are shared between the inverted and regular table
 377                 exp += (scale < 0) ? (-multexp + 1) : multexp;
 378
 379                 ulong multval = s_Pow10MantissaTable[index + ((scale < 0) ? 15 : 0) - 1];
 380                 val = Mul64Lossy(val, multval, ref exp);
 381             }
 382
 383             index = absscale >> 4;
 384             if (index != 0)
 385             {
 386                 int multexp = s_Pow10By16ExponentTable[index - 1];
 387                 // the exponents are shared between the inverted and regular table
 388                 exp += (scale < 0) ? (-multexp + 1) : multexp;
 389
 390                 ulong multval = s_Pow10By16MantissaTable[index + ((scale < 0) ? 21 : 0) - 1];
 391                 val = Mul64Lossy(val, multval, ref exp);
 392             }
 393
 394             // round & scale down
 395             if (((uint)(val) & (1 << 10)) != 0)
 396             {
 397                 // IEEE round to even
 398                 ulong tmp = val + ((1 << 10) - 1) + (((uint)(val) >> 11) & 1);
 399                 if (tmp < val)
 400                 {
 401                     // overflow
 402                     tmp = (tmp >> 1) | 0x8000000000000000;
 403                     exp += 1;
 404                 }
 405                 val = tmp;
 406             }
 407
 408             // return the exponent to a biased state
 409             exp += 0x3FE;
 410
 411             // handle overflow, underflow, "Epsilon - 1/2 Epsilon", denormalized, and the normal case
 412             if (exp <= 0)
 413             {
 414                 if (exp == -52 && (val >= 0x8000000000000058))
 415                 {
 416                     // round X where {Epsilon > X >= 2.470328229206232730000000E-324} up to Epsilon (instead of down to zero)
 417                     val = 0x0000000000000001;
 418                 }
 419                 else if (exp <= -52)
 420                 {
 421                     // underflow
 422                     val = 0;
 423                 }
 424                 else
 425                 {
 426                     // denormalized
 427                     val >>= (-exp + 11 + 1);
 428                 }
 429             }
 430             else if (exp >= 0x7FF)
 431             {
 432                 // overflow
 433                 val = 0x7FF0000000000000;
 434             }
 435             else
 436             {
 437                 // normal postive exponent case
 438                 val = ((ulong)(exp) << 52) + ((val >> 11) & 0x000FFFFFFFFFFFFF);
 439             }
 440
 441             if (number.sign)
 442             {
 443                 val |= 0x8000000000000000;
 444             }
 445
 446             return BitConverter.Int64BitsToDouble((long)(val));
 447         }
 448     }
 449 }