strtodg.c 23 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156
  1. /****************************************************************
  2. The author of this software is David M. Gay.
  3. Copyright (C) 1998-2001 by Lucent Technologies
  4. All Rights Reserved
  5. Permission to use, copy, modify, and distribute this software and
  6. its documentation for any purpose and without fee is hereby
  7. granted, provided that the above copyright notice appear in all
  8. copies and that both that the copyright notice and this
  9. permission notice and warranty disclaimer appear in supporting
  10. documentation, and that the name of Lucent or any of its entities
  11. not be used in advertising or publicity pertaining to
  12. distribution of the software without specific, written prior
  13. permission.
  14. LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
  15. INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
  16. IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
  17. SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  18. WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
  19. IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
  20. ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
  21. THIS SOFTWARE.
  22. ****************************************************************/
  23. /* Please send bug reports to David M. Gay (dmg at acm dot org,
  24. * with " at " changed at "@" and " dot " changed to "."). */
  25. #include "gdtoaimp.h"
  26. #ifdef USE_LOCALE
  27. #include "locale.h"
  28. #endif
  29. static CONST int
  30. fivesbits[] = { 0, 3, 5, 7, 10, 12, 14, 17, 19, 21,
  31. 24, 26, 28, 31, 33, 35, 38, 40, 42, 45,
  32. 47, 49, 52
  33. #ifdef VAX
  34. , 54, 56
  35. #endif
  36. };
  37. Bigint *
  38. #ifdef KR_headers
  39. increment(b) Bigint *b;
  40. #else
  41. increment(Bigint *b)
  42. #endif
  43. {
  44. ULong *x, *xe;
  45. Bigint *b1;
  46. #ifdef Pack_16
  47. ULong carry = 1, y;
  48. #endif
  49. x = b->x;
  50. xe = x + b->wds;
  51. #ifdef Pack_32
  52. do {
  53. if (*x < (ULong)0xffffffffL) {
  54. ++*x;
  55. return b;
  56. }
  57. *x++ = 0;
  58. } while(x < xe);
  59. #else
  60. do {
  61. y = *x + carry;
  62. carry = y >> 16;
  63. *x++ = y & 0xffff;
  64. if (!carry)
  65. return b;
  66. } while(x < xe);
  67. if (carry)
  68. #endif
  69. {
  70. if (b->wds >= b->maxwds) {
  71. b1 = Balloc(b->k+1);
  72. if (b1 == NULL)
  73. return (NULL);
  74. Bcopy(b1,b);
  75. Bfree(b);
  76. b = b1;
  77. }
  78. b->x[b->wds++] = 1;
  79. }
  80. return b;
  81. }
  82. void
  83. #ifdef KR_headers
  84. decrement(b) Bigint *b;
  85. #else
  86. decrement(Bigint *b)
  87. #endif
  88. {
  89. ULong *x, *xe;
  90. #ifdef Pack_16
  91. ULong borrow = 1, y;
  92. #endif
  93. x = b->x;
  94. xe = x + b->wds;
  95. #ifdef Pack_32
  96. do {
  97. if (*x) {
  98. --*x;
  99. break;
  100. }
  101. *x++ = 0xffffffffL;
  102. }
  103. while(x < xe);
  104. #else
  105. do {
  106. y = *x - borrow;
  107. borrow = (y & 0x10000) >> 16;
  108. *x++ = y & 0xffff;
  109. } while(borrow && x < xe);
  110. #endif
  111. }
  112. static int
  113. #ifdef KR_headers
  114. all_on(b, n) Bigint *b; int n;
  115. #else
  116. all_on(Bigint *b, int n)
  117. #endif
  118. {
  119. ULong *x, *xe;
  120. x = b->x;
  121. xe = x + (n >> kshift);
  122. while(x < xe)
  123. if ((*x++ & ALL_ON) != ALL_ON)
  124. return 0;
  125. if (n &= kmask)
  126. return ((*x | (ALL_ON << n)) & ALL_ON) == ALL_ON;
  127. return 1;
  128. }
  129. Bigint *
  130. #ifdef KR_headers
  131. set_ones(b, n) Bigint *b; int n;
  132. #else
  133. set_ones(Bigint *b, int n)
  134. #endif
  135. {
  136. int k;
  137. ULong *x, *xe;
  138. k = (n + ((1 << kshift) - 1)) >> kshift;
  139. if (b->k < k) {
  140. Bfree(b);
  141. b = Balloc(k);
  142. if (b == NULL)
  143. return (NULL);
  144. }
  145. k = n >> kshift;
  146. if (n &= kmask)
  147. k++;
  148. b->wds = k;
  149. x = b->x;
  150. xe = x + k;
  151. while(x < xe)
  152. *x++ = ALL_ON;
  153. if (n)
  154. x[-1] >>= ULbits - n;
  155. return b;
  156. }
  157. static int
  158. rvOK
  159. #ifdef KR_headers
  160. (d, fpi, exp, bits, exact, rd, irv)
  161. U *d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv;
  162. #else
  163. (U *d, FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv)
  164. #endif
  165. {
  166. Bigint *b;
  167. ULong carry, inex, lostbits;
  168. int bdif, e, j, k, k1, nb, rv;
  169. carry = rv = 0;
  170. b = d2b(dval(d), &e, &bdif);
  171. if (b == NULL) {
  172. *irv = STRTOG_NoMemory;
  173. return (1);
  174. }
  175. bdif -= nb = fpi->nbits;
  176. e += bdif;
  177. if (bdif <= 0) {
  178. if (exact)
  179. goto trunc;
  180. goto ret;
  181. }
  182. if (P == nb) {
  183. if (
  184. #ifndef IMPRECISE_INEXACT
  185. exact &&
  186. #endif
  187. fpi->rounding ==
  188. #ifdef RND_PRODQUOT
  189. FPI_Round_near
  190. #else
  191. Flt_Rounds
  192. #endif
  193. ) goto trunc;
  194. goto ret;
  195. }
  196. switch(rd) {
  197. case 1: /* round down (toward -Infinity) */
  198. goto trunc;
  199. case 2: /* round up (toward +Infinity) */
  200. break;
  201. default: /* round near */
  202. k = bdif - 1;
  203. if (k < 0)
  204. goto trunc;
  205. if (!k) {
  206. if (!exact)
  207. goto ret;
  208. if (b->x[0] & 2)
  209. break;
  210. goto trunc;
  211. }
  212. if (b->x[k>>kshift] & ((ULong)1 << (k & kmask)))
  213. break;
  214. goto trunc;
  215. }
  216. /* "break" cases: round up 1 bit, then truncate; bdif > 0 */
  217. carry = 1;
  218. trunc:
  219. inex = lostbits = 0;
  220. if (bdif > 0) {
  221. if ( (lostbits = any_on(b, bdif)) !=0)
  222. inex = STRTOG_Inexlo;
  223. rshift(b, bdif);
  224. if (carry) {
  225. inex = STRTOG_Inexhi;
  226. b = increment(b);
  227. if (b == NULL) {
  228. *irv = STRTOG_NoMemory;
  229. return (1);
  230. }
  231. if ( (j = nb & kmask) !=0)
  232. j = ULbits - j;
  233. if (hi0bits(b->x[b->wds - 1]) != j) {
  234. if (!lostbits)
  235. lostbits = b->x[0] & 1;
  236. rshift(b, 1);
  237. e++;
  238. }
  239. }
  240. }
  241. else if (bdif < 0) {
  242. b = lshift(b, -bdif);
  243. if (b == NULL) {
  244. *irv = STRTOG_NoMemory;
  245. return (1);
  246. }
  247. }
  248. if (e < fpi->emin) {
  249. k = fpi->emin - e;
  250. e = fpi->emin;
  251. if (k > nb || fpi->sudden_underflow) {
  252. b->wds = inex = 0;
  253. *irv = STRTOG_Underflow | STRTOG_Inexlo;
  254. }
  255. else {
  256. k1 = k - 1;
  257. if (k1 > 0 && !lostbits)
  258. lostbits = any_on(b, k1);
  259. if (!lostbits && !exact)
  260. goto ret;
  261. lostbits |=
  262. carry = b->x[k1>>kshift] & (1 << (k1 & kmask));
  263. rshift(b, k);
  264. *irv = STRTOG_Denormal;
  265. if (carry) {
  266. b = increment(b);
  267. if (b == NULL) {
  268. *irv = STRTOG_NoMemory;
  269. return (1);
  270. }
  271. inex = STRTOG_Inexhi | STRTOG_Underflow;
  272. }
  273. else if (lostbits)
  274. inex = STRTOG_Inexlo | STRTOG_Underflow;
  275. }
  276. }
  277. else if (e > fpi->emax) {
  278. e = fpi->emax + 1;
  279. *irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
  280. #ifndef NO_ERRNO
  281. errno = ERANGE;
  282. #endif
  283. b->wds = inex = 0;
  284. }
  285. *exp = e;
  286. copybits(bits, nb, b);
  287. *irv |= inex;
  288. rv = 1;
  289. ret:
  290. Bfree(b);
  291. return rv;
  292. }
  293. static int
  294. #ifdef KR_headers
  295. mantbits(d) U *d;
  296. #else
  297. mantbits(U *d)
  298. #endif
  299. {
  300. ULong L;
  301. #ifdef VAX
  302. L = word1(d) << 16 | word1(d) >> 16;
  303. if (L)
  304. #else
  305. if ( (L = word1(d)) !=0)
  306. #endif
  307. return P - lo0bits(&L);
  308. #ifdef VAX
  309. L = word0(d) << 16 | word0(d) >> 16 | Exp_msk11;
  310. #else
  311. L = word0(d) | Exp_msk1;
  312. #endif
  313. return P - 32 - lo0bits(&L);
  314. }
  315. int
  316. strtodg
  317. #ifdef KR_headers
  318. (s00, se, fpi, exp, bits)
  319. CONST char *s00; char **se; FPI *fpi; Long *exp; ULong *bits;
  320. #else
  321. (CONST char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
  322. #endif
  323. {
  324. int abe, abits, asub;
  325. int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, denorm;
  326. int dsign, e, e1, e2, emin, esign, finished, i, inex, irv;
  327. int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign;
  328. int sudden_underflow;
  329. CONST char *s, *s0, *s1;
  330. double adj0, tol;
  331. Long L;
  332. U adj, rv;
  333. ULong *b, *be, y, z;
  334. Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
  335. #ifdef USE_LOCALE /*{{*/
  336. #ifdef NO_LOCALE_CACHE
  337. char *decimalpoint = localeconv()->decimal_point;
  338. int dplen = strlen(decimalpoint);
  339. #else
  340. char *decimalpoint;
  341. static char *decimalpoint_cache;
  342. static int dplen;
  343. if (!(s0 = decimalpoint_cache)) {
  344. s0 = localeconv()->decimal_point;
  345. if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
  346. strlcpy(decimalpoint_cache, s0, strlen(s0) + 1);
  347. s0 = decimalpoint_cache;
  348. }
  349. dplen = strlen(s0);
  350. }
  351. decimalpoint = (char*)s0;
  352. #endif /*NO_LOCALE_CACHE*/
  353. #else /*USE_LOCALE}{*/
  354. #define dplen 1
  355. #endif /*USE_LOCALE}}*/
  356. irv = STRTOG_Zero;
  357. denorm = sign = nz0 = nz = 0;
  358. dval(&rv) = 0.;
  359. rvb = 0;
  360. nbits = fpi->nbits;
  361. for(s = s00;;s++) switch(*s) {
  362. case '-':
  363. sign = 1;
  364. /* no break */
  365. case '+':
  366. if (*++s)
  367. goto break2;
  368. /* no break */
  369. case 0:
  370. sign = 0;
  371. irv = STRTOG_NoNumber;
  372. s = s00;
  373. goto ret;
  374. case '\t':
  375. case '\n':
  376. case '\v':
  377. case '\f':
  378. case '\r':
  379. case ' ':
  380. continue;
  381. default:
  382. goto break2;
  383. }
  384. break2:
  385. if (*s == '0') {
  386. #ifndef NO_HEX_FP
  387. switch(s[1]) {
  388. case 'x':
  389. case 'X':
  390. irv = gethex(&s, fpi, exp, &rvb, sign);
  391. if (irv == STRTOG_NoMemory)
  392. return (STRTOG_NoMemory);
  393. if (irv == STRTOG_NoNumber) {
  394. s = s00;
  395. sign = 0;
  396. }
  397. goto ret;
  398. }
  399. #endif
  400. nz0 = 1;
  401. while(*++s == '0') ;
  402. if (!*s)
  403. goto ret;
  404. }
  405. sudden_underflow = fpi->sudden_underflow;
  406. s0 = s;
  407. y = z = 0;
  408. for(decpt = nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
  409. if (nd < 9)
  410. y = 10*y + c - '0';
  411. else if (nd < 16)
  412. z = 10*z + c - '0';
  413. nd0 = nd;
  414. #ifdef USE_LOCALE
  415. if (c == *decimalpoint) {
  416. for(i = 1; decimalpoint[i]; ++i)
  417. if (s[i] != decimalpoint[i])
  418. goto dig_done;
  419. s += i;
  420. c = *s;
  421. #else
  422. if (c == '.') {
  423. c = *++s;
  424. #endif
  425. decpt = 1;
  426. if (!nd) {
  427. for(; c == '0'; c = *++s)
  428. nz++;
  429. if (c > '0' && c <= '9') {
  430. s0 = s;
  431. nf += nz;
  432. nz = 0;
  433. goto have_dig;
  434. }
  435. goto dig_done;
  436. }
  437. for(; c >= '0' && c <= '9'; c = *++s) {
  438. have_dig:
  439. nz++;
  440. if (c -= '0') {
  441. nf += nz;
  442. for(i = 1; i < nz; i++)
  443. if (nd++ < 9)
  444. y *= 10;
  445. else if (nd <= DBL_DIG + 1)
  446. z *= 10;
  447. if (nd++ < 9)
  448. y = 10*y + c;
  449. else if (nd <= DBL_DIG + 1)
  450. z = 10*z + c;
  451. nz = 0;
  452. }
  453. }
  454. }/*}*/
  455. dig_done:
  456. e = 0;
  457. if (c == 'e' || c == 'E') {
  458. if (!nd && !nz && !nz0) {
  459. irv = STRTOG_NoNumber;
  460. s = s00;
  461. goto ret;
  462. }
  463. s00 = s;
  464. esign = 0;
  465. switch(c = *++s) {
  466. case '-':
  467. esign = 1;
  468. case '+':
  469. c = *++s;
  470. }
  471. if (c >= '0' && c <= '9') {
  472. while(c == '0')
  473. c = *++s;
  474. if (c > '0' && c <= '9') {
  475. L = c - '0';
  476. s1 = s;
  477. while((c = *++s) >= '0' && c <= '9')
  478. L = 10*L + c - '0';
  479. if (s - s1 > 8 || L > 19999)
  480. /* Avoid confusion from exponents
  481. * so large that e might overflow.
  482. */
  483. e = 19999; /* safe for 16 bit ints */
  484. else
  485. e = (int)L;
  486. if (esign)
  487. e = -e;
  488. }
  489. else
  490. e = 0;
  491. }
  492. else
  493. s = s00;
  494. }
  495. if (!nd) {
  496. if (!nz && !nz0) {
  497. #ifdef INFNAN_CHECK
  498. /* Check for Nan and Infinity */
  499. if (!decpt)
  500. switch(c) {
  501. case 'i':
  502. case 'I':
  503. if (match(&s,"nf")) {
  504. --s;
  505. if (!match(&s,"inity"))
  506. ++s;
  507. irv = STRTOG_Infinite;
  508. goto infnanexp;
  509. }
  510. break;
  511. case 'n':
  512. case 'N':
  513. if (match(&s, "an")) {
  514. irv = STRTOG_NaN;
  515. *exp = fpi->emax + 1;
  516. #ifndef No_Hex_NaN
  517. if (*s == '(') /*)*/
  518. irv = hexnan(&s, fpi, bits);
  519. #endif
  520. goto infnanexp;
  521. }
  522. }
  523. #endif /* INFNAN_CHECK */
  524. irv = STRTOG_NoNumber;
  525. s = s00;
  526. }
  527. goto ret;
  528. }
  529. irv = STRTOG_Normal;
  530. e1 = e -= nf;
  531. rd = 0;
  532. switch(fpi->rounding & 3) {
  533. case FPI_Round_up:
  534. rd = 2 - sign;
  535. break;
  536. case FPI_Round_zero:
  537. rd = 1;
  538. break;
  539. case FPI_Round_down:
  540. rd = 1 + sign;
  541. }
  542. /* Now we have nd0 digits, starting at s0, followed by a
  543. * decimal point, followed by nd-nd0 digits. The number we're
  544. * after is the integer represented by those digits times
  545. * 10**e */
  546. if (!nd0)
  547. nd0 = nd;
  548. k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
  549. dval(&rv) = y;
  550. if (k > 9)
  551. dval(&rv) = tens[k - 9] * dval(&rv) + z;
  552. bd0 = 0;
  553. if (nbits <= P && nd <= DBL_DIG) {
  554. if (!e) {
  555. if (rvOK(&rv, fpi, exp, bits, 1, rd, &irv)) {
  556. if (irv == STRTOG_NoMemory)
  557. return (STRTOG_NoMemory);
  558. goto ret;
  559. }
  560. }
  561. else if (e > 0) {
  562. if (e <= Ten_pmax) {
  563. #ifdef VAX
  564. goto vax_ovfl_check;
  565. #else
  566. i = fivesbits[e] + mantbits(&rv) <= P;
  567. /* rv = */ rounded_product(dval(&rv), tens[e]);
  568. if (rvOK(&rv, fpi, exp, bits, i, rd, &irv)) {
  569. if (irv == STRTOG_NoMemory)
  570. return (STRTOG_NoMemory);
  571. goto ret;
  572. }
  573. e1 -= e;
  574. goto rv_notOK;
  575. #endif
  576. }
  577. i = DBL_DIG - nd;
  578. if (e <= Ten_pmax + i) {
  579. /* A fancier test would sometimes let us do
  580. * this for larger i values.
  581. */
  582. e2 = e - i;
  583. e1 -= i;
  584. dval(&rv) *= tens[i];
  585. #ifdef VAX
  586. /* VAX exponent range is so narrow we must
  587. * worry about overflow here...
  588. */
  589. vax_ovfl_check:
  590. dval(&adj) = dval(&rv);
  591. word0(&adj) -= P*Exp_msk1;
  592. /* adj = */ rounded_product(dval(&adj), tens[e2]);
  593. if ((word0(&adj) & Exp_mask)
  594. > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
  595. goto rv_notOK;
  596. word0(&adj) += P*Exp_msk1;
  597. dval(&rv) = dval(&adj);
  598. #else
  599. /* rv = */ rounded_product(dval(&rv), tens[e2]);
  600. #endif
  601. if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv)) {
  602. if (irv == STRTOG_NoMemory)
  603. return (STRTOG_NoMemory);
  604. goto ret;
  605. }
  606. e1 -= e2;
  607. }
  608. }
  609. #ifndef Inaccurate_Divide
  610. else if (e >= -Ten_pmax) {
  611. /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
  612. if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv)) {
  613. if (irv == STRTOG_NoMemory)
  614. return (STRTOG_NoMemory);
  615. goto ret;
  616. }
  617. e1 -= e;
  618. }
  619. #endif
  620. }
  621. rv_notOK:
  622. e1 += nd - k;
  623. /* Get starting approximation = rv * 10**e1 */
  624. e2 = 0;
  625. if (e1 > 0) {
  626. if ( (i = e1 & 15) !=0)
  627. dval(&rv) *= tens[i];
  628. if (e1 &= ~15) {
  629. e1 >>= 4;
  630. while(e1 >= (1 << (n_bigtens-1))) {
  631. e2 += ((word0(&rv) & Exp_mask)
  632. >> Exp_shift1) - Bias;
  633. word0(&rv) &= ~Exp_mask;
  634. word0(&rv) |= Bias << Exp_shift1;
  635. dval(&rv) *= bigtens[n_bigtens-1];
  636. e1 -= 1 << (n_bigtens-1);
  637. }
  638. e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
  639. word0(&rv) &= ~Exp_mask;
  640. word0(&rv) |= Bias << Exp_shift1;
  641. for(j = 0; e1 > 0; j++, e1 >>= 1)
  642. if (e1 & 1)
  643. dval(&rv) *= bigtens[j];
  644. }
  645. }
  646. else if (e1 < 0) {
  647. e1 = -e1;
  648. if ( (i = e1 & 15) !=0)
  649. dval(&rv) /= tens[i];
  650. if (e1 &= ~15) {
  651. e1 >>= 4;
  652. while(e1 >= (1 << (n_bigtens-1))) {
  653. e2 += ((word0(&rv) & Exp_mask)
  654. >> Exp_shift1) - Bias;
  655. word0(&rv) &= ~Exp_mask;
  656. word0(&rv) |= Bias << Exp_shift1;
  657. dval(&rv) *= tinytens[n_bigtens-1];
  658. e1 -= 1 << (n_bigtens-1);
  659. }
  660. e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
  661. word0(&rv) &= ~Exp_mask;
  662. word0(&rv) |= Bias << Exp_shift1;
  663. for(j = 0; e1 > 0; j++, e1 >>= 1)
  664. if (e1 & 1)
  665. dval(&rv) *= tinytens[j];
  666. }
  667. }
  668. #ifdef IBM
  669. /* e2 is a correction to the (base 2) exponent of the return
  670. * value, reflecting adjustments above to avoid overflow in the
  671. * native arithmetic. For native IBM (base 16) arithmetic, we
  672. * must multiply e2 by 4 to change from base 16 to 2.
  673. */
  674. e2 <<= 2;
  675. #endif
  676. rvb = d2b(dval(&rv), &rve, &rvbits); /* rv = rvb * 2^rve */
  677. if (rvb == NULL)
  678. return (STRTOG_NoMemory);
  679. rve += e2;
  680. if ((j = rvbits - nbits) > 0) {
  681. rshift(rvb, j);
  682. rvbits = nbits;
  683. rve += j;
  684. }
  685. bb0 = 0; /* trailing zero bits in rvb */
  686. e2 = rve + rvbits - nbits;
  687. if (e2 > fpi->emax + 1)
  688. goto huge;
  689. rve1 = rve + rvbits - nbits;
  690. if (e2 < (emin = fpi->emin)) {
  691. denorm = 1;
  692. j = rve - emin;
  693. if (j > 0) {
  694. rvb = lshift(rvb, j);
  695. if (rvb == NULL)
  696. return (STRTOG_NoMemory);
  697. rvbits += j;
  698. }
  699. else if (j < 0) {
  700. rvbits += j;
  701. if (rvbits <= 0) {
  702. if (rvbits < -1) {
  703. ufl:
  704. rvb->wds = 0;
  705. rvb->x[0] = 0;
  706. *exp = emin;
  707. irv = STRTOG_Underflow | STRTOG_Inexlo;
  708. goto ret;
  709. }
  710. rvb->x[0] = rvb->wds = rvbits = 1;
  711. }
  712. else
  713. rshift(rvb, -j);
  714. }
  715. rve = rve1 = emin;
  716. if (sudden_underflow && e2 + 1 < emin)
  717. goto ufl;
  718. }
  719. /* Now the hard part -- adjusting rv to the correct value.*/
  720. /* Put digits into bd: true value = bd * 10^e */
  721. bd0 = s2b(s0, nd0, nd, y, dplen);
  722. if (bd0 == NULL)
  723. return (STRTOG_NoMemory);
  724. for(;;) {
  725. bd = Balloc(bd0->k);
  726. if (bd == NULL)
  727. return (STRTOG_NoMemory);
  728. Bcopy(bd, bd0);
  729. bb = Balloc(rvb->k);
  730. if (bb == NULL)
  731. return (STRTOG_NoMemory);
  732. Bcopy(bb, rvb);
  733. bbbits = rvbits - bb0;
  734. bbe = rve + bb0;
  735. bs = i2b(1);
  736. if (bs == NULL)
  737. return (STRTOG_NoMemory);
  738. if (e >= 0) {
  739. bb2 = bb5 = 0;
  740. bd2 = bd5 = e;
  741. }
  742. else {
  743. bb2 = bb5 = -e;
  744. bd2 = bd5 = 0;
  745. }
  746. if (bbe >= 0)
  747. bb2 += bbe;
  748. else
  749. bd2 -= bbe;
  750. bs2 = bb2;
  751. j = nbits + 1 - bbbits;
  752. i = bbe + bbbits - nbits;
  753. if (i < emin) /* denormal */
  754. j += i - emin;
  755. bb2 += j;
  756. bd2 += j;
  757. i = bb2 < bd2 ? bb2 : bd2;
  758. if (i > bs2)
  759. i = bs2;
  760. if (i > 0) {
  761. bb2 -= i;
  762. bd2 -= i;
  763. bs2 -= i;
  764. }
  765. if (bb5 > 0) {
  766. bs = pow5mult(bs, bb5);
  767. if (bs == NULL)
  768. return (STRTOG_NoMemory);
  769. bb1 = mult(bs, bb);
  770. if (bb1 == NULL)
  771. return (STRTOG_NoMemory);
  772. Bfree(bb);
  773. bb = bb1;
  774. }
  775. bb2 -= bb0;
  776. if (bb2 > 0) {
  777. bb = lshift(bb, bb2);
  778. if (bb == NULL)
  779. return (STRTOG_NoMemory);
  780. }
  781. else if (bb2 < 0)
  782. rshift(bb, -bb2);
  783. if (bd5 > 0) {
  784. bd = pow5mult(bd, bd5);
  785. if (bd == NULL)
  786. return (STRTOG_NoMemory);
  787. }
  788. if (bd2 > 0) {
  789. bd = lshift(bd, bd2);
  790. if (bd == NULL)
  791. return (STRTOG_NoMemory);
  792. }
  793. if (bs2 > 0) {
  794. bs = lshift(bs, bs2);
  795. if (bs == NULL)
  796. return (STRTOG_NoMemory);
  797. }
  798. asub = 1;
  799. inex = STRTOG_Inexhi;
  800. delta = diff(bb, bd);
  801. if (delta == NULL)
  802. return (STRTOG_NoMemory);
  803. if (delta->wds <= 1 && !delta->x[0])
  804. break;
  805. dsign = delta->sign;
  806. delta->sign = finished = 0;
  807. L = 0;
  808. i = cmp(delta, bs);
  809. if (rd && i <= 0) {
  810. irv = STRTOG_Normal;
  811. if ( (finished = dsign ^ (rd&1)) !=0) {
  812. if (dsign != 0) {
  813. irv |= STRTOG_Inexhi;
  814. goto adj1;
  815. }
  816. irv |= STRTOG_Inexlo;
  817. if (rve1 == emin)
  818. goto adj1;
  819. for(i = 0, j = nbits; j >= ULbits;
  820. i++, j -= ULbits) {
  821. if (rvb->x[i] & ALL_ON)
  822. goto adj1;
  823. }
  824. if (j > 1 && lo0bits(rvb->x + i) < j - 1)
  825. goto adj1;
  826. rve = rve1 - 1;
  827. rvb = set_ones(rvb, rvbits = nbits);
  828. if (rvb == NULL)
  829. return (STRTOG_NoMemory);
  830. break;
  831. }
  832. irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi;
  833. break;
  834. }
  835. if (i < 0) {
  836. /* Error is less than half an ulp -- check for
  837. * special case of mantissa a power of two.
  838. */
  839. irv = dsign
  840. ? STRTOG_Normal | STRTOG_Inexlo
  841. : STRTOG_Normal | STRTOG_Inexhi;
  842. if (dsign || bbbits > 1 || denorm || rve1 == emin)
  843. break;
  844. delta = lshift(delta,1);
  845. if (delta == NULL)
  846. return (STRTOG_NoMemory);
  847. if (cmp(delta, bs) > 0) {
  848. irv = STRTOG_Normal | STRTOG_Inexlo;
  849. goto drop_down;
  850. }
  851. break;
  852. }
  853. if (i == 0) {
  854. /* exactly half-way between */
  855. if (dsign) {
  856. if (denorm && all_on(rvb, rvbits)) {
  857. /*boundary case -- increment exponent*/
  858. rvb->wds = 1;
  859. rvb->x[0] = 1;
  860. rve = emin + nbits - (rvbits = 1);
  861. irv = STRTOG_Normal | STRTOG_Inexhi;
  862. denorm = 0;
  863. break;
  864. }
  865. irv = STRTOG_Normal | STRTOG_Inexlo;
  866. }
  867. else if (bbbits == 1) {
  868. irv = STRTOG_Normal;
  869. drop_down:
  870. /* boundary case -- decrement exponent */
  871. if (rve1 == emin) {
  872. irv = STRTOG_Normal | STRTOG_Inexhi;
  873. if (rvb->wds == 1 && rvb->x[0] == 1)
  874. sudden_underflow = 1;
  875. break;
  876. }
  877. rve -= nbits;
  878. rvb = set_ones(rvb, rvbits = nbits);
  879. if (rvb == NULL)
  880. return (STRTOG_NoMemory);
  881. break;
  882. }
  883. else
  884. irv = STRTOG_Normal | STRTOG_Inexhi;
  885. if ((bbbits < nbits && !denorm) || !(rvb->x[0] & 1))
  886. break;
  887. if (dsign) {
  888. rvb = increment(rvb);
  889. if (rvb == NULL)
  890. return (STRTOG_NoMemory);
  891. j = kmask & (ULbits - (rvbits & kmask));
  892. if (hi0bits(rvb->x[rvb->wds - 1]) != j)
  893. rvbits++;
  894. irv = STRTOG_Normal | STRTOG_Inexhi;
  895. }
  896. else {
  897. if (bbbits == 1)
  898. goto undfl;
  899. decrement(rvb);
  900. irv = STRTOG_Normal | STRTOG_Inexlo;
  901. }
  902. break;
  903. }
  904. if ((dval(&adj) = ratio(delta, bs)) <= 2.) {
  905. adj1:
  906. inex = STRTOG_Inexlo;
  907. if (dsign) {
  908. asub = 0;
  909. inex = STRTOG_Inexhi;
  910. }
  911. else if (denorm && bbbits <= 1) {
  912. undfl:
  913. rvb->wds = 0;
  914. rve = emin;
  915. irv = STRTOG_Underflow | STRTOG_Inexlo;
  916. break;
  917. }
  918. adj0 = dval(&adj) = 1.;
  919. }
  920. else {
  921. adj0 = dval(&adj) *= 0.5;
  922. if (dsign) {
  923. asub = 0;
  924. inex = STRTOG_Inexlo;
  925. }
  926. if (dval(&adj) < 2147483647.) {
  927. L = adj0;
  928. adj0 -= L;
  929. switch(rd) {
  930. case 0:
  931. if (adj0 >= .5)
  932. goto inc_L;
  933. break;
  934. case 1:
  935. if (asub && adj0 > 0.)
  936. goto inc_L;
  937. break;
  938. case 2:
  939. if (!asub && adj0 > 0.) {
  940. inc_L:
  941. L++;
  942. inex = STRTOG_Inexact - inex;
  943. }
  944. }
  945. dval(&adj) = L;
  946. }
  947. }
  948. y = rve + rvbits;
  949. /* adj *= ulp(dval(&rv)); */
  950. /* if (asub) rv -= adj; else rv += adj; */
  951. if (!denorm && rvbits < nbits) {
  952. rvb = lshift(rvb, j = nbits - rvbits);
  953. if (rvb == NULL)
  954. return (STRTOG_NoMemory);
  955. rve -= j;
  956. rvbits = nbits;
  957. }
  958. ab = d2b(dval(&adj), &abe, &abits);
  959. if (ab == NULL)
  960. return (STRTOG_NoMemory);
  961. if (abe < 0)
  962. rshift(ab, -abe);
  963. else if (abe > 0) {
  964. ab = lshift(ab, abe);
  965. if (ab == NULL)
  966. return (STRTOG_NoMemory);
  967. }
  968. rvb0 = rvb;
  969. if (asub) {
  970. /* rv -= adj; */
  971. j = hi0bits(rvb->x[rvb->wds-1]);
  972. rvb = diff(rvb, ab);
  973. if (rvb == NULL)
  974. return (STRTOG_NoMemory);
  975. k = rvb0->wds - 1;
  976. if (denorm)
  977. /* do nothing */;
  978. else if (rvb->wds <= k
  979. || hi0bits( rvb->x[k]) >
  980. hi0bits(rvb0->x[k])) {
  981. /* unlikely; can only have lost 1 high bit */
  982. if (rve1 == emin) {
  983. --rvbits;
  984. denorm = 1;
  985. }
  986. else {
  987. rvb = lshift(rvb, 1);
  988. if (rvb == NULL)
  989. return (STRTOG_NoMemory);
  990. --rve;
  991. --rve1;
  992. L = finished = 0;
  993. }
  994. }
  995. }
  996. else {
  997. rvb = sum(rvb, ab);
  998. if (rvb == NULL)
  999. return (STRTOG_NoMemory);
  1000. k = rvb->wds - 1;
  1001. if (k >= rvb0->wds
  1002. || hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) {
  1003. if (denorm) {
  1004. if (++rvbits == nbits)
  1005. denorm = 0;
  1006. }
  1007. else {
  1008. rshift(rvb, 1);
  1009. rve++;
  1010. rve1++;
  1011. L = 0;
  1012. }
  1013. }
  1014. }
  1015. Bfree(ab);
  1016. Bfree(rvb0);
  1017. if (finished)
  1018. break;
  1019. z = rve + rvbits;
  1020. if (y == z && L) {
  1021. /* Can we stop now? */
  1022. tol = dval(&adj) * 5e-16; /* > max rel error */
  1023. dval(&adj) = adj0 - .5;
  1024. if (dval(&adj) < -tol) {
  1025. if (adj0 > tol) {
  1026. irv |= inex;
  1027. break;
  1028. }
  1029. }
  1030. else if (dval(&adj) > tol && adj0 < 1. - tol) {
  1031. irv |= inex;
  1032. break;
  1033. }
  1034. }
  1035. bb0 = denorm ? 0 : trailz(rvb);
  1036. Bfree(bb);
  1037. Bfree(bd);
  1038. Bfree(bs);
  1039. Bfree(delta);
  1040. }
  1041. if (!denorm && (j = nbits - rvbits)) {
  1042. if (j > 0) {
  1043. rvb = lshift(rvb, j);
  1044. if (rvb == NULL)
  1045. return (STRTOG_NoMemory);
  1046. }
  1047. else
  1048. rshift(rvb, -j);
  1049. rve -= j;
  1050. }
  1051. *exp = rve;
  1052. Bfree(bb);
  1053. Bfree(bd);
  1054. Bfree(bs);
  1055. Bfree(bd0);
  1056. Bfree(delta);
  1057. if (rve > fpi->emax) {
  1058. switch(fpi->rounding & 3) {
  1059. case FPI_Round_near:
  1060. goto huge;
  1061. case FPI_Round_up:
  1062. if (!sign)
  1063. goto huge;
  1064. break;
  1065. case FPI_Round_down:
  1066. if (sign)
  1067. goto huge;
  1068. }
  1069. /* Round to largest representable magnitude */
  1070. Bfree(rvb);
  1071. rvb = 0;
  1072. irv = STRTOG_Normal | STRTOG_Inexlo;
  1073. *exp = fpi->emax;
  1074. b = bits;
  1075. be = b + ((fpi->nbits + 31) >> 5);
  1076. while(b < be)
  1077. *b++ = -1;
  1078. if ((j = fpi->nbits & 0x1f))
  1079. *--be >>= (32 - j);
  1080. goto ret;
  1081. huge:
  1082. rvb->wds = 0;
  1083. irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
  1084. #ifndef NO_ERRNO
  1085. errno = ERANGE;
  1086. #endif
  1087. infnanexp:
  1088. *exp = fpi->emax + 1;
  1089. }
  1090. ret:
  1091. if (denorm) {
  1092. if (sudden_underflow) {
  1093. if (!rvb) abort();
  1094. rvb->wds = 0;
  1095. irv = STRTOG_Underflow | STRTOG_Inexlo;
  1096. #ifndef NO_ERRNO
  1097. errno = ERANGE;
  1098. #endif
  1099. }
  1100. else {
  1101. if (!rvb) abort();
  1102. irv = (irv & ~STRTOG_Retmask) |
  1103. (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero);
  1104. if (irv & STRTOG_Inexact) {
  1105. irv |= STRTOG_Underflow;
  1106. #ifndef NO_ERRNO
  1107. errno = ERANGE;
  1108. #endif
  1109. }
  1110. }
  1111. }
  1112. if (se)
  1113. *se = (char *)s;
  1114. if (sign)
  1115. irv |= STRTOG_Neg;
  1116. if (rvb) {
  1117. copybits(bits, nbits, rvb);
  1118. Bfree(rvb);
  1119. }
  1120. return irv;
  1121. }