Line data Source code
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 3 : void _PR_InitDtoa(void)
26 : {
27 3 : dtoa_lock[0] = PR_NewLock();
28 3 : dtoa_lock[1] = PR_NewLock();
29 3 : }
30 :
31 0 : void _PR_CleanupDtoa(void)
32 : {
33 0 : PR_DestroyLock(dtoa_lock[0]);
34 0 : dtoa_lock[0] = NULL;
35 0 : PR_DestroyLock(dtoa_lock[1]);
36 0 : dtoa_lock[1] = NULL;
37 :
38 : /* FIXME: deal with freelist and p5s. */
39 0 : }
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 0 : 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 0 : ACQUIRE_DTOA_LOCK(0);
560 : /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
561 : /* but this case seems very unlikely. */
562 0 : if (k <= Kmax && (rv = freelist[k]))
563 0 : freelist[k] = rv->next;
564 : else {
565 0 : x = 1 << k;
566 : #ifdef Omit_Private_Memory
567 : rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
568 : #else
569 0 : len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
570 0 : /sizeof(double);
571 0 : if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
572 0 : rv = (Bigint*)pmem_next;
573 0 : pmem_next += len;
574 : }
575 : else
576 0 : rv = (Bigint*)MALLOC(len*sizeof(double));
577 : #endif
578 0 : rv->k = k;
579 0 : rv->maxwds = x;
580 : }
581 0 : FREE_DTOA_LOCK(0);
582 0 : rv->sign = rv->wds = 0;
583 0 : return rv;
584 : }
585 :
586 : static void
587 0 : Bfree
588 : #ifdef KR_headers
589 : (v) Bigint *v;
590 : #else
591 : (Bigint *v)
592 : #endif
593 : {
594 0 : if (v) {
595 0 : if (v->k > Kmax)
596 : #ifdef FREE
597 : FREE((void*)v);
598 : #else
599 0 : free((void*)v);
600 : #endif
601 : else {
602 0 : ACQUIRE_DTOA_LOCK(0);
603 0 : v->next = freelist[v->k];
604 0 : freelist[v->k] = v;
605 0 : FREE_DTOA_LOCK(0);
606 : }
607 : }
608 0 : }
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 0 : 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 0 : wds = b->wds;
634 0 : x = b->x;
635 0 : i = 0;
636 0 : 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 0 : xi = *x;
645 0 : y = (xi & 0xffff) * m + carry;
646 0 : z = (xi >> 16) * m + (y >> 16);
647 0 : carry = z >> 16;
648 0 : *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 0 : while(++i < wds);
657 0 : if (carry) {
658 0 : if (wds >= b->maxwds) {
659 0 : b1 = Balloc(b->k+1);
660 0 : Bcopy(b1, b);
661 0 : Bfree(b);
662 0 : b = b1;
663 : }
664 0 : b->x[wds++] = carry;
665 0 : b->wds = wds;
666 : }
667 0 : return b;
668 : }
669 :
670 : static Bigint *
671 0 : 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 0 : x = (nd + 8) / 9;
683 0 : for(k = 0, y = 1; x > y; y <<= 1, k++) ;
684 : #ifdef Pack_32
685 0 : b = Balloc(k);
686 0 : b->x[0] = y9;
687 0 : 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 0 : i = 9;
695 0 : if (9 < nd0) {
696 0 : s += 9;
697 0 : do b = multadd(b, 10, *s++ - '0');
698 0 : while(++i < nd0);
699 0 : s++;
700 : }
701 : else
702 0 : s += 10;
703 0 : for(; i < nd; i++)
704 0 : b = multadd(b, 10, *s++ - '0');
705 0 : return b;
706 : }
707 :
708 : static int
709 0 : 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 0 : 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 0 : 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 0 : ULong x = *y;
757 :
758 0 : if (x>1)
759 0 : *y = ( x >> (k = pr_bitscan_ctz32(x)) );
760 : else
761 0 : 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 0 : return k;
802 : }
803 :
804 : static Bigint *
805 0 : i2b
806 : #ifdef KR_headers
807 : (i) int i;
808 : #else
809 : (int i)
810 : #endif
811 : {
812 : Bigint *b;
813 :
814 0 : b = Balloc(1);
815 0 : b->x[0] = i;
816 0 : b->wds = 1;
817 0 : return b;
818 : }
819 :
820 : static Bigint *
821 0 : 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 0 : if (a->wds < b->wds) {
842 0 : c = a;
843 0 : a = b;
844 0 : b = c;
845 : }
846 0 : k = a->k;
847 0 : wa = a->wds;
848 0 : wb = b->wds;
849 0 : wc = wa + wb;
850 0 : if (wc > a->maxwds)
851 0 : k++;
852 0 : c = Balloc(k);
853 0 : for(x = c->x, xa = x + wc; x < xa; x++)
854 0 : *x = 0;
855 0 : xa = a->x;
856 0 : xae = xa + wa;
857 0 : xb = b->x;
858 0 : xbe = xb + wb;
859 0 : 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 0 : for(; xb < xbe; xb++, xc0++) {
878 0 : if (y = *xb & 0xffff) {
879 0 : x = xa;
880 0 : xc = xc0;
881 0 : carry = 0;
882 : do {
883 0 : z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
884 0 : carry = z >> 16;
885 0 : z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
886 0 : carry = z2 >> 16;
887 0 : Storeinc(xc, z2, z);
888 : }
889 0 : while(x < xae);
890 0 : *xc = carry;
891 : }
892 0 : if (y = *xb >> 16) {
893 0 : x = xa;
894 0 : xc = xc0;
895 0 : carry = 0;
896 0 : z2 = *xc;
897 : do {
898 0 : z = (*x & 0xffff) * y + (*xc >> 16) + carry;
899 0 : carry = z >> 16;
900 0 : Storeinc(xc, z, z2);
901 0 : z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
902 0 : carry = z2 >> 16;
903 : }
904 0 : while(x < xae);
905 0 : *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 0 : for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
926 0 : c->wds = wc;
927 0 : return c;
928 : }
929 :
930 : static Bigint *p5s;
931 :
932 : static Bigint *
933 0 : 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 0 : if (i = k & 3)
945 0 : b = multadd(b, p05[i-1], 0);
946 :
947 0 : if (!(k >>= 2))
948 0 : return b;
949 0 : if (!(p5 = p5s)) {
950 : /* first time */
951 : #ifdef MULTIPLE_THREADS
952 0 : ACQUIRE_DTOA_LOCK(1);
953 0 : if (!(p5 = p5s)) {
954 0 : p5 = p5s = i2b(625);
955 0 : p5->next = 0;
956 : }
957 0 : FREE_DTOA_LOCK(1);
958 : #else
959 : p5 = p5s = i2b(625);
960 : p5->next = 0;
961 : #endif
962 : }
963 : for(;;) {
964 0 : if (k & 1) {
965 0 : b1 = mult(b, p5);
966 0 : Bfree(b);
967 0 : b = b1;
968 : }
969 0 : if (!(k >>= 1))
970 0 : break;
971 0 : if (!(p51 = p5->next)) {
972 : #ifdef MULTIPLE_THREADS
973 0 : ACQUIRE_DTOA_LOCK(1);
974 0 : if (!(p51 = p5->next)) {
975 0 : p51 = p5->next = mult(p5,p5);
976 0 : p51->next = 0;
977 : }
978 0 : FREE_DTOA_LOCK(1);
979 : #else
980 : p51 = p5->next = mult(p5,p5);
981 : p51->next = 0;
982 : #endif
983 : }
984 0 : p5 = p51;
985 : }
986 0 : return b;
987 : }
988 :
989 : static Bigint *
990 0 : 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 0 : n = k >> 5;
1003 : #else
1004 : n = k >> 4;
1005 : #endif
1006 0 : k1 = b->k;
1007 0 : n1 = n + b->wds + 1;
1008 0 : for(i = b->maxwds; n1 > i; i <<= 1)
1009 0 : k1++;
1010 0 : b1 = Balloc(k1);
1011 0 : x1 = b1->x;
1012 0 : for(i = 0; i < n; i++)
1013 0 : *x1++ = 0;
1014 0 : x = b->x;
1015 0 : xe = x + b->wds;
1016 : #ifdef Pack_32
1017 0 : if (k &= 0x1f) {
1018 0 : k1 = 32 - k;
1019 0 : z = 0;
1020 : do {
1021 0 : *x1++ = *x << k | z;
1022 0 : z = *x++ >> k1;
1023 : }
1024 0 : while(x < xe);
1025 0 : if (*x1 = z)
1026 0 : ++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 0 : *x1++ = *x++;
1043 0 : while(x < xe);
1044 0 : b1->wds = n1 - 1;
1045 0 : Bfree(b);
1046 0 : return b1;
1047 : }
1048 :
1049 : static int
1050 0 : 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 0 : i = a->wds;
1061 0 : j = b->wds;
1062 : #ifdef DEBUG
1063 0 : if (i > 1 && !a->x[i-1])
1064 0 : Bug("cmp called with a->x[a->wds-1] == 0");
1065 0 : if (j > 1 && !b->x[j-1])
1066 0 : Bug("cmp called with b->x[b->wds-1] == 0");
1067 : #endif
1068 0 : if (i -= j)
1069 0 : return i;
1070 0 : xa0 = a->x;
1071 0 : xa = xa0 + j;
1072 0 : xb0 = b->x;
1073 0 : xb = xb0 + j;
1074 : for(;;) {
1075 0 : if (*--xa != *--xb)
1076 0 : return *xa < *xb ? -1 : 1;
1077 0 : if (xa <= xa0)
1078 0 : break;
1079 : }
1080 0 : return 0;
1081 : }
1082 :
1083 : static Bigint *
1084 0 : 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 0 : i = cmp(a,b);
1104 0 : if (!i) {
1105 0 : c = Balloc(0);
1106 0 : c->wds = 1;
1107 0 : c->x[0] = 0;
1108 0 : return c;
1109 : }
1110 0 : if (i < 0) {
1111 0 : c = a;
1112 0 : a = b;
1113 0 : b = c;
1114 0 : i = 1;
1115 : }
1116 : else
1117 0 : i = 0;
1118 0 : c = Balloc(a->k);
1119 0 : c->sign = i;
1120 0 : wa = a->wds;
1121 0 : xa = a->x;
1122 0 : xae = xa + wa;
1123 0 : wb = b->wds;
1124 0 : xb = b->x;
1125 0 : xbe = xb + wb;
1126 0 : xc = c->x;
1127 0 : 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 0 : y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
1144 0 : borrow = (y & 0x10000) >> 16;
1145 0 : z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
1146 0 : borrow = (z & 0x10000) >> 16;
1147 0 : Storeinc(xc, z, y);
1148 : }
1149 0 : while(xb < xbe);
1150 0 : while(xa < xae) {
1151 0 : y = (*xa & 0xffff) - borrow;
1152 0 : borrow = (y & 0x10000) >> 16;
1153 0 : z = (*xa++ >> 16) - borrow;
1154 0 : borrow = (z & 0x10000) >> 16;
1155 0 : 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 0 : while(!*--xc)
1172 0 : wa--;
1173 0 : c->wds = wa;
1174 0 : return c;
1175 : }
1176 :
1177 : static double
1178 0 : 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 0 : dval(x) = dx;
1189 0 : 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 0 : word0(a) = L;
1199 0 : 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 0 : return dval(a);
1218 : }
1219 :
1220 : static double
1221 0 : 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 0 : xa0 = a->x;
1239 0 : xa = xa0 + a->wds;
1240 0 : y = *--xa;
1241 : #ifdef DEBUG
1242 0 : if (!y) Bug("zero y in b2d");
1243 : #endif
1244 0 : k = hi0bits(y);
1245 0 : *e = 32 - k;
1246 : #ifdef Pack_32
1247 0 : if (k < Ebits) {
1248 0 : d0 = Exp_1 | y >> Ebits - k;
1249 0 : w = xa > xa0 ? *--xa : 0;
1250 0 : d1 = y << (32-Ebits) + k | w >> Ebits - k;
1251 0 : goto ret_d;
1252 : }
1253 0 : z = xa > xa0 ? *--xa : 0;
1254 0 : if (k -= Ebits) {
1255 0 : d0 = Exp_1 | y << k | z >> 32 - k;
1256 0 : y = xa > xa0 ? *--xa : 0;
1257 0 : d1 = z << k | y >> 32 - k;
1258 : }
1259 : else {
1260 0 : d0 = Exp_1 | y;
1261 0 : 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 0 : return dval(d);
1288 : }
1289 :
1290 : static Bigint *
1291 0 : 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 0 : 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 0 : b = Balloc(1);
1320 : #else
1321 : b = Balloc(2);
1322 : #endif
1323 0 : x = b->x;
1324 :
1325 0 : z = d0 & Frac_mask;
1326 0 : 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 0 : if (de = (int)(d0 >> Exp_shift))
1334 0 : z |= Exp_msk1;
1335 : #endif
1336 : #ifdef Pack_32
1337 0 : if (y = d1) {
1338 0 : if (k = lo0bits(&y)) {
1339 0 : x[0] = y | z << 32 - k;
1340 0 : z >>= k;
1341 : }
1342 : else
1343 0 : x[0] = y;
1344 : #ifndef Sudden_Underflow
1345 0 : i =
1346 : #endif
1347 0 : b->wds = (x[1] = z) ? 2 : 1;
1348 : }
1349 : else {
1350 0 : k = lo0bits(&z);
1351 0 : x[0] = z;
1352 : #ifndef Sudden_Underflow
1353 0 : i =
1354 : #endif
1355 0 : b->wds = 1;
1356 0 : 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 0 : 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 0 : *e = de - Bias - (P-1) + k;
1411 0 : *bits = P - k;
1412 : #endif
1413 : #ifndef Sudden_Underflow
1414 : }
1415 : else {
1416 0 : *e = de - Bias - (P-1) + 1 + k;
1417 : #ifdef Pack_32
1418 0 : *bits = 32*i - hi0bits(x[i-1]);
1419 : #else
1420 : *bits = (i+2)*16 - hi0bits(x[i]);
1421 : #endif
1422 : }
1423 : #endif
1424 0 : return b;
1425 : }
1426 : #undef d0
1427 : #undef d1
1428 :
1429 : static double
1430 0 : 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 0 : dval(da) = b2d(a, &ka);
1441 0 : dval(db) = b2d(b, &kb);
1442 : #ifdef Pack_32
1443 0 : 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 0 : if (k > 0)
1461 0 : word0(da) += k*Exp_msk1;
1462 : else {
1463 0 : k = -k;
1464 0 : word0(db) += k*Exp_msk1;
1465 : }
1466 : #endif
1467 0 : 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 131 : if (!_pr_initialized) _PR_ImplicitInitialization();
1627 :
1628 131 : sign = nz0 = nz = 0;
1629 131 : dval(rv) = 0.;
1630 131 : for(s = s00;;s++) switch(*s) {
1631 : case '-':
1632 8 : sign = 1;
1633 : /* no break */
1634 : case '+':
1635 8 : if (*++s)
1636 8 : goto break2;
1637 : /* no break */
1638 : case 0:
1639 0 : goto ret0;
1640 : case '\t':
1641 : case '\n':
1642 : case '\v':
1643 : case '\f':
1644 : case '\r':
1645 : case ' ':
1646 0 : continue;
1647 : default:
1648 123 : goto break2;
1649 : }
1650 : break2:
1651 131 : if (*s == '0') {
1652 68 : nz0 = 1;
1653 68 : while(*++s == '0') ;
1654 68 : if (!*s)
1655 0 : goto ret;
1656 : }
1657 131 : s0 = s;
1658 131 : y = z = 0;
1659 200 : for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1660 69 : if (nd < 9)
1661 69 : y = 10*y + c - '0';
1662 0 : else if (nd < 16)
1663 0 : z = 10*z + c - '0';
1664 131 : 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 131 : if (c == '.') {
1685 131 : c = *++s;
1686 131 : if (!nd) {
1687 113 : for(; c == '0'; c = *++s)
1688 45 : nz++;
1689 68 : if (c > '0' && c <= '9') {
1690 53 : s0 = s;
1691 53 : nf += nz;
1692 53 : nz = 0;
1693 53 : goto have_dig;
1694 : }
1695 15 : goto dig_done;
1696 : }
1697 272 : for(; c >= '0' && c <= '9'; c = *++s) {
1698 : have_dig:
1699 209 : nz++;
1700 209 : if (c -= '0') {
1701 164 : nf += nz;
1702 167 : for(i = 1; i < nz; i++)
1703 3 : if (nd++ < 9)
1704 3 : y *= 10;
1705 0 : else if (nd <= DBL_DIG + 1)
1706 0 : z *= 10;
1707 164 : if (nd++ < 9)
1708 164 : y = 10*y + c;
1709 0 : else if (nd <= DBL_DIG + 1)
1710 0 : z = 10*z + c;
1711 164 : nz = 0;
1712 : }
1713 : }
1714 : }
1715 : dig_done:
1716 131 : if (nd > 64 * 1024)
1717 0 : goto ret0;
1718 131 : e = 0;
1719 131 : if (c == 'e' || c == 'E') {
1720 0 : if (!nd && !nz && !nz0) {
1721 0 : goto ret0;
1722 : }
1723 0 : s00 = s;
1724 0 : esign = 0;
1725 0 : switch(c = *++s) {
1726 : case '-':
1727 0 : esign = 1;
1728 : case '+':
1729 0 : c = *++s;
1730 : }
1731 0 : if (c >= '0' && c <= '9') {
1732 0 : while(c == '0')
1733 0 : c = *++s;
1734 0 : if (c > '0' && c <= '9') {
1735 0 : L = c - '0';
1736 0 : s1 = s;
1737 0 : while((c = *++s) >= '0' && c <= '9')
1738 0 : L = 10*L + c - '0';
1739 0 : if (s - s1 > 8 || L > 19999)
1740 : /* Avoid confusion from exponents
1741 : * so large that e might overflow.
1742 : */
1743 0 : e = 19999; /* safe for 16 bit ints */
1744 : else
1745 0 : e = (int)L;
1746 0 : if (esign)
1747 0 : e = -e;
1748 : }
1749 : else
1750 0 : e = 0;
1751 : }
1752 : else
1753 0 : s = s00;
1754 : }
1755 131 : if (!nd) {
1756 15 : 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 0 : s = s00;
1786 0 : sign = 0;
1787 : }
1788 15 : goto ret;
1789 : }
1790 116 : 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 116 : if (!nd0)
1798 53 : nd0 = nd;
1799 116 : k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
1800 116 : dval(rv) = y;
1801 116 : if (k > 9) {
1802 : #ifdef SET_INEXACT
1803 : if (k > DBL_DIG)
1804 : oldinexact = get_inexact();
1805 : #endif
1806 0 : dval(rv) = tens[k - 9] * dval(rv) + z;
1807 : }
1808 116 : bd0 = 0;
1809 116 : if (nd <= DBL_DIG
1810 : #ifndef RND_PRODQUOT
1811 : #ifndef Honor_FLT_ROUNDS
1812 : && Flt_Rounds == 1
1813 : #endif
1814 : #endif
1815 : ) {
1816 116 : if (!e)
1817 39 : goto ret;
1818 77 : if (e > 0) {
1819 0 : 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 0 : /* rv = */ rounded_product(dval(rv), tens[e]);
1831 0 : goto ret;
1832 : #endif
1833 : }
1834 0 : i = DBL_DIG - nd;
1835 0 : 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 0 : e -= i;
1847 0 : 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 0 : /* rv = */ rounded_product(dval(rv), tens[e]);
1861 : #endif
1862 0 : goto ret;
1863 : }
1864 : }
1865 : #ifndef Inaccurate_Divide
1866 77 : 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 77 : /* rv = */ rounded_quotient(dval(rv), tens[-e]);
1875 77 : goto ret;
1876 : }
1877 : #endif
1878 : }
1879 0 : 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 0 : 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 0 : if (e1 > 0) {
1904 0 : if (i = e1 & 15)
1905 0 : dval(rv) *= tens[i];
1906 0 : if (e1 &= ~15) {
1907 0 : if (e1 > DBL_MAX_10_EXP) {
1908 : ovfl:
1909 : #ifndef NO_ERRNO
1910 0 : 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 0 : word0(rv) = Exp_mask;
1927 0 : 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 0 : if (bd0)
1939 0 : goto retfree;
1940 0 : goto ret;
1941 : }
1942 0 : e1 >>= 4;
1943 0 : for(j = 0; e1 > 1; j++, e1 >>= 1)
1944 0 : if (e1 & 1)
1945 0 : dval(rv) *= bigtens[j];
1946 : /* The last multiplication could overflow. */
1947 0 : word0(rv) -= P*Exp_msk1;
1948 0 : dval(rv) *= bigtens[j];
1949 0 : if ((z = word0(rv) & Exp_mask)
1950 : > Exp_msk1*(DBL_MAX_EXP+Bias-P))
1951 0 : goto ovfl;
1952 0 : if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
1953 : /* set to largest number */
1954 : /* (Can't trust DBL_MAX) */
1955 0 : word0(rv) = Big0;
1956 0 : word1(rv) = Big1;
1957 : }
1958 : else
1959 0 : word0(rv) += P*Exp_msk1;
1960 : }
1961 : }
1962 0 : else if (e1 < 0) {
1963 0 : e1 = -e1;
1964 0 : if (i = e1 & 15)
1965 0 : dval(rv) /= tens[i];
1966 0 : if (e1 >>= 4) {
1967 0 : if (e1 >= 1 << n_bigtens)
1968 0 : goto undfl;
1969 : #ifdef Avoid_Underflow
1970 0 : if (e1 & Scale_Bit)
1971 0 : scale = 2*P;
1972 0 : for(j = 0; e1 > 0; j++, e1 >>= 1)
1973 0 : if (e1 & 1)
1974 0 : dval(rv) *= tinytens[j];
1975 0 : if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
1976 0 : >> Exp_shift)) > 0) {
1977 : /* scaled rv is denormal; zap j low bits */
1978 0 : if (j >= 32) {
1979 0 : word1(rv) = 0;
1980 0 : if (j >= 53)
1981 0 : word0(rv) = (P+2)*Exp_msk1;
1982 : else
1983 0 : word0(rv) &= 0xffffffff << j-32;
1984 : }
1985 : else
1986 0 : 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 0 : if (!dval(rv)) {
2000 : undfl:
2001 0 : dval(rv) = 0.;
2002 : #ifndef NO_ERRNO
2003 0 : PR_SetError(PR_RANGE_ERROR, 0);
2004 : #endif
2005 0 : if (bd0)
2006 0 : goto retfree;
2007 0 : 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 0 : bd0 = s2b(s0, nd0, nd, y);
2025 :
2026 : for(;;) {
2027 0 : bd = Balloc(bd0->k);
2028 0 : Bcopy(bd, bd0);
2029 0 : bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
2030 0 : bs = i2b(1);
2031 :
2032 0 : if (e >= 0) {
2033 0 : bb2 = bb5 = 0;
2034 0 : bd2 = bd5 = e;
2035 : }
2036 : else {
2037 0 : bb2 = bb5 = -e;
2038 0 : bd2 = bd5 = 0;
2039 : }
2040 0 : if (bbe >= 0)
2041 0 : bb2 += bbe;
2042 : else
2043 0 : bd2 -= bbe;
2044 0 : bs2 = bb2;
2045 : #ifdef Honor_FLT_ROUNDS
2046 : if (rounding != 1)
2047 : bs2++;
2048 : #endif
2049 : #ifdef Avoid_Underflow
2050 0 : j = bbe - scale;
2051 0 : i = j + bbbits - 1; /* logb(rv) */
2052 0 : if (i < Emin) /* denormal */
2053 0 : j += P - Emin;
2054 : else
2055 0 : 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 0 : bb2 += j;
2073 0 : bd2 += j;
2074 : #ifdef Avoid_Underflow
2075 0 : bd2 += scale;
2076 : #endif
2077 0 : i = bb2 < bd2 ? bb2 : bd2;
2078 0 : if (i > bs2)
2079 0 : i = bs2;
2080 0 : if (i > 0) {
2081 0 : bb2 -= i;
2082 0 : bd2 -= i;
2083 0 : bs2 -= i;
2084 : }
2085 0 : if (bb5 > 0) {
2086 0 : bs = pow5mult(bs, bb5);
2087 0 : bb1 = mult(bs, bb);
2088 0 : Bfree(bb);
2089 0 : bb = bb1;
2090 : }
2091 0 : if (bb2 > 0)
2092 0 : bb = lshift(bb, bb2);
2093 0 : if (bd5 > 0)
2094 0 : bd = pow5mult(bd, bd5);
2095 0 : if (bd2 > 0)
2096 0 : bd = lshift(bd, bd2);
2097 0 : if (bs2 > 0)
2098 0 : bs = lshift(bs, bs2);
2099 0 : delta = diff(bb, bd);
2100 0 : dsign = delta->sign;
2101 0 : delta->sign = 0;
2102 0 : 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 0 : if (i < 0) {
2195 : /* Error is less than half an ulp -- check for
2196 : * special case of mantissa a power of two.
2197 : */
2198 0 : if (dsign || word1(rv) || word0(rv) & Bndry_mask
2199 : #ifdef IEEE_Arith
2200 : #ifdef Avoid_Underflow
2201 0 : || (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 0 : if (!delta->x[0] && delta->wds <= 1) {
2214 : /* exact result */
2215 : #ifdef SET_INEXACT
2216 : inexact = 0;
2217 : #endif
2218 0 : break;
2219 : }
2220 0 : delta = lshift(delta,Log2P);
2221 0 : if (cmp(delta, bs) > 0)
2222 0 : goto drop_down;
2223 0 : break;
2224 : }
2225 0 : if (i == 0) {
2226 : /* exactly half-way between */
2227 0 : if (dsign) {
2228 0 : if ((word0(rv) & Bndry_mask1) == Bndry_mask1
2229 0 : && word1(rv) == (
2230 : #ifdef Avoid_Underflow
2231 0 : (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
2232 0 : ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
2233 : #endif
2234 : 0xffffffff)) {
2235 : /*boundary case -- increment exponent*/
2236 0 : word0(rv) = (word0(rv) & Exp_mask)
2237 0 : + Exp_msk1
2238 : #ifdef IBM
2239 : | Exp_msk1 >> 4
2240 : #endif
2241 : ;
2242 0 : word1(rv) = 0;
2243 : #ifdef Avoid_Underflow
2244 0 : dsign = 0;
2245 : #endif
2246 0 : break;
2247 : }
2248 : }
2249 0 : 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 0 : if (scale) {
2268 0 : L = word0(rv) & Exp_mask;
2269 0 : if (L <= (2*P+1)*Exp_msk1) {
2270 0 : if (L > (P+2)*Exp_msk1)
2271 : /* round even ==> */
2272 : /* accept rv */
2273 0 : break;
2274 : /* rv = smallest denormal */
2275 0 : goto undfl;
2276 : }
2277 : }
2278 : #endif /*Avoid_Underflow*/
2279 0 : L = (word0(rv) & Exp_mask) - Exp_msk1;
2280 : #endif /*Sudden_Underflow}}*/
2281 0 : word0(rv) = L | Bndry_mask1;
2282 0 : word1(rv) = 0xffffffff;
2283 : #ifdef IBM
2284 : goto cont;
2285 : #else
2286 0 : break;
2287 : #endif
2288 : }
2289 : #ifndef ROUND_BIASED
2290 0 : if (!(word1(rv) & LSB))
2291 0 : break;
2292 : #endif
2293 0 : if (dsign)
2294 0 : dval(rv) += ulp(dval(rv));
2295 : #ifndef ROUND_BIASED
2296 : else {
2297 0 : dval(rv) -= ulp(dval(rv));
2298 : #ifndef Sudden_Underflow
2299 0 : if (!dval(rv))
2300 0 : goto undfl;
2301 : #endif
2302 : }
2303 : #ifdef Avoid_Underflow
2304 0 : dsign = 1 - dsign;
2305 : #endif
2306 : #endif
2307 0 : break;
2308 : }
2309 0 : if ((aadj = ratio(delta, bs)) <= 2.) {
2310 0 : if (dsign)
2311 0 : aadj = aadj1 = 1.;
2312 0 : else if (word1(rv) || word0(rv) & Bndry_mask) {
2313 : #ifndef Sudden_Underflow
2314 0 : if (word1(rv) == Tiny1 && !word0(rv))
2315 0 : goto undfl;
2316 : #endif
2317 0 : aadj = 1.;
2318 0 : aadj1 = -1.;
2319 : }
2320 : else {
2321 : /* special case -- power of FLT_RADIX to be */
2322 : /* rounded down... */
2323 :
2324 0 : if (aadj < 2./FLT_RADIX)
2325 0 : aadj = 1./FLT_RADIX;
2326 : else
2327 0 : aadj *= 0.5;
2328 0 : aadj1 = -aadj;
2329 : }
2330 : }
2331 : else {
2332 0 : aadj *= 0.5;
2333 0 : 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 0 : y = word0(rv) & Exp_mask;
2349 :
2350 : /* Check for overflow */
2351 :
2352 0 : if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
2353 0 : dval(rv0) = dval(rv);
2354 0 : word0(rv) -= P*Exp_msk1;
2355 0 : adj = aadj1 * ulp(dval(rv));
2356 0 : dval(rv) += adj;
2357 0 : if ((word0(rv) & Exp_mask) >=
2358 : Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
2359 0 : if (word0(rv0) == Big0 && word1(rv0) == Big1)
2360 0 : goto ovfl;
2361 0 : word0(rv) = Big0;
2362 0 : word1(rv) = Big1;
2363 0 : goto cont;
2364 : }
2365 : else
2366 0 : word0(rv) += P*Exp_msk1;
2367 : }
2368 : else {
2369 : #ifdef Avoid_Underflow
2370 0 : if (scale && y <= 2*P*Exp_msk1) {
2371 0 : if (aadj <= 0x7fffffff) {
2372 0 : if ((z = aadj) <= 0)
2373 0 : z = 1;
2374 0 : aadj = z;
2375 0 : aadj1 = dsign ? aadj : -aadj;
2376 : }
2377 0 : dval(aadj2) = aadj1;
2378 0 : word0(aadj2) += (2*P+1)*Exp_msk1 - y;
2379 0 : aadj1 = dval(aadj2);
2380 : }
2381 0 : adj = aadj1 * ulp(dval(rv));
2382 0 : 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 0 : z = word0(rv) & Exp_mask;
2429 : #ifndef SET_INEXACT
2430 : #ifdef Avoid_Underflow
2431 0 : if (!scale)
2432 : #endif
2433 0 : if (y == z) {
2434 : /* Can we stop now? */
2435 0 : L = (Long)aadj;
2436 0 : aadj -= L;
2437 : /* The tolerances below are conservative. */
2438 0 : if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
2439 0 : if (aadj < .4999999 || aadj > .5000001)
2440 : break;
2441 : }
2442 0 : else if (aadj < .4999999/FLT_RADIX)
2443 0 : break;
2444 : }
2445 : #endif
2446 : cont:
2447 0 : Bfree(bb);
2448 0 : Bfree(bd);
2449 0 : Bfree(bs);
2450 0 : 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 0 : if (scale) {
2465 0 : word0(rv0) = Exp_1 - 2*P*Exp_msk1;
2466 0 : word1(rv0) = 0;
2467 0 : dval(rv) *= dval(rv0);
2468 : #ifndef NO_ERRNO
2469 : /* try to avoid the bug of testing an 8087 register value */
2470 0 : if (word0(rv) == 0 && word1(rv) == 0)
2471 0 : 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 0 : Bfree(bb);
2484 0 : Bfree(bd);
2485 0 : Bfree(bs);
2486 0 : Bfree(bd0);
2487 0 : Bfree(delta);
2488 : ret:
2489 131 : if (se)
2490 128 : *se = (char *)s;
2491 131 : return sign ? -dval(rv) : dval(rv);
2492 : }
2493 :
2494 : static int
2495 0 : 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 0 : n = S->wds;
2514 : #ifdef DEBUG
2515 0 : /*debug*/ if (b->wds > n)
2516 0 : /*debug*/ Bug("oversize b in quorem");
2517 : #endif
2518 0 : if (b->wds < n)
2519 0 : return 0;
2520 0 : sx = S->x;
2521 0 : sxe = sx + --n;
2522 0 : bx = b->x;
2523 0 : bxe = bx + n;
2524 0 : q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
2525 : #ifdef DEBUG
2526 0 : /*debug*/ if (q > 9)
2527 0 : /*debug*/ Bug("oversized quotient in quorem");
2528 : #endif
2529 0 : if (q) {
2530 0 : borrow = 0;
2531 0 : 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 0 : si = *sx++;
2542 0 : ys = (si & 0xffff) * q + carry;
2543 0 : zs = (si >> 16) * q + (ys >> 16);
2544 0 : carry = zs >> 16;
2545 0 : y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2546 0 : borrow = (y & 0x10000) >> 16;
2547 0 : z = (*bx >> 16) - (zs & 0xffff) - borrow;
2548 0 : borrow = (z & 0x10000) >> 16;
2549 0 : 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 0 : while(sx <= sxe);
2560 0 : if (!*bxe) {
2561 0 : bx = b->x;
2562 0 : while(--bxe > bx && !*bxe)
2563 0 : --n;
2564 0 : b->wds = n;
2565 : }
2566 : }
2567 0 : if (cmp(b, S) >= 0) {
2568 0 : q++;
2569 0 : borrow = 0;
2570 0 : carry = 0;
2571 0 : bx = b->x;
2572 0 : 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 0 : si = *sx++;
2583 0 : ys = (si & 0xffff) + carry;
2584 0 : zs = (si >> 16) + (ys >> 16);
2585 0 : carry = zs >> 16;
2586 0 : y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
2587 0 : borrow = (y & 0x10000) >> 16;
2588 0 : z = (*bx >> 16) - (zs & 0xffff) - borrow;
2589 0 : borrow = (z & 0x10000) >> 16;
2590 0 : 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 0 : while(sx <= sxe);
2601 0 : bx = b->x;
2602 0 : bxe = bx + n;
2603 0 : if (!*bxe) {
2604 0 : while(--bxe > bx && !*bxe)
2605 0 : --n;
2606 0 : b->wds = n;
2607 : }
2608 : }
2609 0 : 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 0 : rv_alloc(int i)
2621 : #endif
2622 : {
2623 : int j, k, *r;
2624 :
2625 0 : j = sizeof(ULong);
2626 0 : for(k = 0;
2627 0 : sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
2628 0 : j <<= 1)
2629 0 : k++;
2630 0 : r = (int*)Balloc(k);
2631 0 : *r = k;
2632 : return
2633 : #ifndef MULTIPLE_THREADS
2634 : dtoa_result =
2635 : #endif
2636 0 : (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 0 : nrv_alloc(char *s, char **rve, int n)
2644 : #endif
2645 : {
2646 : char *rv, *t;
2647 :
2648 0 : t = rv = rv_alloc(n);
2649 0 : while(*t = *s++) t++;
2650 0 : if (rve)
2651 0 : *rve = t;
2652 0 : 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 0 : freedtoa(char *s)
2666 : #endif
2667 : {
2668 0 : Bigint *b = (Bigint *)((int *)s - 1);
2669 0 : b->maxwds = 1 << (b->k = *(int*)b);
2670 0 : Bfree(b);
2671 : #ifndef MULTIPLE_THREADS
2672 : if (s == dtoa_result)
2673 : dtoa_result = 0;
2674 : #endif
2675 0 : }
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 0 : 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 0 : dval(d) = dd;
2781 0 : if (word0(d) & Sign_bit) {
2782 : /* set sign for everything, including 0's and NaNs */
2783 0 : *sign = 1;
2784 0 : word0(d) &= ~Sign_bit; /* clear sign bit */
2785 : }
2786 : else
2787 0 : *sign = 0;
2788 :
2789 : #if defined(IEEE_Arith) + defined(VAX)
2790 : #ifdef IEEE_Arith
2791 0 : if ((word0(d) & Exp_mask) == Exp_mask)
2792 : #else
2793 : if (word0(d) == 0x8000)
2794 : #endif
2795 : {
2796 : /* Infinity or NaN */
2797 0 : *decpt = 9999;
2798 : #ifdef IEEE_Arith
2799 0 : if (!word1(d) && !(word0(d) & 0xfffff))
2800 0 : return nrv_alloc("Infinity", rve, 8);
2801 : #endif
2802 0 : return nrv_alloc("NaN", rve, 3);
2803 : }
2804 : #endif
2805 : #ifdef IBM
2806 : dval(d) += 0; /* normalize */
2807 : #endif
2808 0 : if (!dval(d)) {
2809 0 : *decpt = 1;
2810 0 : 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 0 : b = d2b(dval(d), &be, &bbits);
2828 : #ifdef Sudden_Underflow
2829 : i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
2830 : #else
2831 0 : if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) {
2832 : #endif
2833 0 : dval(d2) = dval(d);
2834 0 : word0(d2) &= Frac_mask1;
2835 0 : 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 0 : i -= Bias;
2864 : #ifdef IBM
2865 : i <<= 2;
2866 : i += j;
2867 : #endif
2868 : #ifndef Sudden_Underflow
2869 0 : denorm = 0;
2870 : }
2871 : else {
2872 : /* d is denormalized */
2873 :
2874 0 : i = bbits + be + (Bias + (P-1) - 1);
2875 0 : x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32
2876 0 : : word1(d) << 32 - i;
2877 0 : dval(d2) = x;
2878 0 : word0(d2) -= 31*Exp_msk1; /* adjust exponent */
2879 0 : i -= (Bias + (P-1) - 1) + 1;
2880 0 : denorm = 1;
2881 : }
2882 : #endif
2883 0 : ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
2884 0 : k = (int)ds;
2885 0 : if (ds < 0. && ds != k)
2886 0 : k--; /* want k = floor(ds) */
2887 0 : k_check = 1;
2888 0 : if (k >= 0 && k <= Ten_pmax) {
2889 0 : if (dval(d) < tens[k])
2890 0 : k--;
2891 0 : k_check = 0;
2892 : }
2893 0 : j = bbits - i - 1;
2894 0 : if (j >= 0) {
2895 0 : b2 = 0;
2896 0 : s2 = j;
2897 : }
2898 : else {
2899 0 : b2 = -j;
2900 0 : s2 = 0;
2901 : }
2902 0 : if (k >= 0) {
2903 0 : b5 = 0;
2904 0 : s5 = k;
2905 0 : s2 += k;
2906 : }
2907 : else {
2908 0 : b2 -= k;
2909 0 : b5 = -k;
2910 0 : s5 = 0;
2911 : }
2912 0 : if (mode < 0 || mode > 9)
2913 0 : mode = 0;
2914 :
2915 : #ifndef SET_INEXACT
2916 : #ifdef Check_FLT_ROUNDS
2917 : try_quick = Rounding == 1;
2918 : #else
2919 0 : try_quick = 1;
2920 : #endif
2921 : #endif /*SET_INEXACT*/
2922 :
2923 0 : if (mode > 5) {
2924 0 : mode -= 4;
2925 0 : try_quick = 0;
2926 : }
2927 0 : leftright = 1;
2928 0 : switch(mode) {
2929 : case 0:
2930 : case 1:
2931 0 : ilim = ilim1 = -1;
2932 0 : i = 18;
2933 0 : ndigits = 0;
2934 0 : break;
2935 : case 2:
2936 0 : leftright = 0;
2937 : /* no break */
2938 : case 4:
2939 0 : if (ndigits <= 0)
2940 0 : ndigits = 1;
2941 0 : ilim = ilim1 = i = ndigits;
2942 0 : break;
2943 : case 3:
2944 0 : leftright = 0;
2945 : /* no break */
2946 : case 5:
2947 0 : i = ndigits + k + 1;
2948 0 : ilim = i;
2949 0 : ilim1 = i - 1;
2950 0 : if (i <= 0)
2951 0 : i = 1;
2952 : }
2953 0 : s = s0 = rv_alloc(i);
2954 :
2955 : #ifdef Honor_FLT_ROUNDS
2956 : if (mode > 1 && rounding != 1)
2957 : leftright = 0;
2958 : #endif
2959 :
2960 0 : if (ilim >= 0 && ilim <= Quick_max && try_quick) {
2961 :
2962 : /* Try to get by with floating-point arithmetic. */
2963 :
2964 0 : i = 0;
2965 0 : dval(d2) = dval(d);
2966 0 : k0 = k;
2967 0 : ilim0 = ilim;
2968 0 : ieps = 2; /* conservative */
2969 0 : if (k > 0) {
2970 0 : ds = tens[k&0xf];
2971 0 : j = k >> 4;
2972 0 : if (j & Bletch) {
2973 : /* prevent overflows */
2974 0 : j &= Bletch - 1;
2975 0 : dval(d) /= bigtens[n_bigtens-1];
2976 0 : ieps++;
2977 : }
2978 0 : for(; j; j >>= 1, i++)
2979 0 : if (j & 1) {
2980 0 : ieps++;
2981 0 : ds *= bigtens[i];
2982 : }
2983 0 : dval(d) /= ds;
2984 : }
2985 0 : else if (j1 = -k) {
2986 0 : dval(d) *= tens[j1 & 0xf];
2987 0 : for(j = j1 >> 4; j; j >>= 1, i++)
2988 0 : if (j & 1) {
2989 0 : ieps++;
2990 0 : dval(d) *= bigtens[i];
2991 : }
2992 : }
2993 0 : if (k_check && dval(d) < 1. && ilim > 0) {
2994 0 : if (ilim1 <= 0)
2995 0 : goto fast_failed;
2996 0 : ilim = ilim1;
2997 0 : k--;
2998 0 : dval(d) *= 10.;
2999 0 : ieps++;
3000 : }
3001 0 : dval(eps) = ieps*dval(d) + 7.;
3002 0 : word0(eps) -= (P-1)*Exp_msk1;
3003 0 : if (ilim == 0) {
3004 0 : S = mhi = 0;
3005 0 : dval(d) -= 5.;
3006 0 : if (dval(d) > dval(eps))
3007 0 : goto one_digit;
3008 0 : if (dval(d) < -dval(eps))
3009 0 : goto no_digits;
3010 0 : goto fast_failed;
3011 : }
3012 : #ifndef No_leftright
3013 0 : if (leftright) {
3014 : /* Use Steele & White method of only
3015 : * generating digits needed.
3016 : */
3017 0 : dval(eps) = 0.5/tens[ilim-1] - dval(eps);
3018 0 : for(i = 0;;) {
3019 0 : L = dval(d);
3020 0 : dval(d) -= L;
3021 0 : *s++ = '0' + (int)L;
3022 0 : if (dval(d) < dval(eps))
3023 0 : goto ret1;
3024 0 : if (1. - dval(d) < dval(eps))
3025 0 : goto bump_up;
3026 0 : if (++i >= ilim)
3027 0 : break;
3028 0 : dval(eps) *= 10.;
3029 0 : dval(d) *= 10.;
3030 : }
3031 : }
3032 : else {
3033 : #endif
3034 : /* Generate ilim digits, then fix them up. */
3035 0 : dval(eps) *= tens[ilim-1];
3036 0 : for(i = 1;; i++, dval(d) *= 10.) {
3037 0 : L = (Long)(dval(d));
3038 0 : if (!(dval(d) -= L))
3039 0 : ilim = i;
3040 0 : *s++ = '0' + (int)L;
3041 0 : if (i == ilim) {
3042 0 : if (dval(d) > 0.5 + dval(eps))
3043 0 : goto bump_up;
3044 0 : else if (dval(d) < 0.5 - dval(eps)) {
3045 0 : while(*--s == '0');
3046 0 : s++;
3047 0 : goto ret1;
3048 : }
3049 0 : break;
3050 : }
3051 : }
3052 : #ifndef No_leftright
3053 : }
3054 : #endif
3055 : fast_failed:
3056 0 : s = s0;
3057 0 : dval(d) = dval(d2);
3058 0 : k = k0;
3059 0 : ilim = ilim0;
3060 : }
3061 :
3062 : /* Do we have a "small" integer? */
3063 :
3064 0 : if (be >= 0 && k <= Int_max) {
3065 : /* Yes. */
3066 0 : ds = tens[k];
3067 0 : if (ndigits < 0 && ilim <= 0) {
3068 0 : S = mhi = 0;
3069 0 : if (ilim < 0 || dval(d) <= 5*ds)
3070 : goto no_digits;
3071 0 : goto one_digit;
3072 : }
3073 0 : for(i = 1; i <= k+1; i++, dval(d) *= 10.) {
3074 0 : L = (Long)(dval(d) / ds);
3075 0 : 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 0 : *s++ = '0' + (int)L;
3084 0 : if (!dval(d)) {
3085 : #ifdef SET_INEXACT
3086 : inexact = 0;
3087 : #endif
3088 0 : break;
3089 : }
3090 0 : 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 0 : dval(d) += dval(d);
3099 0 : if (dval(d) > ds || dval(d) == ds && L & 1) {
3100 : bump_up:
3101 0 : while(*--s == '9')
3102 0 : if (s == s0) {
3103 0 : k++;
3104 0 : *s = '0';
3105 0 : break;
3106 : }
3107 0 : ++*s++;
3108 : }
3109 0 : break;
3110 : }
3111 : }
3112 0 : goto ret1;
3113 : }
3114 :
3115 0 : m2 = b2;
3116 0 : m5 = b5;
3117 0 : mhi = mlo = 0;
3118 0 : if (leftright) {
3119 0 : i =
3120 : #ifndef Sudden_Underflow
3121 0 : 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 0 : 1 + P - bbits;
3127 : #endif
3128 0 : b2 += i;
3129 0 : s2 += i;
3130 0 : mhi = i2b(1);
3131 : }
3132 0 : if (m2 > 0 && s2 > 0) {
3133 0 : i = m2 < s2 ? m2 : s2;
3134 0 : b2 -= i;
3135 0 : m2 -= i;
3136 0 : s2 -= i;
3137 : }
3138 0 : if (b5 > 0) {
3139 0 : if (leftright) {
3140 0 : if (m5 > 0) {
3141 0 : mhi = pow5mult(mhi, m5);
3142 0 : b1 = mult(mhi, b);
3143 0 : Bfree(b);
3144 0 : b = b1;
3145 : }
3146 0 : if (j = b5 - m5)
3147 0 : b = pow5mult(b, j);
3148 : }
3149 : else
3150 0 : b = pow5mult(b, b5);
3151 : }
3152 0 : S = i2b(1);
3153 0 : if (s5 > 0)
3154 0 : S = pow5mult(S, s5);
3155 :
3156 : /* Check for special case that d is a normalized power of 2. */
3157 :
3158 0 : spec_case = 0;
3159 0 : if ((mode < 2 || leftright)
3160 : #ifdef Honor_FLT_ROUNDS
3161 : && rounding == 1
3162 : #endif
3163 : ) {
3164 0 : if (!word1(d) && !(word0(d) & Bndry_mask)
3165 : #ifndef Sudden_Underflow
3166 0 : && word0(d) & (Exp_mask & ~Exp_msk1)
3167 : #endif
3168 : ) {
3169 : /* The special case */
3170 0 : b2 += Log2P;
3171 0 : s2 += Log2P;
3172 0 : 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 0 : if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
3185 0 : 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 0 : if (i > 4) {
3191 0 : i -= 4;
3192 0 : b2 += i;
3193 0 : m2 += i;
3194 0 : s2 += i;
3195 : }
3196 0 : else if (i < 4) {
3197 0 : i += 28;
3198 0 : b2 += i;
3199 0 : m2 += i;
3200 0 : s2 += i;
3201 : }
3202 0 : if (b2 > 0)
3203 0 : b = lshift(b, b2);
3204 0 : if (s2 > 0)
3205 0 : S = lshift(S, s2);
3206 0 : if (k_check) {
3207 0 : if (cmp(b,S) < 0) {
3208 0 : k--;
3209 0 : b = multadd(b, 10, 0); /* we botched the k estimate */
3210 0 : if (leftright)
3211 0 : mhi = multadd(mhi, 10, 0);
3212 0 : ilim = ilim1;
3213 : }
3214 : }
3215 0 : if (ilim <= 0 && (mode == 3 || mode == 5)) {
3216 0 : if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
3217 : /* no digits, fcvt style */
3218 : no_digits:
3219 0 : k = -1 - ndigits;
3220 0 : goto ret;
3221 : }
3222 : one_digit:
3223 0 : *s++ = '1';
3224 0 : k++;
3225 0 : goto ret;
3226 : }
3227 0 : if (leftright) {
3228 0 : if (m2 > 0)
3229 0 : mhi = lshift(mhi, m2);
3230 :
3231 : /* Compute mlo -- check for special case
3232 : * that d is a normalized power of 2.
3233 : */
3234 :
3235 0 : mlo = mhi;
3236 0 : if (spec_case) {
3237 0 : mhi = Balloc(mhi->k);
3238 0 : Bcopy(mhi, mlo);
3239 0 : mhi = lshift(mhi, Log2P);
3240 : }
3241 :
3242 0 : for(i = 1;;i++) {
3243 0 : dig = quorem(b,S) + '0';
3244 : /* Do we yet have the shortest decimal string
3245 : * that will round to d?
3246 : */
3247 0 : j = cmp(b, mlo);
3248 0 : delta = diff(S, mhi);
3249 0 : j1 = delta->sign ? 1 : cmp(b, delta);
3250 0 : Bfree(delta);
3251 : #ifndef ROUND_BIASED
3252 0 : if (j1 == 0 && mode != 1 && !(word1(d) & 1)
3253 : #ifdef Honor_FLT_ROUNDS
3254 : && rounding >= 1
3255 : #endif
3256 : ) {
3257 0 : if (dig == '9')
3258 0 : goto round_9_up;
3259 0 : if (j > 0)
3260 0 : dig++;
3261 : #ifdef SET_INEXACT
3262 : else if (!b->x[0] && b->wds <= 1)
3263 : inexact = 0;
3264 : #endif
3265 0 : *s++ = dig;
3266 0 : goto ret;
3267 : }
3268 : #endif
3269 0 : if (j < 0 || j == 0 && mode != 1
3270 : #ifndef ROUND_BIASED
3271 0 : && !(word1(d) & 1)
3272 : #endif
3273 : ) {
3274 0 : if (!b->x[0] && b->wds <= 1) {
3275 : #ifdef SET_INEXACT
3276 : inexact = 0;
3277 : #endif
3278 0 : 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 0 : if (j1 > 0) {
3288 0 : b = lshift(b, 1);
3289 0 : j1 = cmp(b, S);
3290 0 : if ((j1 > 0 || j1 == 0 && dig & 1)
3291 0 : && dig++ == '9')
3292 0 : goto round_9_up;
3293 : }
3294 : accept_dig:
3295 0 : *s++ = dig;
3296 0 : goto ret;
3297 : }
3298 0 : if (j1 > 0) {
3299 : #ifdef Honor_FLT_ROUNDS
3300 : if (!rounding)
3301 : goto accept_dig;
3302 : #endif
3303 0 : if (dig == '9') { /* possible if i == 1 */
3304 : round_9_up:
3305 0 : *s++ = '9';
3306 0 : goto roundoff;
3307 : }
3308 0 : *s++ = dig + 1;
3309 0 : goto ret;
3310 : }
3311 : #ifdef Honor_FLT_ROUNDS
3312 : keep_dig:
3313 : #endif
3314 0 : *s++ = dig;
3315 0 : if (i == ilim)
3316 0 : break;
3317 0 : b = multadd(b, 10, 0);
3318 0 : if (mlo == mhi)
3319 0 : mlo = mhi = multadd(mhi, 10, 0);
3320 : else {
3321 0 : mlo = multadd(mlo, 10, 0);
3322 0 : mhi = multadd(mhi, 10, 0);
3323 : }
3324 : }
3325 : }
3326 : else
3327 0 : for(i = 1;; i++) {
3328 0 : *s++ = dig = quorem(b,S) + '0';
3329 0 : if (!b->x[0] && b->wds <= 1) {
3330 : #ifdef SET_INEXACT
3331 : inexact = 0;
3332 : #endif
3333 0 : goto ret;
3334 : }
3335 0 : if (i >= ilim)
3336 0 : break;
3337 0 : 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 0 : b = lshift(b, 1);
3349 0 : j = cmp(b, S);
3350 0 : if (j > 0 || j == 0 && dig & 1) {
3351 : roundoff:
3352 0 : while(*--s == '9')
3353 0 : if (s == s0) {
3354 0 : k++;
3355 0 : *s++ = '1';
3356 0 : goto ret;
3357 : }
3358 0 : ++*s++;
3359 : }
3360 : else {
3361 : #ifdef Honor_FLT_ROUNDS
3362 : trimzeros:
3363 : #endif
3364 0 : while(*--s == '0');
3365 0 : s++;
3366 : }
3367 : ret:
3368 0 : Bfree(S);
3369 0 : if (mhi) {
3370 0 : if (mlo && mlo != mhi)
3371 0 : Bfree(mlo);
3372 0 : 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 0 : Bfree(b);
3387 0 : *s = 0;
3388 0 : *decpt = k + 1;
3389 0 : if (rve)
3390 0 : *rve = s;
3391 0 : 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 0 : PRStatus rv = PR_FAILURE;
3404 :
3405 0 : if (!_pr_initialized) _PR_ImplicitInitialization();
3406 :
3407 0 : if (mode < 0 || mode > 3) {
3408 0 : PR_SetError(PR_INVALID_ARGUMENT_ERROR, 0);
3409 0 : return rv;
3410 : }
3411 0 : result = dtoa(d, mode, ndigits, decpt, sign, rve);
3412 0 : if (!result) {
3413 0 : PR_SetError(PR_OUT_OF_MEMORY_ERROR, 0);
3414 0 : return rv;
3415 : }
3416 0 : resultlen = strlen(result)+1;
3417 0 : if (bufsize < resultlen) {
3418 0 : PR_SetError(PR_BUFFER_OVERFLOW_ERROR, 0);
3419 : } else {
3420 0 : memcpy(buf, result, resultlen);
3421 0 : if (rve) {
3422 0 : *rve = buf + (*rve - result);
3423 : }
3424 0 : rv = PR_SUCCESS;
3425 : }
3426 0 : freedtoa(result);
3427 0 : 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 0 : char *bufp = buf;
3447 : char *endnum;
3448 : U fval;
3449 :
3450 0 : dval(fval) = dfval;
3451 : /* If anything fails, we store an empty string in 'buf' */
3452 0 : num = (char*)PR_MALLOC(bufsz);
3453 0 : if (num == NULL) {
3454 0 : buf[0] = '\0';
3455 0 : return;
3456 : }
3457 : /* XXX Why use mode 1? */
3458 0 : if (PR_dtoa(dval(fval),1,prcsn,&decpt,&sign,&endnum,num,bufsz)
3459 : == PR_FAILURE) {
3460 0 : buf[0] = '\0';
3461 0 : goto done;
3462 : }
3463 0 : numdigits = endnum - num;
3464 0 : nump = num;
3465 :
3466 0 : if (sign &&
3467 0 : !(word0(fval) == Sign_bit && word1(fval) == 0) &&
3468 0 : !((word0(fval) & Exp_mask) == Exp_mask &&
3469 0 : (word1(fval) || (word0(fval) & 0xfffff)))) {
3470 0 : *bufp++ = '-';
3471 : }
3472 :
3473 0 : if (decpt == 9999) {
3474 0 : while ((*bufp++ = *nump++) != 0) {} /* nothing to execute */
3475 0 : goto done;
3476 : }
3477 :
3478 0 : if (decpt > (prcsn+1) || decpt < -(prcsn-1) || decpt < -5) {
3479 0 : *bufp++ = *nump++;
3480 0 : if (numdigits != 1) {
3481 0 : *bufp++ = '.';
3482 : }
3483 :
3484 0 : while (*nump != '\0') {
3485 0 : *bufp++ = *nump++;
3486 : }
3487 0 : *bufp++ = 'e';
3488 0 : PR_snprintf(bufp, bufsz - (bufp - buf), "%+d", decpt-1);
3489 0 : } else if (decpt >= 0) {
3490 0 : if (decpt == 0) {
3491 0 : *bufp++ = '0';
3492 : } else {
3493 0 : while (decpt--) {
3494 0 : if (*nump != '\0') {
3495 0 : *bufp++ = *nump++;
3496 : } else {
3497 0 : *bufp++ = '0';
3498 : }
3499 : }
3500 : }
3501 0 : if (*nump != '\0') {
3502 0 : *bufp++ = '.';
3503 0 : while (*nump != '\0') {
3504 0 : *bufp++ = *nump++;
3505 : }
3506 : }
3507 0 : *bufp++ = '\0';
3508 0 : } else if (decpt < 0) {
3509 0 : *bufp++ = '0';
3510 0 : *bufp++ = '.';
3511 0 : while (decpt++) {
3512 0 : *bufp++ = '0';
3513 : }
3514 :
3515 0 : while (*nump != '\0') {
3516 0 : *bufp++ = *nump++;
3517 : }
3518 0 : *bufp++ = '\0';
3519 : }
3520 : done:
3521 0 : PR_DELETE(num);
3522 : }
|
