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.
5 using System.Diagnostics;
9 internal unsafe partial class Number
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.
15 private static readonly ulong[] s_Pow10MantissaTable = new ulong[]
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
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
52 private static readonly short[] s_Pow10ExponentTable = new short[]
54 // exponents for both powers of 10 and 0.1
72 private static readonly ulong[] s_Pow10By16MantissaTable = new ulong[]
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
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
121 private static readonly short[] s_Pow10By16ExponentTable = new short[]
123 // exponents for both powers of 10^16 and 0.1^16
148 private static bool s_CheckedTables = false;
150 private static void CheckTables()
155 val = 0xA0000000_00000000;
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));
161 val = 0x8E1BC9BF_04000000;
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));
167 val = 0xCCCCCCCC_CCCCCCCD;
170 CheckPow10MantissaTable(val, exp, new Span<ulong>(s_Pow10MantissaTable, 15, 15), nameof(s_Pow10MantissaTable) + " - inv");
172 val = 0xE69594BE_C44DE160;
175 CheckPow10MantissaTable(val, exp, new Span<ulong>(s_Pow10By16MantissaTable, 21, 21), nameof(s_Pow10By16MantissaTable) + " - inv");
178 // debug-only verification of the precomputed tables
179 private static void CheckPow10MantissaTable(ulong val, int exp, Span<ulong> table, string name)
185 for (int i = 0; i < table.Length; i++)
191 Debug.WriteLine(name);
194 Debug.WriteLine($"0x{val:X16}, // {i + 1}");
198 val = Mul64Precise(val, multval, ref exp);
201 Debug.Assert(!isBad, "NumberToDouble table not correct. Correct version dumped to debug output.");
204 // debug-only verification of the precomputed tables
205 private static void CheckPow10ExponentTable(ulong val, int exp, Span<short> table, string name)
211 for (int i = 0; i < table.Length; i++)
217 Debug.WriteLine(name);
220 Debug.WriteLine($"{val}, // {i + 1}");
224 val = Mul64Precise(val, multval, ref exp);
227 Debug.Assert(!isBad, "NumberToDouble table not correct. Correct version dumped to debug output.");
230 // slower high precision version of Mul64 for computation of the tables
231 private static ulong Mul64Precise(ulong a, ulong b, ref int exp)
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;
237 ulong val = Mul32x32To64((uint)(a >> 32), (uint)(b >> 32))
242 if ((val & 0x80000000_00000000) == 0)
252 // get 32-bit integer from at most 9 digits
253 private static uint DigitsToInt(char* p, int count)
255 Debug.Assert((1 <= count) && (count <= 9));
257 char* end = (p + count);
258 uint res = (uint)(p[0] - '0');
260 for (p++; p < end; p++)
262 res = (10 * res) + p[0] - '0';
268 private static int GetLength(char* src)
272 while (src[length] != '\0')
280 // helper function to multiply two 32-bit uints
281 private static ulong Mul32x32To64(uint a, uint b)
283 return (ulong)(a) * b;
286 // multiply two numbers in the internal integer representation
287 private static ulong Mul64Lossy(ulong a, ulong b, ref int exp)
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);
297 if ((val & 0x80000000_00000000) == 0)
306 private static double NumberToDouble(ref NumberBuffer number)
310 if (!s_CheckedTables)
313 s_CheckedTables = true;
317 char* src = number.GetDigitsPointer();
318 int total = GetLength(src);
319 int remaining = total;
321 // skip the leading zeros
322 while (src[0] == '0')
330 return number.sign ? -0.0 : 0.0;
333 int count = Math.Min(remaining, 9);
335 ulong val = DigitsToInt(src, count);
339 count = Math.Min(remaining, 9);
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);
347 int scale = number.scale - (total - remaining);
348 int absscale = Math.Abs(scale);
349 if (absscale >= 22 * 16)
351 // overflow / underflow
354 return number.sign ? double.NegativeInfinity : double.PositiveInfinity;
358 return number.sign ? -0.0 : 0.0;
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; }
372 int index = absscale & 15;
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;
379 ulong multval = s_Pow10MantissaTable[index + ((scale < 0) ? 15 : 0) - 1];
380 val = Mul64Lossy(val, multval, ref exp);
383 index = absscale >> 4;
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;
390 ulong multval = s_Pow10By16MantissaTable[index + ((scale < 0) ? 21 : 0) - 1];
391 val = Mul64Lossy(val, multval, ref exp);
394 // round & scale down
395 if (((uint)(val) & (1 << 10)) != 0)
397 // IEEE round to even
398 ulong tmp = val + ((1 << 10) - 1) + (((uint)(val) >> 11) & 1);
402 tmp = (tmp >> 1) | 0x8000000000000000;
408 // return the exponent to a biased state
411 // handle overflow, underflow, "Epsilon - 1/2 Epsilon", denormalized, and the normal case
414 if (exp == -52 && (val >= 0x8000000000000058))
416 // round X where {Epsilon > X >= 2.470328229206232730000000E-324} up to Epsilon (instead of down to zero)
417 val = 0x0000000000000001;
427 val >>= (-exp + 11 + 1);
430 else if (exp >= 0x7FF)
433 val = 0x7FF0000000000000;
437 // normal postive exponent case
438 val = ((ulong)(exp) << 52) + ((val >> 11) & 0x000FFFFFFFFFFFFF);
443 val |= 0x8000000000000000;
446 return BitConverter.Int64BitsToDouble((long)(val));