FastDtoa.cs 37 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. {
  34. internal sealed class FastDtoa
  35. {
  36. // FastDtoa will produce at most kFastDtoaMaximalLength digits.
  37. public const int KFastDtoaMaximalLength = 17;
  38. // The minimal and maximal target exponent define the range of w's binary
  39. // exponent, where 'w' is the result of multiplying the input by a cached power
  40. // of ten.
  41. //
  42. // A different range might be chosen on a different platform, to optimize digit
  43. // generation, but a smaller range requires more powers of ten to be cached.
  44. private const int MinimalTargetExponent = -60;
  45. private const int MaximalTargetExponent = -32;
  46. // Adjusts the last digit of the generated number, and screens out generated
  47. // solutions that may be inaccurate. A solution may be inaccurate if it is
  48. // outside the safe interval, or if we ctannot prove that it is closer to the
  49. // input than a neighboring representation of the same length.
  50. //
  51. // Input: * buffer containing the digits of too_high / 10^kappa
  52. // * distance_too_high_w == (too_high - w).f() * unit
  53. // * unsafe_interval == (too_high - too_low).f() * unit
  54. // * rest = (too_high - buffer * 10^kappa).f() * unit
  55. // * ten_kappa = 10^kappa * unit
  56. // * unit = the common multiplier
  57. // Output: returns true if the buffer is guaranteed to contain the closest
  58. // representable number to the input.
  59. // Modifies the generated digits in the buffer to approach (round towards) w.
  60. private static bool RoundWeed(
  61. ref DtoaBuilder buffer,
  62. ulong distanceTooHighW,
  63. ulong unsafeInterval,
  64. ulong rest,
  65. ulong tenKappa,
  66. ulong unit)
  67. {
  68. ulong smallDistance = distanceTooHighW - unit;
  69. ulong bigDistance = distanceTooHighW + unit;
  70. // Let w_low = too_high - big_distance, and
  71. // w_high = too_high - small_distance.
  72. // Note: w_low < w < w_high
  73. //
  74. // The real w (* unit) must lie somewhere inside the interval
  75. // ]w_low; w_low[ (often written as "(w_low; w_low)")
  76. // Basically the buffer currently contains a number in the unsafe interval
  77. // ]too_low; too_high[ with too_low < w < too_high
  78. //
  79. // too_high - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  80. // ^v 1 unit ^ ^ ^ ^
  81. // boundary_high --------------------- . . . .
  82. // ^v 1 unit . . . .
  83. // - - - - - - - - - - - - - - - - - - - + - - + - - - - - - . .
  84. // . . ^ . .
  85. // . big_distance . . .
  86. // . . . . rest
  87. // small_distance . . . .
  88. // v . . . .
  89. // w_high - - - - - - - - - - - - - - - - - - . . . .
  90. // ^v 1 unit . . . .
  91. // w ---------------------------------------- . . . .
  92. // ^v 1 unit v . . .
  93. // w_low - - - - - - - - - - - - - - - - - - - - - . . .
  94. // . . v
  95. // buffer --------------------------------------------------+-------+--------
  96. // . .
  97. // safe_interval .
  98. // v .
  99. // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - .
  100. // ^v 1 unit .
  101. // boundary_low ------------------------- unsafe_interval
  102. // ^v 1 unit v
  103. // too_low - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  104. //
  105. //
  106. // Note that the value of buffer could lie anywhere inside the range too_low
  107. // to too_high.
  108. //
  109. // boundary_low, boundary_high and w are approximations of the real boundaries
  110. // and v (the input number). They are guaranteed to be precise up to one unit.
  111. // In fact the error is guaranteed to be strictly less than one unit.
  112. //
  113. // Anything that lies outside the unsafe interval is guaranteed not to round
  114. // to v when read again.
  115. // Anything that lies inside the safe interval is guaranteed to round to v
  116. // when read again.
  117. // If the number inside the buffer lies inside the unsafe interval but not
  118. // inside the safe interval then we simply do not know and bail out (returning
  119. // false).
  120. //
  121. // Similarly we have to take into account the imprecision of 'w' when rounding
  122. // the buffer. If we have two potential representations we need to make sure
  123. // that the chosen one is closer to w_low and w_high since v can be anywhere
  124. // between them.
  125. //
  126. // By generating the digits of too_high we got the largest (closest to
  127. // too_high) buffer that is still in the unsafe interval. In the case where
  128. // w_high < buffer < too_high we try to decrement the buffer.
  129. // This way the buffer approaches (rounds towards) w.
  130. // There are 3 conditions that stop the decrementation process:
  131. // 1) the buffer is already below w_high
  132. // 2) decrementing the buffer would make it leave the unsafe interval
  133. // 3) decrementing the buffer would yield a number below w_high and farther
  134. // away than the current number. In other words:
  135. // (buffer{-1} < w_high) && w_high - buffer{-1} > buffer - w_high
  136. // Instead of using the buffer directly we use its distance to too_high.
  137. // Conceptually rest ~= too_high - buffer
  138. while (rest < smallDistance && // Negated condition 1
  139. unsafeInterval - rest >= tenKappa && // Negated condition 2
  140. (rest + tenKappa < smallDistance || // buffer{-1} > w_high
  141. smallDistance - rest >= rest + tenKappa - smallDistance))
  142. {
  143. buffer.DecreaseLast();
  144. rest += tenKappa;
  145. }
  146. // We have approached w+ as much as possible. We now test if approaching w-
  147. // would require changing the buffer. If yes, then we have two possible
  148. // representations close to w, but we cannot decide which one is closer.
  149. if (rest < bigDistance &&
  150. unsafeInterval - rest >= tenKappa &&
  151. (rest + tenKappa < bigDistance ||
  152. bigDistance - rest > rest + tenKappa - bigDistance))
  153. {
  154. return false;
  155. }
  156. // Weeding test.
  157. // The safe interval is [too_low + 2 ulp; too_high - 2 ulp]
  158. // Since too_low = too_high - unsafe_interval this is equivalent to
  159. // [too_high - unsafe_interval + 4 ulp; too_high - 2 ulp]
  160. // Conceptually we have: rest ~= too_high - buffer
  161. return (2*unit <= rest) && (rest <= unsafeInterval - 4*unit);
  162. }
  163. // Rounds the buffer upwards if the result is closer to v by possibly adding
  164. // 1 to the buffer. If the precision of the calculation is not sufficient to
  165. // round correctly, return false.
  166. // The rounding might shift the whole buffer in which case the kappa is
  167. // adjusted. For example "99", kappa = 3 might become "10", kappa = 4.
  168. //
  169. // If 2*rest > ten_kappa then the buffer needs to be round up.
  170. // rest can have an error of +/- 1 unit. This function accounts for the
  171. // imprecision and returns false, if the rounding direction cannot be
  172. // unambiguously determined.
  173. //
  174. // Precondition: rest < ten_kappa.
  175. static bool RoundWeedCounted(
  176. ref DtoaBuilder buffer,
  177. ulong rest,
  178. ulong ten_kappa,
  179. ulong unit,
  180. ref int kappa)
  181. {
  182. Debug.Assert(rest < ten_kappa);
  183. // The following tests are done in a specific order to avoid overflows. They
  184. // will work correctly with any uint64 values of rest < ten_kappa and unit.
  185. //
  186. // If the unit is too big, then we don't know which way to round. For example
  187. // a unit of 50 means that the real number lies within rest +/- 50. If
  188. // 10^kappa == 40 then there is no way to tell which way to round.
  189. if (unit >= ten_kappa) return false;
  190. // Even if unit is just half the size of 10^kappa we are already completely
  191. // lost. (And after the previous test we know that the expression will not
  192. // over/underflow.)
  193. if (ten_kappa - unit <= unit) return false;
  194. // If 2 * (rest + unit) <= 10^kappa we can safely round down.
  195. if ((ten_kappa - rest > rest) && (ten_kappa - 2 * rest >= 2 * unit))
  196. {
  197. return true;
  198. }
  199. // If 2 * (rest - unit) >= 10^kappa, then we can safely round up.
  200. if ((rest > unit) && (ten_kappa - (rest - unit) <= (rest - unit)))
  201. {
  202. // Increment the last digit recursively until we find a non '9' digit.
  203. buffer._chars[buffer.Length - 1]++;
  204. for (int i = buffer.Length - 1; i > 0; --i)
  205. {
  206. if (buffer._chars[i] != '0' + 10) break;
  207. buffer._chars[i] = '0';
  208. buffer._chars[i - 1]++;
  209. }
  210. // If the first digit is now '0'+ 10 we had a buffer with all '9's. With the
  211. // exception of the first digit all digits are now '0'. Simply switch the
  212. // first digit to '1' and adjust the kappa. Example: "99" becomes "10" and
  213. // the power (the kappa) is increased.
  214. if (buffer._chars[0] == '0' + 10)
  215. {
  216. buffer._chars[0] = '1';
  217. kappa += 1;
  218. }
  219. return true;
  220. }
  221. return false;
  222. }
  223. private const int KTen4 = 10000;
  224. private const int KTen5 = 100000;
  225. private const int KTen6 = 1000000;
  226. private const int KTen7 = 10000000;
  227. private const int KTen8 = 100000000;
  228. private const int KTen9 = 1000000000;
  229. // Returns the biggest power of ten that is less than or equal than the given
  230. // number. We furthermore receive the maximum number of bits 'number' has.
  231. // If number_bits == 0 then 0^-1 is returned
  232. // The number of bits must be <= 32.
  233. // Precondition: (1 << number_bits) <= number < (1 << (number_bits + 1)).
  234. private static void BiggestPowerTen(uint number, int numberBits, out uint power, out int exponent)
  235. {
  236. switch (numberBits)
  237. {
  238. case 32:
  239. case 31:
  240. case 30:
  241. if (KTen9 <= number)
  242. {
  243. power = KTen9;
  244. exponent = 9;
  245. break;
  246. } // else fallthrough
  247. goto case 29;
  248. case 29:
  249. case 28:
  250. case 27:
  251. if (KTen8 <= number)
  252. {
  253. power = KTen8;
  254. exponent = 8;
  255. break;
  256. } // else fallthrough
  257. goto case 26;
  258. case 26:
  259. case 25:
  260. case 24:
  261. if (KTen7 <= number)
  262. {
  263. power = KTen7;
  264. exponent = 7;
  265. break;
  266. } // else fallthrough
  267. goto case 23;
  268. case 23:
  269. case 22:
  270. case 21:
  271. case 20:
  272. if (KTen6 <= number)
  273. {
  274. power = KTen6;
  275. exponent = 6;
  276. break;
  277. } // else fallthrough
  278. goto case 19;
  279. case 19:
  280. case 18:
  281. case 17:
  282. if (KTen5 <= number)
  283. {
  284. power = KTen5;
  285. exponent = 5;
  286. break;
  287. } // else fallthrough
  288. goto case 16;
  289. case 16:
  290. case 15:
  291. case 14:
  292. if (KTen4 <= number)
  293. {
  294. power = KTen4;
  295. exponent = 4;
  296. break;
  297. } // else fallthrough
  298. goto case 13;
  299. case 13:
  300. case 12:
  301. case 11:
  302. case 10:
  303. if (1000 <= number)
  304. {
  305. power = 1000;
  306. exponent = 3;
  307. break;
  308. } // else fallthrough
  309. goto case 9;
  310. case 9:
  311. case 8:
  312. case 7:
  313. if (100 <= number)
  314. {
  315. power = 100;
  316. exponent = 2;
  317. break;
  318. } // else fallthrough
  319. goto case 6;
  320. case 6:
  321. case 5:
  322. case 4:
  323. if (10 <= number)
  324. {
  325. power = 10;
  326. exponent = 1;
  327. break;
  328. } // else fallthrough
  329. goto case 3;
  330. case 3:
  331. case 2:
  332. case 1:
  333. if (1 <= number)
  334. {
  335. power = 1;
  336. exponent = 0;
  337. break;
  338. } // else fallthrough
  339. goto case 0;
  340. case 0:
  341. power = 0;
  342. exponent = -1;
  343. break;
  344. default:
  345. // Following assignments are here to silence compiler warnings.
  346. power = 0;
  347. exponent = 0;
  348. // UNREACHABLE();
  349. break;
  350. }
  351. }
  352. // Generates the digits of input number w.
  353. // w is a floating-point number (DiyFp), consisting of a significand and an
  354. // exponent. Its exponent is bounded by minimal_target_exponent and
  355. // maximal_target_exponent.
  356. // Hence -60 <= w.e() <= -32.
  357. //
  358. // Returns false if it fails, in which case the generated digits in the buffer
  359. // should not be used.
  360. // Preconditions:
  361. // * low, w and high are correct up to 1 ulp (unit in the last place). That
  362. // is, their error must be less that a unit of their last digits.
  363. // * low.e() == w.e() == high.e()
  364. // * low < w < high, and taking into account their error: low~ <= high~
  365. // * minimal_target_exponent <= w.e() <= maximal_target_exponent
  366. // Postconditions: returns false if procedure fails.
  367. // otherwise:
  368. // * buffer is not null-terminated, but len contains the number of digits.
  369. // * buffer contains the shortest possible decimal digit-sequence
  370. // such that LOW < buffer * 10^kappa < HIGH, where LOW and HIGH are the
  371. // correct values of low and high (without their error).
  372. // * if more than one decimal representation gives the minimal number of
  373. // decimal digits then the one closest to W (where W is the correct value
  374. // of w) is chosen.
  375. // Remark: this procedure takes into account the imprecision of its input
  376. // numbers. If the precision is not enough to guarantee all the postconditions
  377. // then false is returned. This usually happens rarely (~0.5%).
  378. //
  379. // Say, for the sake of example, that
  380. // w.e() == -48, and w.f() == 0x1234567890abcdef
  381. // w's value can be computed by w.f() * 2^w.e()
  382. // We can obtain w's integral digits by simply shifting w.f() by -w.e().
  383. // -> w's integral part is 0x1234
  384. // w's fractional part is therefore 0x567890abcdef.
  385. // Printing w's integral part is easy (simply print 0x1234 in decimal).
  386. // In order to print its fraction we repeatedly multiply the fraction by 10 and
  387. // get each digit. Example the first digit after the point would be computed by
  388. // (0x567890abcdef * 10) >> 48. -> 3
  389. // The whole thing becomes slightly more complicated because we want to stop
  390. // once we have enough digits. That is, once the digits inside the buffer
  391. // represent 'w' we can stop. Everything inside the interval low - high
  392. // represents w. However we have to pay attention to low, high and w's
  393. // imprecision.
  394. private static bool DigitGen(
  395. in DiyFp low,
  396. in DiyFp w,
  397. in DiyFp high,
  398. ref DtoaBuilder buffer,
  399. int mk,
  400. out int kappa)
  401. {
  402. // low, w and high are imprecise, but by less than one ulp (unit in the last
  403. // place).
  404. // If we remove (resp. add) 1 ulp from low (resp. high) we are certain that
  405. // the new numbers are outside of the interval we want the final
  406. // representation to lie in.
  407. // Inversely adding (resp. removing) 1 ulp from low (resp. high) would yield
  408. // numbers that are certain to lie in the interval. We will use this fact
  409. // later on.
  410. // We will now start by generating the digits within the uncertain
  411. // interval. Later we will weed out representations that lie outside the safe
  412. // interval and thus _might_ lie outside the correct interval.
  413. ulong unit = 1;
  414. var tooLow = new DiyFp(low.F - unit, low.E);
  415. var tooHigh = new DiyFp(high.F + unit, high.E);
  416. // too_low and too_high are guaranteed to lie outside the interval we want the
  417. // generated number in.
  418. var unsafeInterval = DiyFp.Minus(tooHigh, tooLow);
  419. // We now cut the input number into two parts: the integral digits and the
  420. // fractionals. We will not write any decimal separator though, but adapt
  421. // kappa instead.
  422. // Reminder: we are currently computing the digits (stored inside the buffer)
  423. // such that: too_low < buffer * 10^kappa < too_high
  424. // We use too_high for the digit_generation and stop as soon as possible.
  425. // If we stop early we effectively round down.
  426. var one = new DiyFp(((ulong) 1) << -w.E, w.E);
  427. // Division by one is a shift.
  428. var integrals = (uint) (tooHigh.F.UnsignedShift(-one.E) & 0xffffffffL);
  429. // Modulo by one is an and.
  430. ulong fractionals = tooHigh.F & (one.F - 1);
  431. BiggestPowerTen(
  432. integrals,
  433. DiyFp.KSignificandSize - (-one.E),
  434. out var divider,
  435. out var dividerExponent);
  436. kappa = dividerExponent + 1;
  437. // Loop invariant: buffer = too_high / 10^kappa (integer division)
  438. // The invariant holds for the first iteration: kappa has been initialized
  439. // with the divider exponent + 1. And the divider is the biggest power of ten
  440. // that is smaller than integrals.
  441. while (kappa > 0)
  442. {
  443. int digit = (int) (integrals/divider);
  444. buffer.Append((char) ('0' + digit));
  445. integrals %= divider;
  446. kappa--;
  447. // Note that kappa now equals the exponent of the divider and that the
  448. // invariant thus holds again.
  449. ulong rest = ((ulong) integrals << -one.E) + fractionals;
  450. // Invariant: too_high = buffer * 10^kappa + DiyFp(rest, one.e())
  451. // Reminder: unsafe_interval.e() == one.e()
  452. if (rest < unsafeInterval.F)
  453. {
  454. // Rounding down (by not emitting the remaining digits) yields a number
  455. // that lies within the unsafe interval.
  456. return RoundWeed(
  457. ref buffer,
  458. DiyFp.Minus(tooHigh, w).F,
  459. unsafeInterval.F,
  460. rest,
  461. (ulong) divider << -one.E,
  462. unit);
  463. }
  464. divider /= 10;
  465. }
  466. // The integrals have been generated. We are at the point of the decimal
  467. // separator. In the following loop we simply multiply the remaining digits by
  468. // 10 and divide by one. We just need to pay attention to multiply associated
  469. // data (like the interval or 'unit'), too.
  470. // Instead of multiplying by 10 we multiply by 5 (cheaper operation) and
  471. // increase its (imaginary) exponent. At the same time we decrease the
  472. // divider's (one's) exponent and shift its significand.
  473. // Basically, if fractionals was a DiyFp (with fractionals.e == one.e):
  474. // fractionals.f *= 10;
  475. // fractionals.f >>= 1; fractionals.e++; // value remains unchanged.
  476. // one.f >>= 1; one.e++; // value remains unchanged.
  477. // and we have again fractionals.e == one.e which allows us to divide
  478. // fractionals.f() by one.f()
  479. // We simply combine the *= 10 and the >>= 1.
  480. while (true)
  481. {
  482. fractionals *= 5;
  483. unit *= 5;
  484. unsafeInterval = new DiyFp(unsafeInterval.F*5, unsafeInterval.E + 1); // Will be optimized out.
  485. one = new DiyFp(one.F.UnsignedShift(1), one.E + 1);
  486. // Integer division by one.
  487. var digit = (int) ((fractionals.UnsignedShift(-one.E)) & 0xffffffffL);
  488. buffer.Append((char) ('0' + digit));
  489. fractionals &= one.F - 1; // Modulo by one.
  490. kappa--;
  491. if (fractionals < unsafeInterval.F)
  492. {
  493. return RoundWeed(
  494. ref buffer,
  495. DiyFp.Minus(tooHigh, w).F*unit,
  496. unsafeInterval.F,
  497. fractionals,
  498. one.F,
  499. unit);
  500. }
  501. }
  502. }
  503. // Generates (at most) requested_digits of input number w.
  504. // w is a floating-point number (DiyFp), consisting of a significand and an
  505. // exponent. Its exponent is bounded by kMinimalTargetExponent and
  506. // kMaximalTargetExponent.
  507. // Hence -60 <= w.e() <= -32.
  508. //
  509. // Returns false if it fails, in which case the generated digits in the buffer
  510. // should not be used.
  511. // Preconditions:
  512. // * w is correct up to 1 ulp (unit in the last place). That
  513. // is, its error must be strictly less than a unit of its last digit.
  514. // * kMinimalTargetExponent <= w.e() <= kMaximalTargetExponent
  515. //
  516. // Postconditions: returns false if procedure fails.
  517. // otherwise:
  518. // * buffer is not null-terminated, but length contains the number of
  519. // digits.
  520. // * the representation in buffer is the most precise representation of
  521. // requested_digits digits.
  522. // * buffer contains at most requested_digits digits of w. If there are less
  523. // than requested_digits digits then some trailing '0's have been removed.
  524. // * kappa is such that
  525. // w = buffer * 10^kappa + eps with |eps| < 10^kappa / 2.
  526. //
  527. // Remark: This procedure takes into account the imprecision of its input
  528. // numbers. If the precision is not enough to guarantee all the postconditions
  529. // then false is returned. This usually happens rarely, but the failure-rate
  530. // increases with higher requested_digits.
  531. static bool DigitGenCounted(
  532. in DiyFp w,
  533. int requested_digits,
  534. ref DtoaBuilder buffer,
  535. out int kappa)
  536. {
  537. Debug.Assert(MinimalTargetExponent <= w.E && w.E <= MaximalTargetExponent);
  538. // w is assumed to have an error less than 1 unit. Whenever w is scaled we
  539. // also scale its error.
  540. ulong w_error = 1;
  541. // We cut the input number into two parts: the integral digits and the
  542. // fractional digits. We don't emit any decimal separator, but adapt kappa
  543. // instead. Example: instead of writing "1.2" we put "12" into the buffer and
  544. // increase kappa by 1.
  545. DiyFp one = new DiyFp(((ulong) 1) << -w.E, w.E);
  546. // Division by one is a shift.
  547. uint integrals = (uint) (w.F >> -one.E);
  548. // Modulo by one is an and.
  549. ulong fractionals = w.F & (one.F - 1);
  550. BiggestPowerTen(integrals, DiyFp.KSignificandSize - (-one.E), out var divisor, out var divisor_exponent);
  551. kappa = divisor_exponent + 1;
  552. // Loop invariant: buffer = w / 10^kappa (integer division)
  553. // The invariant holds for the first iteration: kappa has been initialized
  554. // with the divisor exponent + 1. And the divisor is the biggest power of ten
  555. // that is smaller than 'integrals'.
  556. while (kappa > 0)
  557. {
  558. int digit = (int) (integrals / divisor);
  559. buffer.Append((char) ('0' + digit));
  560. requested_digits--;
  561. integrals %= divisor;
  562. kappa--;
  563. // Note that kappa now equals the exponent of the divisor and that the
  564. // invariant thus holds again.
  565. if (requested_digits == 0) break;
  566. divisor /= 10;
  567. }
  568. if (requested_digits == 0)
  569. {
  570. ulong rest = (((ulong) integrals) << -one.E) + fractionals;
  571. return RoundWeedCounted(ref buffer, rest,(ulong) divisor << -one.E, w_error, ref kappa);
  572. }
  573. // The integrals have been generated. We are at the point of the decimal
  574. // separator. In the following loop we simply multiply the remaining digits by
  575. // 10 and divide by one. We just need to pay attention to multiply associated
  576. // data (the 'unit'), too.
  577. // Note that the multiplication by 10 does not overflow, because w.e >= -60
  578. // and thus one.e >= -60.
  579. Debug.Assert(one.E >= -60);
  580. Debug.Assert(fractionals < one.F);
  581. while (requested_digits > 0 && fractionals > w_error) {
  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. }
  725. }