misc.c 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895
  1. /****************************************************************
  2. The author of this software is David M. Gay.
  3. Copyright (C) 1998, 1999 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. static Bigint *freelist[Kmax+1];
  27. #ifndef Omit_Private_Memory
  28. #ifndef PRIVATE_MEM
  29. #define PRIVATE_MEM 2304
  30. #endif
  31. #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
  32. static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
  33. #endif
  34. Bigint *
  35. Balloc
  36. #ifdef KR_headers
  37. (k) int k;
  38. #else
  39. (int k)
  40. #endif
  41. {
  42. int x;
  43. Bigint *rv;
  44. #ifndef Omit_Private_Memory
  45. unsigned int len;
  46. #endif
  47. ACQUIRE_DTOA_LOCK(0);
  48. /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
  49. /* but this case seems very unlikely. */
  50. if (k <= Kmax && (rv = freelist[k]) !=0) {
  51. freelist[k] = rv->next;
  52. }
  53. else {
  54. x = 1 << k;
  55. #ifdef Omit_Private_Memory
  56. rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
  57. if (rv == NULL)
  58. return (NULL);
  59. #else
  60. len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
  61. /sizeof(double);
  62. if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
  63. rv = (Bigint*)pmem_next;
  64. pmem_next += len;
  65. }
  66. else {
  67. rv = (Bigint*)MALLOC(len*sizeof(double));
  68. if (rv == NULL)
  69. return (NULL);
  70. }
  71. #endif
  72. rv->k = k;
  73. rv->maxwds = x;
  74. }
  75. FREE_DTOA_LOCK(0);
  76. rv->sign = rv->wds = 0;
  77. return rv;
  78. }
  79. void
  80. Bfree
  81. #ifdef KR_headers
  82. (v) Bigint *v;
  83. #else
  84. (Bigint *v)
  85. #endif
  86. {
  87. if (v) {
  88. if (v->k > Kmax)
  89. #ifdef FREE
  90. FREE((void*)v);
  91. #else
  92. free((void*)v);
  93. #endif
  94. else {
  95. ACQUIRE_DTOA_LOCK(0);
  96. v->next = freelist[v->k];
  97. freelist[v->k] = v;
  98. FREE_DTOA_LOCK(0);
  99. }
  100. }
  101. }
  102. int
  103. lo0bits
  104. #ifdef KR_headers
  105. (y) ULong *y;
  106. #else
  107. (ULong *y)
  108. #endif
  109. {
  110. int k;
  111. ULong x = *y;
  112. if (x & 7) {
  113. if (x & 1)
  114. return 0;
  115. if (x & 2) {
  116. *y = x >> 1;
  117. return 1;
  118. }
  119. *y = x >> 2;
  120. return 2;
  121. }
  122. k = 0;
  123. if (!(x & 0xffff)) {
  124. k = 16;
  125. x >>= 16;
  126. }
  127. if (!(x & 0xff)) {
  128. k += 8;
  129. x >>= 8;
  130. }
  131. if (!(x & 0xf)) {
  132. k += 4;
  133. x >>= 4;
  134. }
  135. if (!(x & 0x3)) {
  136. k += 2;
  137. x >>= 2;
  138. }
  139. if (!(x & 1)) {
  140. k++;
  141. x >>= 1;
  142. if (!x)
  143. return 32;
  144. }
  145. *y = x;
  146. return k;
  147. }
  148. Bigint *
  149. multadd
  150. #ifdef KR_headers
  151. (b, m, a) Bigint *b; int m, a;
  152. #else
  153. (Bigint *b, int m, int a) /* multiply by m and add a */
  154. #endif
  155. {
  156. int i, wds;
  157. #ifdef ULLong
  158. ULong *x;
  159. ULLong carry, y;
  160. #else
  161. ULong carry, *x, y;
  162. #ifdef Pack_32
  163. ULong xi, z;
  164. #endif
  165. #endif
  166. Bigint *b1;
  167. wds = b->wds;
  168. x = b->x;
  169. i = 0;
  170. carry = a;
  171. do {
  172. #ifdef ULLong
  173. y = *x * (ULLong)m + carry;
  174. carry = y >> 32;
  175. *x++ = y & 0xffffffffUL;
  176. #else
  177. #ifdef Pack_32
  178. xi = *x;
  179. y = (xi & 0xffff) * m + carry;
  180. z = (xi >> 16) * m + (y >> 16);
  181. carry = z >> 16;
  182. *x++ = (z << 16) + (y & 0xffff);
  183. #else
  184. y = *x * m + carry;
  185. carry = y >> 16;
  186. *x++ = y & 0xffff;
  187. #endif
  188. #endif
  189. }
  190. while(++i < wds);
  191. if (carry) {
  192. if (wds >= b->maxwds) {
  193. b1 = Balloc(b->k+1);
  194. if (b1 == NULL)
  195. return (NULL);
  196. Bcopy(b1, b);
  197. Bfree(b);
  198. b = b1;
  199. }
  200. b->x[wds++] = carry;
  201. b->wds = wds;
  202. }
  203. return b;
  204. }
  205. int
  206. hi0bits_D2A
  207. #ifdef KR_headers
  208. (x) ULong x;
  209. #else
  210. (ULong x)
  211. #endif
  212. {
  213. int k = 0;
  214. if (!(x & 0xffff0000)) {
  215. k = 16;
  216. x <<= 16;
  217. }
  218. if (!(x & 0xff000000)) {
  219. k += 8;
  220. x <<= 8;
  221. }
  222. if (!(x & 0xf0000000)) {
  223. k += 4;
  224. x <<= 4;
  225. }
  226. if (!(x & 0xc0000000)) {
  227. k += 2;
  228. x <<= 2;
  229. }
  230. if (!(x & 0x80000000)) {
  231. k++;
  232. if (!(x & 0x40000000))
  233. return 32;
  234. }
  235. return k;
  236. }
  237. Bigint *
  238. i2b
  239. #ifdef KR_headers
  240. (i) int i;
  241. #else
  242. (int i)
  243. #endif
  244. {
  245. Bigint *b;
  246. b = Balloc(1);
  247. if (b == NULL)
  248. return (NULL);
  249. b->x[0] = i;
  250. b->wds = 1;
  251. return b;
  252. }
  253. Bigint *
  254. mult
  255. #ifdef KR_headers
  256. (a, b) Bigint *a, *b;
  257. #else
  258. (Bigint *a, Bigint *b)
  259. #endif
  260. {
  261. Bigint *c;
  262. int k, wa, wb, wc;
  263. ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
  264. ULong y;
  265. #ifdef ULLong
  266. ULLong carry, z;
  267. #else
  268. ULong carry, z;
  269. #ifdef Pack_32
  270. ULong z2;
  271. #endif
  272. #endif
  273. if (a->wds < b->wds) {
  274. c = a;
  275. a = b;
  276. b = c;
  277. }
  278. k = a->k;
  279. wa = a->wds;
  280. wb = b->wds;
  281. wc = wa + wb;
  282. if (wc > a->maxwds)
  283. k++;
  284. c = Balloc(k);
  285. if (c == NULL)
  286. return (NULL);
  287. for(x = c->x, xa = x + wc; x < xa; x++)
  288. *x = 0;
  289. xa = a->x;
  290. xae = xa + wa;
  291. xb = b->x;
  292. xbe = xb + wb;
  293. xc0 = c->x;
  294. #ifdef ULLong
  295. for(; xb < xbe; xc0++) {
  296. if ( (y = *xb++) !=0) {
  297. x = xa;
  298. xc = xc0;
  299. carry = 0;
  300. do {
  301. z = *x++ * (ULLong)y + *xc + carry;
  302. carry = z >> 32;
  303. *xc++ = z & 0xffffffffUL;
  304. }
  305. while(x < xae);
  306. *xc = carry;
  307. }
  308. }
  309. #else
  310. #ifdef Pack_32
  311. for(; xb < xbe; xb++, xc0++) {
  312. if ( (y = *xb & 0xffff) !=0) {
  313. x = xa;
  314. xc = xc0;
  315. carry = 0;
  316. do {
  317. z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
  318. carry = z >> 16;
  319. z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
  320. carry = z2 >> 16;
  321. Storeinc(xc, z2, z);
  322. }
  323. while(x < xae);
  324. *xc = carry;
  325. }
  326. if ( (y = *xb >> 16) !=0) {
  327. x = xa;
  328. xc = xc0;
  329. carry = 0;
  330. z2 = *xc;
  331. do {
  332. z = (*x & 0xffff) * y + (*xc >> 16) + carry;
  333. carry = z >> 16;
  334. Storeinc(xc, z, z2);
  335. z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
  336. carry = z2 >> 16;
  337. }
  338. while(x < xae);
  339. *xc = z2;
  340. }
  341. }
  342. #else
  343. for(; xb < xbe; xc0++) {
  344. if ( (y = *xb++) !=0) {
  345. x = xa;
  346. xc = xc0;
  347. carry = 0;
  348. do {
  349. z = *x++ * y + *xc + carry;
  350. carry = z >> 16;
  351. *xc++ = z & 0xffff;
  352. }
  353. while(x < xae);
  354. *xc = carry;
  355. }
  356. }
  357. #endif
  358. #endif
  359. for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
  360. c->wds = wc;
  361. return c;
  362. }
  363. static Bigint *p5s;
  364. Bigint *
  365. pow5mult
  366. #ifdef KR_headers
  367. (b, k) Bigint *b; int k;
  368. #else
  369. (Bigint *b, int k)
  370. #endif
  371. {
  372. Bigint *b1, *p5, *p51;
  373. int i;
  374. static int p05[3] = { 5, 25, 125 };
  375. if ( (i = k & 3) !=0) {
  376. b = multadd(b, p05[i-1], 0);
  377. if (b == NULL)
  378. return (NULL);
  379. }
  380. if (!(k >>= 2))
  381. return b;
  382. if ((p5 = p5s) == 0) {
  383. /* first time */
  384. #ifdef MULTIPLE_THREADS
  385. ACQUIRE_DTOA_LOCK(1);
  386. if (!(p5 = p5s)) {
  387. p5 = p5s = i2b(625);
  388. if (p5 == NULL)
  389. return (NULL);
  390. p5->next = 0;
  391. }
  392. FREE_DTOA_LOCK(1);
  393. #else
  394. p5 = p5s = i2b(625);
  395. if (p5 == NULL)
  396. return (NULL);
  397. p5->next = 0;
  398. #endif
  399. }
  400. for(;;) {
  401. if (k & 1) {
  402. b1 = mult(b, p5);
  403. if (b1 == NULL)
  404. return (NULL);
  405. Bfree(b);
  406. b = b1;
  407. }
  408. if (!(k >>= 1))
  409. break;
  410. if ((p51 = p5->next) == 0) {
  411. #ifdef MULTIPLE_THREADS
  412. ACQUIRE_DTOA_LOCK(1);
  413. if (!(p51 = p5->next)) {
  414. p51 = p5->next = mult(p5,p5);
  415. if (p51 == NULL)
  416. return (NULL);
  417. p51->next = 0;
  418. }
  419. FREE_DTOA_LOCK(1);
  420. #else
  421. p51 = p5->next = mult(p5,p5);
  422. if (p51 == NULL)
  423. return (NULL);
  424. p51->next = 0;
  425. #endif
  426. }
  427. p5 = p51;
  428. }
  429. return b;
  430. }
  431. Bigint *
  432. lshift
  433. #ifdef KR_headers
  434. (b, k) Bigint *b; int k;
  435. #else
  436. (Bigint *b, int k)
  437. #endif
  438. {
  439. int i, k1, n, n1;
  440. Bigint *b1;
  441. ULong *x, *x1, *xe, z;
  442. n = k >> kshift;
  443. k1 = b->k;
  444. n1 = n + b->wds + 1;
  445. for(i = b->maxwds; n1 > i; i <<= 1)
  446. k1++;
  447. b1 = Balloc(k1);
  448. if (b1 == NULL)
  449. return (NULL);
  450. x1 = b1->x;
  451. for(i = 0; i < n; i++)
  452. *x1++ = 0;
  453. x = b->x;
  454. xe = x + b->wds;
  455. if (k &= kmask) {
  456. #ifdef Pack_32
  457. k1 = 32 - k;
  458. z = 0;
  459. do {
  460. *x1++ = *x << k | z;
  461. z = *x++ >> k1;
  462. }
  463. while(x < xe);
  464. if ((*x1 = z) !=0)
  465. ++n1;
  466. #else
  467. k1 = 16 - k;
  468. z = 0;
  469. do {
  470. *x1++ = *x << k & 0xffff | z;
  471. z = *x++ >> k1;
  472. }
  473. while(x < xe);
  474. if (*x1 = z)
  475. ++n1;
  476. #endif
  477. }
  478. else do
  479. *x1++ = *x++;
  480. while(x < xe);
  481. b1->wds = n1 - 1;
  482. Bfree(b);
  483. return b1;
  484. }
  485. int
  486. cmp
  487. #ifdef KR_headers
  488. (a, b) Bigint *a, *b;
  489. #else
  490. (Bigint *a, Bigint *b)
  491. #endif
  492. {
  493. ULong *xa, *xa0, *xb, *xb0;
  494. int i, j;
  495. i = a->wds;
  496. j = b->wds;
  497. #ifdef DEBUG
  498. if (i > 1 && !a->x[i-1])
  499. Bug("cmp called with a->x[a->wds-1] == 0");
  500. if (j > 1 && !b->x[j-1])
  501. Bug("cmp called with b->x[b->wds-1] == 0");
  502. #endif
  503. if (i -= j)
  504. return i;
  505. xa0 = a->x;
  506. xa = xa0 + j;
  507. xb0 = b->x;
  508. xb = xb0 + j;
  509. for(;;) {
  510. if (*--xa != *--xb)
  511. return *xa < *xb ? -1 : 1;
  512. if (xa <= xa0)
  513. break;
  514. }
  515. return 0;
  516. }
  517. Bigint *
  518. diff
  519. #ifdef KR_headers
  520. (a, b) Bigint *a, *b;
  521. #else
  522. (Bigint *a, Bigint *b)
  523. #endif
  524. {
  525. Bigint *c;
  526. int i, wa, wb;
  527. ULong *xa, *xae, *xb, *xbe, *xc;
  528. #ifdef ULLong
  529. ULLong borrow, y;
  530. #else
  531. ULong borrow, y;
  532. #ifdef Pack_32
  533. ULong z;
  534. #endif
  535. #endif
  536. i = cmp(a,b);
  537. if (!i) {
  538. c = Balloc(0);
  539. if (c == NULL)
  540. return (NULL);
  541. c->wds = 1;
  542. c->x[0] = 0;
  543. return c;
  544. }
  545. if (i < 0) {
  546. c = a;
  547. a = b;
  548. b = c;
  549. i = 1;
  550. }
  551. else
  552. i = 0;
  553. c = Balloc(a->k);
  554. if (c == NULL)
  555. return (NULL);
  556. c->sign = i;
  557. wa = a->wds;
  558. xa = a->x;
  559. xae = xa + wa;
  560. wb = b->wds;
  561. xb = b->x;
  562. xbe = xb + wb;
  563. xc = c->x;
  564. borrow = 0;
  565. #ifdef ULLong
  566. do {
  567. y = (ULLong)*xa++ - *xb++ - borrow;
  568. borrow = y >> 32 & 1UL;
  569. *xc++ = y & 0xffffffffUL;
  570. }
  571. while(xb < xbe);
  572. while(xa < xae) {
  573. y = *xa++ - borrow;
  574. borrow = y >> 32 & 1UL;
  575. *xc++ = y & 0xffffffffUL;
  576. }
  577. #else
  578. #ifdef Pack_32
  579. do {
  580. y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
  581. borrow = (y & 0x10000) >> 16;
  582. z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
  583. borrow = (z & 0x10000) >> 16;
  584. Storeinc(xc, z, y);
  585. }
  586. while(xb < xbe);
  587. while(xa < xae) {
  588. y = (*xa & 0xffff) - borrow;
  589. borrow = (y & 0x10000) >> 16;
  590. z = (*xa++ >> 16) - borrow;
  591. borrow = (z & 0x10000) >> 16;
  592. Storeinc(xc, z, y);
  593. }
  594. #else
  595. do {
  596. y = *xa++ - *xb++ - borrow;
  597. borrow = (y & 0x10000) >> 16;
  598. *xc++ = y & 0xffff;
  599. }
  600. while(xb < xbe);
  601. while(xa < xae) {
  602. y = *xa++ - borrow;
  603. borrow = (y & 0x10000) >> 16;
  604. *xc++ = y & 0xffff;
  605. }
  606. #endif
  607. #endif
  608. while(!*--xc)
  609. wa--;
  610. c->wds = wa;
  611. return c;
  612. }
  613. double
  614. b2d
  615. #ifdef KR_headers
  616. (a, e) Bigint *a; int *e;
  617. #else
  618. (Bigint *a, int *e)
  619. #endif
  620. {
  621. ULong *xa, *xa0, w, y, z;
  622. int k;
  623. U d;
  624. #ifdef VAX
  625. ULong d0, d1;
  626. #else
  627. #define d0 word0(&d)
  628. #define d1 word1(&d)
  629. #endif
  630. xa0 = a->x;
  631. xa = xa0 + a->wds;
  632. y = *--xa;
  633. #ifdef DEBUG
  634. if (!y) Bug("zero y in b2d");
  635. #endif
  636. k = hi0bits(y);
  637. *e = 32 - k;
  638. #ifdef Pack_32
  639. if (k < Ebits) {
  640. d0 = Exp_1 | y >> (Ebits - k);
  641. w = xa > xa0 ? *--xa : 0;
  642. d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
  643. goto ret_d;
  644. }
  645. z = xa > xa0 ? *--xa : 0;
  646. if (k -= Ebits) {
  647. d0 = Exp_1 | y << k | z >> (32 - k);
  648. y = xa > xa0 ? *--xa : 0;
  649. d1 = z << k | y >> (32 - k);
  650. }
  651. else {
  652. d0 = Exp_1 | y;
  653. d1 = z;
  654. }
  655. #else
  656. if (k < Ebits + 16) {
  657. z = xa > xa0 ? *--xa : 0;
  658. d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
  659. w = xa > xa0 ? *--xa : 0;
  660. y = xa > xa0 ? *--xa : 0;
  661. d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
  662. goto ret_d;
  663. }
  664. z = xa > xa0 ? *--xa : 0;
  665. w = xa > xa0 ? *--xa : 0;
  666. k -= Ebits + 16;
  667. d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
  668. y = xa > xa0 ? *--xa : 0;
  669. d1 = w << k + 16 | y << k;
  670. #endif
  671. ret_d:
  672. #ifdef VAX
  673. word0(&d) = d0 >> 16 | d0 << 16;
  674. word1(&d) = d1 >> 16 | d1 << 16;
  675. #endif
  676. return dval(&d);
  677. }
  678. #undef d0
  679. #undef d1
  680. Bigint *
  681. d2b
  682. #ifdef KR_headers
  683. (dd, e, bits) double dd; int *e, *bits;
  684. #else
  685. (double dd, int *e, int *bits)
  686. #endif
  687. {
  688. Bigint *b;
  689. U d;
  690. #ifndef Sudden_Underflow
  691. int i;
  692. #endif
  693. int de, k;
  694. ULong *x, y, z;
  695. #ifdef VAX
  696. ULong d0, d1;
  697. #else
  698. #define d0 word0(&d)
  699. #define d1 word1(&d)
  700. #endif
  701. d.d = dd;
  702. #ifdef VAX
  703. d0 = word0(&d) >> 16 | word0(&d) << 16;
  704. d1 = word1(&d) >> 16 | word1(&d) << 16;
  705. #endif
  706. #ifdef Pack_32
  707. b = Balloc(1);
  708. #else
  709. b = Balloc(2);
  710. #endif
  711. if (b == NULL)
  712. return (NULL);
  713. x = b->x;
  714. z = d0 & Frac_mask;
  715. d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
  716. #ifdef Sudden_Underflow
  717. de = (int)(d0 >> Exp_shift);
  718. #ifndef IBM
  719. z |= Exp_msk11;
  720. #endif
  721. #else
  722. if ( (de = (int)(d0 >> Exp_shift)) !=0)
  723. z |= Exp_msk1;
  724. #endif
  725. #ifdef Pack_32
  726. if ( (y = d1) !=0) {
  727. if ( (k = lo0bits(&y)) !=0) {
  728. x[0] = y | z << (32 - k);
  729. z >>= k;
  730. }
  731. else
  732. x[0] = y;
  733. #ifndef Sudden_Underflow
  734. i =
  735. #endif
  736. b->wds = (x[1] = z) !=0 ? 2 : 1;
  737. }
  738. else {
  739. k = lo0bits(&z);
  740. x[0] = z;
  741. #ifndef Sudden_Underflow
  742. i =
  743. #endif
  744. b->wds = 1;
  745. k += 32;
  746. }
  747. #else
  748. if ( (y = d1) !=0) {
  749. if ( (k = lo0bits(&y)) !=0)
  750. if (k >= 16) {
  751. x[0] = y | z << 32 - k & 0xffff;
  752. x[1] = z >> k - 16 & 0xffff;
  753. x[2] = z >> k;
  754. i = 2;
  755. }
  756. else {
  757. x[0] = y & 0xffff;
  758. x[1] = y >> 16 | z << 16 - k & 0xffff;
  759. x[2] = z >> k & 0xffff;
  760. x[3] = z >> k+16;
  761. i = 3;
  762. }
  763. else {
  764. x[0] = y & 0xffff;
  765. x[1] = y >> 16;
  766. x[2] = z & 0xffff;
  767. x[3] = z >> 16;
  768. i = 3;
  769. }
  770. }
  771. else {
  772. #ifdef DEBUG
  773. if (!z)
  774. Bug("Zero passed to d2b");
  775. #endif
  776. k = lo0bits(&z);
  777. if (k >= 16) {
  778. x[0] = z;
  779. i = 0;
  780. }
  781. else {
  782. x[0] = z & 0xffff;
  783. x[1] = z >> 16;
  784. i = 1;
  785. }
  786. k += 32;
  787. }
  788. while(!x[i])
  789. --i;
  790. b->wds = i + 1;
  791. #endif
  792. #ifndef Sudden_Underflow
  793. if (de) {
  794. #endif
  795. #ifdef IBM
  796. *e = (de - Bias - (P-1) << 2) + k;
  797. *bits = 4*P + 8 - k - hi0bits(word0(&d) & Frac_mask);
  798. #else
  799. *e = de - Bias - (P-1) + k;
  800. *bits = P - k;
  801. #endif
  802. #ifndef Sudden_Underflow
  803. }
  804. else {
  805. *e = de - Bias - (P-1) + 1 + k;
  806. #ifdef Pack_32
  807. *bits = 32*i - hi0bits(x[i-1]);
  808. #else
  809. *bits = (i+2)*16 - hi0bits(x[i]);
  810. #endif
  811. }
  812. #endif
  813. return b;
  814. }
  815. #undef d0
  816. #undef d1
  817. CONST double
  818. #ifdef IEEE_Arith
  819. bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
  820. CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256
  821. };
  822. #else
  823. #ifdef IBM
  824. bigtens[] = { 1e16, 1e32, 1e64 };
  825. CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
  826. #else
  827. bigtens[] = { 1e16, 1e32 };
  828. CONST double tinytens[] = { 1e-16, 1e-32 };
  829. #endif
  830. #endif
  831. CONST double
  832. tens[] = {
  833. 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
  834. 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
  835. 1e20, 1e21, 1e22
  836. #ifdef VAX
  837. , 1e23, 1e24
  838. #endif
  839. };
  840. #ifdef NO_STRING_H
  841. Char *
  842. #ifdef KR_headers
  843. memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
  844. #else
  845. memcpy_D2A(void *a1, void *b1, size_t len)
  846. #endif
  847. {
  848. char *a = (char*)a1, *ae = a + len;
  849. char *b = (char*)b1, *a0 = a;
  850. while(a < ae)
  851. *a++ = *b++;
  852. return a0;
  853. }
  854. #endif /* NO_STRING_H */