FastDtoa.cs 35 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769
  1. #nullable disable
  2. // Copyright 2010 the V8 project authors. All rights reserved.
  3. // Redistribution and use in source and binary forms, with or without
  4. // modification, are permitted provided that the following conditions are
  5. // met:
  6. //
  7. // * Redistributions of source code must retain the above copyright
  8. // notice, this list of conditions and the following disclaimer.
  9. // * Redistributions in binary form must reproduce the above
  10. // copyright notice, this list of conditions and the following
  11. // disclaimer in the documentation and/or other materials provided
  12. // with the distribution.
  13. // * Neither the name of Google Inc. nor the names of its
  14. // contributors may be used to endorse or promote products derived
  15. // from this software without specific prior written permission.
  16. //
  17. // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  18. // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  19. // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  20. // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
  21. // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  22. // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
  23. // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  24. // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  25. // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  26. // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  27. // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  28. // Ported to Java from Mozilla's version of V8-dtoa by Hannes Wallnoefer.
  29. // The original revision was 67d1049b0bf9 from the mozilla-central tree.
  30. using System.Diagnostics;
  31. using Jint.Runtime;
  32. namespace Jint.Native.Number.Dtoa;
  33. internal sealed class FastDtoa
  34. {
  35. // FastDtoa will produce at most kFastDtoaMaximalLength digits.
  36. public const int KFastDtoaMaximalLength = 17;
  37. // The minimal and maximal target exponent define the range of w's binary
  38. // exponent, where 'w' is the result of multiplying the input by a cached power
  39. // of ten.
  40. //
  41. // A different range might be chosen on a different platform, to optimize digit
  42. // generation, but a smaller range requires more powers of ten to be cached.
  43. private const int MinimalTargetExponent = -60;
  44. private const int MaximalTargetExponent = -32;
  45. // Adjusts the last digit of the generated number, and screens out generated
  46. // solutions that may be inaccurate. A solution may be inaccurate if it is
  47. // outside the safe interval, or if we ctannot prove that it is closer to the
  48. // input than a neighboring representation of the same length.
  49. //
  50. // Input: * buffer containing the digits of too_high / 10^kappa
  51. // * distance_too_high_w == (too_high - w).f() * unit
  52. // * unsafe_interval == (too_high - too_low).f() * unit
  53. // * rest = (too_high - buffer * 10^kappa).f() * unit
  54. // * ten_kappa = 10^kappa * unit
  55. // * unit = the common multiplier
  56. // Output: returns true if the buffer is guaranteed to contain the closest
  57. // representable number to the input.
  58. // Modifies the generated digits in the buffer to approach (round towards) w.
  59. private static bool RoundWeed(
  60. ref DtoaBuilder buffer,
  61. ulong distanceTooHighW,
  62. ulong unsafeInterval,
  63. ulong rest,
  64. ulong tenKappa,
  65. ulong unit)
  66. {
  67. ulong smallDistance = distanceTooHighW - unit;
  68. ulong bigDistance = distanceTooHighW + unit;
  69. // Let w_low = too_high - big_distance, and
  70. // w_high = too_high - small_distance.
  71. // Note: w_low < w < w_high
  72. //
  73. // The real w (* unit) must lie somewhere inside the interval
  74. // ]w_low; w_low[ (often written as "(w_low; w_low)")
  75. // Basically the buffer currently contains a number in the unsafe interval
  76. // ]too_low; too_high[ with too_low < w < too_high
  77. //
  78. // too_high - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  79. // ^v 1 unit ^ ^ ^ ^
  80. // boundary_high --------------------- . . . .
  81. // ^v 1 unit . . . .
  82. // - - - - - - - - - - - - - - - - - - - + - - + - - - - - - . .
  83. // . . ^ . .
  84. // . big_distance . . .
  85. // . . . . rest
  86. // small_distance . . . .
  87. // v . . . .
  88. // w_high - - - - - - - - - - - - - - - - - - . . . .
  89. // ^v 1 unit . . . .
  90. // w ---------------------------------------- . . . .
  91. // ^v 1 unit v . . .
  92. // w_low - - - - - - - - - - - - - - - - - - - - - . . .
  93. // . . v
  94. // buffer --------------------------------------------------+-------+--------
  95. // . .
  96. // safe_interval .
  97. // v .
  98. // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - .
  99. // ^v 1 unit .
  100. // boundary_low ------------------------- unsafe_interval
  101. // ^v 1 unit v
  102. // too_low - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  103. //
  104. //
  105. // Note that the value of buffer could lie anywhere inside the range too_low
  106. // to too_high.
  107. //
  108. // boundary_low, boundary_high and w are approximations of the real boundaries
  109. // and v (the input number). They are guaranteed to be precise up to one unit.
  110. // In fact the error is guaranteed to be strictly less than one unit.
  111. //
  112. // Anything that lies outside the unsafe interval is guaranteed not to round
  113. // to v when read again.
  114. // Anything that lies inside the safe interval is guaranteed to round to v
  115. // when read again.
  116. // If the number inside the buffer lies inside the unsafe interval but not
  117. // inside the safe interval then we simply do not know and bail out (returning
  118. // false).
  119. //
  120. // Similarly we have to take into account the imprecision of 'w' when rounding
  121. // the buffer. If we have two potential representations we need to make sure
  122. // that the chosen one is closer to w_low and w_high since v can be anywhere
  123. // between them.
  124. //
  125. // By generating the digits of too_high we got the largest (closest to
  126. // too_high) buffer that is still in the unsafe interval. In the case where
  127. // w_high < buffer < too_high we try to decrement the buffer.
  128. // This way the buffer approaches (rounds towards) w.
  129. // There are 3 conditions that stop the decrementation process:
  130. // 1) the buffer is already below w_high
  131. // 2) decrementing the buffer would make it leave the unsafe interval
  132. // 3) decrementing the buffer would yield a number below w_high and farther
  133. // away than the current number. In other words:
  134. // (buffer{-1} < w_high) && w_high - buffer{-1} > buffer - w_high
  135. // Instead of using the buffer directly we use its distance to too_high.
  136. // Conceptually rest ~= too_high - buffer
  137. while (rest < smallDistance && // Negated condition 1
  138. unsafeInterval - rest >= tenKappa && // Negated condition 2
  139. (rest + tenKappa < smallDistance || // buffer{-1} > w_high
  140. smallDistance - rest >= rest + tenKappa - smallDistance))
  141. {
  142. buffer.DecreaseLast();
  143. rest += tenKappa;
  144. }
  145. // We have approached w+ as much as possible. We now test if approaching w-
  146. // would require changing the buffer. If yes, then we have two possible
  147. // representations close to w, but we cannot decide which one is closer.
  148. if (rest < bigDistance &&
  149. unsafeInterval - rest >= tenKappa &&
  150. (rest + tenKappa < bigDistance ||
  151. bigDistance - rest > rest + tenKappa - bigDistance))
  152. {
  153. return false;
  154. }
  155. // Weeding test.
  156. // The safe interval is [too_low + 2 ulp; too_high - 2 ulp]
  157. // Since too_low = too_high - unsafe_interval this is equivalent to
  158. // [too_high - unsafe_interval + 4 ulp; too_high - 2 ulp]
  159. // Conceptually we have: rest ~= too_high - buffer
  160. return (2 * unit <= rest) && (rest <= unsafeInterval - 4 * unit);
  161. }
  162. // Rounds the buffer upwards if the result is closer to v by possibly adding
  163. // 1 to the buffer. If the precision of the calculation is not sufficient to
  164. // round correctly, return false.
  165. // The rounding might shift the whole buffer in which case the kappa is
  166. // adjusted. For example "99", kappa = 3 might become "10", kappa = 4.
  167. //
  168. // If 2*rest > ten_kappa then the buffer needs to be round up.
  169. // rest can have an error of +/- 1 unit. This function accounts for the
  170. // imprecision and returns false, if the rounding direction cannot be
  171. // unambiguously determined.
  172. //
  173. // Precondition: rest < ten_kappa.
  174. static bool RoundWeedCounted(
  175. ref DtoaBuilder buffer,
  176. ulong rest,
  177. ulong ten_kappa,
  178. ulong unit,
  179. ref int kappa)
  180. {
  181. Debug.Assert(rest < ten_kappa);
  182. // The following tests are done in a specific order to avoid overflows. They
  183. // will work correctly with any uint64 values of rest < ten_kappa and unit.
  184. //
  185. // If the unit is too big, then we don't know which way to round. For example
  186. // a unit of 50 means that the real number lies within rest +/- 50. If
  187. // 10^kappa == 40 then there is no way to tell which way to round.
  188. if (unit >= ten_kappa) return false;
  189. // Even if unit is just half the size of 10^kappa we are already completely
  190. // lost. (And after the previous test we know that the expression will not
  191. // over/underflow.)
  192. if (ten_kappa - unit <= unit) return false;
  193. // If 2 * (rest + unit) <= 10^kappa we can safely round down.
  194. if ((ten_kappa - rest > rest) && (ten_kappa - 2 * rest >= 2 * unit))
  195. {
  196. return true;
  197. }
  198. // If 2 * (rest - unit) >= 10^kappa, then we can safely round up.
  199. if ((rest > unit) && (ten_kappa - (rest - unit) <= (rest - unit)))
  200. {
  201. // Increment the last digit recursively until we find a non '9' digit.
  202. buffer._chars[buffer.Length - 1]++;
  203. for (int i = buffer.Length - 1; i > 0; --i)
  204. {
  205. if (buffer._chars[i] != '0' + 10) break;
  206. buffer._chars[i] = '0';
  207. buffer._chars[i - 1]++;
  208. }
  209. // If the first digit is now '0'+ 10 we had a buffer with all '9's. With the
  210. // exception of the first digit all digits are now '0'. Simply switch the
  211. // first digit to '1' and adjust the kappa. Example: "99" becomes "10" and
  212. // the power (the kappa) is increased.
  213. if (buffer._chars[0] == '0' + 10)
  214. {
  215. buffer._chars[0] = '1';
  216. kappa += 1;
  217. }
  218. return true;
  219. }
  220. return false;
  221. }
  222. private const int KTen4 = 10000;
  223. private const int KTen5 = 100000;
  224. private const int KTen6 = 1000000;
  225. private const int KTen7 = 10000000;
  226. private const int KTen8 = 100000000;
  227. private const int KTen9 = 1000000000;
  228. // Returns the biggest power of ten that is less than or equal than the given
  229. // number. We furthermore receive the maximum number of bits 'number' has.
  230. // If number_bits == 0 then 0^-1 is returned
  231. // The number of bits must be <= 32.
  232. // Precondition: (1 << number_bits) <= number < (1 << (number_bits + 1)).
  233. private static void BiggestPowerTen(uint number, int numberBits, out uint power, out int exponent)
  234. {
  235. switch (numberBits)
  236. {
  237. case 32:
  238. case 31:
  239. case 30:
  240. if (KTen9 <= number)
  241. {
  242. power = KTen9;
  243. exponent = 9;
  244. break;
  245. } // else fallthrough
  246. goto case 29;
  247. case 29:
  248. case 28:
  249. case 27:
  250. if (KTen8 <= number)
  251. {
  252. power = KTen8;
  253. exponent = 8;
  254. break;
  255. } // else fallthrough
  256. goto case 26;
  257. case 26:
  258. case 25:
  259. case 24:
  260. if (KTen7 <= number)
  261. {
  262. power = KTen7;
  263. exponent = 7;
  264. break;
  265. } // else fallthrough
  266. goto case 23;
  267. case 23:
  268. case 22:
  269. case 21:
  270. case 20:
  271. if (KTen6 <= number)
  272. {
  273. power = KTen6;
  274. exponent = 6;
  275. break;
  276. } // else fallthrough
  277. goto case 19;
  278. case 19:
  279. case 18:
  280. case 17:
  281. if (KTen5 <= number)
  282. {
  283. power = KTen5;
  284. exponent = 5;
  285. break;
  286. } // else fallthrough
  287. goto case 16;
  288. case 16:
  289. case 15:
  290. case 14:
  291. if (KTen4 <= number)
  292. {
  293. power = KTen4;
  294. exponent = 4;
  295. break;
  296. } // else fallthrough
  297. goto case 13;
  298. case 13:
  299. case 12:
  300. case 11:
  301. case 10:
  302. if (1000 <= number)
  303. {
  304. power = 1000;
  305. exponent = 3;
  306. break;
  307. } // else fallthrough
  308. goto case 9;
  309. case 9:
  310. case 8:
  311. case 7:
  312. if (100 <= number)
  313. {
  314. power = 100;
  315. exponent = 2;
  316. break;
  317. } // else fallthrough
  318. goto case 6;
  319. case 6:
  320. case 5:
  321. case 4:
  322. if (10 <= number)
  323. {
  324. power = 10;
  325. exponent = 1;
  326. break;
  327. } // else fallthrough
  328. goto case 3;
  329. case 3:
  330. case 2:
  331. case 1:
  332. if (1 <= number)
  333. {
  334. power = 1;
  335. exponent = 0;
  336. break;
  337. } // else fallthrough
  338. goto case 0;
  339. case 0:
  340. power = 0;
  341. exponent = -1;
  342. break;
  343. default:
  344. // Following assignments are here to silence compiler warnings.
  345. power = 0;
  346. exponent = 0;
  347. // UNREACHABLE();
  348. break;
  349. }
  350. }
  351. // Generates the digits of input number w.
  352. // w is a floating-point number (DiyFp), consisting of a significand and an
  353. // exponent. Its exponent is bounded by minimal_target_exponent and
  354. // maximal_target_exponent.
  355. // Hence -60 <= w.e() <= -32.
  356. //
  357. // Returns false if it fails, in which case the generated digits in the buffer
  358. // should not be used.
  359. // Preconditions:
  360. // * low, w and high are correct up to 1 ulp (unit in the last place). That
  361. // is, their error must be less that a unit of their last digits.
  362. // * low.e() == w.e() == high.e()
  363. // * low < w < high, and taking into account their error: low~ <= high~
  364. // * minimal_target_exponent <= w.e() <= maximal_target_exponent
  365. // Postconditions: returns false if procedure fails.
  366. // otherwise:
  367. // * buffer is not null-terminated, but len contains the number of digits.
  368. // * buffer contains the shortest possible decimal digit-sequence
  369. // such that LOW < buffer * 10^kappa < HIGH, where LOW and HIGH are the
  370. // correct values of low and high (without their error).
  371. // * if more than one decimal representation gives the minimal number of
  372. // decimal digits then the one closest to W (where W is the correct value
  373. // of w) is chosen.
  374. // Remark: this procedure takes into account the imprecision of its input
  375. // numbers. If the precision is not enough to guarantee all the postconditions
  376. // then false is returned. This usually happens rarely (~0.5%).
  377. //
  378. // Say, for the sake of example, that
  379. // w.e() == -48, and w.f() == 0x1234567890abcdef
  380. // w's value can be computed by w.f() * 2^w.e()
  381. // We can obtain w's integral digits by simply shifting w.f() by -w.e().
  382. // -> w's integral part is 0x1234
  383. // w's fractional part is therefore 0x567890abcdef.
  384. // Printing w's integral part is easy (simply print 0x1234 in decimal).
  385. // In order to print its fraction we repeatedly multiply the fraction by 10 and
  386. // get each digit. Example the first digit after the point would be computed by
  387. // (0x567890abcdef * 10) >> 48. -> 3
  388. // The whole thing becomes slightly more complicated because we want to stop
  389. // once we have enough digits. That is, once the digits inside the buffer
  390. // represent 'w' we can stop. Everything inside the interval low - high
  391. // represents w. However we have to pay attention to low, high and w's
  392. // imprecision.
  393. private static bool DigitGen(
  394. in DiyFp low,
  395. in DiyFp w,
  396. in DiyFp high,
  397. ref DtoaBuilder buffer,
  398. int mk,
  399. out int kappa)
  400. {
  401. // low, w and high are imprecise, but by less than one ulp (unit in the last
  402. // place).
  403. // If we remove (resp. add) 1 ulp from low (resp. high) we are certain that
  404. // the new numbers are outside of the interval we want the final
  405. // representation to lie in.
  406. // Inversely adding (resp. removing) 1 ulp from low (resp. high) would yield
  407. // numbers that are certain to lie in the interval. We will use this fact
  408. // later on.
  409. // We will now start by generating the digits within the uncertain
  410. // interval. Later we will weed out representations that lie outside the safe
  411. // interval and thus _might_ lie outside the correct interval.
  412. ulong unit = 1;
  413. var tooLow = new DiyFp(low.F - unit, low.E);
  414. var tooHigh = new DiyFp(high.F + unit, high.E);
  415. // too_low and too_high are guaranteed to lie outside the interval we want the
  416. // generated number in.
  417. var unsafeInterval = DiyFp.Minus(tooHigh, tooLow);
  418. // We now cut the input number into two parts: the integral digits and the
  419. // fractionals. We will not write any decimal separator though, but adapt
  420. // kappa instead.
  421. // Reminder: we are currently computing the digits (stored inside the buffer)
  422. // such that: too_low < buffer * 10^kappa < too_high
  423. // We use too_high for the digit_generation and stop as soon as possible.
  424. // If we stop early we effectively round down.
  425. var one = new DiyFp(((ulong) 1) << -w.E, w.E);
  426. // Division by one is a shift.
  427. var integrals = (uint) (tooHigh.F.UnsignedShift(-one.E) & 0xffffffffL);
  428. // Modulo by one is an and.
  429. ulong fractionals = tooHigh.F & (one.F - 1);
  430. BiggestPowerTen(
  431. integrals,
  432. DiyFp.KSignificandSize - (-one.E),
  433. out var divider,
  434. out var dividerExponent);
  435. kappa = dividerExponent + 1;
  436. // Loop invariant: buffer = too_high / 10^kappa (integer division)
  437. // The invariant holds for the first iteration: kappa has been initialized
  438. // with the divider exponent + 1. And the divider is the biggest power of ten
  439. // that is smaller than integrals.
  440. while (kappa > 0)
  441. {
  442. int digit = (int) (integrals / divider);
  443. buffer.Append((char) ('0' + digit));
  444. integrals %= divider;
  445. kappa--;
  446. // Note that kappa now equals the exponent of the divider and that the
  447. // invariant thus holds again.
  448. ulong rest = ((ulong) integrals << -one.E) + fractionals;
  449. // Invariant: too_high = buffer * 10^kappa + DiyFp(rest, one.e())
  450. // Reminder: unsafe_interval.e() == one.e()
  451. if (rest < unsafeInterval.F)
  452. {
  453. // Rounding down (by not emitting the remaining digits) yields a number
  454. // that lies within the unsafe interval.
  455. return RoundWeed(
  456. ref buffer,
  457. DiyFp.Minus(tooHigh, w).F,
  458. unsafeInterval.F,
  459. rest,
  460. (ulong) divider << -one.E,
  461. unit);
  462. }
  463. divider /= 10;
  464. }
  465. // The integrals have been generated. We are at the point of the decimal
  466. // separator. In the following loop we simply multiply the remaining digits by
  467. // 10 and divide by one. We just need to pay attention to multiply associated
  468. // data (like the interval or 'unit'), too.
  469. // Instead of multiplying by 10 we multiply by 5 (cheaper operation) and
  470. // increase its (imaginary) exponent. At the same time we decrease the
  471. // divider's (one's) exponent and shift its significand.
  472. // Basically, if fractionals was a DiyFp (with fractionals.e == one.e):
  473. // fractionals.f *= 10;
  474. // fractionals.f >>= 1; fractionals.e++; // value remains unchanged.
  475. // one.f >>= 1; one.e++; // value remains unchanged.
  476. // and we have again fractionals.e == one.e which allows us to divide
  477. // fractionals.f() by one.f()
  478. // We simply combine the *= 10 and the >>= 1.
  479. while (true)
  480. {
  481. fractionals *= 5;
  482. unit *= 5;
  483. unsafeInterval = new DiyFp(unsafeInterval.F * 5, unsafeInterval.E + 1); // Will be optimized out.
  484. one = new DiyFp(one.F.UnsignedShift(1), one.E + 1);
  485. // Integer division by one.
  486. var digit = (int) ((fractionals.UnsignedShift(-one.E)) & 0xffffffffL);
  487. buffer.Append((char) ('0' + digit));
  488. fractionals &= one.F - 1; // Modulo by one.
  489. kappa--;
  490. if (fractionals < unsafeInterval.F)
  491. {
  492. return RoundWeed(
  493. ref buffer,
  494. DiyFp.Minus(tooHigh, w).F * unit,
  495. unsafeInterval.F,
  496. fractionals,
  497. one.F,
  498. unit);
  499. }
  500. }
  501. }
  502. // Generates (at most) requested_digits of input number w.
  503. // w is a floating-point number (DiyFp), consisting of a significand and an
  504. // exponent. Its exponent is bounded by kMinimalTargetExponent and
  505. // kMaximalTargetExponent.
  506. // Hence -60 <= w.e() <= -32.
  507. //
  508. // Returns false if it fails, in which case the generated digits in the buffer
  509. // should not be used.
  510. // Preconditions:
  511. // * w is correct up to 1 ulp (unit in the last place). That
  512. // is, its error must be strictly less than a unit of its last digit.
  513. // * kMinimalTargetExponent <= w.e() <= kMaximalTargetExponent
  514. //
  515. // Postconditions: returns false if procedure fails.
  516. // otherwise:
  517. // * buffer is not null-terminated, but length contains the number of
  518. // digits.
  519. // * the representation in buffer is the most precise representation of
  520. // requested_digits digits.
  521. // * buffer contains at most requested_digits digits of w. If there are less
  522. // than requested_digits digits then some trailing '0's have been removed.
  523. // * kappa is such that
  524. // w = buffer * 10^kappa + eps with |eps| < 10^kappa / 2.
  525. //
  526. // Remark: This procedure takes into account the imprecision of its input
  527. // numbers. If the precision is not enough to guarantee all the postconditions
  528. // then false is returned. This usually happens rarely, but the failure-rate
  529. // increases with higher requested_digits.
  530. static bool DigitGenCounted(
  531. in DiyFp w,
  532. int requested_digits,
  533. ref DtoaBuilder buffer,
  534. out int kappa)
  535. {
  536. Debug.Assert(MinimalTargetExponent <= w.E && w.E <= MaximalTargetExponent);
  537. // w is assumed to have an error less than 1 unit. Whenever w is scaled we
  538. // also scale its error.
  539. ulong w_error = 1;
  540. // We cut the input number into two parts: the integral digits and the
  541. // fractional digits. We don't emit any decimal separator, but adapt kappa
  542. // instead. Example: instead of writing "1.2" we put "12" into the buffer and
  543. // increase kappa by 1.
  544. DiyFp one = new DiyFp(((ulong) 1) << -w.E, w.E);
  545. // Division by one is a shift.
  546. uint integrals = (uint) (w.F >> -one.E);
  547. // Modulo by one is an and.
  548. ulong fractionals = w.F & (one.F - 1);
  549. BiggestPowerTen(integrals, DiyFp.KSignificandSize - (-one.E), out var divisor, out var divisor_exponent);
  550. kappa = divisor_exponent + 1;
  551. // Loop invariant: buffer = w / 10^kappa (integer division)
  552. // The invariant holds for the first iteration: kappa has been initialized
  553. // with the divisor exponent + 1. And the divisor is the biggest power of ten
  554. // that is smaller than 'integrals'.
  555. while (kappa > 0)
  556. {
  557. int digit = (int) (integrals / divisor);
  558. buffer.Append((char) ('0' + digit));
  559. requested_digits--;
  560. integrals %= divisor;
  561. kappa--;
  562. // Note that kappa now equals the exponent of the divisor and that the
  563. // invariant thus holds again.
  564. if (requested_digits == 0) break;
  565. divisor /= 10;
  566. }
  567. if (requested_digits == 0)
  568. {
  569. ulong rest = (((ulong) integrals) << -one.E) + fractionals;
  570. return RoundWeedCounted(ref buffer, rest, (ulong) divisor << -one.E, w_error, ref kappa);
  571. }
  572. // The integrals have been generated. We are at the point of the decimal
  573. // separator. In the following loop we simply multiply the remaining digits by
  574. // 10 and divide by one. We just need to pay attention to multiply associated
  575. // data (the 'unit'), too.
  576. // Note that the multiplication by 10 does not overflow, because w.e >= -60
  577. // and thus one.e >= -60.
  578. Debug.Assert(one.E >= -60);
  579. Debug.Assert(fractionals < one.F);
  580. while (requested_digits > 0 && fractionals > w_error)
  581. {
  582. fractionals *= 10;
  583. w_error *= 10;
  584. // Integer division by one.
  585. int digit = (int) (fractionals >> -one.E);
  586. buffer.Append((char) ('0' + digit));
  587. requested_digits--;
  588. fractionals &= one.F - 1; // Modulo by one.
  589. (kappa)--;
  590. }
  591. if (requested_digits != 0) return false;
  592. return RoundWeedCounted(ref buffer, fractionals, one.F, w_error, ref kappa);
  593. }
  594. // Provides a decimal representation of v.
  595. // Returns true if it succeeds, otherwise the result cannot be trusted.
  596. // There will be *length digits inside the buffer (not null-terminated).
  597. // If the function returns true then
  598. // v == (double) (buffer * 10^decimal_exponent).
  599. // The digits in the buffer are the shortest representation possible: no
  600. // 0.09999999999999999 instead of 0.1. The shorter representation will even be
  601. // chosen even if the longer one would be closer to v.
  602. // The last digit will be closest to the actual v. That is, even if several
  603. // digits might correctly yield 'v' when read again, the closest will be
  604. // computed.
  605. private static bool Grisu3(double v, ref DtoaBuilder buffer, out int decimal_exponent)
  606. {
  607. ulong bits = (ulong) BitConverter.DoubleToInt64Bits(v);
  608. DiyFp w = DoubleHelper.AsNormalizedDiyFp(bits);
  609. // boundary_minus and boundary_plus are the boundaries between v and its
  610. // closest floating-point neighbors. Any number strictly between
  611. // boundary_minus and boundary_plus will round to v when convert to a double.
  612. // Grisu3 will never output representations that lie exactly on a boundary.
  613. var boundaries = DoubleHelper.NormalizedBoundaries(bits);
  614. var boundaryMinus = boundaries.Minus;
  615. var boundaryPlus = boundaries.Plus;
  616. Debug.Assert(boundaryPlus.E == w.E);
  617. var result = CachedPowers.GetCachedPowerForBinaryExponentRange(
  618. MinimalTargetExponent - (w.E + DiyFp.KSignificandSize),
  619. MaximalTargetExponent - (w.E + DiyFp.KSignificandSize));
  620. var mk = result.decimalExponent;
  621. var tenMk = result.cMk;
  622. Debug.Assert(MinimalTargetExponent <= w.E + tenMk.E +
  623. DiyFp.KSignificandSize &&
  624. MaximalTargetExponent >= w.E + tenMk.E +
  625. DiyFp.KSignificandSize);
  626. // Note that ten_mk is only an approximation of 10^-k. A DiyFp only contains a
  627. // 64 bit significand and ten_mk is thus only precise up to 64 bits.
  628. // The DiyFp::Times procedure rounds its result, and ten_mk is approximated
  629. // too. The variable scaled_w (as well as scaled_boundary_minus/plus) are now
  630. // off by a small amount.
  631. // In fact: scaled_w - w*10^k < 1ulp (unit in the last place) of scaled_w.
  632. // In other words: let f = scaled_w.f() and e = scaled_w.e(), then
  633. // (f-1) * 2^e < w*10^k < (f+1) * 2^e
  634. DiyFp scaledW = DiyFp.Times(w, tenMk);
  635. Debug.Assert(scaledW.E ==
  636. boundaryPlus.E + tenMk.E + DiyFp.KSignificandSize);
  637. // In theory it would be possible to avoid some recomputations by computing
  638. // the difference between w and boundary_minus/plus (a power of 2) and to
  639. // compute scaled_boundary_minus/plus by subtracting/adding from
  640. // scaled_w. However the code becomes much less readable and the speed
  641. // enhancements are not terriffic.
  642. DiyFp scaledBoundaryMinus = DiyFp.Times(boundaryMinus, tenMk);
  643. DiyFp scaledBoundaryPlus = DiyFp.Times(boundaryPlus, tenMk);
  644. // DigitGen will generate the digits of scaled_w. Therefore we have
  645. // v == (double) (scaled_w * 10^-mk).
  646. // Set decimal_exponent == -mk and pass it to DigitGen. If scaled_w is not an
  647. // integer than it will be updated. For instance if scaled_w == 1.23 then
  648. // the buffer will be filled with "123" und the decimal_exponent will be
  649. // decreased by 2.
  650. int kappa;
  651. var digitGen = DigitGen(scaledBoundaryMinus, scaledW, scaledBoundaryPlus, ref buffer, mk, out kappa);
  652. decimal_exponent = -mk + kappa;
  653. return digitGen;
  654. }
  655. // The "counted" version of grisu3 (see above) only generates requested_digits
  656. // number of digits. This version does not generate the shortest representation,
  657. // and with enough requested digits 0.1 will at some point print as 0.9999999...
  658. // Grisu3 is too imprecise for real halfway cases (1.5 will not work) and
  659. // therefore the rounding strategy for halfway cases is irrelevant.
  660. static bool Grisu3Counted(
  661. double v,
  662. int requested_digits,
  663. ref DtoaBuilder buffer,
  664. out int decimal_exponent)
  665. {
  666. ulong bits = (ulong) BitConverter.DoubleToInt64Bits(v);
  667. DiyFp w = DoubleHelper.AsNormalizedDiyFp(bits);
  668. var powerResult = CachedPowers.GetCachedPowerForBinaryExponentRange(
  669. MinimalTargetExponent - (w.E + DiyFp.KSignificandSize),
  670. MaximalTargetExponent - (w.E + DiyFp.KSignificandSize));
  671. var mk = powerResult.decimalExponent;
  672. var ten_mk = powerResult.cMk;
  673. Debug.Assert((MinimalTargetExponent <= w.E + ten_mk.E + DiyFp.KSignificandSize) && (MaximalTargetExponent >= w.E + ten_mk.E + DiyFp.KSignificandSize));
  674. // Note that ten_mk is only an approximation of 10^-k. A DiyFp only contains a
  675. // 64 bit significand and ten_mk is thus only precise up to 64 bits.
  676. // The DiyFp::Times procedure rounds its result, and ten_mk is approximated
  677. // too. The variable scaled_w (as well as scaled_boundary_minus/plus) are now
  678. // off by a small amount.
  679. // In fact: scaled_w - w*10^k < 1ulp (unit in the last place) of scaled_w.
  680. // In other words: let f = scaled_w.f() and e = scaled_w.e(), then
  681. // (f-1) * 2^e < w*10^k < (f+1) * 2^e
  682. DiyFp scaled_w = DiyFp.Times(w, ten_mk);
  683. // We now have (double) (scaled_w * 10^-mk).
  684. // DigitGen will generate the first requested_digits digits of scaled_w and
  685. // return together with a kappa such that scaled_w ~= buffer * 10^kappa. (It
  686. // will not always be exactly the same since DigitGenCounted only produces a
  687. // limited number of digits.)
  688. bool result = DigitGenCounted(scaled_w, requested_digits, ref buffer, out var kappa);
  689. decimal_exponent = -mk + kappa;
  690. return result;
  691. }
  692. public static bool NumberToString(
  693. double v,
  694. DtoaMode mode,
  695. int requested_digits,
  696. out int decimal_point,
  697. ref DtoaBuilder buffer)
  698. {
  699. Debug.Assert(v > 0);
  700. Debug.Assert(!double.IsNaN(v));
  701. Debug.Assert(!double.IsInfinity(v));
  702. var result = false;
  703. var decimal_exponent = 0;
  704. switch (mode)
  705. {
  706. case DtoaMode.Shortest:
  707. result = Grisu3(v, ref buffer, out decimal_exponent);
  708. break;
  709. case DtoaMode.Precision:
  710. result = Grisu3Counted(v, requested_digits, ref buffer, out decimal_exponent);
  711. break;
  712. default:
  713. ExceptionHelper.ThrowArgumentOutOfRangeException();
  714. break;
  715. }
  716. if (result)
  717. {
  718. decimal_point = buffer.Length + decimal_exponent;
  719. return true;
  720. }
  721. decimal_point = -1;
  722. return false;
  723. }
  724. }