comparison nspr/pr/src/misc/prdtoa.c @ 0:1e5118fa0cb1

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