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