1 /*	$NetBSD: strtod.c,v 1.45.2.1 2005/04/19 13:35:54 tron Exp $	*/
2 
3 /****************************************************************
4  *
5  * The author of this software is David M. Gay.
6  *
7  * Copyright (c) 1991 by AT&T.
8  *
9  * Permission to use, copy, modify, and distribute this software for any
10  * purpose without fee is hereby granted, provided that this entire notice
11  * is included in all copies of any software which is or includes a copy
12  * or modification of this software and in all copies of the supporting
13  * documentation for such software.
14  *
15  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
16  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
17  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
18  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
19  *
20  ***************************************************************/
21 
22 /* Please send bug reports to
23 	David M. Gay
24 	AT&T Bell Laboratories, Room 2C-463
25 	600 Mountain Avenue
26 	Murray Hill, NJ 07974-2070
27 	U.S.A.
28 	dmg@research.att.com or research!dmg
29  */
30 
31 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
32  *
33  * This strtod returns a nearest machine number to the input decimal
34  * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
35  * broken by the IEEE round-even rule.  Otherwise ties are broken by
36  * biased rounding (add half and chop).
37  *
38  * Inspired loosely by William D. Clinger's paper "How to Read Floating
39  * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
40  *
41  * Modifications:
42  *
43  *	1. We only require IEEE, IBM, or VAX double-precision
44  *		arithmetic (not IEEE double-extended).
45  *	2. We get by with floating-point arithmetic in a case that
46  *		Clinger missed -- when we're computing d * 10^n
47  *		for a small integer d and the integer n is not too
48  *		much larger than 22 (the maximum integer k for which
49  *		we can represent 10^k exactly), we may be able to
50  *		compute (d*10^k) * 10^(e-k) with just one roundoff.
51  *	3. Rather than a bit-at-a-time adjustment of the binary
52  *		result in the hard case, we use floating-point
53  *		arithmetic to determine the adjustment to within
54  *		one bit; only in really hard cases do we need to
55  *		compute a second residual.
56  *	4. Because of 3., we don't need a large table of powers of 10
57  *		for ten-to-e (just some small tables, e.g. of 10^k
58  *		for 0 <= k <= 22).
59  */
60 
61 /*
62  * #define IEEE_LITTLE_ENDIAN for IEEE-arithmetic machines where the least
63  *	significant byte has the lowest address.
64  * #define IEEE_BIG_ENDIAN for IEEE-arithmetic machines where the most
65  *	significant byte has the lowest address.
66  * #define Long int on machines with 32-bit ints and 64-bit longs.
67  * #define Sudden_Underflow for IEEE-format machines without gradual
68  *	underflow (i.e., that flush to zero on underflow).
69  * #define IBM for IBM mainframe-style floating-point arithmetic.
70  * #define VAX for VAX-style floating-point arithmetic.
71  * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
72  * #define No_leftright to omit left-right logic in fast floating-point
73  *	computation of dtoa.
74  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
75  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
76  *	that use extended-precision instructions to compute rounded
77  *	products and quotients) with IBM.
78  * #define ROUND_BIASED for IEEE-format with biased rounding.
79  * #define Inaccurate_Divide for IEEE-format with correctly rounded
80  *	products but inaccurate quotients, e.g., for Intel i860.
81  * #define Just_16 to store 16 bits per 32-bit Long when doing high-precision
82  *	integer arithmetic.  Whether this speeds things up or slows things
83  *	down depends on the machine and the number being converted.
84  * #define KR_headers for old-style C function headers.
85  * #define Bad_float_h if your system lacks a float.h or if it does not
86  *	define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
87  *	FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
88  * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
89  *	if memory is available and otherwise does something you deem
90  *	appropriate.  If MALLOC is undefined, malloc will be invoked
91  *	directly -- and assumed always to succeed.
92  */
93 
94 #define ANDROID_CHANGES
95 
96 #ifdef ANDROID_CHANGES
97 /* Needs to be above math.h include below */
98 #include "fpmath.h"
99 
100 #include <pthread.h>
101 #define mutex_lock(x) pthread_mutex_lock(x)
102 #define mutex_unlock(x) pthread_mutex_unlock(x)
103 #endif
104 
105 #include <sys/cdefs.h>
106 #if defined(LIBC_SCCS) && !defined(lint)
107 __RCSID("$NetBSD: strtod.c,v 1.45.2.1 2005/04/19 13:35:54 tron Exp $");
108 #endif /* LIBC_SCCS and not lint */
109 
110 #define Unsigned_Shifts
111 #if defined(__m68k__) || defined(__sparc__) || defined(__i386__) || \
112     defined(__mips__) || defined(__ns32k__) || defined(__alpha__) || \
113     defined(__powerpc__) || defined(__sh__) || defined(__x86_64__) || \
114     defined(__hppa__) || \
115     defined(__arm__) || defined(__aarch64__)
116 #include <endian.h>
117 #if BYTE_ORDER == BIG_ENDIAN
118 #define IEEE_BIG_ENDIAN
119 #else
120 #define IEEE_LITTLE_ENDIAN
121 #endif
122 #endif
123 
124 #ifdef __vax__
125 #define VAX
126 #endif
127 
128 #if defined(__hppa__) || defined(__mips__) || defined(__sh__)
129 #define	NAN_WORD0	0x7ff40000
130 #else
131 #define	NAN_WORD0	0x7ff80000
132 #endif
133 #define	NAN_WORD1	0
134 
135 #define Long	int32_t
136 #define ULong	u_int32_t
137 
138 #ifdef DEBUG
139 #include "stdio.h"
140 #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
141 #define BugPrintf(x, v) {fprintf(stderr, x, v); exit(1);}
142 #endif
143 
144 #ifdef __cplusplus
145 #include "malloc.h"
146 #include "memory.h"
147 #else
148 #ifndef KR_headers
149 #include "stdlib.h"
150 #include "string.h"
151 #ifndef ANDROID_CHANGES
152 #include "locale.h"
153 #endif /* ANDROID_CHANGES */
154 #else
155 #include "malloc.h"
156 #include "memory.h"
157 #endif
158 #endif
159 #ifndef ANDROID_CHANGES
160 #include "extern.h"
161 #include "reentrant.h"
162 #endif /* ANDROID_CHANGES */
163 
164 #ifdef MALLOC
165 #ifdef KR_headers
166 extern char *MALLOC();
167 #else
168 extern void *MALLOC(size_t);
169 #endif
170 #else
171 #define MALLOC malloc
172 #endif
173 
174 #include "ctype.h"
175 #include "errno.h"
176 #include "float.h"
177 
178 #ifndef __MATH_H__
179 #include "math.h"
180 #endif
181 
182 #ifdef __cplusplus
183 extern "C" {
184 #endif
185 
186 #ifndef CONST
187 #ifdef KR_headers
188 #define CONST /* blank */
189 #else
190 #define CONST const
191 #endif
192 #endif
193 
194 #ifdef Unsigned_Shifts
195 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
196 #else
197 #define Sign_Extend(a,b) /*no-op*/
198 #endif
199 
200 #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN) + defined(VAX) + \
201     defined(IBM) != 1
202 Exactly one of IEEE_LITTLE_ENDIAN IEEE_BIG_ENDIAN, VAX, or
203 IBM should be defined.
204 #endif
205 
206 typedef union {
207 	double d;
208 	ULong ul[2];
209 } _double;
210 #define value(x) ((x).d)
211 #ifdef IEEE_LITTLE_ENDIAN
212 #define word0(x) ((x).ul[1])
213 #define word1(x) ((x).ul[0])
214 #else
215 #define word0(x) ((x).ul[0])
216 #define word1(x) ((x).ul[1])
217 #endif
218 
219 /* The following definition of Storeinc is appropriate for MIPS processors.
220  * An alternative that might be better on some machines is
221  * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
222  */
223 #if defined(IEEE_LITTLE_ENDIAN) + defined(VAX) + defined(__arm__)
224 #define Storeinc(a,b,c) \
225     (((u_short *)(void *)a)[1] = \
226 	(u_short)b, ((u_short *)(void *)a)[0] = (u_short)c, a++)
227 #else
228 #define Storeinc(a,b,c) \
229     (((u_short *)(void *)a)[0] = \
230 	(u_short)b, ((u_short *)(void *)a)[1] = (u_short)c, a++)
231 #endif
232 
233 /* #define P DBL_MANT_DIG */
234 /* Ten_pmax = floor(P*log(2)/log(5)) */
235 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
236 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
237 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
238 
239 #if defined(IEEE_LITTLE_ENDIAN) + defined(IEEE_BIG_ENDIAN)
240 #define Exp_shift  20
241 #define Exp_shift1 20
242 #define Exp_msk1    0x100000
243 #define Exp_msk11   0x100000
244 #define Exp_mask  0x7ff00000
245 #define P 53
246 #define Bias 1023
247 #define IEEE_Arith
248 #define Emin (-1022)
249 #define Exp_1  0x3ff00000
250 #define Exp_11 0x3ff00000
251 #define Ebits 11
252 #define Frac_mask  0xfffff
253 #define Frac_mask1 0xfffff
254 #define Ten_pmax 22
255 #define Bletch 0x10
256 #define Bndry_mask  0xfffff
257 #define Bndry_mask1 0xfffff
258 #define LSB 1
259 #define Sign_bit 0x80000000
260 #define Log2P 1
261 #define Tiny0 0
262 #define Tiny1 1
263 #define Quick_max 14
264 #define Int_max 14
265 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
266 #else
267 #undef  Sudden_Underflow
268 #define Sudden_Underflow
269 #ifdef IBM
270 #define Exp_shift  24
271 #define Exp_shift1 24
272 #define Exp_msk1   0x1000000
273 #define Exp_msk11  0x1000000
274 #define Exp_mask  0x7f000000
275 #define P 14
276 #define Bias 65
277 #define Exp_1  0x41000000
278 #define Exp_11 0x41000000
279 #define Ebits 8	/* exponent has 7 bits, but 8 is the right value in b2d */
280 #define Frac_mask  0xffffff
281 #define Frac_mask1 0xffffff
282 #define Bletch 4
283 #define Ten_pmax 22
284 #define Bndry_mask  0xefffff
285 #define Bndry_mask1 0xffffff
286 #define LSB 1
287 #define Sign_bit 0x80000000
288 #define Log2P 4
289 #define Tiny0 0x100000
290 #define Tiny1 0
291 #define Quick_max 14
292 #define Int_max 15
293 #else /* VAX */
294 #define Exp_shift  23
295 #define Exp_shift1 7
296 #define Exp_msk1    0x80
297 #define Exp_msk11   0x800000
298 #define Exp_mask  0x7f80
299 #define P 56
300 #define Bias 129
301 #define Exp_1  0x40800000
302 #define Exp_11 0x4080
303 #define Ebits 8
304 #define Frac_mask  0x7fffff
305 #define Frac_mask1 0xffff007f
306 #define Ten_pmax 24
307 #define Bletch 2
308 #define Bndry_mask  0xffff007f
309 #define Bndry_mask1 0xffff007f
310 #define LSB 0x10000
311 #define Sign_bit 0x8000
312 #define Log2P 1
313 #define Tiny0 0x80
314 #define Tiny1 0
315 #define Quick_max 15
316 #define Int_max 15
317 #endif
318 #endif
319 
320 #ifndef IEEE_Arith
321 #define ROUND_BIASED
322 #endif
323 
324 #ifdef RND_PRODQUOT
325 #define rounded_product(a,b) a = rnd_prod(a, b)
326 #define rounded_quotient(a,b) a = rnd_quot(a, b)
327 #ifdef KR_headers
328 extern double rnd_prod(), rnd_quot();
329 #else
330 extern double rnd_prod(double, double), rnd_quot(double, double);
331 #endif
332 #else
333 #define rounded_product(a,b) a *= b
334 #define rounded_quotient(a,b) a /= b
335 #endif
336 
337 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
338 #define Big1 0xffffffff
339 
340 #ifndef Just_16
341 /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
342  * This makes some inner loops simpler and sometimes saves work
343  * during multiplications, but it often seems to make things slightly
344  * slower.  Hence the default is now to store 32 bits per Long.
345  */
346 #ifndef Pack_32
347 #define Pack_32
348 #endif
349 #endif
350 
351 #define Kmax 15
352 
353 #ifdef Pack_32
354 #define ULbits 32
355 #define kshift 5
356 #define kmask 31
357 #define ALL_ON 0xffffffff
358 #else
359 #define ULbits 16
360 #define kshift 4
361 #define kmask 15
362 #define ALL_ON 0xffff
363 #endif
364 
365 #define Kmax 15
366 
367  enum {	/* return values from strtodg */
368 	STRTOG_Zero	= 0,
369 	STRTOG_Normal	= 1,
370 	STRTOG_Denormal	= 2,
371 	STRTOG_Infinite	= 3,
372 	STRTOG_NaN	= 4,
373 	STRTOG_NaNbits	= 5,
374 	STRTOG_NoNumber	= 6,
375 	STRTOG_Retmask	= 7,
376 
377 	/* The following may be or-ed into one of the above values. */
378 
379 	STRTOG_Neg	= 0x08, /* does not affect STRTOG_Inexlo or STRTOG_Inexhi */
380 	STRTOG_Inexlo	= 0x10,	/* returned result rounded toward zero */
381 	STRTOG_Inexhi	= 0x20, /* returned result rounded away from zero */
382 	STRTOG_Inexact	= 0x30,
383 	STRTOG_Underflow= 0x40,
384 	STRTOG_Overflow	= 0x80
385 	};
386 
387  typedef struct
388 FPI {
389 	int nbits;
390 	int emin;
391 	int emax;
392 	int rounding;
393 	int sudden_underflow;
394 	} FPI;
395 
396 enum {	/* FPI.rounding values: same as FLT_ROUNDS */
397 	FPI_Round_zero = 0,
398 	FPI_Round_near = 1,
399 	FPI_Round_up = 2,
400 	FPI_Round_down = 3
401 	};
402 
403 #undef SI
404 #ifdef Sudden_Underflow
405 #define SI 1
406 #else
407 #define SI 0
408 #endif
409 
410 #ifdef __cplusplus
411 extern "C" double strtod(const char *s00, char **se);
412 extern "C" char *__dtoa(double d, int mode, int ndigits,
413 			int *decpt, int *sign, char **rve);
414 #endif
415 
416  struct
417 Bigint {
418 	struct Bigint *next;
419 	int k, maxwds, sign, wds;
420 	ULong x[1];
421 };
422 
423  typedef struct Bigint Bigint;
424 
425 CONST unsigned char hexdig[256] = {
426 	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
427 	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
428 	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
429 	0x10, 0x11, 0x12, 0x13, 0x14, 0x15, 0x16, 0x17, 0x18, 0x19, 0, 0, 0, 0, 0, 0,
430 	0, 0x1a, 0x1b, 0x1c, 0x1d, 0x1e, 0x1f, 0, 0, 0, 0, 0, 0, 0, 0, 0,
431 	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
432 	0, 0x1a, 0x1b, 0x1c, 0x1d, 0x1e, 0x1f, 0, 0, 0, 0, 0, 0, 0, 0, 0,
433 	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
434 	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
435 	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
436 	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
437 	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
438 	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
439 	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
440 	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
441 	0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
442 };
443 
444 static int
445 gethex(CONST char **, CONST FPI *, Long *, Bigint **, int, locale_t);
446 
447 
448  static Bigint *freelist[Kmax+1];
449 
450 #ifdef ANDROID_CHANGES
451  static pthread_mutex_t freelist_mutex = PTHREAD_MUTEX_INITIALIZER;
452 #else
453 #ifdef _REENTRANT
454  static mutex_t freelist_mutex = MUTEX_INITIALIZER;
455 #endif
456 #endif
457 
458 /* Special value used to indicate an invalid Bigint value,
459  * e.g. when a memory allocation fails. The idea is that we
460  * want to avoid introducing NULL checks everytime a bigint
461  * computation is performed. Also the NULL value can also be
462  * already used to indicate "value not initialized yet" and
463  * returning NULL might alter the execution code path in
464  * case of OOM.
465  */
466 #define  BIGINT_INVALID   ((Bigint *)&bigint_invalid_value)
467 
468 static const Bigint bigint_invalid_value;
469 
470 
471 static void
copybits(ULong * c,int n,Bigint * b)472 copybits(ULong *c, int n, Bigint *b)
473 {
474 	ULong *ce, *x, *xe;
475 #ifdef Pack_16
476 	int nw, nw1;
477 #endif
478 
479 	ce = c + ((n-1) >> kshift) + 1;
480 	x = b->x;
481 #ifdef Pack_32
482 	xe = x + b->wds;
483 	while(x < xe)
484 		*c++ = *x++;
485 #else
486 	nw = b->wds;
487 	nw1 = nw & 1;
488 	for(xe = x + (nw - nw1); x < xe; x += 2)
489 		Storeinc(c, x[1], x[0]);
490 	if (nw1)
491 		*c++ = *x;
492 #endif
493 	while(c < ce)
494 		*c++ = 0;
495 	}
496 
497  ULong
any_on(Bigint * b,int k)498 any_on(Bigint *b, int k)
499 {
500 	int n, nwds;
501 	ULong *x, *x0, x1, x2;
502 
503 	x = b->x;
504 	nwds = b->wds;
505 	n = k >> kshift;
506 	if (n > nwds)
507 		n = nwds;
508 	else if (n < nwds && (k &= kmask)) {
509 		x1 = x2 = x[n];
510 		x1 >>= k;
511 		x1 <<= k;
512 		if (x1 != x2)
513 			return 1;
514 		}
515 	x0 = x;
516 	x += n;
517 	while(x > x0)
518 		if (*--x)
519 			return 1;
520 	return 0;
521 	}
522 
523  void
rshift(Bigint * b,int k)524 rshift(Bigint *b, int k)
525 {
526 	ULong *x, *x1, *xe, y;
527 	int n;
528 
529 	x = x1 = b->x;
530 	n = k >> kshift;
531 	if (n < b->wds) {
532 		xe = x + b->wds;
533 		x += n;
534 		if (k &= kmask) {
535 			n = ULbits - k;
536 			y = *x++ >> k;
537 			while(x < xe) {
538 				*x1++ = (y | (*x << n)) & ALL_ON;
539 				y = *x++ >> k;
540 				}
541 			if ((*x1 = y) !=0)
542 				x1++;
543 			}
544 		else
545 			while(x < xe)
546 				*x1++ = *x++;
547 		}
548 	if ((b->wds = x1 - b->x) == 0)
549 		b->x[0] = 0;
550 	}
551 
552 
553 typedef union { double d; ULong L[2]; } U;
554 
555 static void
ULtod(ULong * L,ULong * bits,Long exp,int k)556 ULtod(ULong *L, ULong *bits, Long exp, int k)
557 {
558 #  define _0 1
559 #  define _1 0
560 
561 	switch(k & STRTOG_Retmask) {
562 	  case STRTOG_NoNumber:
563 	  case STRTOG_Zero:
564 		L[0] = L[1] = 0;
565 		break;
566 
567 	  case STRTOG_Denormal:
568 		L[_1] = bits[0];
569 		L[_0] = bits[1];
570 		break;
571 
572 	  case STRTOG_Normal:
573 	  case STRTOG_NaNbits:
574 		L[_1] = bits[0];
575 		L[_0] = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20);
576 		break;
577 
578 	  case STRTOG_Infinite:
579 		L[_0] = 0x7ff00000;
580 		L[_1] = 0;
581 		break;
582 
583 #define d_QNAN0 0x7ff80000
584 #define d_QNAN1 0x0
585 	  case STRTOG_NaN:
586 		L[0] = d_QNAN0;
587 		L[1] = d_QNAN1;
588 	  }
589 	if (k & STRTOG_Neg)
590 		L[_0] |= 0x80000000L;
591 	}
592 
593 
594 
595 /* Return BIGINT_INVALID on allocation failure.
596  *
597  * Most of the code here depends on the fact that this function
598  * never returns NULL.
599  */
600  static Bigint *
Balloc(k)601 Balloc
602 #ifdef KR_headers
603 	(k) int k;
604 #else
605 	(int k)
606 #endif
607 {
608 	int x;
609 	Bigint *rv;
610 
611 	mutex_lock(&freelist_mutex);
612 
613 	if ((rv = freelist[k]) != NULL) {
614 		freelist[k] = rv->next;
615 	}
616 	else {
617 		x = 1 << k;
618 		rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(Long));
619 		if (rv == NULL) {
620 		        rv = BIGINT_INVALID;
621 			goto EXIT;
622 		}
623 		rv->k = k;
624 		rv->maxwds = x;
625 	}
626 	rv->sign = rv->wds = 0;
627 EXIT:
628 	mutex_unlock(&freelist_mutex);
629 
630 	return rv;
631 }
632 
633  static void
Bfree(v)634 Bfree
635 #ifdef KR_headers
636 	(v) Bigint *v;
637 #else
638 	(Bigint *v)
639 #endif
640 {
641 	if (v && v != BIGINT_INVALID) {
642 		mutex_lock(&freelist_mutex);
643 
644 		v->next = freelist[v->k];
645 		freelist[v->k] = v;
646 
647 		mutex_unlock(&freelist_mutex);
648 	}
649 }
650 
651 #define Bcopy_valid(x,y) memcpy(&(x)->sign, &(y)->sign, \
652     (y)->wds*sizeof(Long) + 2*sizeof(int))
653 
654 #define Bcopy(x,y)  Bcopy_ptr(&(x),(y))
655 
656  static void
Bcopy_ptr(Bigint ** px,Bigint * y)657 Bcopy_ptr(Bigint **px, Bigint *y)
658 {
659 	if (*px == BIGINT_INVALID)
660 		return; /* no space to store copy */
661 	if (y == BIGINT_INVALID) {
662 		Bfree(*px); /* invalid input */
663 		*px = BIGINT_INVALID;
664 	} else {
665 		Bcopy_valid(*px,y);
666 	}
667 }
668 
669  static Bigint *
multadd(b,m,a)670 multadd
671 #ifdef KR_headers
672 	(b, m, a) Bigint *b; int m, a;
673 #else
674 	(Bigint *b, int m, int a)	/* multiply by m and add a */
675 #endif
676 {
677 	int i, wds;
678 	ULong *x, y;
679 #ifdef Pack_32
680 	ULong xi, z;
681 #endif
682 	Bigint *b1;
683 
684 	if (b == BIGINT_INVALID)
685 		return b;
686 
687 	wds = b->wds;
688 	x = b->x;
689 	i = 0;
690 	do {
691 #ifdef Pack_32
692 		xi = *x;
693 		y = (xi & 0xffff) * m + a;
694 		z = (xi >> 16) * m + (y >> 16);
695 		a = (int)(z >> 16);
696 		*x++ = (z << 16) + (y & 0xffff);
697 #else
698 		y = *x * m + a;
699 		a = (int)(y >> 16);
700 		*x++ = y & 0xffff;
701 #endif
702 	}
703 	while(++i < wds);
704 	if (a) {
705 		if (wds >= b->maxwds) {
706 			b1 = Balloc(b->k+1);
707 			if (b1 == BIGINT_INVALID) {
708 				Bfree(b);
709 				return b1;
710 			}
711 			Bcopy_valid(b1, b);
712 			Bfree(b);
713 			b = b1;
714 			}
715 		b->x[wds++] = a;
716 		b->wds = wds;
717 	}
718 	return b;
719 }
720 
721  Bigint *
increment(Bigint * b)722 increment(Bigint *b)
723 {
724 	ULong *x, *xe;
725 	Bigint *b1;
726 #ifdef Pack_16
727 	ULong carry = 1, y;
728 #endif
729 
730 	x = b->x;
731 	xe = x + b->wds;
732 #ifdef Pack_32
733 	do {
734 		if (*x < (ULong)0xffffffffL) {
735 			++*x;
736 			return b;
737 			}
738 		*x++ = 0;
739 		} while(x < xe);
740 #else
741 	do {
742 		y = *x + carry;
743 		carry = y >> 16;
744 		*x++ = y & 0xffff;
745 		if (!carry)
746 			return b;
747 		} while(x < xe);
748 	if (carry)
749 #endif
750 	{
751 		if (b->wds >= b->maxwds) {
752 			b1 = Balloc(b->k+1);
753 			Bcopy(b1,b);
754 			Bfree(b);
755 			b = b1;
756 			}
757 		b->x[b->wds++] = 1;
758 		}
759 	return b;
760 	}
761 
762 
763  static Bigint *
s2b(s,nd0,nd,y9)764 s2b
765 #ifdef KR_headers
766 	(s, nd0, nd, y9) CONST char *s; int nd0, nd; ULong y9;
767 #else
768 	(CONST char *s, int nd0, int nd, ULong y9)
769 #endif
770 {
771 	Bigint *b;
772 	int i, k;
773 	Long x, y;
774 
775 	x = (nd + 8) / 9;
776 	for(k = 0, y = 1; x > y; y <<= 1, k++) ;
777 #ifdef Pack_32
778 	b = Balloc(k);
779 	if (b == BIGINT_INVALID)
780 		return b;
781 	b->x[0] = y9;
782 	b->wds = 1;
783 #else
784 	b = Balloc(k+1);
785 	if (b == BIGINT_INVALID)
786 		return b;
787 
788 	b->x[0] = y9 & 0xffff;
789 	b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
790 #endif
791 
792 	i = 9;
793 	if (9 < nd0) {
794 		s += 9;
795 		do b = multadd(b, 10, *s++ - '0');
796 			while(++i < nd0);
797 		s++;
798 	}
799 	else
800 		s += 10;
801 	for(; i < nd; i++)
802 		b = multadd(b, 10, *s++ - '0');
803 	return b;
804 }
805 
806  static int
hi0bits(x)807 hi0bits
808 #ifdef KR_headers
809 	(x) ULong x;
810 #else
811 	(ULong x)
812 #endif
813 {
814 	int k = 0;
815 
816 	if (!(x & 0xffff0000)) {
817 		k = 16;
818 		x <<= 16;
819 	}
820 	if (!(x & 0xff000000)) {
821 		k += 8;
822 		x <<= 8;
823 	}
824 	if (!(x & 0xf0000000)) {
825 		k += 4;
826 		x <<= 4;
827 	}
828 	if (!(x & 0xc0000000)) {
829 		k += 2;
830 		x <<= 2;
831 	}
832 	if (!(x & 0x80000000)) {
833 		k++;
834 		if (!(x & 0x40000000))
835 			return 32;
836 	}
837 	return k;
838 }
839 
840  static int
lo0bits(y)841 lo0bits
842 #ifdef KR_headers
843 	(y) ULong *y;
844 #else
845 	(ULong *y)
846 #endif
847 {
848 	int k;
849 	ULong x = *y;
850 
851 	if (x & 7) {
852 		if (x & 1)
853 			return 0;
854 		if (x & 2) {
855 			*y = x >> 1;
856 			return 1;
857 			}
858 		*y = x >> 2;
859 		return 2;
860 	}
861 	k = 0;
862 	if (!(x & 0xffff)) {
863 		k = 16;
864 		x >>= 16;
865 	}
866 	if (!(x & 0xff)) {
867 		k += 8;
868 		x >>= 8;
869 	}
870 	if (!(x & 0xf)) {
871 		k += 4;
872 		x >>= 4;
873 	}
874 	if (!(x & 0x3)) {
875 		k += 2;
876 		x >>= 2;
877 	}
878 	if (!(x & 1)) {
879 		k++;
880 		x >>= 1;
881 		if (!x & 1)
882 			return 32;
883 	}
884 	*y = x;
885 	return k;
886 }
887 
888  static Bigint *
i2b(i)889 i2b
890 #ifdef KR_headers
891 	(i) int i;
892 #else
893 	(int i)
894 #endif
895 {
896 	Bigint *b;
897 
898 	b = Balloc(1);
899 	if (b != BIGINT_INVALID) {
900 		b->x[0] = i;
901 		b->wds = 1;
902 		}
903 	return b;
904 }
905 
906  static Bigint *
mult(a,b)907 mult
908 #ifdef KR_headers
909 	(a, b) Bigint *a, *b;
910 #else
911 	(Bigint *a, Bigint *b)
912 #endif
913 {
914 	Bigint *c;
915 	int k, wa, wb, wc;
916 	ULong carry, y, z;
917 	ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
918 #ifdef Pack_32
919 	ULong z2;
920 #endif
921 
922 	if (a == BIGINT_INVALID || b == BIGINT_INVALID)
923 		return BIGINT_INVALID;
924 
925 	if (a->wds < b->wds) {
926 		c = a;
927 		a = b;
928 		b = c;
929 	}
930 	k = a->k;
931 	wa = a->wds;
932 	wb = b->wds;
933 	wc = wa + wb;
934 	if (wc > a->maxwds)
935 		k++;
936 	c = Balloc(k);
937 	if (c == BIGINT_INVALID)
938 		return c;
939 	for(x = c->x, xa = x + wc; x < xa; x++)
940 		*x = 0;
941 	xa = a->x;
942 	xae = xa + wa;
943 	xb = b->x;
944 	xbe = xb + wb;
945 	xc0 = c->x;
946 #ifdef Pack_32
947 	for(; xb < xbe; xb++, xc0++) {
948 		if ((y = *xb & 0xffff) != 0) {
949 			x = xa;
950 			xc = xc0;
951 			carry = 0;
952 			do {
953 				z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
954 				carry = z >> 16;
955 				z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
956 				carry = z2 >> 16;
957 				Storeinc(xc, z2, z);
958 			}
959 			while(x < xae);
960 			*xc = carry;
961 		}
962 		if ((y = *xb >> 16) != 0) {
963 			x = xa;
964 			xc = xc0;
965 			carry = 0;
966 			z2 = *xc;
967 			do {
968 				z = (*x & 0xffff) * y + (*xc >> 16) + carry;
969 				carry = z >> 16;
970 				Storeinc(xc, z, z2);
971 				z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
972 				carry = z2 >> 16;
973 			}
974 			while(x < xae);
975 			*xc = z2;
976 		}
977 	}
978 #else
979 	for(; xb < xbe; xc0++) {
980 		if (y = *xb++) {
981 			x = xa;
982 			xc = xc0;
983 			carry = 0;
984 			do {
985 				z = *x++ * y + *xc + carry;
986 				carry = z >> 16;
987 				*xc++ = z & 0xffff;
988 			}
989 			while(x < xae);
990 			*xc = carry;
991 		}
992 	}
993 #endif
994 	for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
995 	c->wds = wc;
996 	return c;
997 }
998 
999  static Bigint *p5s;
1000  static pthread_mutex_t p5s_mutex = PTHREAD_MUTEX_INITIALIZER;
1001 
1002  static Bigint *
pow5mult(b,k)1003 pow5mult
1004 #ifdef KR_headers
1005 	(b, k) Bigint *b; int k;
1006 #else
1007 	(Bigint *b, int k)
1008 #endif
1009 {
1010 	Bigint *b1, *p5, *p51;
1011 	int i;
1012 	static const int p05[3] = { 5, 25, 125 };
1013 
1014 	if (b == BIGINT_INVALID)
1015 		return b;
1016 
1017 	if ((i = k & 3) != 0)
1018 		b = multadd(b, p05[i-1], 0);
1019 
1020 	if (!(k = (unsigned int) k >> 2))
1021 		return b;
1022 	mutex_lock(&p5s_mutex);
1023 	if (!(p5 = p5s)) {
1024 		/* first time */
1025 		p5 = i2b(625);
1026 		if (p5 == BIGINT_INVALID) {
1027 			Bfree(b);
1028 			mutex_unlock(&p5s_mutex);
1029 			return p5;
1030 		}
1031 		p5s = p5;
1032 		p5->next = 0;
1033 	}
1034 	for(;;) {
1035 		if (k & 1) {
1036 			b1 = mult(b, p5);
1037 			Bfree(b);
1038 			b = b1;
1039 		}
1040 		if (!(k = (unsigned int) k >> 1))
1041 			break;
1042 		if (!(p51 = p5->next)) {
1043 			p51 = mult(p5,p5);
1044 			if (p51 == BIGINT_INVALID) {
1045 				Bfree(b);
1046 				mutex_unlock(&p5s_mutex);
1047 				return p51;
1048 			}
1049 			p5->next = p51;
1050 			p51->next = 0;
1051 		}
1052 		p5 = p51;
1053 	}
1054 	mutex_unlock(&p5s_mutex);
1055 	return b;
1056 }
1057 
1058  static Bigint *
lshift(b,k)1059 lshift
1060 #ifdef KR_headers
1061 	(b, k) Bigint *b; int k;
1062 #else
1063 	(Bigint *b, int k)
1064 #endif
1065 {
1066 	int i, k1, n, n1;
1067 	Bigint *b1;
1068 	ULong *x, *x1, *xe, z;
1069 
1070 	if (b == BIGINT_INVALID)
1071 		return b;
1072 
1073 #ifdef Pack_32
1074 	n = (unsigned int)k >> 5;
1075 #else
1076 	n = (unsigned int)k >> 4;
1077 #endif
1078 	k1 = b->k;
1079 	n1 = n + b->wds + 1;
1080 	for(i = b->maxwds; n1 > i; i <<= 1)
1081 		k1++;
1082 	b1 = Balloc(k1);
1083 	if (b1 == BIGINT_INVALID) {
1084 		Bfree(b);
1085 		return b1;
1086 	}
1087 	x1 = b1->x;
1088 	for(i = 0; i < n; i++)
1089 		*x1++ = 0;
1090 	x = b->x;
1091 	xe = x + b->wds;
1092 #ifdef Pack_32
1093 	if (k &= 0x1f) {
1094 		k1 = 32 - k;
1095 		z = 0;
1096 		do {
1097 			*x1++ = *x << k | z;
1098 			z = *x++ >> k1;
1099 		}
1100 		while(x < xe);
1101 		if ((*x1 = z) != 0)
1102 			++n1;
1103 	}
1104 #else
1105 	if (k &= 0xf) {
1106 		k1 = 16 - k;
1107 		z = 0;
1108 		do {
1109 			*x1++ = *x << k  & 0xffff | z;
1110 			z = *x++ >> k1;
1111 		}
1112 		while(x < xe);
1113 		if (*x1 = z)
1114 			++n1;
1115 	}
1116 #endif
1117 	else do
1118 		*x1++ = *x++;
1119 		while(x < xe);
1120 	b1->wds = n1 - 1;
1121 	Bfree(b);
1122 	return b1;
1123 }
1124 
1125  static int
cmp(a,b)1126 cmp
1127 #ifdef KR_headers
1128 	(a, b) Bigint *a, *b;
1129 #else
1130 	(Bigint *a, Bigint *b)
1131 #endif
1132 {
1133 	ULong *xa, *xa0, *xb, *xb0;
1134 	int i, j;
1135 
1136 	if (a == BIGINT_INVALID || b == BIGINT_INVALID)
1137 #ifdef DEBUG
1138 		Bug("cmp called with a or b invalid");
1139 #else
1140 		return 0; /* equal - the best we can do right now */
1141 #endif
1142 
1143 	i = a->wds;
1144 	j = b->wds;
1145 #ifdef DEBUG
1146 	if (i > 1 && !a->x[i-1])
1147 		Bug("cmp called with a->x[a->wds-1] == 0");
1148 	if (j > 1 && !b->x[j-1])
1149 		Bug("cmp called with b->x[b->wds-1] == 0");
1150 #endif
1151 	if (i -= j)
1152 		return i;
1153 	xa0 = a->x;
1154 	xa = xa0 + j;
1155 	xb0 = b->x;
1156 	xb = xb0 + j;
1157 	for(;;) {
1158 		if (*--xa != *--xb)
1159 			return *xa < *xb ? -1 : 1;
1160 		if (xa <= xa0)
1161 			break;
1162 	}
1163 	return 0;
1164 }
1165 
1166  static Bigint *
diff(a,b)1167 diff
1168 #ifdef KR_headers
1169 	(a, b) Bigint *a, *b;
1170 #else
1171 	(Bigint *a, Bigint *b)
1172 #endif
1173 {
1174 	Bigint *c;
1175 	int i, wa, wb;
1176 	Long borrow, y;	/* We need signed shifts here. */
1177 	ULong *xa, *xae, *xb, *xbe, *xc;
1178 #ifdef Pack_32
1179 	Long z;
1180 #endif
1181 
1182 	if (a == BIGINT_INVALID || b == BIGINT_INVALID)
1183 		return BIGINT_INVALID;
1184 
1185 	i = cmp(a,b);
1186 	if (!i) {
1187 		c = Balloc(0);
1188 		if (c != BIGINT_INVALID) {
1189 			c->wds = 1;
1190 			c->x[0] = 0;
1191 			}
1192 		return c;
1193 	}
1194 	if (i < 0) {
1195 		c = a;
1196 		a = b;
1197 		b = c;
1198 		i = 1;
1199 	}
1200 	else
1201 		i = 0;
1202 	c = Balloc(a->k);
1203 	if (c == BIGINT_INVALID)
1204 		return c;
1205 	c->sign = i;
1206 	wa = a->wds;
1207 	xa = a->x;
1208 	xae = xa + wa;
1209 	wb = b->wds;
1210 	xb = b->x;
1211 	xbe = xb + wb;
1212 	xc = c->x;
1213 	borrow = 0;
1214 #ifdef Pack_32
1215 	do {
1216 		y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
1217 		borrow = (ULong)y >> 16;
1218 		Sign_Extend(borrow, y);
1219 		z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
1220 		borrow = (ULong)z >> 16;
1221 		Sign_Extend(borrow, z);
1222 		Storeinc(xc, z, y);
1223 	}
1224 	while(xb < xbe);
1225 	while(xa < xae) {
1226 		y = (*xa & 0xffff) + borrow;
1227 		borrow = (ULong)y >> 16;
1228 		Sign_Extend(borrow, y);
1229 		z = (*xa++ >> 16) + borrow;
1230 		borrow = (ULong)z >> 16;
1231 		Sign_Extend(borrow, z);
1232 		Storeinc(xc, z, y);
1233 	}
1234 #else
1235 	do {
1236 		y = *xa++ - *xb++ + borrow;
1237 		borrow = y >> 16;
1238 		Sign_Extend(borrow, y);
1239 		*xc++ = y & 0xffff;
1240 	}
1241 	while(xb < xbe);
1242 	while(xa < xae) {
1243 		y = *xa++ + borrow;
1244 		borrow = y >> 16;
1245 		Sign_Extend(borrow, y);
1246 		*xc++ = y & 0xffff;
1247 	}
1248 #endif
1249 	while(!*--xc)
1250 		wa--;
1251 	c->wds = wa;
1252 	return c;
1253 }
1254 
1255  static double
ulp(_x)1256 ulp
1257 #ifdef KR_headers
1258 	(_x) double _x;
1259 #else
1260 	(double _x)
1261 #endif
1262 {
1263 	_double x;
1264 	Long L;
1265 	_double a;
1266 
1267 	value(x) = _x;
1268 	L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
1269 #ifndef Sudden_Underflow
1270 	if (L > 0) {
1271 #endif
1272 #ifdef IBM
1273 		L |= Exp_msk1 >> 4;
1274 #endif
1275 		word0(a) = L;
1276 		word1(a) = 0;
1277 #ifndef Sudden_Underflow
1278 	}
1279 	else {
1280 		L = (ULong)-L >> Exp_shift;
1281 		if (L < Exp_shift) {
1282 			word0(a) = 0x80000 >> L;
1283 			word1(a) = 0;
1284 		}
1285 		else {
1286 			word0(a) = 0;
1287 			L -= Exp_shift;
1288 			word1(a) = L >= 31 ? 1 : 1 << (31 - L);
1289 		}
1290 	}
1291 #endif
1292 	return value(a);
1293 }
1294 
1295  static double
b2d(a,e)1296 b2d
1297 #ifdef KR_headers
1298 	(a, e) Bigint *a; int *e;
1299 #else
1300 	(Bigint *a, int *e)
1301 #endif
1302 {
1303 	ULong *xa, *xa0, w, y, z;
1304 	int k;
1305 	_double d;
1306 #ifdef VAX
1307 	ULong d0, d1;
1308 #else
1309 #define d0 word0(d)
1310 #define d1 word1(d)
1311 #endif
1312 
1313 	if (a == BIGINT_INVALID)
1314 		return NAN;
1315 
1316 	xa0 = a->x;
1317 	xa = xa0 + a->wds;
1318 	y = *--xa;
1319 #ifdef DEBUG
1320 	if (!y) Bug("zero y in b2d");
1321 #endif
1322 	k = hi0bits(y);
1323 	*e = 32 - k;
1324 #ifdef Pack_32
1325 	if (k < Ebits) {
1326 		d0 = Exp_1 | y >> (Ebits - k);
1327 		w = xa > xa0 ? *--xa : 0;
1328 		d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
1329 		goto ret_d;
1330 	}
1331 	z = xa > xa0 ? *--xa : 0;
1332 	if (k -= Ebits) {
1333 		d0 = Exp_1 | y << k | z >> (32 - k);
1334 		y = xa > xa0 ? *--xa : 0;
1335 		d1 = z << k | y >> (32 - k);
1336 	}
1337 	else {
1338 		d0 = Exp_1 | y;
1339 		d1 = z;
1340 	}
1341 #else
1342 	if (k < Ebits + 16) {
1343 		z = xa > xa0 ? *--xa : 0;
1344 		d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
1345 		w = xa > xa0 ? *--xa : 0;
1346 		y = xa > xa0 ? *--xa : 0;
1347 		d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
1348 		goto ret_d;
1349 	}
1350 	z = xa > xa0 ? *--xa : 0;
1351 	w = xa > xa0 ? *--xa : 0;
1352 	k -= Ebits + 16;
1353 	d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
1354 	y = xa > xa0 ? *--xa : 0;
1355 	d1 = w << k + 16 | y << k;
1356 #endif
1357  ret_d:
1358 #ifdef VAX
1359 	word0(d) = d0 >> 16 | d0 << 16;
1360 	word1(d) = d1 >> 16 | d1 << 16;
1361 #else
1362 #undef d0
1363 #undef d1
1364 #endif
1365 	return value(d);
1366 }
1367 
1368  static Bigint *
d2b(_d,e,bits)1369 d2b
1370 #ifdef KR_headers
1371 	(_d, e, bits) double d; int *e, *bits;
1372 #else
1373 	(double _d, int *e, int *bits)
1374 #endif
1375 {
1376 	Bigint *b;
1377 	int de, i, k;
1378 	ULong *x, y, z;
1379 	_double d;
1380 #ifdef VAX
1381 	ULong d0, d1;
1382 #endif
1383 
1384 	value(d) = _d;
1385 #ifdef VAX
1386 	d0 = word0(d) >> 16 | word0(d) << 16;
1387 	d1 = word1(d) >> 16 | word1(d) << 16;
1388 #else
1389 #define d0 word0(d)
1390 #define d1 word1(d)
1391 #endif
1392 
1393 #ifdef Pack_32
1394 	b = Balloc(1);
1395 #else
1396 	b = Balloc(2);
1397 #endif
1398 	if (b == BIGINT_INVALID)
1399 		return b;
1400 	x = b->x;
1401 
1402 	z = d0 & Frac_mask;
1403 	d0 &= 0x7fffffff;	/* clear sign bit, which we ignore */
1404 #ifdef Sudden_Underflow
1405 	de = (int)(d0 >> Exp_shift);
1406 #ifndef IBM
1407 	z |= Exp_msk11;
1408 #endif
1409 #else
1410 	if ((de = (int)(d0 >> Exp_shift)) != 0)
1411 		z |= Exp_msk1;
1412 #endif
1413 #ifdef Pack_32
1414 	if ((y = d1) != 0) {
1415 		if ((k = lo0bits(&y)) != 0) {
1416 			x[0] = y | z << (32 - k);
1417 			z >>= k;
1418 		}
1419 		else
1420 			x[0] = y;
1421 		i = b->wds = (x[1] = z) ? 2 : 1;
1422 	}
1423 	else {
1424 #ifdef DEBUG
1425 		if (!z)
1426 			Bug("Zero passed to d2b");
1427 #endif
1428 		k = lo0bits(&z);
1429 		x[0] = z;
1430 		i = b->wds = 1;
1431 		k += 32;
1432 	}
1433 #else
1434 	if (y = d1) {
1435 		if (k = lo0bits(&y))
1436 			if (k >= 16) {
1437 				x[0] = y | z << 32 - k & 0xffff;
1438 				x[1] = z >> k - 16 & 0xffff;
1439 				x[2] = z >> k;
1440 				i = 2;
1441 			}
1442 			else {
1443 				x[0] = y & 0xffff;
1444 				x[1] = y >> 16 | z << 16 - k & 0xffff;
1445 				x[2] = z >> k & 0xffff;
1446 				x[3] = z >> k+16;
1447 				i = 3;
1448 			}
1449 		else {
1450 			x[0] = y & 0xffff;
1451 			x[1] = y >> 16;
1452 			x[2] = z & 0xffff;
1453 			x[3] = z >> 16;
1454 			i = 3;
1455 		}
1456 	}
1457 	else {
1458 #ifdef DEBUG
1459 		if (!z)
1460 			Bug("Zero passed to d2b");
1461 #endif
1462 		k = lo0bits(&z);
1463 		if (k >= 16) {
1464 			x[0] = z;
1465 			i = 0;
1466 		}
1467 		else {
1468 			x[0] = z & 0xffff;
1469 			x[1] = z >> 16;
1470 			i = 1;
1471 		}
1472 		k += 32;
1473 	}
1474 	while(!x[i])
1475 		--i;
1476 	b->wds = i + 1;
1477 #endif
1478 #ifndef Sudden_Underflow
1479 	if (de) {
1480 #endif
1481 #ifdef IBM
1482 		*e = (de - Bias - (P-1) << 2) + k;
1483 		*bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
1484 #else
1485 		*e = de - Bias - (P-1) + k;
1486 		*bits = P - k;
1487 #endif
1488 #ifndef Sudden_Underflow
1489 	}
1490 	else {
1491 		*e = de - Bias - (P-1) + 1 + k;
1492 #ifdef Pack_32
1493 		*bits = 32*i - hi0bits(x[i-1]);
1494 #else
1495 		*bits = (i+2)*16 - hi0bits(x[i]);
1496 #endif
1497 		}
1498 #endif
1499 	return b;
1500 }
1501 #undef d0
1502 #undef d1
1503 
1504  static double
ratio(a,b)1505 ratio
1506 #ifdef KR_headers
1507 	(a, b) Bigint *a, *b;
1508 #else
1509 	(Bigint *a, Bigint *b)
1510 #endif
1511 {
1512 	_double da, db;
1513 	int k, ka, kb;
1514 
1515 	if (a == BIGINT_INVALID || b == BIGINT_INVALID)
1516 		return NAN; /* for lack of better value ? */
1517 
1518 	value(da) = b2d(a, &ka);
1519 	value(db) = b2d(b, &kb);
1520 #ifdef Pack_32
1521 	k = ka - kb + 32*(a->wds - b->wds);
1522 #else
1523 	k = ka - kb + 16*(a->wds - b->wds);
1524 #endif
1525 #ifdef IBM
1526 	if (k > 0) {
1527 		word0(da) += (k >> 2)*Exp_msk1;
1528 		if (k &= 3)
1529 			da *= 1 << k;
1530 	}
1531 	else {
1532 		k = -k;
1533 		word0(db) += (k >> 2)*Exp_msk1;
1534 		if (k &= 3)
1535 			db *= 1 << k;
1536 	}
1537 #else
1538 	if (k > 0)
1539 		word0(da) += k*Exp_msk1;
1540 	else {
1541 		k = -k;
1542 		word0(db) += k*Exp_msk1;
1543 	}
1544 #endif
1545 	return value(da) / value(db);
1546 }
1547 
1548 static CONST double
1549 tens[] = {
1550 		1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1551 		1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1552 		1e20, 1e21, 1e22
1553 #ifdef VAX
1554 		, 1e23, 1e24
1555 #endif
1556 };
1557 
1558 #ifdef IEEE_Arith
1559 static CONST double bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
1560 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
1561 #define n_bigtens 5
1562 #else
1563 #ifdef IBM
1564 static CONST double bigtens[] = { 1e16, 1e32, 1e64 };
1565 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
1566 #define n_bigtens 3
1567 #else
1568 static CONST double bigtens[] = { 1e16, 1e32 };
1569 static CONST double tinytens[] = { 1e-16, 1e-32 };
1570 #define n_bigtens 2
1571 #endif
1572 #endif
1573 
1574  double
strtod(s00,se)1575 strtod
1576 #ifdef KR_headers
1577 	(s00, se) CONST char *s00; char **se;
1578 #else
1579 	(CONST char *s00, char **se)
1580 #endif
1581 {
1582 	int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
1583 		 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
1584 	CONST char *s, *s0, *s1;
1585 	double aadj, aadj1, adj;
1586 	_double rv, rv0;
1587 	Long L;
1588 	ULong y, z;
1589 	Bigint *bb1, *bd0;
1590 	Bigint *bb = NULL, *bd = NULL, *bs = NULL, *delta = NULL;/* pacify gcc */
1591 
1592 #ifdef ANDROID_CHANGES
1593 	CONST char decimal_point = '.';
1594 #else /* ANDROID_CHANGES */
1595 #ifndef KR_headers
1596 	CONST char decimal_point = localeconv()->decimal_point[0];
1597 #else
1598 	CONST char decimal_point = '.';
1599 #endif
1600 
1601 #endif /* ANDROID_CHANGES */
1602 
1603 	sign = nz0 = nz = 0;
1604 	value(rv) = 0.;
1605 
1606 
1607 	for(s = s00; isspace((unsigned char) *s); s++)
1608 		;
1609 
1610 	if (*s == '-') {
1611 		sign = 1;
1612 		s++;
1613 	} else if (*s == '+') {
1614 		s++;
1615 	}
1616 
1617 	if (*s == '\0') {
1618 		s = s00;
1619 		goto ret;
1620 	}
1621 
1622 	/* "INF" or "INFINITY" */
1623 	if (tolower((unsigned char)*s) == 'i' && strncasecmp(s, "inf", 3) == 0) {
1624 		if (strncasecmp(s + 3, "inity", 5) == 0)
1625 			s += 8;
1626 		else
1627 			s += 3;
1628 
1629 		value(rv) = HUGE_VAL;
1630 		goto ret;
1631 	}
1632 
1633 #ifdef IEEE_Arith
1634 	/* "NAN" or "NAN(n-char-sequence-opt)" */
1635 	if (tolower((unsigned char)*s) == 'n' && strncasecmp(s, "nan", 3) == 0) {
1636 		/* Build a quiet NaN. */
1637 		word0(rv) = NAN_WORD0;
1638 		word1(rv) = NAN_WORD1;
1639 		s+= 3;
1640 
1641 		/* Don't interpret (n-char-sequence-opt), for now. */
1642 		if (*s == '(') {
1643 			s0 = s;
1644 			for (s++; *s != ')' && *s != '\0'; s++)
1645 				;
1646 			if (*s == ')')
1647 				s++;	/* Skip over closing paren ... */
1648 			else
1649 				s = s0;	/* ... otherwise go back. */
1650 		}
1651 
1652 		goto ret;
1653 	}
1654 #endif
1655 
1656 	if (*s == '0') {
1657 #ifndef NO_HEX_FP /*{*/
1658 		{
1659 		static CONST FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
1660 		Long exp;
1661 		ULong bits[2];
1662 		switch(s[1]) {
1663 		  case 'x':
1664 		  case 'X':
1665 			{
1666 #ifdef Honor_FLT_ROUNDS
1667 			FPI fpi1 = fpi;
1668 			fpi1.rounding = Rounding;
1669 #else
1670 #define fpi1 fpi
1671 #endif
1672 			switch((i = gethex(&s, &fpi1, &exp, &bb, sign, 0)) & STRTOG_Retmask) {
1673 			  case STRTOG_NoNumber:
1674 				s = s00;
1675 				sign = 0;
1676 			  case STRTOG_Zero:
1677 				break;
1678 			  default:
1679 				if (bb) {
1680 					copybits(bits, fpi.nbits, bb);
1681 					Bfree(bb);
1682 					}
1683 				ULtod(((U*)&rv)->L, bits, exp, i);
1684 			  }}
1685 			goto ret;
1686 		  }
1687 		}
1688 #endif /*}*/
1689 		nz0 = 1;
1690 		while(*++s == '0') ;
1691 		if (!*s)
1692 			goto ret;
1693 	}
1694 	s0 = s;
1695 	y = z = 0;
1696 	for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1697 		if (nd < 9)
1698 			y = 10*y + c - '0';
1699 		else if (nd < 16)
1700 			z = 10*z + c - '0';
1701 	nd0 = nd;
1702 	if (c == decimal_point) {
1703 		c = *++s;
1704 		if (!nd) {
1705 			for(; c == '0'; c = *++s)
1706 				nz++;
1707 			if (c > '0' && c <= '9') {
1708 				s0 = s;
1709 				nf += nz;
1710 				nz = 0;
1711 				goto have_dig;
1712 				}
1713 			goto dig_done;
1714 		}
1715 		for(; c >= '0' && c <= '9'; c = *++s) {
1716  have_dig:
1717 			nz++;
1718 			if (c -= '0') {
1719 				nf += nz;
1720 				for(i = 1; i < nz; i++)
1721 					if (nd++ < 9)
1722 						y *= 10;
1723 					else if (nd <= DBL_DIG + 1)
1724 						z *= 10;
1725 				if (nd++ < 9)
1726 					y = 10*y + c;
1727 				else if (nd <= DBL_DIG + 1)
1728 					z = 10*z + c;
1729 				nz = 0;
1730 			}
1731 		}
1732 	}
1733  dig_done:
1734 	e = 0;
1735 	if (c == 'e' || c == 'E') {
1736 		if (!nd && !nz && !nz0) {
1737 			s = s00;
1738 			goto ret;
1739 		}
1740 		s00 = s;
1741 		esign = 0;
1742 		switch(c = *++s) {
1743 			case '-':
1744 				esign = 1;
1745 				/* FALLTHROUGH */
1746 			case '+':
1747 				c = *++s;
1748 		}
1749 		if (c >= '0' && c <= '9') {
1750 			while(c == '0')
1751 				c = *++s;
1752 			if (c > '0' && c <= '9') {
1753 				L = c - '0';
1754 				s1 = s;
1755 				while((c = *++s) >= '0' && c <= '9')
1756 					L = 10*L + c - '0';
1757 				if (s - s1 > 8 || L > 19999)
1758 					/* Avoid confusion from exponents
1759 					 * so large that e might overflow.
1760 					 */
1761 					e = 19999; /* safe for 16 bit ints */
1762 				else
1763 					e = (int)L;
1764 				if (esign)
1765 					e = -e;
1766 			}
1767 			else
1768 				e = 0;
1769 		}
1770 		else
1771 			s = s00;
1772 	}
1773 	if (!nd) {
1774 		if (!nz && !nz0)
1775 			s = s00;
1776 		goto ret;
1777 	}
1778 	e1 = e -= nf;
1779 
1780 	/* Now we have nd0 digits, starting at s0, followed by a
1781 	 * decimal point, followed by nd-nd0 digits.  The number we're
1782 	 * after is the integer represented by those digits times
1783 	 * 10**e */
1784 
1785 	if (!nd0)
1786 		nd0 = nd;
1787 	k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1788 	value(rv) = y;
1789 	if (k > 9)
1790 		value(rv) = tens[k - 9] * value(rv) + z;
1791 	bd0 = 0;
1792 	if (nd <= DBL_DIG
1793 #ifndef RND_PRODQUOT
1794 		&& FLT_ROUNDS == 1
1795 #endif
1796 		) {
1797 		if (!e)
1798 			goto ret;
1799 		if (e > 0) {
1800 			if (e <= Ten_pmax) {
1801 #ifdef VAX
1802 				goto vax_ovfl_check;
1803 #else
1804 				/* value(rv) = */ rounded_product(value(rv),
1805 				    tens[e]);
1806 				goto ret;
1807 #endif
1808 			}
1809 			i = DBL_DIG - nd;
1810 			if (e <= Ten_pmax + i) {
1811 				/* A fancier test would sometimes let us do
1812 				 * this for larger i values.
1813 				 */
1814 				e -= i;
1815 				value(rv) *= tens[i];
1816 #ifdef VAX
1817 				/* VAX exponent range is so narrow we must
1818 				 * worry about overflow here...
1819 				 */
1820  vax_ovfl_check:
1821 				word0(rv) -= P*Exp_msk1;
1822 				/* value(rv) = */ rounded_product(value(rv),
1823 				    tens[e]);
1824 				if ((word0(rv) & Exp_mask)
1825 				 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
1826 					goto ovfl;
1827 				word0(rv) += P*Exp_msk1;
1828 #else
1829 				/* value(rv) = */ rounded_product(value(rv),
1830 				    tens[e]);
1831 #endif
1832 				goto ret;
1833 			}
1834 		}
1835 #ifndef Inaccurate_Divide
1836 		else if (e >= -Ten_pmax) {
1837 			/* value(rv) = */ rounded_quotient(value(rv),
1838 			    tens[-e]);
1839 			goto ret;
1840 		}
1841 #endif
1842 	}
1843 	e1 += nd - k;
1844 
1845 	/* Get starting approximation = rv * 10**e1 */
1846 
1847 	if (e1 > 0) {
1848 		if ((i = e1 & 15) != 0)
1849 			value(rv) *= tens[i];
1850 		if (e1 &= ~15) {
1851 			if (e1 > DBL_MAX_10_EXP) {
1852  ovfl:
1853 				errno = ERANGE;
1854 				value(rv) = HUGE_VAL;
1855 				if (bd0)
1856 					goto retfree;
1857 				goto ret;
1858 			}
1859 			if ((e1 = (unsigned int)e1 >> 4) != 0) {
1860 				for(j = 0; e1 > 1; j++,
1861 				    e1 = (unsigned int)e1 >> 1)
1862 					if (e1 & 1)
1863 						value(rv) *= bigtens[j];
1864 			/* The last multiplication could overflow. */
1865 				word0(rv) -= P*Exp_msk1;
1866 				value(rv) *= bigtens[j];
1867 				if ((z = word0(rv) & Exp_mask)
1868 				 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1869 					goto ovfl;
1870 				if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1871 					/* set to largest number */
1872 					/* (Can't trust DBL_MAX) */
1873 					word0(rv) = Big0;
1874 					word1(rv) = Big1;
1875 					}
1876 				else
1877 					word0(rv) += P*Exp_msk1;
1878 			}
1879 		}
1880 	}
1881 	else if (e1 < 0) {
1882 		e1 = -e1;
1883 		if ((i = e1 & 15) != 0)
1884 			value(rv) /= tens[i];
1885 		if (e1 &= ~15) {
1886 			e1 = (unsigned int)e1 >> 4;
1887 			if (e1 >= 1 << n_bigtens)
1888 				goto undfl;
1889 			for(j = 0; e1 > 1; j++,
1890 			    e1 = (unsigned int)e1 >> 1)
1891 				if (e1 & 1)
1892 					value(rv) *= tinytens[j];
1893 			/* The last multiplication could underflow. */
1894 			value(rv0) = value(rv);
1895 			value(rv) *= tinytens[j];
1896 			if (!value(rv)) {
1897 				value(rv) = 2.*value(rv0);
1898 				value(rv) *= tinytens[j];
1899 				if (!value(rv)) {
1900  undfl:
1901 					value(rv) = 0.;
1902 					errno = ERANGE;
1903 					if (bd0)
1904 						goto retfree;
1905 					goto ret;
1906 				}
1907 				word0(rv) = Tiny0;
1908 				word1(rv) = Tiny1;
1909 				/* The refinement below will clean
1910 				 * this approximation up.
1911 				 */
1912 			}
1913 		}
1914 	}
1915 
1916 	/* Now the hard part -- adjusting rv to the correct value.*/
1917 
1918 	/* Put digits into bd: true value = bd * 10^e */
1919 
1920 	bd0 = s2b(s0, nd0, nd, y);
1921 
1922 	for(;;) {
1923 		bd = Balloc(bd0->k);
1924 		Bcopy(bd, bd0);
1925 		bb = d2b(value(rv), &bbe, &bbbits);	/* rv = bb * 2^bbe */
1926 		bs = i2b(1);
1927 
1928 		if (e >= 0) {
1929 			bb2 = bb5 = 0;
1930 			bd2 = bd5 = e;
1931 		}
1932 		else {
1933 			bb2 = bb5 = -e;
1934 			bd2 = bd5 = 0;
1935 		}
1936 		if (bbe >= 0)
1937 			bb2 += bbe;
1938 		else
1939 			bd2 -= bbe;
1940 		bs2 = bb2;
1941 #ifdef Sudden_Underflow
1942 #ifdef IBM
1943 		j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
1944 #else
1945 		j = P + 1 - bbbits;
1946 #endif
1947 #else
1948 		i = bbe + bbbits - 1;	/* logb(rv) */
1949 		if (i < Emin)	/* denormal */
1950 			j = bbe + (P-Emin);
1951 		else
1952 			j = P + 1 - bbbits;
1953 #endif
1954 		bb2 += j;
1955 		bd2 += j;
1956 		i = bb2 < bd2 ? bb2 : bd2;
1957 		if (i > bs2)
1958 			i = bs2;
1959 		if (i > 0) {
1960 			bb2 -= i;
1961 			bd2 -= i;
1962 			bs2 -= i;
1963 		}
1964 		if (bb5 > 0) {
1965 			bs = pow5mult(bs, bb5);
1966 			bb1 = mult(bs, bb);
1967 			Bfree(bb);
1968 			bb = bb1;
1969 		}
1970 		if (bb2 > 0)
1971 			bb = lshift(bb, bb2);
1972 		if (bd5 > 0)
1973 			bd = pow5mult(bd, bd5);
1974 		if (bd2 > 0)
1975 			bd = lshift(bd, bd2);
1976 		if (bs2 > 0)
1977 			bs = lshift(bs, bs2);
1978 		delta = diff(bb, bd);
1979 		dsign = delta->sign;
1980 		delta->sign = 0;
1981 		i = cmp(delta, bs);
1982 		if (i < 0) {
1983 			/* Error is less than half an ulp -- check for
1984 			 * special case of mantissa a power of two.
1985 			 */
1986 			if (dsign || word1(rv) || word0(rv) & Bndry_mask)
1987 				break;
1988 			delta = lshift(delta,Log2P);
1989 			if (cmp(delta, bs) > 0)
1990 				goto drop_down;
1991 			break;
1992 		}
1993 		if (i == 0) {
1994 			/* exactly half-way between */
1995 			if (dsign) {
1996 				if ((word0(rv) & Bndry_mask1) == Bndry_mask1
1997 				 &&  word1(rv) == 0xffffffff) {
1998 					/*boundary case -- increment exponent*/
1999 					word0(rv) = (word0(rv) & Exp_mask)
2000 						+ Exp_msk1
2001 #ifdef IBM
2002 						| Exp_msk1 >> 4
2003 #endif
2004 						;
2005 					word1(rv) = 0;
2006 					break;
2007 				}
2008 			}
2009 			else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
2010  drop_down:
2011 				/* boundary case -- decrement exponent */
2012 #ifdef Sudden_Underflow
2013 				L = word0(rv) & Exp_mask;
2014 #ifdef IBM
2015 				if (L <  Exp_msk1)
2016 #else
2017 				if (L <= Exp_msk1)
2018 #endif
2019 					goto undfl;
2020 				L -= Exp_msk1;
2021 #else
2022 				L = (word0(rv) & Exp_mask) - Exp_msk1;
2023 #endif
2024 				word0(rv) = L | Bndry_mask1;
2025 				word1(rv) = 0xffffffff;
2026 #ifdef IBM
2027 				goto cont;
2028 #else
2029 				break;
2030 #endif
2031 			}
2032 #ifndef ROUND_BIASED
2033 			if (!(word1(rv) & LSB))
2034 				break;
2035 #endif
2036 			if (dsign)
2037 				value(rv) += ulp(value(rv));
2038 #ifndef ROUND_BIASED
2039 			else {
2040 				value(rv) -= ulp(value(rv));
2041 #ifndef Sudden_Underflow
2042 				if (!value(rv))
2043 					goto undfl;
2044 #endif
2045 			}
2046 #endif
2047 			break;
2048 		}
2049 		if ((aadj = ratio(delta, bs)) <= 2.) {
2050 			if (dsign)
2051 				aadj = aadj1 = 1.;
2052 			else if (word1(rv) || word0(rv) & Bndry_mask) {
2053 #ifndef Sudden_Underflow
2054 				if (word1(rv) == Tiny1 && !word0(rv))
2055 					goto undfl;
2056 #endif
2057 				aadj = 1.;
2058 				aadj1 = -1.;
2059 			}
2060 			else {
2061 				/* special case -- power of FLT_RADIX to be */
2062 				/* rounded down... */
2063 
2064 				if (aadj < 2./FLT_RADIX)
2065 					aadj = 1./FLT_RADIX;
2066 				else
2067 					aadj *= 0.5;
2068 				aadj1 = -aadj;
2069 				}
2070 		}
2071 		else {
2072 			aadj *= 0.5;
2073 			aadj1 = dsign ? aadj : -aadj;
2074 #ifdef Check_FLT_ROUNDS
2075 			switch(FLT_ROUNDS) {
2076 				case 2: /* towards +infinity */
2077 					aadj1 -= 0.5;
2078 					break;
2079 				case 0: /* towards 0 */
2080 				case 3: /* towards -infinity */
2081 					aadj1 += 0.5;
2082 			}
2083 #else
2084 			if (FLT_ROUNDS == 0)
2085 				aadj1 += 0.5;
2086 #endif
2087 		}
2088 		y = word0(rv) & Exp_mask;
2089 
2090 		/* Check for overflow */
2091 
2092 		if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2093 			value(rv0) = value(rv);
2094 			word0(rv) -= P*Exp_msk1;
2095 			adj = aadj1 * ulp(value(rv));
2096 			value(rv) += adj;
2097 			if ((word0(rv) & Exp_mask) >=
2098 					Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2099 				if (word0(rv0) == Big0 && word1(rv0) == Big1)
2100 					goto ovfl;
2101 				word0(rv) = Big0;
2102 				word1(rv) = Big1;
2103 				goto cont;
2104 			}
2105 			else
2106 				word0(rv) += P*Exp_msk1;
2107 		}
2108 		else {
2109 #ifdef Sudden_Underflow
2110 			if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
2111 				value(rv0) = value(rv);
2112 				word0(rv) += P*Exp_msk1;
2113 				adj = aadj1 * ulp(value(rv));
2114 				value(rv) += adj;
2115 #ifdef IBM
2116 				if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
2117 #else
2118 				if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
2119 #endif
2120 				{
2121 					if (word0(rv0) == Tiny0
2122 					 && word1(rv0) == Tiny1)
2123 						goto undfl;
2124 					word0(rv) = Tiny0;
2125 					word1(rv) = Tiny1;
2126 					goto cont;
2127 				}
2128 				else
2129 					word0(rv) -= P*Exp_msk1;
2130 				}
2131 			else {
2132 				adj = aadj1 * ulp(value(rv));
2133 				value(rv) += adj;
2134 			}
2135 #else
2136 			/* Compute adj so that the IEEE rounding rules will
2137 			 * correctly round rv + adj in some half-way cases.
2138 			 * If rv * ulp(rv) is denormalized (i.e.,
2139 			 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
2140 			 * trouble from bits lost to denormalization;
2141 			 * example: 1.2e-307 .
2142 			 */
2143 			if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
2144 				aadj1 = (double)(int)(aadj + 0.5);
2145 				if (!dsign)
2146 					aadj1 = -aadj1;
2147 			}
2148 			adj = aadj1 * ulp(value(rv));
2149 			value(rv) += adj;
2150 #endif
2151 		}
2152 		z = word0(rv) & Exp_mask;
2153 		if (y == z) {
2154 			/* Can we stop now? */
2155 			L = aadj;
2156 			aadj -= L;
2157 			/* The tolerances below are conservative. */
2158 			if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
2159 				if (aadj < .4999999 || aadj > .5000001)
2160 					break;
2161 			}
2162 			else if (aadj < .4999999/FLT_RADIX)
2163 				break;
2164 		}
2165  cont:
2166 		Bfree(bb);
2167 		Bfree(bd);
2168 		Bfree(bs);
2169 		Bfree(delta);
2170 	}
2171  retfree:
2172 	Bfree(bb);
2173 	Bfree(bd);
2174 	Bfree(bs);
2175 	Bfree(bd0);
2176 	Bfree(delta);
2177  ret:
2178 	if (se)
2179 		/* LINTED interface specification */
2180 		*se = (char *)s;
2181 	return sign ? -value(rv) : value(rv);
2182 }
2183 
2184  static int
quorem(b,S)2185 quorem
2186 #ifdef KR_headers
2187 	(b, S) Bigint *b, *S;
2188 #else
2189 	(Bigint *b, Bigint *S)
2190 #endif
2191 {
2192 	int n;
2193 	Long borrow, y;
2194 	ULong carry, q, ys;
2195 	ULong *bx, *bxe, *sx, *sxe;
2196 #ifdef Pack_32
2197 	Long z;
2198 	ULong si, zs;
2199 #endif
2200 
2201 	if (b == BIGINT_INVALID || S == BIGINT_INVALID)
2202 		return 0;
2203 
2204 	n = S->wds;
2205 #ifdef DEBUG
2206 	/*debug*/ if (b->wds > n)
2207 	/*debug*/	Bug("oversize b in quorem");
2208 #endif
2209 	if (b->wds < n)
2210 		return 0;
2211 	sx = S->x;
2212 	sxe = sx + --n;
2213 	bx = b->x;
2214 	bxe = bx + n;
2215 	q = *bxe / (*sxe + 1);	/* ensure q <= true quotient */
2216 #ifdef DEBUG
2217 	/*debug*/ if (q > 9)
2218 	/*debug*/	Bug("oversized quotient in quorem");
2219 #endif
2220 	if (q) {
2221 		borrow = 0;
2222 		carry = 0;
2223 		do {
2224 #ifdef Pack_32
2225 			si = *sx++;
2226 			ys = (si & 0xffff) * q + carry;
2227 			zs = (si >> 16) * q + (ys >> 16);
2228 			carry = zs >> 16;
2229 			y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
2230 			borrow = (ULong)y >> 16;
2231 			Sign_Extend(borrow, y);
2232 			z = (*bx >> 16) - (zs & 0xffff) + borrow;
2233 			borrow = (ULong)z >> 16;
2234 			Sign_Extend(borrow, z);
2235 			Storeinc(bx, z, y);
2236 #else
2237 			ys = *sx++ * q + carry;
2238 			carry = ys >> 16;
2239 			y = *bx - (ys & 0xffff) + borrow;
2240 			borrow = y >> 16;
2241 			Sign_Extend(borrow, y);
2242 			*bx++ = y & 0xffff;
2243 #endif
2244 		}
2245 		while(sx <= sxe);
2246 		if (!*bxe) {
2247 			bx = b->x;
2248 			while(--bxe > bx && !*bxe)
2249 				--n;
2250 			b->wds = n;
2251 		}
2252 	}
2253 	if (cmp(b, S) >= 0) {
2254 		q++;
2255 		borrow = 0;
2256 		carry = 0;
2257 		bx = b->x;
2258 		sx = S->x;
2259 		do {
2260 #ifdef Pack_32
2261 			si = *sx++;
2262 			ys = (si & 0xffff) + carry;
2263 			zs = (si >> 16) + (ys >> 16);
2264 			carry = zs >> 16;
2265 			y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
2266 			borrow = (ULong)y >> 16;
2267 			Sign_Extend(borrow, y);
2268 			z = (*bx >> 16) - (zs & 0xffff) + borrow;
2269 			borrow = (ULong)z >> 16;
2270 			Sign_Extend(borrow, z);
2271 			Storeinc(bx, z, y);
2272 #else
2273 			ys = *sx++ + carry;
2274 			carry = ys >> 16;
2275 			y = *bx - (ys & 0xffff) + borrow;
2276 			borrow = y >> 16;
2277 			Sign_Extend(borrow, y);
2278 			*bx++ = y & 0xffff;
2279 #endif
2280 		}
2281 		while(sx <= sxe);
2282 		bx = b->x;
2283 		bxe = bx + n;
2284 		if (!*bxe) {
2285 			while(--bxe > bx && !*bxe)
2286 				--n;
2287 			b->wds = n;
2288 		}
2289 	}
2290 	return q;
2291 }
2292 
2293 /* freedtoa(s) must be used to free values s returned by dtoa
2294  * when MULTIPLE_THREADS is #defined.  It should be used in all cases,
2295  * but for consistency with earlier versions of dtoa, it is optional
2296  * when MULTIPLE_THREADS is not defined.
2297  */
2298 
2299 void
2300 #ifdef KR_headers
freedtoa(s)2301 freedtoa(s) char *s;
2302 #else
2303 freedtoa(char *s)
2304 #endif
2305 {
2306 	free(s);
2307 }
2308 
2309 
2310 
2311 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
2312  *
2313  * Inspired by "How to Print Floating-Point Numbers Accurately" by
2314  * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
2315  *
2316  * Modifications:
2317  *	1. Rather than iterating, we use a simple numeric overestimate
2318  *	   to determine k = floor(log10(d)).  We scale relevant
2319  *	   quantities using O(log2(k)) rather than O(k) multiplications.
2320  *	2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
2321  *	   try to generate digits strictly left to right.  Instead, we
2322  *	   compute with fewer bits and propagate the carry if necessary
2323  *	   when rounding the final digit up.  This is often faster.
2324  *	3. Under the assumption that input will be rounded nearest,
2325  *	   mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
2326  *	   That is, we allow equality in stopping tests when the
2327  *	   round-nearest rule will give the same floating-point value
2328  *	   as would satisfaction of the stopping test with strict
2329  *	   inequality.
2330  *	4. We remove common factors of powers of 2 from relevant
2331  *	   quantities.
2332  *	5. When converting floating-point integers less than 1e16,
2333  *	   we use floating-point arithmetic rather than resorting
2334  *	   to multiple-precision integers.
2335  *	6. When asked to produce fewer than 15 digits, we first try
2336  *	   to get by with floating-point arithmetic; we resort to
2337  *	   multiple-precision integer arithmetic only if we cannot
2338  *	   guarantee that the floating-point calculation has given
2339  *	   the correctly rounded result.  For k requested digits and
2340  *	   "uniformly" distributed input, the probability is
2341  *	   something like 10^(k-15) that we must resort to the Long
2342  *	   calculation.
2343  */
2344 
2345 __LIBC_HIDDEN__  char *
__dtoa(_d,mode,ndigits,decpt,sign,rve)2346 __dtoa
2347 #ifdef KR_headers
2348 	(_d, mode, ndigits, decpt, sign, rve)
2349 	double _d; int mode, ndigits, *decpt, *sign; char **rve;
2350 #else
2351 	(double _d, int mode, int ndigits, int *decpt, int *sign, char **rve)
2352 #endif
2353 {
2354  /*	Arguments ndigits, decpt, sign are similar to those
2355 	of ecvt and fcvt; trailing zeros are suppressed from
2356 	the returned string.  If not null, *rve is set to point
2357 	to the end of the return value.  If d is +-Infinity or NaN,
2358 	then *decpt is set to 9999.
2359 
2360 	mode:
2361 		0 ==> shortest string that yields d when read in
2362 			and rounded to nearest.
2363 		1 ==> like 0, but with Steele & White stopping rule;
2364 			e.g. with IEEE P754 arithmetic , mode 0 gives
2365 			1e23 whereas mode 1 gives 9.999999999999999e22.
2366 		2 ==> max(1,ndigits) significant digits.  This gives a
2367 			return value similar to that of ecvt, except
2368 			that trailing zeros are suppressed.
2369 		3 ==> through ndigits past the decimal point.  This
2370 			gives a return value similar to that from fcvt,
2371 			except that trailing zeros are suppressed, and
2372 			ndigits can be negative.
2373 		4-9 should give the same return values as 2-3, i.e.,
2374 			4 <= mode <= 9 ==> same return as mode
2375 			2 + (mode & 1).  These modes are mainly for
2376 			debugging; often they run slower but sometimes
2377 			faster than modes 2-3.
2378 		4,5,8,9 ==> left-to-right digit generation.
2379 		6-9 ==> don't try fast floating-point estimate
2380 			(if applicable).
2381 
2382 		Values of mode other than 0-9 are treated as mode 0.
2383 
2384 		Sufficient space is allocated to the return value
2385 		to hold the suppressed trailing zeros.
2386 	*/
2387 
2388 	int bbits, b2, b5, be, dig, i, ieps, ilim0,
2389 		j, jj1, k, k0, k_check, leftright, m2, m5, s2, s5,
2390 		try_quick;
2391 	int ilim = 0, ilim1 = 0, spec_case = 0;	/* pacify gcc */
2392 	Long L;
2393 #ifndef Sudden_Underflow
2394 	int denorm;
2395 	ULong x;
2396 #endif
2397 	Bigint *b, *b1, *delta, *mhi, *S;
2398 	Bigint *mlo = NULL; /* pacify gcc */
2399 	double ds;
2400 	char *s, *s0;
2401 	Bigint *result = NULL;
2402 	int result_k = 0;
2403 	_double d, d2, eps;
2404 
2405 	value(d) = _d;
2406 
2407 	if (word0(d) & Sign_bit) {
2408 		/* set sign for everything, including 0's and NaNs */
2409 		*sign = 1;
2410 		word0(d) &= ~Sign_bit;	/* clear sign bit */
2411 	}
2412 	else
2413 		*sign = 0;
2414 
2415 #if defined(IEEE_Arith) + defined(VAX)
2416 #ifdef IEEE_Arith
2417 	if ((word0(d) & Exp_mask) == Exp_mask)
2418 #else
2419 	if (word0(d)  == 0x8000)
2420 #endif
2421 	{
2422 		/* Infinity or NaN */
2423 		*decpt = 9999;
2424 		s =
2425 #ifdef IEEE_Arith
2426 			!word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
2427 #endif
2428 				"NaN";
2429 		result = Balloc(strlen(s)+1);
2430 		if (result == BIGINT_INVALID)
2431 			return NULL;
2432 		s0 = (char *)(void *)result;
2433 		strcpy(s0, s);
2434 		if (rve)
2435 			*rve =
2436 #ifdef IEEE_Arith
2437 				s0[3] ? s0 + 8 :
2438 #endif
2439 				s0 + 3;
2440 		return s0;
2441 	}
2442 #endif
2443 #ifdef IBM
2444 	value(d) += 0; /* normalize */
2445 #endif
2446 	if (!value(d)) {
2447 		*decpt = 1;
2448 		result = Balloc(2);
2449 		if (result == BIGINT_INVALID)
2450 			return NULL;
2451 		s0 = (char *)(void *)result;
2452 		strcpy(s0, "0");
2453 		if (rve)
2454 			*rve = s0 + 1;
2455 		return s0;
2456 	}
2457 
2458 	b = d2b(value(d), &be, &bbits);
2459 #ifdef Sudden_Underflow
2460 	i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
2461 #else
2462 	if ((i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) != 0) {
2463 #endif
2464 		value(d2) = value(d);
2465 		word0(d2) &= Frac_mask1;
2466 		word0(d2) |= Exp_11;
2467 #ifdef IBM
2468 		if (j = 11 - hi0bits(word0(d2) & Frac_mask))
2469 			value(d2) /= 1 << j;
2470 #endif
2471 
2472 		/* log(x)	~=~ log(1.5) + (x-1.5)/1.5
2473 		 * log10(x)	 =  log(x) / log(10)
2474 		 *		~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
2475 		 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
2476 		 *
2477 		 * This suggests computing an approximation k to log10(d) by
2478 		 *
2479 		 * k = (i - Bias)*0.301029995663981
2480 		 *	+ ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
2481 		 *
2482 		 * We want k to be too large rather than too small.
2483 		 * The error in the first-order Taylor series approximation
2484 		 * is in our favor, so we just round up the constant enough
2485 		 * to compensate for any error in the multiplication of
2486 		 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
2487 		 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
2488 		 * adding 1e-13 to the constant term more than suffices.
2489 		 * Hence we adjust the constant term to 0.1760912590558.
2490 		 * (We could get a more accurate k by invoking log10,
2491 		 *  but this is probably not worthwhile.)
2492 		 */
2493 
2494 		i -= Bias;
2495 #ifdef IBM
2496 		i <<= 2;
2497 		i += j;
2498 #endif
2499 #ifndef Sudden_Underflow
2500 		denorm = 0;
2501 	}
2502 	else {
2503 		/* d is denormalized */
2504 
2505 		i = bbits + be + (Bias + (P-1) - 1);
2506 		x = i > 32  ? word0(d) << (64 - i) | word1(d) >> (i - 32)
2507 			    : word1(d) << (32 - i);
2508 		value(d2) = x;
2509 		word0(d2) -= 31*Exp_msk1; /* adjust exponent */
2510 		i -= (Bias + (P-1) - 1) + 1;
2511 		denorm = 1;
2512 	}
2513 #endif
2514 	ds = (value(d2)-1.5)*0.289529654602168 + 0.1760912590558 +
2515 	    i*0.301029995663981;
2516 	k = (int)ds;
2517 	if (ds < 0. && ds != k)
2518 		k--;	/* want k = floor(ds) */
2519 	k_check = 1;
2520 	if (k >= 0 && k <= Ten_pmax) {
2521 		if (value(d) < tens[k])
2522 			k--;
2523 		k_check = 0;
2524 	}
2525 	j = bbits - i - 1;
2526 	if (j >= 0) {
2527 		b2 = 0;
2528 		s2 = j;
2529 	}
2530 	else {
2531 		b2 = -j;
2532 		s2 = 0;
2533 	}
2534 	if (k >= 0) {
2535 		b5 = 0;
2536 		s5 = k;
2537 		s2 += k;
2538 	}
2539 	else {
2540 		b2 -= k;
2541 		b5 = -k;
2542 		s5 = 0;
2543 	}
2544 	if (mode < 0 || mode > 9)
2545 		mode = 0;
2546 	try_quick = 1;
2547 	if (mode > 5) {
2548 		mode -= 4;
2549 		try_quick = 0;
2550 	}
2551 	leftright = 1;
2552 	switch(mode) {
2553 		case 0:
2554 		case 1:
2555 			ilim = ilim1 = -1;
2556 			i = 18;
2557 			ndigits = 0;
2558 			break;
2559 		case 2:
2560 			leftright = 0;
2561 			/* FALLTHROUGH */
2562 		case 4:
2563 			if (ndigits <= 0)
2564 				ndigits = 1;
2565 			ilim = ilim1 = i = ndigits;
2566 			break;
2567 		case 3:
2568 			leftright = 0;
2569 			/* FALLTHROUGH */
2570 		case 5:
2571 			i = ndigits + k + 1;
2572 			ilim = i;
2573 			ilim1 = i - 1;
2574 			if (i <= 0)
2575 				i = 1;
2576 	}
2577 	j = sizeof(ULong);
2578         for(result_k = 0; (int)(sizeof(Bigint) - sizeof(ULong)) + j <= i;
2579 		j <<= 1) result_k++;
2580         // this is really a ugly hack, the code uses Balloc
2581         // instead of malloc, but casts the result into a char*
2582         // it seems the only reason to do that is due to the
2583         // complicated way the block size need to be computed
2584         // buuurk....
2585 	result = Balloc(result_k);
2586 	if (result == BIGINT_INVALID) {
2587 		Bfree(b);
2588 		return NULL;
2589 	}
2590 	s = s0 = (char *)(void *)result;
2591 
2592 	if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2593 
2594 		/* Try to get by with floating-point arithmetic. */
2595 
2596 		i = 0;
2597 		value(d2) = value(d);
2598 		k0 = k;
2599 		ilim0 = ilim;
2600 		ieps = 2; /* conservative */
2601 		if (k > 0) {
2602 			ds = tens[k&0xf];
2603 			j = (unsigned int)k >> 4;
2604 			if (j & Bletch) {
2605 				/* prevent overflows */
2606 				j &= Bletch - 1;
2607 				value(d) /= bigtens[n_bigtens-1];
2608 				ieps++;
2609 				}
2610 			for(; j; j = (unsigned int)j >> 1, i++)
2611 				if (j & 1) {
2612 					ieps++;
2613 					ds *= bigtens[i];
2614 					}
2615 			value(d) /= ds;
2616 		}
2617 		else if ((jj1 = -k) != 0) {
2618 			value(d) *= tens[jj1 & 0xf];
2619 			for(j = (unsigned int)jj1 >> 4; j;
2620 			    j = (unsigned int)j >> 1, i++)
2621 				if (j & 1) {
2622 					ieps++;
2623 					value(d) *= bigtens[i];
2624 				}
2625 		}
2626 		if (k_check && value(d) < 1. && ilim > 0) {
2627 			if (ilim1 <= 0)
2628 				goto fast_failed;
2629 			ilim = ilim1;
2630 			k--;
2631 			value(d) *= 10.;
2632 			ieps++;
2633 		}
2634 		value(eps) = ieps*value(d) + 7.;
2635 		word0(eps) -= (P-1)*Exp_msk1;
2636 		if (ilim == 0) {
2637 			S = mhi = 0;
2638 			value(d) -= 5.;
2639 			if (value(d) > value(eps))
2640 				goto one_digit;
2641 			if (value(d) < -value(eps))
2642 				goto no_digits;
2643 			goto fast_failed;
2644 		}
2645 #ifndef No_leftright
2646 		if (leftright) {
2647 			/* Use Steele & White method of only
2648 			 * generating digits needed.
2649 			 */
2650 			value(eps) = 0.5/tens[ilim-1] - value(eps);
2651 			for(i = 0;;) {
2652 				L = value(d);
2653 				value(d) -= L;
2654 				*s++ = '0' + (int)L;
2655 				if (value(d) < value(eps))
2656 					goto ret1;
2657 				if (1. - value(d) < value(eps))
2658 					goto bump_up;
2659 				if (++i >= ilim)
2660 					break;
2661 				value(eps) *= 10.;
2662 				value(d) *= 10.;
2663 				}
2664 		}
2665 		else {
2666 #endif
2667 			/* Generate ilim digits, then fix them up. */
2668 			value(eps) *= tens[ilim-1];
2669 			for(i = 1;; i++, value(d) *= 10.) {
2670 				L = value(d);
2671 				value(d) -= L;
2672 				*s++ = '0' + (int)L;
2673 				if (i == ilim) {
2674 					if (value(d) > 0.5 + value(eps))
2675 						goto bump_up;
2676 					else if (value(d) < 0.5 - value(eps)) {
2677 						while(*--s == '0');
2678 						s++;
2679 						goto ret1;
2680 						}
2681 					break;
2682 				}
2683 			}
2684 #ifndef No_leftright
2685 		}
2686 #endif
2687  fast_failed:
2688 		s = s0;
2689 		value(d) = value(d2);
2690 		k = k0;
2691 		ilim = ilim0;
2692 	}
2693 
2694 	/* Do we have a "small" integer? */
2695 
2696 	if (be >= 0 && k <= Int_max) {
2697 		/* Yes. */
2698 		ds = tens[k];
2699 		if (ndigits < 0 && ilim <= 0) {
2700 			S = mhi = 0;
2701 			if (ilim < 0 || value(d) <= 5*ds)
2702 				goto no_digits;
2703 			goto one_digit;
2704 		}
2705 		for(i = 1;; i++) {
2706 			L = value(d) / ds;
2707 			value(d) -= L*ds;
2708 #ifdef Check_FLT_ROUNDS
2709 			/* If FLT_ROUNDS == 2, L will usually be high by 1 */
2710 			if (value(d) < 0) {
2711 				L--;
2712 				value(d) += ds;
2713 			}
2714 #endif
2715 			*s++ = '0' + (int)L;
2716 			if (i == ilim) {
2717 				value(d) += value(d);
2718 				if (value(d) > ds || (value(d) == ds && L & 1)) {
2719  bump_up:
2720 					while(*--s == '9')
2721 						if (s == s0) {
2722 							k++;
2723 							*s = '0';
2724 							break;
2725 						}
2726 					++*s++;
2727 				}
2728 				break;
2729 			}
2730 			if (!(value(d) *= 10.))
2731 				break;
2732 			}
2733 		goto ret1;
2734 	}
2735 
2736 	m2 = b2;
2737 	m5 = b5;
2738 	mhi = mlo = 0;
2739 	if (leftright) {
2740 		if (mode < 2) {
2741 			i =
2742 #ifndef Sudden_Underflow
2743 				denorm ? be + (Bias + (P-1) - 1 + 1) :
2744 #endif
2745 #ifdef IBM
2746 				1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
2747 #else
2748 				1 + P - bbits;
2749 #endif
2750 		}
2751 		else {
2752 			j = ilim - 1;
2753 			if (m5 >= j)
2754 				m5 -= j;
2755 			else {
2756 				s5 += j -= m5;
2757 				b5 += j;
2758 				m5 = 0;
2759 			}
2760 			if ((i = ilim) < 0) {
2761 				m2 -= i;
2762 				i = 0;
2763 			}
2764 		}
2765 		b2 += i;
2766 		s2 += i;
2767 		mhi = i2b(1);
2768 	}
2769 	if (m2 > 0 && s2 > 0) {
2770 		i = m2 < s2 ? m2 : s2;
2771 		b2 -= i;
2772 		m2 -= i;
2773 		s2 -= i;
2774 	}
2775 	if (b5 > 0) {
2776 		if (leftright) {
2777 			if (m5 > 0) {
2778 				mhi = pow5mult(mhi, m5);
2779 				b1 = mult(mhi, b);
2780 				Bfree(b);
2781 				b = b1;
2782 			}
2783 			if ((j = b5 - m5) != 0)
2784 				b = pow5mult(b, j);
2785 			}
2786 		else
2787 			b = pow5mult(b, b5);
2788 	}
2789 	S = i2b(1);
2790 	if (s5 > 0)
2791 		S = pow5mult(S, s5);
2792 
2793 	/* Check for special case that d is a normalized power of 2. */
2794 
2795 	if (mode < 2) {
2796 		if (!word1(d) && !(word0(d) & Bndry_mask)
2797 #ifndef Sudden_Underflow
2798 		 && word0(d) & Exp_mask
2799 #endif
2800 				) {
2801 			/* The special case */
2802 			b2 += Log2P;
2803 			s2 += Log2P;
2804 			spec_case = 1;
2805 			}
2806 		else
2807 			spec_case = 0;
2808 	}
2809 
2810 	/* Arrange for convenient computation of quotients:
2811 	 * shift left if necessary so divisor has 4 leading 0 bits.
2812 	 *
2813 	 * Perhaps we should just compute leading 28 bits of S once
2814 	 * and for all and pass them and a shift to quorem, so it
2815 	 * can do shifts and ors to compute the numerator for q.
2816 	 */
2817 	if (S == BIGINT_INVALID) {
2818 		i = 0;
2819 	} else {
2820 #ifdef Pack_32
2821 		if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) != 0)
2822 			i = 32 - i;
2823 #else
2824 		if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf)
2825 			i = 16 - i;
2826 #endif
2827 	}
2828 
2829 	if (i > 4) {
2830 		i -= 4;
2831 		b2 += i;
2832 		m2 += i;
2833 		s2 += i;
2834 	}
2835 	else if (i < 4) {
2836 		i += 28;
2837 		b2 += i;
2838 		m2 += i;
2839 		s2 += i;
2840 	}
2841 	if (b2 > 0)
2842 		b = lshift(b, b2);
2843 	if (s2 > 0)
2844 		S = lshift(S, s2);
2845 	if (k_check) {
2846 		if (cmp(b,S) < 0) {
2847 			k--;
2848 			b = multadd(b, 10, 0);	/* we botched the k estimate */
2849 			if (leftright)
2850 				mhi = multadd(mhi, 10, 0);
2851 			ilim = ilim1;
2852 			}
2853 	}
2854 	if (ilim <= 0 && mode > 2) {
2855 		if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
2856 			/* no digits, fcvt style */
2857  no_digits:
2858 			k = -1 - ndigits;
2859 			goto ret;
2860 		}
2861  one_digit:
2862 		*s++ = '1';
2863 		k++;
2864 		goto ret;
2865 	}
2866 	if (leftright) {
2867 		if (m2 > 0)
2868 			mhi = lshift(mhi, m2);
2869 
2870 		/* Compute mlo -- check for special case
2871 		 * that d is a normalized power of 2.
2872 		 */
2873 
2874 		mlo = mhi;
2875 		if (spec_case) {
2876 			mhi = Balloc(mhi->k);
2877 			Bcopy(mhi, mlo);
2878 			mhi = lshift(mhi, Log2P);
2879 		}
2880 
2881 		for(i = 1;;i++) {
2882 			dig = quorem(b,S) + '0';
2883 			/* Do we yet have the shortest decimal string
2884 			 * that will round to d?
2885 			 */
2886 			j = cmp(b, mlo);
2887 			delta = diff(S, mhi);
2888 			jj1 = delta->sign ? 1 : cmp(b, delta);
2889 			Bfree(delta);
2890 #ifndef ROUND_BIASED
2891 			if (jj1 == 0 && !mode && !(word1(d) & 1)) {
2892 				if (dig == '9')
2893 					goto round_9_up;
2894 				if (j > 0)
2895 					dig++;
2896 				*s++ = dig;
2897 				goto ret;
2898 			}
2899 #endif
2900 			if (j < 0 || (j == 0 && !mode
2901 #ifndef ROUND_BIASED
2902 							&& !(word1(d) & 1)
2903 #endif
2904 					)) {
2905 				if (jj1 > 0) {
2906 					b = lshift(b, 1);
2907 					jj1 = cmp(b, S);
2908 					if ((jj1 > 0 || (jj1 == 0 && dig & 1))
2909 					&& dig++ == '9')
2910 						goto round_9_up;
2911 					}
2912 				*s++ = dig;
2913 				goto ret;
2914 			}
2915 			if (jj1 > 0) {
2916 				if (dig == '9') { /* possible if i == 1 */
2917  round_9_up:
2918 					*s++ = '9';
2919 					goto roundoff;
2920 					}
2921 				*s++ = dig + 1;
2922 				goto ret;
2923 			}
2924 			*s++ = dig;
2925 			if (i == ilim)
2926 				break;
2927 			b = multadd(b, 10, 0);
2928 			if (mlo == mhi)
2929 				mlo = mhi = multadd(mhi, 10, 0);
2930 			else {
2931 				mlo = multadd(mlo, 10, 0);
2932 				mhi = multadd(mhi, 10, 0);
2933 			}
2934 		}
2935 	}
2936 	else
2937 		for(i = 1;; i++) {
2938 			*s++ = dig = quorem(b,S) + '0';
2939 			if (i >= ilim)
2940 				break;
2941 			b = multadd(b, 10, 0);
2942 		}
2943 
2944 	/* Round off last digit */
2945 
2946 	b = lshift(b, 1);
2947 	j = cmp(b, S);
2948 	if (j > 0 || (j == 0 && dig & 1)) {
2949  roundoff:
2950 		while(*--s == '9')
2951 			if (s == s0) {
2952 				k++;
2953 				*s++ = '1';
2954 				goto ret;
2955 				}
2956 		++*s++;
2957 	}
2958 	else {
2959 		while(*--s == '0');
2960 		s++;
2961 	}
2962  ret:
2963 	Bfree(S);
2964 	if (mhi) {
2965 		if (mlo && mlo != mhi)
2966 			Bfree(mlo);
2967 		Bfree(mhi);
2968 	}
2969  ret1:
2970 	Bfree(b);
2971 	if (s == s0) {				/* don't return empty string */
2972 		*s++ = '0';
2973 		k = 0;
2974 	}
2975 	*s = 0;
2976 	*decpt = k + 1;
2977 	if (rve)
2978 		*rve = s;
2979 	return s0;
2980 }
2981 
2982 #include <limits.h>
2983 
2984  char *
rv_alloc(int i)2985 rv_alloc(int i)
2986 {
2987 	int j, k, *r;
2988 
2989 	j = sizeof(ULong);
2990 	for(k = 0;
2991 		sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
2992 		j <<= 1)
2993 			k++;
2994 	r = (int*)Balloc(k);
2995 	*r = k;
2996 	return (char *)(r+1);
2997 	}
2998 
2999  char *
nrv_alloc(char * s,char ** rve,int n)3000 nrv_alloc(char *s, char **rve, int n)
3001 {
3002 	char *rv, *t;
3003 
3004 	t = rv = rv_alloc(n);
3005 	while((*t = *s++) !=0)
3006 		t++;
3007 	if (rve)
3008 		*rve = t;
3009 	return rv;
3010 	}
3011 
3012 
3013 /* Strings values used by dtoa() */
3014 #define	INFSTR	"Infinity"
3015 #define	NANSTR	"NaN"
3016 
3017 #define	DBL_ADJ		(DBL_MAX_EXP - 2 + ((DBL_MANT_DIG - 1) % 4))
3018 #define	LDBL_ADJ	(LDBL_MAX_EXP - 2 + ((LDBL_MANT_DIG - 1) % 4))
3019 
3020 /*
3021  * Round up the given digit string.  If the digit string is fff...f,
3022  * this procedure sets it to 100...0 and returns 1 to indicate that
3023  * the exponent needs to be bumped.  Otherwise, 0 is returned.
3024  */
3025 static int
roundup(char * s0,int ndigits)3026 roundup(char *s0, int ndigits)
3027 {
3028 	char *s;
3029 
3030 	for (s = s0 + ndigits - 1; *s == 0xf; s--) {
3031 		if (s == s0) {
3032 			*s = 1;
3033 			return (1);
3034 		}
3035 		*s = 0;
3036 	}
3037 	++*s;
3038 	return (0);
3039 }
3040 
3041 /*
3042  * Round the given digit string to ndigits digits according to the
3043  * current rounding mode.  Note that this could produce a string whose
3044  * value is not representable in the corresponding floating-point
3045  * type.  The exponent pointed to by decpt is adjusted if necessary.
3046  */
3047 static void
dorounding(char * s0,int ndigits,int sign,int * decpt)3048 dorounding(char *s0, int ndigits, int sign, int *decpt)
3049 {
3050 	int adjust = 0;	/* do we need to adjust the exponent? */
3051 
3052 	switch (FLT_ROUNDS) {
3053 	case 0:		/* toward zero */
3054 	default:	/* implementation-defined */
3055 		break;
3056 	case 1:		/* to nearest, halfway rounds to even */
3057 		if ((s0[ndigits] > 8) ||
3058 		    (s0[ndigits] == 8 && s0[ndigits + 1] & 1))
3059 			adjust = roundup(s0, ndigits);
3060 		break;
3061 	case 2:		/* toward +inf */
3062 		if (sign == 0)
3063 			adjust = roundup(s0, ndigits);
3064 		break;
3065 	case 3:		/* toward -inf */
3066 		if (sign != 0)
3067 			adjust = roundup(s0, ndigits);
3068 		break;
3069 	}
3070 
3071 	if (adjust)
3072 		*decpt += 4;
3073 }
3074 
3075 /*
3076  * This procedure converts a double-precision number in IEEE format
3077  * into a string of hexadecimal digits and an exponent of 2.  Its
3078  * behavior is bug-for-bug compatible with dtoa() in mode 2, with the
3079  * following exceptions:
3080  *
3081  * - An ndigits < 0 causes it to use as many digits as necessary to
3082  *   represent the number exactly.
3083  * - The additional xdigs argument should point to either the string
3084  *   "0123456789ABCDEF" or the string "0123456789abcdef", depending on
3085  *   which case is desired.
3086  * - This routine does not repeat dtoa's mistake of setting decpt
3087  *   to 9999 in the case of an infinity or NaN.  INT_MAX is used
3088  *   for this purpose instead.
3089  *
3090  * Note that the C99 standard does not specify what the leading digit
3091  * should be for non-zero numbers.  For instance, 0x1.3p3 is the same
3092  * as 0x2.6p2 is the same as 0x4.cp3.  This implementation chooses the
3093  * first digit so that subsequent digits are aligned on nibble
3094  * boundaries (before rounding).
3095  *
3096  * Inputs:	d, xdigs, ndigits
3097  * Outputs:	decpt, sign, rve
3098  */
3099 char *
__hdtoa(double d,const char * xdigs,int ndigits,int * decpt,int * sign,char ** rve)3100 __hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign,
3101     char **rve)
3102 {
3103 	static const int sigfigs = (DBL_MANT_DIG + 3) / 4;
3104 	union IEEEd2bits u;
3105 	char *s, *s0;
3106 	int bufsize, f;
3107 
3108 	u.d = d;
3109 	*sign = u.bits.sign;
3110 
3111 	switch (f = fpclassify(d)) {
3112 	case FP_NORMAL:
3113 		*decpt = u.bits.exp - DBL_ADJ;
3114 		break;
3115 	case FP_ZERO:
3116 return_zero:
3117 		*decpt = 1;
3118 		return (nrv_alloc("0", rve, 1));
3119 	case FP_SUBNORMAL:
3120 		/*
3121 		 * For processors that treat subnormals as zero, comparison
3122 		 * with zero will be equal, so we jump to the FP_ZERO case.
3123 		 */
3124 		if(u.d == 0.0) goto return_zero;
3125 		u.d *= 0x1p514;
3126 		*decpt = u.bits.exp - (514 + DBL_ADJ);
3127 		break;
3128 	case FP_INFINITE:
3129 		*decpt = INT_MAX;
3130 		return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
3131 	case FP_NAN:
3132 		*decpt = INT_MAX;
3133 		return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
3134 	default:
3135 #ifdef DEBUG
3136 		BugPrintf("fpclassify returned %d\n", f);
3137 #endif
3138 		return 0; // FIXME??
3139 	}
3140 
3141 	/* FP_NORMAL or FP_SUBNORMAL */
3142 
3143 	if (ndigits == 0)		/* dtoa() compatibility */
3144 		ndigits = 1;
3145 
3146 	/*
3147 	 * For simplicity, we generate all the digits even if the
3148 	 * caller has requested fewer.
3149 	 */
3150 	bufsize = (sigfigs > ndigits) ? sigfigs : ndigits;
3151 	s0 = rv_alloc(bufsize);
3152 
3153 	/*
3154 	 * We work from right to left, first adding any requested zero
3155 	 * padding, then the least significant portion of the
3156 	 * mantissa, followed by the most significant.  The buffer is
3157 	 * filled with the byte values 0x0 through 0xf, which are
3158 	 * converted to xdigs[0x0] through xdigs[0xf] after the
3159 	 * rounding phase.
3160 	 */
3161 	for (s = s0 + bufsize - 1; s > s0 + sigfigs - 1; s--)
3162 		*s = 0;
3163 	for (; s > s0 + sigfigs - (DBL_MANL_SIZE / 4) - 1 && s > s0; s--) {
3164 		*s = u.bits.manl & 0xf;
3165 		u.bits.manl >>= 4;
3166 	}
3167 	for (; s > s0; s--) {
3168 		*s = u.bits.manh & 0xf;
3169 		u.bits.manh >>= 4;
3170 	}
3171 
3172 	/*
3173 	 * At this point, we have snarfed all the bits in the
3174 	 * mantissa, with the possible exception of the highest-order
3175 	 * (partial) nibble, which is dealt with by the next
3176 	 * statement.  We also tack on the implicit normalization bit.
3177 	 */
3178 	*s = u.bits.manh | (1U << ((DBL_MANT_DIG - 1) % 4));
3179 
3180 	/* If ndigits < 0, we are expected to auto-size the precision. */
3181 	if (ndigits < 0) {
3182 		for (ndigits = sigfigs; s0[ndigits - 1] == 0; ndigits--)
3183 			;
3184 	}
3185 
3186 	if (sigfigs > ndigits && s0[ndigits] != 0)
3187 		dorounding(s0, ndigits, u.bits.sign, decpt);
3188 
3189 	s = s0 + ndigits;
3190 	if (rve != NULL)
3191 		*rve = s;
3192 	*s-- = '\0';
3193 	for (; s >= s0; s--)
3194 		*s = xdigs[(unsigned int)*s];
3195 
3196 	return (s0);
3197 }
3198 
3199 #ifndef NO_HEX_FP /*{*/
3200 
3201 static int
gethex(CONST char ** sp,CONST FPI * fpi,Long * exp,Bigint ** bp,int sign,locale_t loc)3202 gethex( CONST char **sp, CONST FPI *fpi, Long *exp, Bigint **bp, int sign, locale_t loc)
3203 {
3204 	Bigint *b;
3205 	CONST unsigned char *decpt, *s0, *s, *s1;
3206 	unsigned char *strunc;
3207 	int big, esign, havedig, irv, j, k, n, n0, nbits, up, zret;
3208 	ULong L, lostbits, *x;
3209 	Long e, e1;
3210 #ifdef USE_LOCALE
3211 	int i;
3212 	NORMALIZE_LOCALE(loc);
3213 #ifdef NO_LOCALE_CACHE
3214 	const unsigned char *decimalpoint = (unsigned char*)localeconv_l(loc)->decimal_point;
3215 #else
3216 	const unsigned char *decimalpoint;
3217 	static unsigned char *decimalpoint_cache;
3218 	if (!(s0 = decimalpoint_cache)) {
3219 		s0 = (unsigned char*)localeconv_l(loc)->decimal_point;
3220 		if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
3221 			strcpy(decimalpoint_cache, s0);
3222 			s0 = decimalpoint_cache;
3223 			}
3224 		}
3225 	decimalpoint = s0;
3226 #endif
3227 #endif
3228 
3229 #ifndef ANDROID_CHANGES
3230 	if (!hexdig['0'])
3231 		hexdig_init_D2A();
3232 #endif
3233 
3234 	*bp = 0;
3235 	havedig = 0;
3236 	s0 = *(CONST unsigned char **)sp + 2;
3237 	while(s0[havedig] == '0')
3238 		havedig++;
3239 	s0 += havedig;
3240 	s = s0;
3241 	decpt = 0;
3242 	zret = 0;
3243 	e = 0;
3244 	if (hexdig[*s])
3245 		havedig++;
3246 	else {
3247 		zret = 1;
3248 #ifdef USE_LOCALE
3249 		for(i = 0; decimalpoint[i]; ++i) {
3250 			if (s[i] != decimalpoint[i])
3251 				goto pcheck;
3252 			}
3253 		decpt = s += i;
3254 #else
3255 		if (*s != '.')
3256 			goto pcheck;
3257 		decpt = ++s;
3258 #endif
3259 		if (!hexdig[*s])
3260 			goto pcheck;
3261 		while(*s == '0')
3262 			s++;
3263 		if (hexdig[*s])
3264 			zret = 0;
3265 		havedig = 1;
3266 		s0 = s;
3267 		}
3268 	while(hexdig[*s])
3269 		s++;
3270 #ifdef USE_LOCALE
3271 	if (*s == *decimalpoint && !decpt) {
3272 		for(i = 1; decimalpoint[i]; ++i) {
3273 			if (s[i] != decimalpoint[i])
3274 				goto pcheck;
3275 			}
3276 		decpt = s += i;
3277 #else
3278 	if (*s == '.' && !decpt) {
3279 		decpt = ++s;
3280 #endif
3281 		while(hexdig[*s])
3282 			s++;
3283 		}/*}*/
3284 	if (decpt)
3285 		e = -(((Long)(s-decpt)) << 2);
3286  pcheck:
3287 	s1 = s;
3288 	big = esign = 0;
3289 	switch(*s) {
3290 	  case 'p':
3291 	  case 'P':
3292 		switch(*++s) {
3293 		  case '-':
3294 			esign = 1;
3295 			/* no break */
3296 		  case '+':
3297 			s++;
3298 		  }
3299 		if ((n = hexdig[*s]) == 0 || n > 0x19) {
3300 			s = s1;
3301 			break;
3302 			}
3303 		e1 = n - 0x10;
3304 		while((n = hexdig[*++s]) !=0 && n <= 0x19) {
3305 			if (e1 & 0xf8000000)
3306 				big = 1;
3307 			e1 = 10*e1 + n - 0x10;
3308 			}
3309 		if (esign)
3310 			e1 = -e1;
3311 		e += e1;
3312 	  }
3313 	*sp = (char*)s;
3314 	if (!havedig)
3315 		*sp = (char*)s0 - 1;
3316 	if (zret)
3317 		return STRTOG_Zero;
3318 	if (big) {
3319 		if (esign) {
3320 			switch(fpi->rounding) {
3321 			  case FPI_Round_up:
3322 				if (sign)
3323 					break;
3324 				goto ret_tiny;
3325 			  case FPI_Round_down:
3326 				if (!sign)
3327 					break;
3328 				goto ret_tiny;
3329 			  }
3330 			goto retz;
3331  ret_tiny:
3332 			b = Balloc(0);
3333 			b->wds = 1;
3334 			b->x[0] = 1;
3335 			goto dret;
3336 			}
3337 		switch(fpi->rounding) {
3338 		  case FPI_Round_near:
3339 			goto ovfl1;
3340 		  case FPI_Round_up:
3341 			if (!sign)
3342 				goto ovfl1;
3343 			goto ret_big;
3344 		  case FPI_Round_down:
3345 			if (sign)
3346 				goto ovfl1;
3347 			goto ret_big;
3348 		  }
3349  ret_big:
3350 		nbits = fpi->nbits;
3351 		n0 = n = nbits >> kshift;
3352 		if (nbits & kmask)
3353 			++n;
3354 		for(j = n, k = 0; j >>= 1; ++k);
3355 		*bp = b = Balloc(k);
3356 		b->wds = n;
3357 		for(j = 0; j < n0; ++j)
3358 			b->x[j] = ALL_ON;
3359 		if (n > n0)
3360 			b->x[j] = ULbits >> (ULbits - (nbits & kmask));
3361 		*exp = fpi->emin;
3362 		return STRTOG_Normal | STRTOG_Inexlo;
3363 		}
3364 	/*
3365 	 * Truncate the hex string if it is longer than the precision needed,
3366 	 * to avoid denial-of-service issues with very large strings.  Use
3367 	 * additional digits to insure precision.  Scan to-be-truncated digits
3368 	 * and replace with either '1' or '0' to ensure proper rounding.
3369 	 */
3370 	{
3371 		int maxdigits = ((fpi->nbits + 3) >> 2) + 2;
3372 		size_t nd = s1 - s0;
3373 #ifdef USE_LOCALE
3374 		int dplen = strlen((const char *)decimalpoint);
3375 #else
3376 		int dplen = 1;
3377 #endif
3378 
3379 		if (decpt && s0 < decpt)
3380 			nd -= dplen;
3381 		if (nd > maxdigits && (strunc = alloca(maxdigits + dplen + 2)) != NULL) {
3382 			ssize_t nd0 = decpt ? decpt - s0 - dplen : nd;
3383 			unsigned char *tp = strunc + maxdigits;
3384 			int found = 0;
3385 			if ((nd0 -= maxdigits) >= 0 || s0 >= decpt)
3386 				memcpy(strunc, s0, maxdigits);
3387 			else {
3388 				memcpy(strunc, s0, maxdigits + dplen);
3389 				tp += dplen;
3390 				}
3391 			s0 += maxdigits;
3392 			e += (nd - (maxdigits + 1)) << 2;
3393 			if (nd0 > 0) {
3394 				while(nd0-- > 0)
3395 					if (*s0++ != '0') {
3396 						found++;
3397 						break;
3398 						}
3399 				s0 += dplen;
3400 				}
3401 			if (!found && decpt) {
3402 				while(s0 < s1)
3403 					if(*s0++ != '0') {
3404 						found++;
3405 						break;
3406 						}
3407 				}
3408 			*tp++ = found ? '1' : '0';
3409 			*tp = 0;
3410 			s0 = strunc;
3411 			s1 = tp;
3412 			}
3413 		}
3414 
3415 	n = s1 - s0 - 1;
3416 	for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
3417 		k++;
3418 	b = Balloc(k);
3419 	x = b->x;
3420 	n = 0;
3421 	L = 0;
3422 #ifdef USE_LOCALE
3423 	for(i = 0; decimalpoint[i+1]; ++i);
3424 #endif
3425 	while(s1 > s0) {
3426 #ifdef USE_LOCALE
3427 		if (*--s1 == decimalpoint[i]) {
3428 			s1 -= i;
3429 			continue;
3430 			}
3431 #else
3432 		if (*--s1 == '.')
3433 			continue;
3434 #endif
3435 		if (n == ULbits) {
3436 			*x++ = L;
3437 			L = 0;
3438 			n = 0;
3439 			}
3440 		L |= (hexdig[*s1] & 0x0f) << n;
3441 		n += 4;
3442 		}
3443 	*x++ = L;
3444 	b->wds = n = x - b->x;
3445 	n = ULbits*n - hi0bits(L);
3446 	nbits = fpi->nbits;
3447 	lostbits = 0;
3448 	x = b->x;
3449 	if (n > nbits) {
3450 		n -= nbits;
3451 		if (any_on(b,n)) {
3452 			lostbits = 1;
3453 			k = n - 1;
3454 			if (x[k>>kshift] & 1 << (k & kmask)) {
3455 				lostbits = 2;
3456 				if (k > 0 && any_on(b,k))
3457 					lostbits = 3;
3458 				}
3459 			}
3460 		rshift(b, n);
3461 		e += n;
3462 		}
3463 	else if (n < nbits) {
3464 		n = nbits - n;
3465 		b = lshift(b, n);
3466 		e -= n;
3467 		x = b->x;
3468 		}
3469 	if (e > fpi->emax) {
3470  ovfl:
3471 		Bfree(b);
3472  ovfl1:
3473 #ifndef NO_ERRNO
3474 		errno = ERANGE;
3475 #endif
3476 		return STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
3477 		}
3478 	irv = STRTOG_Normal;
3479 	if (e < fpi->emin) {
3480 		irv = STRTOG_Denormal;
3481 		n = fpi->emin - e;
3482 		if (n >= nbits) {
3483 			switch (fpi->rounding) {
3484 			  case FPI_Round_near:
3485 				if (n == nbits && (n < 2 || any_on(b,n-1)))
3486 					goto one_bit;
3487 				break;
3488 			  case FPI_Round_up:
3489 				if (!sign)
3490 					goto one_bit;
3491 				break;
3492 			  case FPI_Round_down:
3493 				if (sign) {
3494  one_bit:
3495 					x[0] = b->wds = 1;
3496  dret:
3497 					*bp = b;
3498 					*exp = fpi->emin;
3499 #ifndef NO_ERRNO
3500 					errno = ERANGE;
3501 #endif
3502 					return STRTOG_Denormal | STRTOG_Inexhi
3503 						| STRTOG_Underflow;
3504 					}
3505 			  }
3506 			Bfree(b);
3507  retz:
3508 #ifndef NO_ERRNO
3509 			errno = ERANGE;
3510 #endif
3511 			return STRTOG_Zero | STRTOG_Inexlo | STRTOG_Underflow;
3512 			}
3513 		k = n - 1;
3514 		if (lostbits)
3515 			lostbits = 1;
3516 		else if (k > 0)
3517 			lostbits = any_on(b,k);
3518 		if (x[k>>kshift] & 1 << (k & kmask))
3519 			lostbits |= 2;
3520 		nbits -= n;
3521 		rshift(b,n);
3522 		e = fpi->emin;
3523 		}
3524 	if (lostbits) {
3525 		up = 0;
3526 		switch(fpi->rounding) {
3527 		  case FPI_Round_zero:
3528 			break;
3529 		  case FPI_Round_near:
3530 			if (lostbits & 2
3531 			 && (lostbits | x[0]) & 1)
3532 				up = 1;
3533 			break;
3534 		  case FPI_Round_up:
3535 			up = 1 - sign;
3536 			break;
3537 		  case FPI_Round_down:
3538 			up = sign;
3539 		  }
3540 		if (up) {
3541 			k = b->wds;
3542 			b = increment(b);
3543 			x = b->x;
3544 			if (irv == STRTOG_Denormal) {
3545 				if (nbits == fpi->nbits - 1
3546 				 && x[nbits >> kshift] & 1 << (nbits & kmask))
3547 					irv =  STRTOG_Normal;
3548 				}
3549 			else if (b->wds > k
3550 			 || ((n = nbits & kmask) !=0
3551 			      && hi0bits(x[k-1]) < 32-n)) {
3552 				rshift(b,1);
3553 				if (++e > fpi->emax)
3554 					goto ovfl;
3555 				}
3556 			irv |= STRTOG_Inexhi;
3557 			}
3558 		else
3559 			irv |= STRTOG_Inexlo;
3560 		}
3561 	*bp = b;
3562 	*exp = e;
3563 	return irv;
3564 	}
3565 
3566 #endif /*}*/
3567 
3568 #ifdef __cplusplus
3569 }
3570 #endif
3571