summaryrefslogtreecommitdiff
path: root/src/lib/libc/stdlib/strtod.c
diff options
context:
space:
mode:
authorderaadt <>1995-10-18 08:42:23 +0000
committerderaadt <>1995-10-18 08:42:23 +0000
commit0527d29da443886d92e9a418180c5b25a5f8d270 (patch)
tree86b3a64928451a669cefa27900e5884036b4e349 /src/lib/libc/stdlib/strtod.c
downloadopenbsd-0527d29da443886d92e9a418180c5b25a5f8d270.tar.gz
openbsd-0527d29da443886d92e9a418180c5b25a5f8d270.tar.bz2
openbsd-0527d29da443886d92e9a418180c5b25a5f8d270.zip
initial import of NetBSD tree
Diffstat (limited to 'src/lib/libc/stdlib/strtod.c')
-rw-r--r--src/lib/libc/stdlib/strtod.c2499
1 files changed, 2499 insertions, 0 deletions
diff --git a/src/lib/libc/stdlib/strtod.c b/src/lib/libc/stdlib/strtod.c
new file mode 100644
index 0000000000..b13fa128f5
--- /dev/null
+++ b/src/lib/libc/stdlib/strtod.c
@@ -0,0 +1,2499 @@
1/****************************************************************
2 *
3 * The author of this software is David M. Gay.
4 *
5 * Copyright (c) 1991 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
20/* Please send bug reports to
21 David M. Gay
22 AT&T Bell Laboratories, Room 2C-463
23 600 Mountain Avenue
24 Murray Hill, NJ 07974-2070
25 U.S.A.
26 dmg@research.att.com or research!dmg
27 */
28
29/* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
30 *
31 * This strtod returns a nearest machine number to the input decimal
32 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
33 * broken by the IEEE round-even rule. Otherwise ties are broken by
34 * biased rounding (add half and chop).
35 *
36 * Inspired loosely by William D. Clinger's paper "How to Read Floating
37 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
38 *
39 * Modifications:
40 *
41 * 1. We only require IEEE, IBM, or VAX double-precision
42 * arithmetic (not IEEE double-extended).
43 * 2. We get by with floating-point arithmetic in a case that
44 * Clinger missed -- when we're computing d * 10^n
45 * for a small integer d and the integer n is not too
46 * much larger than 22 (the maximum integer k for which
47 * we can represent 10^k exactly), we may be able to
48 * compute (d*10^k) * 10^(e-k) with just one roundoff.
49 * 3. Rather than a bit-at-a-time adjustment of the binary
50 * result in the hard case, we use floating-point
51 * arithmetic to determine the adjustment to within
52 * one bit; only in really hard cases do we need to
53 * compute a second residual.
54 * 4. Because of 3., we don't need a large table of powers of 10
55 * for ten-to-e (just some small tables, e.g. of 10^k
56 * for 0 <= k <= 22).
57 */
58
59/*
60 * #define IEEE_LITTLE_ENDIAN for IEEE-arithmetic machines where the least
61 * significant byte has the lowest address.
62 * #define IEEE_BIG_ENDIAN for IEEE-arithmetic machines where the most
63 * significant byte has the lowest address.
64 * #define Long int on machines with 32-bit ints and 64-bit longs.
65 * #define Sudden_Underflow for IEEE-format machines without gradual
66 * underflow (i.e., that flush to zero on underflow).
67 * #define IBM for IBM mainframe-style floating-point arithmetic.
68 * #define VAX for VAX-style floating-point arithmetic.
69 * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
70 * #define No_leftright to omit left-right logic in fast floating-point
71 * computation of dtoa.
72 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
73 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
74 * that use extended-precision instructions to compute rounded
75 * products and quotients) with IBM.
76 * #define ROUND_BIASED for IEEE-format with biased rounding.
77 * #define Inaccurate_Divide for IEEE-format with correctly rounded
78 * products but inaccurate quotients, e.g., for Intel i860.
79 * #define Just_16 to store 16 bits per 32-bit Long when doing high-precision
80 * integer arithmetic. Whether this speeds things up or slows things
81 * down depends on the machine and the number being converted.
82 * #define KR_headers for old-style C function headers.
83 * #define Bad_float_h if your system lacks a float.h or if it does not
84 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
85 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
86 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
87 * if memory is available and otherwise does something you deem
88 * appropriate. If MALLOC is undefined, malloc will be invoked
89 * directly -- and assumed always to succeed.
90 */
91
92#if defined(LIBC_SCCS) && !defined(lint)
93static char *rcsid = "$Id: strtod.c,v 1.1.1.1 1995/10/18 08:42:19 deraadt Exp $";
94#endif /* LIBC_SCCS and not lint */
95
96#if defined(__m68k__) || defined(__sparc__) || defined(__i386__) || \
97 defined(__mips__) || defined(__ns32k__) || defined(__alpha__)
98#include <machine/endian.h>
99#if BYTE_ORDER == BIG_ENDIAN
100#define IEEE_BIG_ENDIAN
101#else
102#define IEEE_LITTLE_ENDIAN
103#endif
104#endif
105
106#ifdef vax
107#define VAX
108#endif
109
110#define Long int32_t
111#define ULong u_int32_t
112
113#ifdef DEBUG
114#include "stdio.h"
115#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
116#endif
117
118#ifdef __cplusplus
119#include "malloc.h"
120#include "memory.h"
121#else
122#ifndef KR_headers
123#include "stdlib.h"
124#include "string.h"
125#include "locale.h"
126#else
127#include "malloc.h"
128#include "memory.h"
129#endif
130#endif
131
132#ifdef MALLOC
133#ifdef KR_headers
134extern char *MALLOC();
135#else
136extern void *MALLOC(size_t);
137#endif
138#else
139#define MALLOC malloc
140#endif
141
142#include "ctype.h"
143#include "errno.h"
144
145#ifdef Bad_float_h
146#undef __STDC__
147#ifdef IEEE_BIG_ENDIAN
148#define IEEE_ARITHMETIC
149#endif
150#ifdef IEEE_LITTLE_ENDIAN
151#define IEEE_ARITHMETIC
152#endif
153
154#ifdef IEEE_ARITHMETIC
155#define DBL_DIG 15
156#define DBL_MAX_10_EXP 308
157#define DBL_MAX_EXP 1024
158#define FLT_RADIX 2
159#define FLT_ROUNDS 1
160#define DBL_MAX 1.7976931348623157e+308
161#endif
162
163#ifdef IBM
164#define DBL_DIG 16
165#define DBL_MAX_10_EXP 75
166#define DBL_MAX_EXP 63
167#define FLT_RADIX 16
168#define FLT_ROUNDS 0
169#define DBL_MAX 7.2370055773322621e+75
170#endif
171
172#ifdef VAX
173#define DBL_DIG 16
174#define DBL_MAX_10_EXP 38
175#define DBL_MAX_EXP 127
176#define FLT_RADIX 2
177#define FLT_ROUNDS 1
178#define DBL_MAX 1.7014118346046923e+38
179#endif
180
181#ifndef LONG_MAX
182#define LONG_MAX 2147483647
183#endif
184#else
185#include "float.h"
186#endif
187#ifndef __MATH_H__
188#include "math.h"
189#endif
190
191#ifdef __cplusplus
192extern "C" {
193#endif
194
195#ifndef CONST
196#ifdef KR_headers
197#define CONST /* blank */
198#else
199#define CONST const
200#endif
201#endif
202
203#ifdef Unsigned_Shifts
204#define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
205#else
206#define Sign_Extend(a,b) /*no-op*/
207#endif
208
209#if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN) + defined(VAX) + \
210 defined(IBM) != 1
211Exactly one of IEEE_LITTLE_ENDIAN IEEE_BIG_ENDIAN, VAX, or
212IBM should be defined.
213#endif
214
215#ifdef IEEE_LITTLE_ENDIAN
216#define word0(x) ((ULong *)&x)[1]
217#define word1(x) ((ULong *)&x)[0]
218#else
219#define word0(x) ((ULong *)&x)[0]
220#define word1(x) ((ULong *)&x)[1]
221#endif
222
223/* The following definition of Storeinc is appropriate for MIPS processors.
224 * An alternative that might be better on some machines is
225 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
226 */
227#if defined(IEEE_LITTLE_ENDIAN) + defined(VAX)
228#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
229((unsigned short *)a)[0] = (unsigned short)c, a++)
230#else
231#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
232((unsigned short *)a)[1] = (unsigned short)c, a++)
233#endif
234
235/* #define P DBL_MANT_DIG */
236/* Ten_pmax = floor(P*log(2)/log(5)) */
237/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
238/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
239/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
240
241#if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN)
242#define Exp_shift 20
243#define Exp_shift1 20
244#define Exp_msk1 0x100000
245#define Exp_msk11 0x100000
246#define Exp_mask 0x7ff00000
247#define P 53
248#define Bias 1023
249#define IEEE_Arith
250#define Emin (-1022)
251#define Exp_1 0x3ff00000
252#define Exp_11 0x3ff00000
253#define Ebits 11
254#define Frac_mask 0xfffff
255#define Frac_mask1 0xfffff
256#define Ten_pmax 22
257#define Bletch 0x10
258#define Bndry_mask 0xfffff
259#define Bndry_mask1 0xfffff
260#define LSB 1
261#define Sign_bit 0x80000000
262#define Log2P 1
263#define Tiny0 0
264#define Tiny1 1
265#define Quick_max 14
266#define Int_max 14
267#define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
268#else
269#undef Sudden_Underflow
270#define Sudden_Underflow
271#ifdef IBM
272#define Exp_shift 24
273#define Exp_shift1 24
274#define Exp_msk1 0x1000000
275#define Exp_msk11 0x1000000
276#define Exp_mask 0x7f000000
277#define P 14
278#define Bias 65
279#define Exp_1 0x41000000
280#define Exp_11 0x41000000
281#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
282#define Frac_mask 0xffffff
283#define Frac_mask1 0xffffff
284#define Bletch 4
285#define Ten_pmax 22
286#define Bndry_mask 0xefffff
287#define Bndry_mask1 0xffffff
288#define LSB 1
289#define Sign_bit 0x80000000
290#define Log2P 4
291#define Tiny0 0x100000
292#define Tiny1 0
293#define Quick_max 14
294#define Int_max 15
295#else /* VAX */
296#define Exp_shift 23
297#define Exp_shift1 7
298#define Exp_msk1 0x80
299#define Exp_msk11 0x800000
300#define Exp_mask 0x7f80
301#define P 56
302#define Bias 129
303#define Exp_1 0x40800000
304#define Exp_11 0x4080
305#define Ebits 8
306#define Frac_mask 0x7fffff
307#define Frac_mask1 0xffff007f
308#define Ten_pmax 24
309#define Bletch 2
310#define Bndry_mask 0xffff007f
311#define Bndry_mask1 0xffff007f
312#define LSB 0x10000
313#define Sign_bit 0x8000
314#define Log2P 1
315#define Tiny0 0x80
316#define Tiny1 0
317#define Quick_max 15
318#define Int_max 15
319#endif
320#endif
321
322#ifndef IEEE_Arith
323#define ROUND_BIASED
324#endif
325
326#ifdef RND_PRODQUOT
327#define rounded_product(a,b) a = rnd_prod(a, b)
328#define rounded_quotient(a,b) a = rnd_quot(a, b)
329#ifdef KR_headers
330extern double rnd_prod(), rnd_quot();
331#else
332extern double rnd_prod(double, double), rnd_quot(double, double);
333#endif
334#else
335#define rounded_product(a,b) a *= b
336#define rounded_quotient(a,b) a /= b
337#endif
338
339#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
340#define Big1 0xffffffff
341
342#ifndef Just_16
343/* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
344 * This makes some inner loops simpler and sometimes saves work
345 * during multiplications, but it often seems to make things slightly
346 * slower. Hence the default is now to store 32 bits per Long.
347 */
348#ifndef Pack_32
349#define Pack_32
350#endif
351#endif
352
353#define Kmax 15
354
355#ifdef __cplusplus
356extern "C" double strtod(const char *s00, char **se);
357extern "C" char *__dtoa(double d, int mode, int ndigits,
358 int *decpt, int *sign, char **rve);
359#endif
360
361 struct
362Bigint {
363 struct Bigint *next;
364 int k, maxwds, sign, wds;
365 ULong x[1];
366 };
367
368 typedef struct Bigint Bigint;
369
370 static Bigint *freelist[Kmax+1];
371
372 static Bigint *
373Balloc
374#ifdef KR_headers
375 (k) int k;
376#else
377 (int k)
378#endif
379{
380 int x;
381 Bigint *rv;
382
383 if (rv = freelist[k]) {
384 freelist[k] = rv->next;
385 }
386 else {
387 x = 1 << k;
388 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(Long));
389 rv->k = k;
390 rv->maxwds = x;
391 }
392 rv->sign = rv->wds = 0;
393 return rv;
394 }
395
396 static void
397Bfree
398#ifdef KR_headers
399 (v) Bigint *v;
400#else
401 (Bigint *v)
402#endif
403{
404 if (v) {
405 v->next = freelist[v->k];
406 freelist[v->k] = v;
407 }
408 }
409
410#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
411y->wds*sizeof(Long) + 2*sizeof(int))
412
413 static Bigint *
414multadd
415#ifdef KR_headers
416 (b, m, a) Bigint *b; int m, a;
417#else
418 (Bigint *b, int m, int a) /* multiply by m and add a */
419#endif
420{
421 int i, wds;
422 ULong *x, y;
423#ifdef Pack_32
424 ULong xi, z;
425#endif
426 Bigint *b1;
427
428 wds = b->wds;
429 x = b->x;
430 i = 0;
431 do {
432#ifdef Pack_32
433 xi = *x;
434 y = (xi & 0xffff) * m + a;
435 z = (xi >> 16) * m + (y >> 16);
436 a = (int)(z >> 16);
437 *x++ = (z << 16) + (y & 0xffff);
438#else
439 y = *x * m + a;
440 a = (int)(y >> 16);
441 *x++ = y & 0xffff;
442#endif
443 }
444 while(++i < wds);
445 if (a) {
446 if (wds >= b->maxwds) {
447 b1 = Balloc(b->k+1);
448 Bcopy(b1, b);
449 Bfree(b);
450 b = b1;
451 }
452 b->x[wds++] = a;
453 b->wds = wds;
454 }
455 return b;
456 }
457
458 static Bigint *
459s2b
460#ifdef KR_headers
461 (s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
462#else
463 (CONST char *s, int nd0, int nd, ULong y9)
464#endif
465{
466 Bigint *b;
467 int i, k;
468 Long x, y;
469
470 x = (nd + 8) / 9;
471 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
472#ifdef Pack_32
473 b = Balloc(k);
474 b->x[0] = y9;
475 b->wds = 1;
476#else
477 b = Balloc(k+1);
478 b->x[0] = y9 & 0xffff;
479 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
480#endif
481
482 i = 9;
483 if (9 < nd0) {
484 s += 9;
485 do b = multadd(b, 10, *s++ - '0');
486 while(++i < nd0);
487 s++;
488 }
489 else
490 s += 10;
491 for(; i < nd; i++)
492 b = multadd(b, 10, *s++ - '0');
493 return b;
494 }
495
496 static int
497hi0bits
498#ifdef KR_headers
499 (x) register ULong x;
500#else
501 (register ULong x)
502#endif
503{
504 register int k = 0;
505
506 if (!(x & 0xffff0000)) {
507 k = 16;
508 x <<= 16;
509 }
510 if (!(x & 0xff000000)) {
511 k += 8;
512 x <<= 8;
513 }
514 if (!(x & 0xf0000000)) {
515 k += 4;
516 x <<= 4;
517 }
518 if (!(x & 0xc0000000)) {
519 k += 2;
520 x <<= 2;
521 }
522 if (!(x & 0x80000000)) {
523 k++;
524 if (!(x & 0x40000000))
525 return 32;
526 }
527 return k;
528 }
529
530 static int
531lo0bits
532#ifdef KR_headers
533 (y) ULong *y;
534#else
535 (ULong *y)
536#endif
537{
538 register int k;
539 register ULong x = *y;
540
541 if (x & 7) {
542 if (x & 1)
543 return 0;
544 if (x & 2) {
545 *y = x >> 1;
546 return 1;
547 }
548 *y = x >> 2;
549 return 2;
550 }
551 k = 0;
552 if (!(x & 0xffff)) {
553 k = 16;
554 x >>= 16;
555 }
556 if (!(x & 0xff)) {
557 k += 8;
558 x >>= 8;
559 }
560 if (!(x & 0xf)) {
561 k += 4;
562 x >>= 4;
563 }
564 if (!(x & 0x3)) {
565 k += 2;
566 x >>= 2;
567 }
568 if (!(x & 1)) {
569 k++;
570 x >>= 1;
571 if (!x & 1)
572 return 32;
573 }
574 *y = x;
575 return k;
576 }
577
578 static Bigint *
579i2b
580#ifdef KR_headers
581 (i) int i;
582#else
583 (int i)
584#endif
585{
586 Bigint *b;
587
588 b = Balloc(1);
589 b->x[0] = i;
590 b->wds = 1;
591 return b;
592 }
593
594 static Bigint *
595mult
596#ifdef KR_headers
597 (a, b) Bigint *a, *b;
598#else
599 (Bigint *a, Bigint *b)
600#endif
601{
602 Bigint *c;
603 int k, wa, wb, wc;
604 ULong carry, y, z;
605 ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
606#ifdef Pack_32
607 ULong z2;
608#endif
609
610 if (a->wds < b->wds) {
611 c = a;
612 a = b;
613 b = c;
614 }
615 k = a->k;
616 wa = a->wds;
617 wb = b->wds;
618 wc = wa + wb;
619 if (wc > a->maxwds)
620 k++;
621 c = Balloc(k);
622 for(x = c->x, xa = x + wc; x < xa; x++)
623 *x = 0;
624 xa = a->x;
625 xae = xa + wa;
626 xb = b->x;
627 xbe = xb + wb;
628 xc0 = c->x;
629#ifdef Pack_32
630 for(; xb < xbe; xb++, xc0++) {
631 if (y = *xb & 0xffff) {
632 x = xa;
633 xc = xc0;
634 carry = 0;
635 do {
636 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
637 carry = z >> 16;
638 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
639 carry = z2 >> 16;
640 Storeinc(xc, z2, z);
641 }
642 while(x < xae);
643 *xc = carry;
644 }
645 if (y = *xb >> 16) {
646 x = xa;
647 xc = xc0;
648 carry = 0;
649 z2 = *xc;
650 do {
651 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
652 carry = z >> 16;
653 Storeinc(xc, z, z2);
654 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
655 carry = z2 >> 16;
656 }
657 while(x < xae);
658 *xc = z2;
659 }
660 }
661#else
662 for(; xb < xbe; xc0++) {
663 if (y = *xb++) {
664 x = xa;
665 xc = xc0;
666 carry = 0;
667 do {
668 z = *x++ * y + *xc + carry;
669 carry = z >> 16;
670 *xc++ = z & 0xffff;
671 }
672 while(x < xae);
673 *xc = carry;
674 }
675 }
676#endif
677 for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
678 c->wds = wc;
679 return c;
680 }
681
682 static Bigint *p5s;
683
684 static Bigint *
685pow5mult
686#ifdef KR_headers
687 (b, k) Bigint *b; int k;
688#else
689 (Bigint *b, int k)
690#endif
691{
692 Bigint *b1, *p5, *p51;
693 int i;
694 static int p05[3] = { 5, 25, 125 };
695
696 if (i = k & 3)
697 b = multadd(b, p05[i-1], 0);
698
699 if (!(k >>= 2))
700 return b;
701 if (!(p5 = p5s)) {
702 /* first time */
703 p5 = p5s = i2b(625);
704 p5->next = 0;
705 }
706 for(;;) {
707 if (k & 1) {
708 b1 = mult(b, p5);
709 Bfree(b);
710 b = b1;
711 }
712 if (!(k >>= 1))
713 break;
714 if (!(p51 = p5->next)) {
715 p51 = p5->next = mult(p5,p5);
716 p51->next = 0;
717 }
718 p5 = p51;
719 }
720 return b;
721 }
722
723 static Bigint *
724lshift
725#ifdef KR_headers
726 (b, k) Bigint *b; int k;
727#else
728 (Bigint *b, int k)
729#endif
730{
731 int i, k1, n, n1;
732 Bigint *b1;
733 ULong *x, *x1, *xe, z;
734
735#ifdef Pack_32
736 n = k >> 5;
737#else
738 n = k >> 4;
739#endif
740 k1 = b->k;
741 n1 = n + b->wds + 1;
742 for(i = b->maxwds; n1 > i; i <<= 1)
743 k1++;
744 b1 = Balloc(k1);
745 x1 = b1->x;
746 for(i = 0; i < n; i++)
747 *x1++ = 0;
748 x = b->x;
749 xe = x + b->wds;
750#ifdef Pack_32
751 if (k &= 0x1f) {
752 k1 = 32 - k;
753 z = 0;
754 do {
755 *x1++ = *x << k | z;
756 z = *x++ >> k1;
757 }
758 while(x < xe);
759 if (*x1 = z)
760 ++n1;
761 }
762#else
763 if (k &= 0xf) {
764 k1 = 16 - k;
765 z = 0;
766 do {
767 *x1++ = *x << k & 0xffff | z;
768 z = *x++ >> k1;
769 }
770 while(x < xe);
771 if (*x1 = z)
772 ++n1;
773 }
774#endif
775 else do
776 *x1++ = *x++;
777 while(x < xe);
778 b1->wds = n1 - 1;
779 Bfree(b);
780 return b1;
781 }
782
783 static int
784cmp
785#ifdef KR_headers
786 (a, b) Bigint *a, *b;
787#else
788 (Bigint *a, Bigint *b)
789#endif
790{
791 ULong *xa, *xa0, *xb, *xb0;
792 int i, j;
793
794 i = a->wds;
795 j = b->wds;
796#ifdef DEBUG
797 if (i > 1 && !a->x[i-1])
798 Bug("cmp called with a->x[a->wds-1] == 0");
799 if (j > 1 && !b->x[j-1])
800 Bug("cmp called with b->x[b->wds-1] == 0");
801#endif
802 if (i -= j)
803 return i;
804 xa0 = a->x;
805 xa = xa0 + j;
806 xb0 = b->x;
807 xb = xb0 + j;
808 for(;;) {
809 if (*--xa != *--xb)
810 return *xa < *xb ? -1 : 1;
811 if (xa <= xa0)
812 break;
813 }
814 return 0;
815 }
816
817 static Bigint *
818diff
819#ifdef KR_headers
820 (a, b) Bigint *a, *b;
821#else
822 (Bigint *a, Bigint *b)
823#endif
824{
825 Bigint *c;
826 int i, wa, wb;
827 Long borrow, y; /* We need signed shifts here. */
828 ULong *xa, *xae, *xb, *xbe, *xc;
829#ifdef Pack_32
830 Long z;
831#endif
832
833 i = cmp(a,b);
834 if (!i) {
835 c = Balloc(0);
836 c->wds = 1;
837 c->x[0] = 0;
838 return c;
839 }
840 if (i < 0) {
841 c = a;
842 a = b;
843 b = c;
844 i = 1;
845 }
846 else
847 i = 0;
848 c = Balloc(a->k);
849 c->sign = i;
850 wa = a->wds;
851 xa = a->x;
852 xae = xa + wa;
853 wb = b->wds;
854 xb = b->x;
855 xbe = xb + wb;
856 xc = c->x;
857 borrow = 0;
858#ifdef Pack_32
859 do {
860 y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
861 borrow = y >> 16;
862 Sign_Extend(borrow, y);
863 z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
864 borrow = z >> 16;
865 Sign_Extend(borrow, z);
866 Storeinc(xc, z, y);
867 }
868 while(xb < xbe);
869 while(xa < xae) {
870 y = (*xa & 0xffff) + borrow;
871 borrow = y >> 16;
872 Sign_Extend(borrow, y);
873 z = (*xa++ >> 16) + borrow;
874 borrow = z >> 16;
875 Sign_Extend(borrow, z);
876 Storeinc(xc, z, y);
877 }
878#else
879 do {
880 y = *xa++ - *xb++ + borrow;
881 borrow = y >> 16;
882 Sign_Extend(borrow, y);
883 *xc++ = y & 0xffff;
884 }
885 while(xb < xbe);
886 while(xa < xae) {
887 y = *xa++ + borrow;
888 borrow = y >> 16;
889 Sign_Extend(borrow, y);
890 *xc++ = y & 0xffff;
891 }
892#endif
893 while(!*--xc)
894 wa--;
895 c->wds = wa;
896 return c;
897 }
898
899 static double
900ulp
901#ifdef KR_headers
902 (x) double x;
903#else
904 (double x)
905#endif
906{
907 register Long L;
908 double a;
909
910 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
911#ifndef Sudden_Underflow
912 if (L > 0) {
913#endif
914#ifdef IBM
915 L |= Exp_msk1 >> 4;
916#endif
917 word0(a) = L;
918 word1(a) = 0;
919#ifndef Sudden_Underflow
920 }
921 else {
922 L = -L >> Exp_shift;
923 if (L < Exp_shift) {
924 word0(a) = 0x80000 >> L;
925 word1(a) = 0;
926 }
927 else {
928 word0(a) = 0;
929 L -= Exp_shift;
930 word1(a) = L >= 31 ? 1 : 1 << 31 - L;
931 }
932 }
933#endif
934 return a;
935 }
936
937 static double
938b2d
939#ifdef KR_headers
940 (a, e) Bigint *a; int *e;
941#else
942 (Bigint *a, int *e)
943#endif
944{
945 ULong *xa, *xa0, w, y, z;
946 int k;
947 double d;
948#ifdef VAX
949 ULong d0, d1;
950#else
951#define d0 word0(d)
952#define d1 word1(d)
953#endif
954
955 xa0 = a->x;
956 xa = xa0 + a->wds;
957 y = *--xa;
958#ifdef DEBUG
959 if (!y) Bug("zero y in b2d");
960#endif
961 k = hi0bits(y);
962 *e = 32 - k;
963#ifdef Pack_32
964 if (k < Ebits) {
965 d0 = Exp_1 | y >> Ebits - k;
966 w = xa > xa0 ? *--xa : 0;
967 d1 = y << (32-Ebits) + k | w >> Ebits - k;
968 goto ret_d;
969 }
970 z = xa > xa0 ? *--xa : 0;
971 if (k -= Ebits) {
972 d0 = Exp_1 | y << k | z >> 32 - k;
973 y = xa > xa0 ? *--xa : 0;
974 d1 = z << k | y >> 32 - k;
975 }
976 else {
977 d0 = Exp_1 | y;
978 d1 = z;
979 }
980#else
981 if (k < Ebits + 16) {
982 z = xa > xa0 ? *--xa : 0;
983 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
984 w = xa > xa0 ? *--xa : 0;
985 y = xa > xa0 ? *--xa : 0;
986 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
987 goto ret_d;
988 }
989 z = xa > xa0 ? *--xa : 0;
990 w = xa > xa0 ? *--xa : 0;
991 k -= Ebits + 16;
992 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
993 y = xa > xa0 ? *--xa : 0;
994 d1 = w << k + 16 | y << k;
995#endif
996 ret_d:
997#ifdef VAX
998 word0(d) = d0 >> 16 | d0 << 16;
999 word1(d) = d1 >> 16 | d1 << 16;
1000#else
1001#undef d0
1002#undef d1
1003#endif
1004 return d;
1005 }
1006
1007 static Bigint *
1008d2b
1009#ifdef KR_headers
1010 (d, e, bits) double d; int *e, *bits;
1011#else
1012 (double d, int *e, int *bits)
1013#endif
1014{
1015 Bigint *b;
1016 int de, i, k;
1017 ULong *x, y, z;
1018#ifdef VAX
1019 ULong d0, d1;
1020 d0 = word0(d) >> 16 | word0(d) << 16;
1021 d1 = word1(d) >> 16 | word1(d) << 16;
1022#else
1023#define d0 word0(d)
1024#define d1 word1(d)
1025#endif
1026
1027#ifdef Pack_32
1028 b = Balloc(1);
1029#else
1030 b = Balloc(2);
1031#endif
1032 x = b->x;
1033
1034 z = d0 & Frac_mask;
1035 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
1036#ifdef Sudden_Underflow
1037 de = (int)(d0 >> Exp_shift);
1038#ifndef IBM
1039 z |= Exp_msk11;
1040#endif
1041#else
1042 if (de = (int)(d0 >> Exp_shift))
1043 z |= Exp_msk1;
1044#endif
1045#ifdef Pack_32
1046 if (y = d1) {
1047 if (k = lo0bits(&y)) {
1048 x[0] = y | z << 32 - k;
1049 z >>= k;
1050 }
1051 else
1052 x[0] = y;
1053 i = b->wds = (x[1] = z) ? 2 : 1;
1054 }
1055 else {
1056#ifdef DEBUG
1057 if (!z)
1058 Bug("Zero passed to d2b");
1059#endif
1060 k = lo0bits(&z);
1061 x[0] = z;
1062 i = b->wds = 1;
1063 k += 32;
1064 }
1065#else
1066 if (y = d1) {
1067 if (k = lo0bits(&y))
1068 if (k >= 16) {
1069 x[0] = y | z << 32 - k & 0xffff;
1070 x[1] = z >> k - 16 & 0xffff;
1071 x[2] = z >> k;
1072 i = 2;
1073 }
1074 else {
1075 x[0] = y & 0xffff;
1076 x[1] = y >> 16 | z << 16 - k & 0xffff;
1077 x[2] = z >> k & 0xffff;
1078 x[3] = z >> k+16;
1079 i = 3;
1080 }
1081 else {
1082 x[0] = y & 0xffff;
1083 x[1] = y >> 16;
1084 x[2] = z & 0xffff;
1085 x[3] = z >> 16;
1086 i = 3;
1087 }
1088 }
1089 else {
1090#ifdef DEBUG
1091 if (!z)
1092 Bug("Zero passed to d2b");
1093#endif
1094 k = lo0bits(&z);
1095 if (k >= 16) {
1096 x[0] = z;
1097 i = 0;
1098 }
1099 else {
1100 x[0] = z & 0xffff;
1101 x[1] = z >> 16;
1102 i = 1;
1103 }
1104 k += 32;
1105 }
1106 while(!x[i])
1107 --i;
1108 b->wds = i + 1;
1109#endif
1110#ifndef Sudden_Underflow
1111 if (de) {
1112#endif
1113#ifdef IBM
1114 *e = (de - Bias - (P-1) << 2) + k;
1115 *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1116#else
1117 *e = de - Bias - (P-1) + k;
1118 *bits = P - k;
1119#endif
1120#ifndef Sudden_Underflow
1121 }
1122 else {
1123 *e = de - Bias - (P-1) + 1 + k;
1124#ifdef Pack_32
1125 *bits = 32*i - hi0bits(x[i-1]);
1126#else
1127 *bits = (i+2)*16 - hi0bits(x[i]);
1128#endif
1129 }
1130#endif
1131 return b;
1132 }
1133#undef d0
1134#undef d1
1135
1136 static double
1137ratio
1138#ifdef KR_headers
1139 (a, b) Bigint *a, *b;
1140#else
1141 (Bigint *a, Bigint *b)
1142#endif
1143{
1144 double da, db;
1145 int k, ka, kb;
1146
1147 da = b2d(a, &ka);
1148 db = b2d(b, &kb);
1149#ifdef Pack_32
1150 k = ka - kb + 32*(a->wds - b->wds);
1151#else
1152 k = ka - kb + 16*(a->wds - b->wds);
1153#endif
1154#ifdef IBM
1155 if (k > 0) {
1156 word0(da) += (k >> 2)*Exp_msk1;
1157 if (k &= 3)
1158 da *= 1 << k;
1159 }
1160 else {
1161 k = -k;
1162 word0(db) += (k >> 2)*Exp_msk1;
1163 if (k &= 3)
1164 db *= 1 << k;
1165 }
1166#else
1167 if (k > 0)
1168 word0(da) += k*Exp_msk1;
1169 else {
1170 k = -k;
1171 word0(db) += k*Exp_msk1;
1172 }
1173#endif
1174 return da / db;
1175 }
1176
1177static CONST double
1178tens[] = {
1179 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1180 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1181 1e20, 1e21, 1e22
1182#ifdef VAX
1183 , 1e23, 1e24
1184#endif
1185 };
1186
1187#ifdef IEEE_Arith
1188static CONST double bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1189static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1190#define n_bigtens 5
1191#else
1192#ifdef IBM
1193static CONST double bigtens[] = { 1e16, 1e32, 1e64 };
1194static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1195#define n_bigtens 3
1196#else
1197static CONST double bigtens[] = { 1e16, 1e32 };
1198static CONST double tinytens[] = { 1e-16, 1e-32 };
1199#define n_bigtens 2
1200#endif
1201#endif
1202
1203 double
1204strtod
1205#ifdef KR_headers
1206 (s00, se) CONST char *s00; char **se;
1207#else
1208 (CONST char *s00, char **se)
1209#endif
1210{
1211 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1212 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1213 CONST char *s, *s0, *s1;
1214 double aadj, aadj1, adj, rv, rv0;
1215 Long L;
1216 ULong y, z;
1217 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
1218
1219#ifndef KR_headers
1220 CONST char decimal_point = localeconv()->decimal_point[0];
1221#else
1222 CONST char decimal_point = '.';
1223#endif
1224
1225 sign = nz0 = nz = 0;
1226 rv = 0.;
1227
1228
1229 for(s = s00; isspace(*s); s++)
1230 ;
1231
1232 if (*s == '-') {
1233 sign = 1;
1234 s++;
1235 } else if (*s == '+') {
1236 s++;
1237 }
1238
1239 if (*s == '\0') {
1240 s = s00;
1241 goto ret;
1242 }
1243
1244 if (*s == '0') {
1245 nz0 = 1;
1246 while(*++s == '0') ;
1247 if (!*s)
1248 goto ret;
1249 }
1250 s0 = s;
1251 y = z = 0;
1252 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1253 if (nd < 9)
1254 y = 10*y + c - '0';
1255 else if (nd < 16)
1256 z = 10*z + c - '0';
1257 nd0 = nd;
1258 if (c == decimal_point) {
1259 c = *++s;
1260 if (!nd) {
1261 for(; c == '0'; c = *++s)
1262 nz++;
1263 if (c > '0' && c <= '9') {
1264 s0 = s;
1265 nf += nz;
1266 nz = 0;
1267 goto have_dig;
1268 }
1269 goto dig_done;
1270 }
1271 for(; c >= '0' && c <= '9'; c = *++s) {
1272 have_dig:
1273 nz++;
1274 if (c -= '0') {
1275 nf += nz;
1276 for(i = 1; i < nz; i++)
1277 if (nd++ < 9)
1278 y *= 10;
1279 else if (nd <= DBL_DIG + 1)
1280 z *= 10;
1281 if (nd++ < 9)
1282 y = 10*y + c;
1283 else if (nd <= DBL_DIG + 1)
1284 z = 10*z + c;
1285 nz = 0;
1286 }
1287 }
1288 }
1289 dig_done:
1290 e = 0;
1291 if (c == 'e' || c == 'E') {
1292 if (!nd && !nz && !nz0) {
1293 s = s00;
1294 goto ret;
1295 }
1296 s00 = s;
1297 esign = 0;
1298 switch(c = *++s) {
1299 case '-':
1300 esign = 1;
1301 case '+':
1302 c = *++s;
1303 }
1304 if (c >= '0' && c <= '9') {
1305 while(c == '0')
1306 c = *++s;
1307 if (c > '0' && c <= '9') {
1308 L = c - '0';
1309 s1 = s;
1310 while((c = *++s) >= '0' && c <= '9')
1311 L = 10*L + c - '0';
1312 if (s - s1 > 8 || L > 19999)
1313 /* Avoid confusion from exponents
1314 * so large that e might overflow.
1315 */
1316 e = 19999; /* safe for 16 bit ints */
1317 else
1318 e = (int)L;
1319 if (esign)
1320 e = -e;
1321 }
1322 else
1323 e = 0;
1324 }
1325 else
1326 s = s00;
1327 }
1328 if (!nd) {
1329 if (!nz && !nz0)
1330 s = s00;
1331 goto ret;
1332 }
1333 e1 = e -= nf;
1334
1335 /* Now we have nd0 digits, starting at s0, followed by a
1336 * decimal point, followed by nd-nd0 digits. The number we're
1337 * after is the integer represented by those digits times
1338 * 10**e */
1339
1340 if (!nd0)
1341 nd0 = nd;
1342 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1343 rv = y;
1344 if (k > 9)
1345 rv = tens[k - 9] * rv + z;
1346 bd0 = 0;
1347 if (nd <= DBL_DIG
1348#ifndef RND_PRODQUOT
1349 && FLT_ROUNDS == 1
1350#endif
1351 ) {
1352 if (!e)
1353 goto ret;
1354 if (e > 0) {
1355 if (e <= Ten_pmax) {
1356#ifdef VAX
1357 goto vax_ovfl_check;
1358#else
1359 /* rv = */ rounded_product(rv, tens[e]);
1360 goto ret;
1361#endif
1362 }
1363 i = DBL_DIG - nd;
1364 if (e <= Ten_pmax + i) {
1365 /* A fancier test would sometimes let us do
1366 * this for larger i values.
1367 */
1368 e -= i;
1369 rv *= tens[i];
1370#ifdef VAX
1371 /* VAX exponent range is so narrow we must
1372 * worry about overflow here...
1373 */
1374 vax_ovfl_check:
1375 word0(rv) -= P*Exp_msk1;
1376 /* rv = */ rounded_product(rv, tens[e]);
1377 if ((word0(rv) & Exp_mask)
1378 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1379 goto ovfl;
1380 word0(rv) += P*Exp_msk1;
1381#else
1382 /* rv = */ rounded_product(rv, tens[e]);
1383#endif
1384 goto ret;
1385 }
1386 }
1387#ifndef Inaccurate_Divide
1388 else if (e >= -Ten_pmax) {
1389 /* rv = */ rounded_quotient(rv, tens[-e]);
1390 goto ret;
1391 }
1392#endif
1393 }
1394 e1 += nd - k;
1395
1396 /* Get starting approximation = rv * 10**e1 */
1397
1398 if (e1 > 0) {
1399 if (i = e1 & 15)
1400 rv *= tens[i];
1401 if (e1 &= ~15) {
1402 if (e1 > DBL_MAX_10_EXP) {
1403 ovfl:
1404 errno = ERANGE;
1405#ifdef __STDC__
1406 rv = HUGE_VAL;
1407#else
1408 /* Can't trust HUGE_VAL */
1409#ifdef IEEE_Arith
1410 word0(rv) = Exp_mask;
1411 word1(rv) = 0;
1412#else
1413 word0(rv) = Big0;
1414 word1(rv) = Big1;
1415#endif
1416#endif
1417 if (bd0)
1418 goto retfree;
1419 goto ret;
1420 }
1421 if (e1 >>= 4) {
1422 for(j = 0; e1 > 1; j++, e1 >>= 1)
1423 if (e1 & 1)
1424 rv *= bigtens[j];
1425 /* The last multiplication could overflow. */
1426 word0(rv) -= P*Exp_msk1;
1427 rv *= bigtens[j];
1428 if ((z = word0(rv) & Exp_mask)
1429 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1430 goto ovfl;
1431 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1432 /* set to largest number */
1433 /* (Can't trust DBL_MAX) */
1434 word0(rv) = Big0;
1435 word1(rv) = Big1;
1436 }
1437 else
1438 word0(rv) += P*Exp_msk1;
1439 }
1440
1441 }
1442 }
1443 else if (e1 < 0) {
1444 e1 = -e1;
1445 if (i = e1 & 15)
1446 rv /= tens[i];
1447 if (e1 &= ~15) {
1448 e1 >>= 4;
1449 if (e1 >= 1 << n_bigtens)
1450 goto undfl;
1451 for(j = 0; e1 > 1; j++, e1 >>= 1)
1452 if (e1 & 1)
1453 rv *= tinytens[j];
1454 /* The last multiplication could underflow. */
1455 rv0 = rv;
1456 rv *= tinytens[j];
1457 if (!rv) {
1458 rv = 2.*rv0;
1459 rv *= tinytens[j];
1460 if (!rv) {
1461 undfl:
1462 rv = 0.;
1463 errno = ERANGE;
1464 if (bd0)
1465 goto retfree;
1466 goto ret;
1467 }
1468 word0(rv) = Tiny0;
1469 word1(rv) = Tiny1;
1470 /* The refinement below will clean
1471 * this approximation up.
1472 */
1473 }
1474 }
1475 }
1476
1477 /* Now the hard part -- adjusting rv to the correct value.*/
1478
1479 /* Put digits into bd: true value = bd * 10^e */
1480
1481 bd0 = s2b(s0, nd0, nd, y);
1482
1483 for(;;) {
1484 bd = Balloc(bd0->k);
1485 Bcopy(bd, bd0);
1486 bb = d2b(rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
1487 bs = i2b(1);
1488
1489 if (e >= 0) {
1490 bb2 = bb5 = 0;
1491 bd2 = bd5 = e;
1492 }
1493 else {
1494 bb2 = bb5 = -e;
1495 bd2 = bd5 = 0;
1496 }
1497 if (bbe >= 0)
1498 bb2 += bbe;
1499 else
1500 bd2 -= bbe;
1501 bs2 = bb2;
1502#ifdef Sudden_Underflow
1503#ifdef IBM
1504 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1505#else
1506 j = P + 1 - bbbits;
1507#endif
1508#else
1509 i = bbe + bbbits - 1; /* logb(rv) */
1510 if (i < Emin) /* denormal */
1511 j = bbe + (P-Emin);
1512 else
1513 j = P + 1 - bbbits;
1514#endif
1515 bb2 += j;
1516 bd2 += j;
1517 i = bb2 < bd2 ? bb2 : bd2;
1518 if (i > bs2)
1519 i = bs2;
1520 if (i > 0) {
1521 bb2 -= i;
1522 bd2 -= i;
1523 bs2 -= i;
1524 }
1525 if (bb5 > 0) {
1526 bs = pow5mult(bs, bb5);
1527 bb1 = mult(bs, bb);
1528 Bfree(bb);
1529 bb = bb1;
1530 }
1531 if (bb2 > 0)
1532 bb = lshift(bb, bb2);
1533 if (bd5 > 0)
1534 bd = pow5mult(bd, bd5);
1535 if (bd2 > 0)
1536 bd = lshift(bd, bd2);
1537 if (bs2 > 0)
1538 bs = lshift(bs, bs2);
1539 delta = diff(bb, bd);
1540 dsign = delta->sign;
1541 delta->sign = 0;
1542 i = cmp(delta, bs);
1543 if (i < 0) {
1544 /* Error is less than half an ulp -- check for
1545 * special case of mantissa a power of two.
1546 */
1547 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1548 break;
1549 delta = lshift(delta,Log2P);
1550 if (cmp(delta, bs) > 0)
1551 goto drop_down;
1552 break;
1553 }
1554 if (i == 0) {
1555 /* exactly half-way between */
1556 if (dsign) {
1557 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1558 && word1(rv) == 0xffffffff) {
1559 /*boundary case -- increment exponent*/
1560 word0(rv) = (word0(rv) & Exp_mask)
1561 + Exp_msk1
1562#ifdef IBM
1563 | Exp_msk1 >> 4
1564#endif
1565 ;
1566 word1(rv) = 0;
1567 break;
1568 }
1569 }
1570 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
1571 drop_down:
1572 /* boundary case -- decrement exponent */
1573#ifdef Sudden_Underflow
1574 L = word0(rv) & Exp_mask;
1575#ifdef IBM
1576 if (L < Exp_msk1)
1577#else
1578 if (L <= Exp_msk1)
1579#endif
1580 goto undfl;
1581 L -= Exp_msk1;
1582#else
1583 L = (word0(rv) & Exp_mask) - Exp_msk1;
1584#endif
1585 word0(rv) = L | Bndry_mask1;
1586 word1(rv) = 0xffffffff;
1587#ifdef IBM
1588 goto cont;
1589#else
1590 break;
1591#endif
1592 }
1593#ifndef ROUND_BIASED
1594 if (!(word1(rv) & LSB))
1595 break;
1596#endif
1597 if (dsign)
1598 rv += ulp(rv);
1599#ifndef ROUND_BIASED
1600 else {
1601 rv -= ulp(rv);
1602#ifndef Sudden_Underflow
1603 if (!rv)
1604 goto undfl;
1605#endif
1606 }
1607#endif
1608 break;
1609 }
1610 if ((aadj = ratio(delta, bs)) <= 2.) {
1611 if (dsign)
1612 aadj = aadj1 = 1.;
1613 else if (word1(rv) || word0(rv) & Bndry_mask) {
1614#ifndef Sudden_Underflow
1615 if (word1(rv) == Tiny1 && !word0(rv))
1616 goto undfl;
1617#endif
1618 aadj = 1.;
1619 aadj1 = -1.;
1620 }
1621 else {
1622 /* special case -- power of FLT_RADIX to be */
1623 /* rounded down... */
1624
1625 if (aadj < 2./FLT_RADIX)
1626 aadj = 1./FLT_RADIX;
1627 else
1628 aadj *= 0.5;
1629 aadj1 = -aadj;
1630 }
1631 }
1632 else {
1633 aadj *= 0.5;
1634 aadj1 = dsign ? aadj : -aadj;
1635#ifdef Check_FLT_ROUNDS
1636 switch(FLT_ROUNDS) {
1637 case 2: /* towards +infinity */
1638 aadj1 -= 0.5;
1639 break;
1640 case 0: /* towards 0 */
1641 case 3: /* towards -infinity */
1642 aadj1 += 0.5;
1643 }
1644#else
1645 if (FLT_ROUNDS == 0)
1646 aadj1 += 0.5;
1647#endif
1648 }
1649 y = word0(rv) & Exp_mask;
1650
1651 /* Check for overflow */
1652
1653 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1654 rv0 = rv;
1655 word0(rv) -= P*Exp_msk1;
1656 adj = aadj1 * ulp(rv);
1657 rv += adj;
1658 if ((word0(rv) & Exp_mask) >=
1659 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1660 if (word0(rv0) == Big0 && word1(rv0) == Big1)
1661 goto ovfl;
1662 word0(rv) = Big0;
1663 word1(rv) = Big1;
1664 goto cont;
1665 }
1666 else
1667 word0(rv) += P*Exp_msk1;
1668 }
1669 else {
1670#ifdef Sudden_Underflow
1671 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
1672 rv0 = rv;
1673 word0(rv) += P*Exp_msk1;
1674 adj = aadj1 * ulp(rv);
1675 rv += adj;
1676#ifdef IBM
1677 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
1678#else
1679 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
1680#endif
1681 {
1682 if (word0(rv0) == Tiny0
1683 && word1(rv0) == Tiny1)
1684 goto undfl;
1685 word0(rv) = Tiny0;
1686 word1(rv) = Tiny1;
1687 goto cont;
1688 }
1689 else
1690 word0(rv) -= P*Exp_msk1;
1691 }
1692 else {
1693 adj = aadj1 * ulp(rv);
1694 rv += adj;
1695 }
1696#else
1697 /* Compute adj so that the IEEE rounding rules will
1698 * correctly round rv + adj in some half-way cases.
1699 * If rv * ulp(rv) is denormalized (i.e.,
1700 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1701 * trouble from bits lost to denormalization;
1702 * example: 1.2e-307 .
1703 */
1704 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
1705 aadj1 = (double)(int)(aadj + 0.5);
1706 if (!dsign)
1707 aadj1 = -aadj1;
1708 }
1709 adj = aadj1 * ulp(rv);
1710 rv += adj;
1711#endif
1712 }
1713 z = word0(rv) & Exp_mask;
1714 if (y == z) {
1715 /* Can we stop now? */
1716 L = aadj;
1717 aadj -= L;
1718 /* The tolerances below are conservative. */
1719 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
1720 if (aadj < .4999999 || aadj > .5000001)
1721 break;
1722 }
1723 else if (aadj < .4999999/FLT_RADIX)
1724 break;
1725 }
1726 cont:
1727 Bfree(bb);
1728 Bfree(bd);
1729 Bfree(bs);
1730 Bfree(delta);
1731 }
1732 retfree:
1733 Bfree(bb);
1734 Bfree(bd);
1735 Bfree(bs);
1736 Bfree(bd0);
1737 Bfree(delta);
1738 ret:
1739 if (se)
1740 *se = (char *)s;
1741 return sign ? -rv : rv;
1742 }
1743
1744 static int
1745quorem
1746#ifdef KR_headers
1747 (b, S) Bigint *b, *S;
1748#else
1749 (Bigint *b, Bigint *S)
1750#endif
1751{
1752 int n;
1753 Long borrow, y;
1754 ULong carry, q, ys;
1755 ULong *bx, *bxe, *sx, *sxe;
1756#ifdef Pack_32
1757 Long z;
1758 ULong si, zs;
1759#endif
1760
1761 n = S->wds;
1762#ifdef DEBUG
1763 /*debug*/ if (b->wds > n)
1764 /*debug*/ Bug("oversize b in quorem");
1765#endif
1766 if (b->wds < n)
1767 return 0;
1768 sx = S->x;
1769 sxe = sx + --n;
1770 bx = b->x;
1771 bxe = bx + n;
1772 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
1773#ifdef DEBUG
1774 /*debug*/ if (q > 9)
1775 /*debug*/ Bug("oversized quotient in quorem");
1776#endif
1777 if (q) {
1778 borrow = 0;
1779 carry = 0;
1780 do {
1781#ifdef Pack_32
1782 si = *sx++;
1783 ys = (si & 0xffff) * q + carry;
1784 zs = (si >> 16) * q + (ys >> 16);
1785 carry = zs >> 16;
1786 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1787 borrow = y >> 16;
1788 Sign_Extend(borrow, y);
1789 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1790 borrow = z >> 16;
1791 Sign_Extend(borrow, z);
1792 Storeinc(bx, z, y);
1793#else
1794 ys = *sx++ * q + carry;
1795 carry = ys >> 16;
1796 y = *bx - (ys & 0xffff) + borrow;
1797 borrow = y >> 16;
1798 Sign_Extend(borrow, y);
1799 *bx++ = y & 0xffff;
1800#endif
1801 }
1802 while(sx <= sxe);
1803 if (!*bxe) {
1804 bx = b->x;
1805 while(--bxe > bx && !*bxe)
1806 --n;
1807 b->wds = n;
1808 }
1809 }
1810 if (cmp(b, S) >= 0) {
1811 q++;
1812 borrow = 0;
1813 carry = 0;
1814 bx = b->x;
1815 sx = S->x;
1816 do {
1817#ifdef Pack_32
1818 si = *sx++;
1819 ys = (si & 0xffff) + carry;
1820 zs = (si >> 16) + (ys >> 16);
1821 carry = zs >> 16;
1822 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
1823 borrow = y >> 16;
1824 Sign_Extend(borrow, y);
1825 z = (*bx >> 16) - (zs & 0xffff) + borrow;
1826 borrow = z >> 16;
1827 Sign_Extend(borrow, z);
1828 Storeinc(bx, z, y);
1829#else
1830 ys = *sx++ + carry;
1831 carry = ys >> 16;
1832 y = *bx - (ys & 0xffff) + borrow;
1833 borrow = y >> 16;
1834 Sign_Extend(borrow, y);
1835 *bx++ = y & 0xffff;
1836#endif
1837 }
1838 while(sx <= sxe);
1839 bx = b->x;
1840 bxe = bx + n;
1841 if (!*bxe) {
1842 while(--bxe > bx && !*bxe)
1843 --n;
1844 b->wds = n;
1845 }
1846 }
1847 return q;
1848 }
1849
1850/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
1851 *
1852 * Inspired by "How to Print Floating-Point Numbers Accurately" by
1853 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
1854 *
1855 * Modifications:
1856 * 1. Rather than iterating, we use a simple numeric overestimate
1857 * to determine k = floor(log10(d)). We scale relevant
1858 * quantities using O(log2(k)) rather than O(k) multiplications.
1859 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
1860 * try to generate digits strictly left to right. Instead, we
1861 * compute with fewer bits and propagate the carry if necessary
1862 * when rounding the final digit up. This is often faster.
1863 * 3. Under the assumption that input will be rounded nearest,
1864 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
1865 * That is, we allow equality in stopping tests when the
1866 * round-nearest rule will give the same floating-point value
1867 * as would satisfaction of the stopping test with strict
1868 * inequality.
1869 * 4. We remove common factors of powers of 2 from relevant
1870 * quantities.
1871 * 5. When converting floating-point integers less than 1e16,
1872 * we use floating-point arithmetic rather than resorting
1873 * to multiple-precision integers.
1874 * 6. When asked to produce fewer than 15 digits, we first try
1875 * to get by with floating-point arithmetic; we resort to
1876 * multiple-precision integer arithmetic only if we cannot
1877 * guarantee that the floating-point calculation has given
1878 * the correctly rounded result. For k requested digits and
1879 * "uniformly" distributed input, the probability is
1880 * something like 10^(k-15) that we must resort to the Long
1881 * calculation.
1882 */
1883
1884 char *
1885__dtoa
1886#ifdef KR_headers
1887 (d, mode, ndigits, decpt, sign, rve)
1888 double d; int mode, ndigits, *decpt, *sign; char **rve;
1889#else
1890 (double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
1891#endif
1892{
1893 /* Arguments ndigits, decpt, sign are similar to those
1894 of ecvt and fcvt; trailing zeros are suppressed from
1895 the returned string. If not null, *rve is set to point
1896 to the end of the return value. If d is +-Infinity or NaN,
1897 then *decpt is set to 9999.
1898
1899 mode:
1900 0 ==> shortest string that yields d when read in
1901 and rounded to nearest.
1902 1 ==> like 0, but with Steele & White stopping rule;
1903 e.g. with IEEE P754 arithmetic , mode 0 gives
1904 1e23 whereas mode 1 gives 9.999999999999999e22.
1905 2 ==> max(1,ndigits) significant digits. This gives a
1906 return value similar to that of ecvt, except
1907 that trailing zeros are suppressed.
1908 3 ==> through ndigits past the decimal point. This
1909 gives a return value similar to that from fcvt,
1910 except that trailing zeros are suppressed, and
1911 ndigits can be negative.
1912 4-9 should give the same return values as 2-3, i.e.,
1913 4 <= mode <= 9 ==> same return as mode
1914 2 + (mode & 1). These modes are mainly for
1915 debugging; often they run slower but sometimes
1916 faster than modes 2-3.
1917 4,5,8,9 ==> left-to-right digit generation.
1918 6-9 ==> don't try fast floating-point estimate
1919 (if applicable).
1920
1921 Values of mode other than 0-9 are treated as mode 0.
1922
1923 Sufficient space is allocated to the return value
1924 to hold the suppressed trailing zeros.
1925 */
1926
1927 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
1928 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
1929 spec_case, try_quick;
1930 Long L;
1931#ifndef Sudden_Underflow
1932 int denorm;
1933 ULong x;
1934#endif
1935 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
1936 double d2, ds, eps;
1937 char *s, *s0;
1938 static Bigint *result;
1939 static int result_k;
1940
1941 if (result) {
1942 result->k = result_k;
1943 result->maxwds = 1 << result_k;
1944 Bfree(result);
1945 result = 0;
1946 }
1947
1948 if (word0(d) & Sign_bit) {
1949 /* set sign for everything, including 0's and NaNs */
1950 *sign = 1;
1951 word0(d) &= ~Sign_bit; /* clear sign bit */
1952 }
1953 else
1954 *sign = 0;
1955
1956#if defined(IEEE_Arith) + defined(VAX)
1957#ifdef IEEE_Arith
1958 if ((word0(d) & Exp_mask) == Exp_mask)
1959#else
1960 if (word0(d) == 0x8000)
1961#endif
1962 {
1963 /* Infinity or NaN */
1964 *decpt = 9999;
1965 s =
1966#ifdef IEEE_Arith
1967 !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
1968#endif
1969 "NaN";
1970 if (rve)
1971 *rve =
1972#ifdef IEEE_Arith
1973 s[3] ? s + 8 :
1974#endif
1975 s + 3;
1976 return s;
1977 }
1978#endif
1979#ifdef IBM
1980 d += 0; /* normalize */
1981#endif
1982 if (!d) {
1983 *decpt = 1;
1984 s = "0";
1985 if (rve)
1986 *rve = s + 1;
1987 return s;
1988 }
1989
1990 b = d2b(d, &be, &bbits);
1991#ifdef Sudden_Underflow
1992 i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
1993#else
1994 if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) {
1995#endif
1996 d2 = d;
1997 word0(d2) &= Frac_mask1;
1998 word0(d2) |= Exp_11;
1999#ifdef IBM
2000 if (j = 11 - hi0bits(word0(d2) & Frac_mask))
2001 d2 /= 1 << j;
2002#endif
2003
2004 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
2005 * log10(x) = log(x) / log(10)
2006 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2007 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2008 *
2009 * This suggests computing an approximation k to log10(d) by
2010 *
2011 * k = (i - Bias)*0.301029995663981
2012 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2013 *
2014 * We want k to be too large rather than too small.
2015 * The error in the first-order Taylor series approximation
2016 * is in our favor, so we just round up the constant enough
2017 * to compensate for any error in the multiplication of
2018 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2019 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2020 * adding 1e-13 to the constant term more than suffices.
2021 * Hence we adjust the constant term to 0.1760912590558.
2022 * (We could get a more accurate k by invoking log10,
2023 * but this is probably not worthwhile.)
2024 */
2025
2026 i -= Bias;
2027#ifdef IBM
2028 i <<= 2;
2029 i += j;
2030#endif
2031#ifndef Sudden_Underflow
2032 denorm = 0;
2033 }
2034 else {
2035 /* d is denormalized */
2036
2037 i = bbits + be + (Bias + (P-1) - 1);
2038 x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32
2039 : word1(d) << 32 - i;
2040 d2 = x;
2041 word0(d2) -= 31*Exp_msk1; /* adjust exponent */
2042 i -= (Bias + (P-1) - 1) + 1;
2043 denorm = 1;
2044 }
2045#endif
2046 ds = (d2-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2047 k = (int)ds;
2048 if (ds < 0. && ds != k)
2049 k--; /* want k = floor(ds) */
2050 k_check = 1;
2051 if (k >= 0 && k <= Ten_pmax) {
2052 if (d < tens[k])
2053 k--;
2054 k_check = 0;
2055 }
2056 j = bbits - i - 1;
2057 if (j >= 0) {
2058 b2 = 0;
2059 s2 = j;
2060 }
2061 else {
2062 b2 = -j;
2063 s2 = 0;
2064 }
2065 if (k >= 0) {
2066 b5 = 0;
2067 s5 = k;
2068 s2 += k;
2069 }
2070 else {
2071 b2 -= k;
2072 b5 = -k;
2073 s5 = 0;
2074 }
2075 if (mode < 0 || mode > 9)
2076 mode = 0;
2077 try_quick = 1;
2078 if (mode > 5) {
2079 mode -= 4;
2080 try_quick = 0;
2081 }
2082 leftright = 1;
2083 switch(mode) {
2084 case 0:
2085 case 1:
2086 ilim = ilim1 = -1;
2087 i = 18;
2088 ndigits = 0;
2089 break;
2090 case 2:
2091 leftright = 0;
2092 /* no break */
2093 case 4:
2094 if (ndigits <= 0)
2095 ndigits = 1;
2096 ilim = ilim1 = i = ndigits;
2097 break;
2098 case 3:
2099 leftright = 0;
2100 /* no break */
2101 case 5:
2102 i = ndigits + k + 1;
2103 ilim = i;
2104 ilim1 = i - 1;
2105 if (i <= 0)
2106 i = 1;
2107 }
2108 j = sizeof(ULong);
2109 for(result_k = 0; sizeof(Bigint) - sizeof(ULong) + j <= i;
2110 j <<= 1) result_k++;
2111 result = Balloc(result_k);
2112 s = s0 = (char *)result;
2113
2114 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2115
2116 /* Try to get by with floating-point arithmetic. */
2117
2118 i = 0;
2119 d2 = d;
2120 k0 = k;
2121 ilim0 = ilim;
2122 ieps = 2; /* conservative */
2123 if (k > 0) {
2124 ds = tens[k&0xf];
2125 j = k >> 4;
2126 if (j & Bletch) {
2127 /* prevent overflows */
2128 j &= Bletch - 1;
2129 d /= bigtens[n_bigtens-1];
2130 ieps++;
2131 }
2132 for(; j; j >>= 1, i++)
2133 if (j & 1) {
2134 ieps++;
2135 ds *= bigtens[i];
2136 }
2137 d /= ds;
2138 }
2139 else if (j1 = -k) {
2140 d *= tens[j1 & 0xf];
2141 for(j = j1 >> 4; j; j >>= 1, i++)
2142 if (j & 1) {
2143 ieps++;
2144 d *= bigtens[i];
2145 }
2146 }
2147 if (k_check && d < 1. && ilim > 0) {
2148 if (ilim1 <= 0)
2149 goto fast_failed;
2150 ilim = ilim1;
2151 k--;
2152 d *= 10.;
2153 ieps++;
2154 }
2155 eps = ieps*d + 7.;
2156 word0(eps) -= (P-1)*Exp_msk1;
2157 if (ilim == 0) {
2158 S = mhi = 0;
2159 d -= 5.;
2160 if (d > eps)
2161 goto one_digit;
2162 if (d < -eps)
2163 goto no_digits;
2164 goto fast_failed;
2165 }
2166#ifndef No_leftright
2167 if (leftright) {
2168 /* Use Steele & White method of only
2169 * generating digits needed.
2170 */
2171 eps = 0.5/tens[ilim-1] - eps;
2172 for(i = 0;;) {
2173 L = d;
2174 d -= L;
2175 *s++ = '0' + (int)L;
2176 if (d < eps)
2177 goto ret1;
2178 if (1. - d < eps)
2179 goto bump_up;
2180 if (++i >= ilim)
2181 break;
2182 eps *= 10.;
2183 d *= 10.;
2184 }
2185 }
2186 else {
2187#endif
2188 /* Generate ilim digits, then fix them up. */
2189 eps *= tens[ilim-1];
2190 for(i = 1;; i++, d *= 10.) {
2191 L = d;
2192 d -= L;
2193 *s++ = '0' + (int)L;
2194 if (i == ilim) {
2195 if (d > 0.5 + eps)
2196 goto bump_up;
2197 else if (d < 0.5 - eps) {
2198 while(*--s == '0');
2199 s++;
2200 goto ret1;
2201 }
2202 break;
2203 }
2204 }
2205#ifndef No_leftright
2206 }
2207#endif
2208 fast_failed:
2209 s = s0;
2210 d = d2;
2211 k = k0;
2212 ilim = ilim0;
2213 }
2214
2215 /* Do we have a "small" integer? */
2216
2217 if (be >= 0 && k <= Int_max) {
2218 /* Yes. */
2219 ds = tens[k];
2220 if (ndigits < 0 && ilim <= 0) {
2221 S = mhi = 0;
2222 if (ilim < 0 || d <= 5*ds)
2223 goto no_digits;
2224 goto one_digit;
2225 }
2226 for(i = 1;; i++) {
2227 L = d / ds;
2228 d -= L*ds;
2229#ifdef Check_FLT_ROUNDS
2230 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
2231 if (d < 0) {
2232 L--;
2233 d += ds;
2234 }
2235#endif
2236 *s++ = '0' + (int)L;
2237 if (i == ilim) {
2238 d += d;
2239 if (d > ds || d == ds && L & 1) {
2240 bump_up:
2241 while(*--s == '9')
2242 if (s == s0) {
2243 k++;
2244 *s = '0';
2245 break;
2246 }
2247 ++*s++;
2248 }
2249 break;
2250 }
2251 if (!(d *= 10.))
2252 break;
2253 }
2254 goto ret1;
2255 }
2256
2257 m2 = b2;
2258 m5 = b5;
2259 mhi = mlo = 0;
2260 if (leftright) {
2261 if (mode < 2) {
2262 i =
2263#ifndef Sudden_Underflow
2264 denorm ? be + (Bias + (P-1) - 1 + 1) :
2265#endif
2266#ifdef IBM
2267 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2268#else
2269 1 + P - bbits;
2270#endif
2271 }
2272 else {
2273 j = ilim - 1;
2274 if (m5 >= j)
2275 m5 -= j;
2276 else {
2277 s5 += j -= m5;
2278 b5 += j;
2279 m5 = 0;
2280 }
2281 if ((i = ilim) < 0) {
2282 m2 -= i;
2283 i = 0;
2284 }
2285 }
2286 b2 += i;
2287 s2 += i;
2288 mhi = i2b(1);
2289 }
2290 if (m2 > 0 && s2 > 0) {
2291 i = m2 < s2 ? m2 : s2;
2292 b2 -= i;
2293 m2 -= i;
2294 s2 -= i;
2295 }
2296 if (b5 > 0) {
2297 if (leftright) {
2298 if (m5 > 0) {
2299 mhi = pow5mult(mhi, m5);
2300 b1 = mult(mhi, b);
2301 Bfree(b);
2302 b = b1;
2303 }
2304 if (j = b5 - m5)
2305 b = pow5mult(b, j);
2306 }
2307 else
2308 b = pow5mult(b, b5);
2309 }
2310 S = i2b(1);
2311 if (s5 > 0)
2312 S = pow5mult(S, s5);
2313
2314 /* Check for special case that d is a normalized power of 2. */
2315
2316 if (mode < 2) {
2317 if (!word1(d) && !(word0(d) & Bndry_mask)
2318#ifndef Sudden_Underflow
2319 && word0(d) & Exp_mask
2320#endif
2321 ) {
2322 /* The special case */
2323 b2 += Log2P;
2324 s2 += Log2P;
2325 spec_case = 1;
2326 }
2327 else
2328 spec_case = 0;
2329 }
2330
2331 /* Arrange for convenient computation of quotients:
2332 * shift left if necessary so divisor has 4 leading 0 bits.
2333 *
2334 * Perhaps we should just compute leading 28 bits of S once
2335 * and for all and pass them and a shift to quorem, so it
2336 * can do shifts and ors to compute the numerator for q.
2337 */
2338#ifdef Pack_32
2339 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
2340 i = 32 - i;
2341#else
2342 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
2343 i = 16 - i;
2344#endif
2345 if (i > 4) {
2346 i -= 4;
2347 b2 += i;
2348 m2 += i;
2349 s2 += i;
2350 }
2351 else if (i < 4) {
2352 i += 28;
2353 b2 += i;
2354 m2 += i;
2355 s2 += i;
2356 }
2357 if (b2 > 0)
2358 b = lshift(b, b2);
2359 if (s2 > 0)
2360 S = lshift(S, s2);
2361 if (k_check) {
2362 if (cmp(b,S) < 0) {
2363 k--;
2364 b = multadd(b, 10, 0); /* we botched the k estimate */
2365 if (leftright)
2366 mhi = multadd(mhi, 10, 0);
2367 ilim = ilim1;
2368 }
2369 }
2370 if (ilim <= 0 && mode > 2) {
2371 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
2372 /* no digits, fcvt style */
2373 no_digits:
2374 k = -1 - ndigits;
2375 goto ret;
2376 }
2377 one_digit:
2378 *s++ = '1';
2379 k++;
2380 goto ret;
2381 }
2382 if (leftright) {
2383 if (m2 > 0)
2384 mhi = lshift(mhi, m2);
2385
2386 /* Compute mlo -- check for special case
2387 * that d is a normalized power of 2.
2388 */
2389
2390 mlo = mhi;
2391 if (spec_case) {
2392 mhi = Balloc(mhi->k);
2393 Bcopy(mhi, mlo);
2394 mhi = lshift(mhi, Log2P);
2395 }
2396
2397 for(i = 1;;i++) {
2398 dig = quorem(b,S) + '0';
2399 /* Do we yet have the shortest decimal string
2400 * that will round to d?
2401 */
2402 j = cmp(b, mlo);
2403 delta = diff(S, mhi);
2404 j1 = delta->sign ? 1 : cmp(b, delta);
2405 Bfree(delta);
2406#ifndef ROUND_BIASED
2407 if (j1 == 0 && !mode && !(word1(d) & 1)) {
2408 if (dig == '9')
2409 goto round_9_up;
2410 if (j > 0)
2411 dig++;
2412 *s++ = dig;
2413 goto ret;
2414 }
2415#endif
2416 if (j < 0 || j == 0 && !mode
2417#ifndef ROUND_BIASED
2418 && !(word1(d) & 1)
2419#endif
2420 ) {
2421 if (j1 > 0) {
2422 b = lshift(b, 1);
2423 j1 = cmp(b, S);
2424 if ((j1 > 0 || j1 == 0 && dig & 1)
2425 && dig++ == '9')
2426 goto round_9_up;
2427 }
2428 *s++ = dig;
2429 goto ret;
2430 }
2431 if (j1 > 0) {
2432 if (dig == '9') { /* possible if i == 1 */
2433 round_9_up:
2434 *s++ = '9';
2435 goto roundoff;
2436 }
2437 *s++ = dig + 1;
2438 goto ret;
2439 }
2440 *s++ = dig;
2441 if (i == ilim)
2442 break;
2443 b = multadd(b, 10, 0);
2444 if (mlo == mhi)
2445 mlo = mhi = multadd(mhi, 10, 0);
2446 else {
2447 mlo = multadd(mlo, 10, 0);
2448 mhi = multadd(mhi, 10, 0);
2449 }
2450 }
2451 }
2452 else
2453 for(i = 1;; i++) {
2454 *s++ = dig = quorem(b,S) + '0';
2455 if (i >= ilim)
2456 break;
2457 b = multadd(b, 10, 0);
2458 }
2459
2460 /* Round off last digit */
2461
2462 b = lshift(b, 1);
2463 j = cmp(b, S);
2464 if (j > 0 || j == 0 && dig & 1) {
2465 roundoff:
2466 while(*--s == '9')
2467 if (s == s0) {
2468 k++;
2469 *s++ = '1';
2470 goto ret;
2471 }
2472 ++*s++;
2473 }
2474 else {
2475 while(*--s == '0');
2476 s++;
2477 }
2478 ret:
2479 Bfree(S);
2480 if (mhi) {
2481 if (mlo && mlo != mhi)
2482 Bfree(mlo);
2483 Bfree(mhi);
2484 }
2485 ret1:
2486 Bfree(b);
2487 if (s == s0) { /* don't return empty string */
2488 *s++ = '0';
2489 k = 0;
2490 }
2491 *s = 0;
2492 *decpt = k + 1;
2493 if (rve)
2494 *rve = s;
2495 return s0;
2496 }
2497#ifdef __cplusplus
2498}
2499#endif
generated by cgit v1.2.3-55-g6feb (git 2.39.1) at 2025-04-26 11:40:34 +0000