Bug Summary

File:pr/Linux6.7_x86_64_gcc_glibc_PTH_64_DBG.OBJ/pr/src/misc/../../../../pr/src/misc/prdtoa.c
Warning:line 901, column 16
Assigned value is garbage or undefined

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/Linux6.7_x86_64_gcc_glibc_PTH_64_DBG.OBJ/pr/src/misc -fcoverage-compilation-dir=/var/lib/jenkins/workspace/nss-scan-build/nspr/Linux6.7_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/Linux6.7_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-07-17-052436-2337001-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])) {
54
←
Assuming 'k' is > Kmax
→
555 freelist[k] = rv->next;
556 }
557 else {
558 x = 1 << k;
55
←
Assuming right operand of bit shift is less than 32
→
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
55.1
'k' is > Kmax
 <= 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));
56
←
Storing uninitialized value
→
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
49.1
'a->wds' is >= 'b->wds'
 < b->wds) {
50
←
Taking false branch
→
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) {
51
←
Assuming 'wc' is <= field 'maxwds'
→
52
←
Taking false branch
→
854 k++;
855 }
856 c = Balloc(k);
53
←
Calling 'Balloc'
→
57
←
Returning from 'Balloc'
→
857 for(x = c->x, xa = x + wc; x < xa; x++) {
58
←
Assuming 'x' is >= 'xa'
→
59
←
Loop condition is false. Execution continues on line 860
→
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++) {
60
←
Assuming 'xb' is < 'xbe'
→
61
←
Loop condition is true. Entering loop body
→
883 if (y = *xb & 0xffff) {
62
←
Assuming 'y' is 0
→
63
←
Taking false branch
→
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) {
64
←
Assuming 'y' is not equal to 0
→
65
←
Taking true branch
→
898 x = xa;
899 xc = xc0;
900 carry = 0;
901 z2 = *xc;
66
←
Assigned value is garbage or undefined
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) {
34
←
Assuming 'i' is 0
→
35
←
Taking false branch
→
950 b = multadd(b, p05[i-1], 0);
951 }
952
953 if (!(k >>= 2)) {
36
←
Assuming the condition is false
→
37
←
Taking false branch
→
954 return b;
955 }
956 if (!(p5 = p5s)) {
38
←
Assuming 'p5' is non-null
→
39
←
Taking false branch
→
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(;;) {
40
←
Loop condition is true. Entering loop body
→
971 if (k & 1) {
41
←
Assuming the condition is false
→
42
←
Taking false branch
→
972 b1 = mult(b, p5);
973 Bfree(b);
974 b = b1;
975 }
976 if (!(k >>= 1)) {
43
←
Assuming the condition is false
→
44
←
Taking false branch
→
977 break;
978 }
979 if (!(p51 = p5->next)) {
45
←
Assuming 'p51' is null
→
46
←
Taking true branch
→
980#ifdef MULTIPLE_THREADS
981 ACQUIRE_DTOA_LOCK(1)PR_Lock(dtoa_lock[1]);
982 if (!(p51 = p5->next)) {
47
←
Assuming 'p51' is null
→
48
←
Taking true branch
→
983 p51 = p5->next = mult(p5,p5);
49
←
Calling 'mult'
→
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 > 32 ? word0(d)(d).L[1] << 64 - i | word1(d)(d).L[0] >> i - 32
9
←
Assuming 'i' is <= 32
→
10
←
'?' condition is false
→
3014 : word1(d)(d).L[0] << 32 - i;
11
←
Assuming right operand of bit shift is non-negative but less than 32
→
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) {
12
←
Assuming the condition is false
→
3024 k--; /* want k = floor(ds) */
3025}
3026k_check = 1;
3027if (k >= 0 && k <= Ten_pmax22) {
13
←
Assuming 'k' is < 0
→
3028 if (dval(d)(d).d < tens[k]) {
3029 k--;
3030 }
3031 k_check = 0;
3032}
3033j = bbits - i - 1;
3034if (j >= 0) {
14
←
Assuming 'j' is < 0
→
15
←
Taking false branch
→
3035 b2 = 0;
3036 s2 = j;
3037}
3038else {
3039 b2 = -j;
3040 s2 = 0;
3041}
3042if (k
15.1
'k' is < 0
 >= 0) {
16
←
Taking false branch
→
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) {
17
←
Assuming 'mode' is >= 0
→
18
←
Assuming 'mode' is <= 9
→
19
←
Taking false branch
→
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) {
20
←
Assuming 'mode' is <= 5
→
21
←
Taking false branch
→
3065 mode -= 4;
3066 try_quick = 0;
3067}
3068leftright = 1;
3069switch(mode) {
22
←
Control jumps to 'case 5:' at line 3088
→
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) {
23
←
Assuming 'i' is > 0
→
24
←
Taking false branch
→
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
24.1
'ilim' is >= 0
 >= 0 && ilim <= Quick_max14 && try_quick) {
25
←
Assuming 'ilim' is > Quick_max
→
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
25.1
'be' is < 0
 >= 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
25.2
'leftright' is 1
) {
3272 i =
26
←
Taking true branch
→
3273#ifndef Sudden_Underflow
3274 denorm
26.1
'denorm' is 1
 ? be + (Bias1023 + (P53-1) - 1 + 1) :
27
←
'?' condition is true
→
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) {
28
←
Assuming 'm2' is <= 0
→
3286 i = m2 < s2 ? m2 : s2;
3287 b2 -= i;
3288 m2 -= i;
3289 s2 -= i;
3290}
3291if (b5 > 0) {
29
←
Assuming 'b5' is > 0
→
30
←
Taking true branch
→
3292 if (leftright
30.1
'leftright' is 1
) {
31
←
Taking true branch
→
3293 if (m5
31.1
'm5' is > 0
 > 0) {
32
←
Taking true branch
→
3294 mhi = pow5mult(mhi, m5);
33
←
Calling 'pow5mult'
→
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}