Bug Summary

File:pr/Linux4.19_x86_64_gcc_glibc_PTH_64_DBG.OBJ/pr/src/misc/../../../../pr/src/misc/prdtoa.c
Warning:line 3013, column 49
The result of right shift is undefined because the right operand '32' is not smaller than 32, the capacity of 'PRUint32'

Annotated Source Code

Press '?' to see keyboard shortcuts

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