PageRenderTime 51ms CodeModel.GetById 14ms RepoModel.GetById 0ms app.codeStats 0ms

/libjava/classpath/native/fdlibm/dtoa.c

https://bitbucket.org/danchr/llvm-gcc
C | 919 lines | 731 code | 38 blank | 150 comment | 172 complexity | a16977b0e84cc4118e568aa73bc1d5ab MD5 | raw file
Possible License(s): GPL-2.0, LGPL-2.1, LGPL-2.0, BSD-3-Clause, CC-BY-SA-3.0
  1. /****************************************************************
  2. *
  3. * The author of this software is David M. Gay.
  4. *
  5. * Copyright (c) 1991, 2006 by AT&T.
  6. *
  7. * Permission to use, copy, modify, and distribute this software for any
  8. * purpose without fee is hereby granted, provided that this entire notice
  9. * is included in all copies of any software which is or includes a copy
  10. * or modification of this software and in all copies of the supporting
  11. * documentation for such software.
  12. *
  13. * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
  14. * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
  15. * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
  16. * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
  17. *
  18. ***************************************************************/
  19. /* Please send bug reports to
  20. David M. Gay
  21. AT&T Bell Laboratories, Room 2C-463
  22. 600 Mountain Avenue
  23. Murray Hill, NJ 07974-2070
  24. U.S.A.
  25. dmg@research.att.com or research!dmg
  26. */
  27. #include "mprec.h"
  28. #include <string.h>
  29. static int
  30. _DEFUN (quorem,
  31. (b, S),
  32. _Jv_Bigint * b _AND _Jv_Bigint * S)
  33. {
  34. int n;
  35. long borrow, y;
  36. unsigned long carry, q, ys;
  37. unsigned long *bx, *bxe, *sx, *sxe;
  38. #ifdef Pack_32
  39. long z;
  40. unsigned long si, zs;
  41. #endif
  42. n = S->_wds;
  43. #ifdef DEBUG
  44. /*debug*/ if (b->_wds > n)
  45. /*debug*/ Bug ("oversize b in quorem");
  46. #endif
  47. if (b->_wds < n)
  48. return 0;
  49. sx = S->_x;
  50. sxe = sx + --n;
  51. bx = b->_x;
  52. bxe = bx + n;
  53. q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
  54. #ifdef DEBUG
  55. /*debug*/ if (q > 9)
  56. /*debug*/ Bug ("oversized quotient in quorem");
  57. #endif
  58. if (q)
  59. {
  60. borrow = 0;
  61. carry = 0;
  62. do
  63. {
  64. #ifdef Pack_32
  65. si = *sx++;
  66. ys = (si & 0xffff) * q + carry;
  67. zs = (si >> 16) * q + (ys >> 16);
  68. carry = zs >> 16;
  69. y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
  70. borrow = y >> 16;
  71. Sign_Extend (borrow, y);
  72. z = (*bx >> 16) - (zs & 0xffff) + borrow;
  73. borrow = z >> 16;
  74. Sign_Extend (borrow, z);
  75. Storeinc (bx, z, y);
  76. #else
  77. ys = *sx++ * q + carry;
  78. carry = ys >> 16;
  79. y = *bx - (ys & 0xffff) + borrow;
  80. borrow = y >> 16;
  81. Sign_Extend (borrow, y);
  82. *bx++ = y & 0xffff;
  83. #endif
  84. }
  85. while (sx <= sxe);
  86. if (!*bxe)
  87. {
  88. bx = b->_x;
  89. while (--bxe > bx && !*bxe)
  90. --n;
  91. b->_wds = n;
  92. }
  93. }
  94. if (cmp (b, S) >= 0)
  95. {
  96. q++;
  97. borrow = 0;
  98. carry = 0;
  99. bx = b->_x;
  100. sx = S->_x;
  101. do
  102. {
  103. #ifdef Pack_32
  104. si = *sx++;
  105. ys = (si & 0xffff) + carry;
  106. zs = (si >> 16) + (ys >> 16);
  107. carry = zs >> 16;
  108. y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
  109. borrow = y >> 16;
  110. Sign_Extend (borrow, y);
  111. z = (*bx >> 16) - (zs & 0xffff) + borrow;
  112. borrow = z >> 16;
  113. Sign_Extend (borrow, z);
  114. Storeinc (bx, z, y);
  115. #else
  116. ys = *sx++ + carry;
  117. carry = ys >> 16;
  118. y = *bx - (ys & 0xffff) + borrow;
  119. borrow = y >> 16;
  120. Sign_Extend (borrow, y);
  121. *bx++ = y & 0xffff;
  122. #endif
  123. }
  124. while (sx <= sxe);
  125. bx = b->_x;
  126. bxe = bx + n;
  127. if (!*bxe)
  128. {
  129. while (--bxe > bx && !*bxe)
  130. --n;
  131. b->_wds = n;
  132. }
  133. }
  134. return q;
  135. }
  136. #ifdef DEBUG
  137. #include <stdio.h>
  138. void
  139. print (_Jv_Bigint * b)
  140. {
  141. int i, wds;
  142. unsigned long *x, y;
  143. wds = b->_wds;
  144. x = b->_x+wds;
  145. i = 0;
  146. do
  147. {
  148. x--;
  149. fprintf (stderr, "%08x", *x);
  150. }
  151. while (++i < wds);
  152. fprintf (stderr, "\n");
  153. }
  154. #endif
  155. /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
  156. *
  157. * Inspired by "How to Print Floating-Point Numbers Accurately" by
  158. * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
  159. *
  160. * Modifications:
  161. * 1. Rather than iterating, we use a simple numeric overestimate
  162. * to determine k = floor(log10(d)). We scale relevant
  163. * quantities using O(log2(k)) rather than O(k) multiplications.
  164. * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
  165. * try to generate digits strictly left to right. Instead, we
  166. * compute with fewer bits and propagate the carry if necessary
  167. * when rounding the final digit up. This is often faster.
  168. * 3. Under the assumption that input will be rounded nearest,
  169. * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
  170. * That is, we allow equality in stopping tests when the
  171. * round-nearest rule will give the same floating-point value
  172. * as would satisfaction of the stopping test with strict
  173. * inequality.
  174. * 4. We remove common factors of powers of 2 from relevant
  175. * quantities.
  176. * 5. When converting floating-point integers less than 1e16,
  177. * we use floating-point arithmetic rather than resorting
  178. * to multiple-precision integers.
  179. * 6. When asked to produce fewer than 15 digits, we first try
  180. * to get by with floating-point arithmetic; we resort to
  181. * multiple-precision integer arithmetic only if we cannot
  182. * guarantee that the floating-point calculation has given
  183. * the correctly rounded result. For k requested digits and
  184. * "uniformly" distributed input, the probability is
  185. * something like 10^(k-15) that we must resort to the long
  186. * calculation.
  187. */
  188. char *
  189. _DEFUN (_dtoa_r,
  190. (ptr, _d, mode, ndigits, decpt, sign, rve, float_type),
  191. struct _Jv_reent *ptr _AND
  192. double _d _AND
  193. int mode _AND
  194. int ndigits _AND
  195. int *decpt _AND
  196. int *sign _AND
  197. char **rve _AND
  198. int float_type)
  199. {
  200. /*
  201. float_type == 0 for double precision, 1 for float.
  202. Arguments ndigits, decpt, sign are similar to those
  203. of ecvt and fcvt; trailing zeros are suppressed from
  204. the returned string. If not null, *rve is set to point
  205. to the end of the return value. If d is +-Infinity or NaN,
  206. then *decpt is set to 9999.
  207. mode:
  208. 0 ==> shortest string that yields d when read in
  209. and rounded to nearest.
  210. 1 ==> like 0, but with Steele & White stopping rule;
  211. e.g. with IEEE P754 arithmetic , mode 0 gives
  212. 1e23 whereas mode 1 gives 9.999999999999999e22.
  213. 2 ==> max(1,ndigits) significant digits. This gives a
  214. return value similar to that of ecvt, except
  215. that trailing zeros are suppressed.
  216. 3 ==> through ndigits past the decimal point. This
  217. gives a return value similar to that from fcvt,
  218. except that trailing zeros are suppressed, and
  219. ndigits can be negative.
  220. 4-9 should give the same return values as 2-3, i.e.,
  221. 4 <= mode <= 9 ==> same return as mode
  222. 2 + (mode & 1). These modes are mainly for
  223. debugging; often they run slower but sometimes
  224. faster than modes 2-3.
  225. 4,5,8,9 ==> left-to-right digit generation.
  226. 6-9 ==> don't try fast floating-point estimate
  227. (if applicable).
  228. > 16 ==> Floating-point arg is treated as single precision.
  229. Values of mode other than 0-9 are treated as mode 0.
  230. Sufficient space is allocated to the return value
  231. to hold the suppressed trailing zeros.
  232. */
  233. int bbits, b2, b5, be, dig, i, ieps, ilim0, j, j1, k, k0,
  234. k_check, leftright, m2, m5, s2, s5, try_quick;
  235. int ilim = 0, ilim1 = 0, spec_case = 0;
  236. union double_union d, d2, eps;
  237. long L;
  238. #ifndef Sudden_Underflow
  239. int denorm;
  240. unsigned long x;
  241. #endif
  242. _Jv_Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
  243. double ds;
  244. char *s, *s0;
  245. d.d = _d;
  246. if (ptr->_result)
  247. {
  248. ptr->_result->_k = ptr->_result_k;
  249. ptr->_result->_maxwds = 1 << ptr->_result_k;
  250. Bfree (ptr, ptr->_result);
  251. ptr->_result = 0;
  252. }
  253. if (word0 (d) & Sign_bit)
  254. {
  255. /* set sign for everything, including 0's and NaNs */
  256. *sign = 1;
  257. word0 (d) &= ~Sign_bit; /* clear sign bit */
  258. }
  259. else
  260. *sign = 0;
  261. #if defined(IEEE_Arith) + defined(VAX)
  262. #ifdef IEEE_Arith
  263. if ((word0 (d) & Exp_mask) == Exp_mask)
  264. #else
  265. if (word0 (d) == 0x8000)
  266. #endif
  267. {
  268. /* Infinity or NaN */
  269. *decpt = 9999;
  270. s =
  271. #ifdef IEEE_Arith
  272. !word1 (d) && !(word0 (d) & 0xfffff) ? "Infinity" :
  273. #endif
  274. "NaN";
  275. if (rve)
  276. *rve =
  277. #ifdef IEEE_Arith
  278. s[3] ? s + 8 :
  279. #endif
  280. s + 3;
  281. return s;
  282. }
  283. #endif
  284. #ifdef IBM
  285. d.d += 0; /* normalize */
  286. #endif
  287. if (!d.d)
  288. {
  289. *decpt = 1;
  290. s = "0";
  291. if (rve)
  292. *rve = s + 1;
  293. return s;
  294. }
  295. b = d2b (ptr, d.d, &be, &bbits);
  296. #ifdef Sudden_Underflow
  297. i = (int) (word0 (d) >> Exp_shift1 & (Exp_mask >> Exp_shift1));
  298. #else
  299. if ((i = (int) (word0 (d) >> Exp_shift1 & (Exp_mask >> Exp_shift1))))
  300. {
  301. #endif
  302. d2.d = d.d;
  303. word0 (d2) &= Frac_mask1;
  304. word0 (d2) |= Exp_11;
  305. #ifdef IBM
  306. if (j = 11 - hi0bits (word0 (d2) & Frac_mask))
  307. d2.d /= 1 << j;
  308. #endif
  309. /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
  310. * log10(x) = log(x) / log(10)
  311. * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
  312. * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
  313. *
  314. * This suggests computing an approximation k to log10(d) by
  315. *
  316. * k = (i - Bias)*0.301029995663981
  317. * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
  318. *
  319. * We want k to be too large rather than too small.
  320. * The error in the first-order Taylor series approximation
  321. * is in our favor, so we just round up the constant enough
  322. * to compensate for any error in the multiplication of
  323. * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
  324. * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
  325. * adding 1e-13 to the constant term more than suffices.
  326. * Hence we adjust the constant term to 0.1760912590558.
  327. * (We could get a more accurate k by invoking log10,
  328. * but this is probably not worthwhile.)
  329. */
  330. i -= Bias;
  331. #ifdef IBM
  332. i <<= 2;
  333. i += j;
  334. #endif
  335. #ifndef Sudden_Underflow
  336. denorm = 0;
  337. }
  338. else
  339. {
  340. /* d is denormalized */
  341. i = bbits + be + (Bias + (P - 1) - 1);
  342. x = i > 32 ? word0 (d) << (64 - i) | word1 (d) >> (i - 32)
  343. : word1 (d) << (32 - i);
  344. d2.d = x;
  345. word0 (d2) -= 31 * Exp_msk1; /* adjust exponent */
  346. i -= (Bias + (P - 1) - 1) + 1;
  347. denorm = 1;
  348. }
  349. #endif
  350. ds = (d2.d - 1.5) * 0.289529654602168 + 0.1760912590558 + i * 0.301029995663981;
  351. k = (int) ds;
  352. if (ds < 0. && ds != k)
  353. k--; /* want k = floor(ds) */
  354. k_check = 1;
  355. if (k >= 0 && k <= Ten_pmax)
  356. {
  357. if (d.d < tens[k])
  358. k--;
  359. k_check = 0;
  360. }
  361. j = bbits - i - 1;
  362. if (j >= 0)
  363. {
  364. b2 = 0;
  365. s2 = j;
  366. }
  367. else
  368. {
  369. b2 = -j;
  370. s2 = 0;
  371. }
  372. if (k >= 0)
  373. {
  374. b5 = 0;
  375. s5 = k;
  376. s2 += k;
  377. }
  378. else
  379. {
  380. b2 -= k;
  381. b5 = -k;
  382. s5 = 0;
  383. }
  384. if (mode < 0 || mode > 9)
  385. mode = 0;
  386. try_quick = 1;
  387. if (mode > 5)
  388. {
  389. mode -= 4;
  390. try_quick = 0;
  391. }
  392. leftright = 1;
  393. switch (mode)
  394. {
  395. case 0:
  396. case 1:
  397. ilim = ilim1 = -1;
  398. i = 18;
  399. ndigits = 0;
  400. break;
  401. case 2:
  402. leftright = 0;
  403. /* no break */
  404. case 4:
  405. if (ndigits <= 0)
  406. ndigits = 1;
  407. ilim = ilim1 = i = ndigits;
  408. break;
  409. case 3:
  410. leftright = 0;
  411. /* no break */
  412. case 5:
  413. i = ndigits + k + 1;
  414. ilim = i;
  415. ilim1 = i - 1;
  416. if (i <= 0)
  417. i = 1;
  418. }
  419. j = sizeof (unsigned long);
  420. for (ptr->_result_k = 0; (int) (sizeof (_Jv_Bigint) - sizeof (unsigned long)) + j <= i;
  421. j <<= 1)
  422. ptr->_result_k++;
  423. ptr->_result = Balloc (ptr, ptr->_result_k);
  424. s = s0 = (char *) ptr->_result;
  425. if (ilim >= 0 && ilim <= Quick_max && try_quick)
  426. {
  427. /* Try to get by with floating-point arithmetic. */
  428. i = 0;
  429. d2.d = d.d;
  430. k0 = k;
  431. ilim0 = ilim;
  432. ieps = 2; /* conservative */
  433. if (k > 0)
  434. {
  435. ds = tens[k & 0xf];
  436. j = k >> 4;
  437. if (j & Bletch)
  438. {
  439. /* prevent overflows */
  440. j &= Bletch - 1;
  441. d.d /= bigtens[n_bigtens - 1];
  442. ieps++;
  443. }
  444. for (; j; j >>= 1, i++)
  445. if (j & 1)
  446. {
  447. ieps++;
  448. ds *= bigtens[i];
  449. }
  450. d.d /= ds;
  451. }
  452. else if ((j1 = -k))
  453. {
  454. d.d *= tens[j1 & 0xf];
  455. for (j = j1 >> 4; j; j >>= 1, i++)
  456. if (j & 1)
  457. {
  458. ieps++;
  459. d.d *= bigtens[i];
  460. }
  461. }
  462. if (k_check && d.d < 1. && ilim > 0)
  463. {
  464. if (ilim1 <= 0)
  465. goto fast_failed;
  466. ilim = ilim1;
  467. k--;
  468. d.d *= 10.;
  469. ieps++;
  470. }
  471. eps.d = ieps * d.d + 7.;
  472. word0 (eps) -= (P - 1) * Exp_msk1;
  473. if (ilim == 0)
  474. {
  475. S = mhi = 0;
  476. d.d -= 5.;
  477. if (d.d > eps.d)
  478. goto one_digit;
  479. if (d.d < -eps.d)
  480. goto no_digits;
  481. goto fast_failed;
  482. }
  483. #ifndef No_leftright
  484. if (leftright)
  485. {
  486. /* Use Steele & White method of only
  487. * generating digits needed.
  488. */
  489. eps.d = 0.5 / tens[ilim - 1] - eps.d;
  490. for (i = 0;;)
  491. {
  492. L = d.d;
  493. d.d -= L;
  494. *s++ = '0' + (int) L;
  495. if (d.d < eps.d)
  496. goto ret1;
  497. if (1. - d.d < eps.d)
  498. goto bump_up;
  499. if (++i >= ilim)
  500. break;
  501. eps.d *= 10.;
  502. d.d *= 10.;
  503. }
  504. }
  505. else
  506. {
  507. #endif
  508. /* Generate ilim digits, then fix them up. */
  509. eps.d *= tens[ilim - 1];
  510. for (i = 1;; i++, d.d *= 10.)
  511. {
  512. L = d.d;
  513. d.d -= L;
  514. *s++ = '0' + (int) L;
  515. if (i == ilim)
  516. {
  517. if (d.d > 0.5 + eps.d)
  518. goto bump_up;
  519. else if (d.d < 0.5 - eps.d)
  520. {
  521. while (*--s == '0');
  522. s++;
  523. goto ret1;
  524. }
  525. break;
  526. }
  527. }
  528. #ifndef No_leftright
  529. }
  530. #endif
  531. fast_failed:
  532. s = s0;
  533. d.d = d2.d;
  534. k = k0;
  535. ilim = ilim0;
  536. }
  537. /* Do we have a "small" integer? */
  538. if (be >= 0 && k <= Int_max)
  539. {
  540. /* Yes. */
  541. ds = tens[k];
  542. if (ndigits < 0 && ilim <= 0)
  543. {
  544. S = mhi = 0;
  545. if (ilim < 0 || d.d <= 5 * ds)
  546. goto no_digits;
  547. goto one_digit;
  548. }
  549. for (i = 1;; i++)
  550. {
  551. L = d.d / ds;
  552. d.d -= L * ds;
  553. #ifdef Check_FLT_ROUNDS
  554. /* If FLT_ROUNDS == 2, L will usually be high by 1 */
  555. if (d.d < 0)
  556. {
  557. L--;
  558. d.d += ds;
  559. }
  560. #endif
  561. *s++ = '0' + (int) L;
  562. if (i == ilim)
  563. {
  564. d.d += d.d;
  565. if (d.d > ds || (d.d == ds && L & 1))
  566. {
  567. bump_up:
  568. while (*--s == '9')
  569. if (s == s0)
  570. {
  571. k++;
  572. *s = '0';
  573. break;
  574. }
  575. ++*s++;
  576. }
  577. break;
  578. }
  579. if (!(d.d *= 10.))
  580. break;
  581. }
  582. goto ret1;
  583. }
  584. m2 = b2;
  585. m5 = b5;
  586. mhi = mlo = 0;
  587. if (leftright)
  588. {
  589. if (mode < 2)
  590. {
  591. i =
  592. #ifndef Sudden_Underflow
  593. denorm ? be + (Bias + (P - 1) - 1 + 1) :
  594. #endif
  595. #ifdef IBM
  596. 1 + 4 * P - 3 - bbits + ((bbits + be - 1) & 3);
  597. #else
  598. 1 + P - bbits;
  599. #endif
  600. }
  601. else
  602. {
  603. j = ilim - 1;
  604. if (m5 >= j)
  605. m5 -= j;
  606. else
  607. {
  608. s5 += j -= m5;
  609. b5 += j;
  610. m5 = 0;
  611. }
  612. if ((i = ilim) < 0)
  613. {
  614. m2 -= i;
  615. i = 0;
  616. }
  617. }
  618. b2 += i;
  619. s2 += i;
  620. mhi = i2b (ptr, 1);
  621. }
  622. if (m2 > 0 && s2 > 0)
  623. {
  624. i = m2 < s2 ? m2 : s2;
  625. b2 -= i;
  626. m2 -= i;
  627. s2 -= i;
  628. }
  629. if (b5 > 0)
  630. {
  631. if (leftright)
  632. {
  633. if (m5 > 0)
  634. {
  635. mhi = pow5mult (ptr, mhi, m5);
  636. b1 = mult (ptr, mhi, b);
  637. Bfree (ptr, b);
  638. b = b1;
  639. }
  640. if ((j = b5 - m5))
  641. b = pow5mult (ptr, b, j);
  642. }
  643. else
  644. b = pow5mult (ptr, b, b5);
  645. }
  646. S = i2b (ptr, 1);
  647. if (s5 > 0)
  648. S = pow5mult (ptr, S, s5);
  649. /* Check for special case that d is a normalized power of 2. */
  650. if (mode < 2)
  651. {
  652. if (!word1 (d) && !(word0 (d) & Bndry_mask)
  653. #ifndef Sudden_Underflow
  654. && word0(d) & Exp_mask
  655. #endif
  656. )
  657. {
  658. /* The special case */
  659. b2 += Log2P;
  660. s2 += Log2P;
  661. spec_case = 1;
  662. }
  663. else
  664. spec_case = 0;
  665. }
  666. /* Arrange for convenient computation of quotients:
  667. * shift left if necessary so divisor has 4 leading 0 bits.
  668. *
  669. * Perhaps we should just compute leading 28 bits of S once
  670. * and for all and pass them and a shift to quorem, so it
  671. * can do shifts and ors to compute the numerator for q.
  672. */
  673. #ifdef Pack_32
  674. if ((i = ((s5 ? 32 - hi0bits (S->_x[S->_wds - 1]) : 1) + s2) & 0x1f))
  675. i = 32 - i;
  676. #else
  677. if ((i = ((s5 ? 32 - hi0bits (S->_x[S->_wds - 1]) : 1) + s2) & 0xf))
  678. i = 16 - i;
  679. #endif
  680. if (i > 4)
  681. {
  682. i -= 4;
  683. b2 += i;
  684. m2 += i;
  685. s2 += i;
  686. }
  687. else if (i < 4)
  688. {
  689. i += 28;
  690. b2 += i;
  691. m2 += i;
  692. s2 += i;
  693. }
  694. if (b2 > 0)
  695. b = lshift (ptr, b, b2);
  696. if (s2 > 0)
  697. S = lshift (ptr, S, s2);
  698. if (k_check)
  699. {
  700. if (cmp (b, S) < 0)
  701. {
  702. k--;
  703. b = multadd (ptr, b, 10, 0); /* we botched the k estimate */
  704. if (leftright)
  705. mhi = multadd (ptr, mhi, 10, 0);
  706. ilim = ilim1;
  707. }
  708. }
  709. if (ilim <= 0 && mode > 2)
  710. {
  711. if (ilim < 0 || cmp (b, S = multadd (ptr, S, 5, 0)) <= 0)
  712. {
  713. /* no digits, fcvt style */
  714. no_digits:
  715. k = -1 - ndigits;
  716. goto ret;
  717. }
  718. one_digit:
  719. *s++ = '1';
  720. k++;
  721. goto ret;
  722. }
  723. if (leftright)
  724. {
  725. if (m2 > 0)
  726. mhi = lshift (ptr, mhi, m2);
  727. /* Single precision case, */
  728. if (float_type)
  729. mhi = lshift (ptr, mhi, 29);
  730. /* Compute mlo -- check for special case
  731. * that d is a normalized power of 2.
  732. */
  733. mlo = mhi;
  734. if (spec_case)
  735. {
  736. mhi = Balloc (ptr, mhi->_k);
  737. Bcopy (mhi, mlo);
  738. mhi = lshift (ptr, mhi, Log2P);
  739. }
  740. for (i = 1;; i++)
  741. {
  742. dig = quorem (b, S) + '0';
  743. /* Do we yet have the shortest decimal string
  744. * that will round to d?
  745. */
  746. j = cmp (b, mlo);
  747. delta = diff (ptr, S, mhi);
  748. j1 = delta->_sign ? 1 : cmp (b, delta);
  749. Bfree (ptr, delta);
  750. #ifndef ROUND_BIASED
  751. if (j1 == 0 && !mode && !(word1 (d) & 1))
  752. {
  753. if (dig == '9')
  754. goto round_9_up;
  755. if (j > 0)
  756. dig++;
  757. *s++ = dig;
  758. goto ret;
  759. }
  760. #endif
  761. if (j < 0 || (j == 0 && !mode
  762. #ifndef ROUND_BIASED
  763. && !(word1 (d) & 1)
  764. #endif
  765. ))
  766. {
  767. if (j1 > 0)
  768. {
  769. b = lshift (ptr, b, 1);
  770. j1 = cmp (b, S);
  771. if ((j1 > 0 || (j1 == 0 && dig & 1))
  772. && dig++ == '9')
  773. goto round_9_up;
  774. }
  775. *s++ = dig;
  776. goto ret;
  777. }
  778. if (j1 > 0)
  779. {
  780. if (dig == '9')
  781. { /* possible if i == 1 */
  782. round_9_up:
  783. *s++ = '9';
  784. goto roundoff;
  785. }
  786. *s++ = dig + 1;
  787. goto ret;
  788. }
  789. *s++ = dig;
  790. if (i == ilim)
  791. break;
  792. b = multadd (ptr, b, 10, 0);
  793. if (mlo == mhi)
  794. mlo = mhi = multadd (ptr, mhi, 10, 0);
  795. else
  796. {
  797. mlo = multadd (ptr, mlo, 10, 0);
  798. mhi = multadd (ptr, mhi, 10, 0);
  799. }
  800. }
  801. }
  802. else
  803. for (i = 1;; i++)
  804. {
  805. *s++ = dig = quorem (b, S) + '0';
  806. if (i >= ilim)
  807. break;
  808. b = multadd (ptr, b, 10, 0);
  809. }
  810. /* Round off last digit */
  811. b = lshift (ptr, b, 1);
  812. j = cmp (b, S);
  813. if (j > 0 || (j == 0 && dig & 1))
  814. {
  815. roundoff:
  816. while (*--s == '9')
  817. if (s == s0)
  818. {
  819. k++;
  820. *s++ = '1';
  821. goto ret;
  822. }
  823. ++*s++;
  824. }
  825. else
  826. {
  827. while (*--s == '0');
  828. s++;
  829. }
  830. ret:
  831. Bfree (ptr, S);
  832. if (mhi)
  833. {
  834. if (mlo && mlo != mhi)
  835. Bfree (ptr, mlo);
  836. Bfree (ptr, mhi);
  837. }
  838. ret1:
  839. Bfree (ptr, b);
  840. *s = 0;
  841. *decpt = k + 1;
  842. if (rve)
  843. *rve = s;
  844. return s0;
  845. }
  846. _VOID
  847. _DEFUN (_dtoa,
  848. (_d, mode, ndigits, decpt, sign, rve, buf, float_type),
  849. double _d _AND
  850. int mode _AND
  851. int ndigits _AND
  852. int *decpt _AND
  853. int *sign _AND
  854. char **rve _AND
  855. char *buf _AND
  856. int float_type)
  857. {
  858. struct _Jv_reent reent;
  859. char *p;
  860. int i;
  861. memset (&reent, 0, sizeof reent);
  862. p = _dtoa_r (&reent, _d, mode, ndigits, decpt, sign, rve, float_type);
  863. strcpy (buf, p);
  864. for (i = 0; i < reent._result_k; ++i)
  865. {
  866. struct _Jv_Bigint *l = reent._freelist[i];
  867. while (l)
  868. {
  869. struct _Jv_Bigint *next = l->_next;
  870. free (l);
  871. l = next;
  872. }
  873. }
  874. if (reent._freelist)
  875. free (reent._freelist);
  876. }