File: | pr/Linux4.19_x86_64_gcc_glibc_PTH_64_DBG.OBJ/pr/src/misc/../../../../pr/src/misc/prdtoa.c |
Warning: | line 888, column 46 The left operand of '&' is a garbage value |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
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 | ||||
23 | static PRLock *dtoa_lock[2]; | |||
24 | ||||
25 | void _PR_InitDtoa(void) | |||
26 | { | |||
27 | dtoa_lock[0] = PR_NewLock(); | |||
28 | dtoa_lock[1] = PR_NewLock(); | |||
29 | } | |||
30 | ||||
31 | void _PR_CleanupDtoa(void) | |||
32 | { | |||
33 | PR_DestroyLock(dtoa_lock[0]); | |||
34 | dtoa_lock[0] = NULL((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 | |||
229 | typedef 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 | |||
246 | extern char *MALLOCmalloc(); | |||
247 | #else | |||
248 | extern 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)) | |||
259 | static 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 | |||
314 | extern "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 | |||
326 | Exactly one of IEEE_8087, IEEE_MC68k, IEEE_ARM, VAX, or IBM should be defined. | |||
327 | #endif | |||
328 | ||||
329 | typedef 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 | |||
478 | extern double rnd_prod(), rnd_quot(); | |||
479 | #else | |||
480 | extern 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 | ||||
526 | struct | |||
527 | Bigint { | |||
528 | struct Bigint *next; | |||
529 | int k, maxwds, sign, wds; | |||
530 | ULongPRUint32 x[1]; | |||
531 | }; | |||
532 | ||||
533 | typedef struct Bigint Bigint; | |||
534 | ||||
535 | static Bigint *freelist[Kmax7+1]; | |||
536 | ||||
537 | static Bigint * | |||
538 | Balloc | |||
539 | #ifdef KR_headers | |||
540 | (k) int k; | |||
541 | #else | |||
542 | (int k) | |||
543 | #endif | |||
544 | { | |||
545 | int x; | |||
546 | Bigint *rv; | |||
547 | #ifndef Omit_Private_Memory | |||
548 | unsigned int len; | |||
549 | #endif | |||
550 | ||||
551 | ACQUIRE_DTOA_LOCK(0)PR_Lock(dtoa_lock[0]); | |||
552 | /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */ | |||
553 | /* but this case seems very unlikely. */ | |||
554 | if (k <= Kmax7 && (rv = freelist[k])) { | |||
555 | freelist[k] = rv->next; | |||
556 | } | |||
557 | else { | |||
558 | x = 1 << k; | |||
559 | #ifdef Omit_Private_Memory | |||
560 | rv = (Bigint *)MALLOCmalloc(sizeof(Bigint) + (x-1)*sizeof(ULongPRUint32)); | |||
561 | #else | |||
562 | len = (sizeof(Bigint) + (x-1)*sizeof(ULongPRUint32) + sizeof(double) - 1) | |||
563 | /sizeof(double); | |||
564 | if (k <= Kmax7 && pmem_next - private_mem + len <= PRIVATE_mem((2304 +sizeof(double)-1)/sizeof(double))) { | |||
565 | rv = (Bigint*)pmem_next; | |||
566 | pmem_next += len; | |||
567 | } | |||
568 | else { | |||
569 | rv = (Bigint*)MALLOCmalloc(len*sizeof(double)); | |||
570 | } | |||
571 | #endif | |||
572 | rv->k = k; | |||
573 | rv->maxwds = x; | |||
574 | } | |||
575 | FREE_DTOA_LOCK(0)PR_Unlock(dtoa_lock[0]); | |||
576 | rv->sign = rv->wds = 0; | |||
577 | return rv; | |||
578 | } | |||
579 | ||||
580 | static void | |||
581 | Bfree | |||
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, \ | |||
605 | y->wds*sizeof(LongPRInt32) + 2*sizeof(int)) | |||
606 | ||||
607 | static Bigint * | |||
608 | multadd | |||
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 | ||||
664 | static Bigint * | |||
665 | s2b | |||
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 | ||||
706 | static int | |||
707 | hi0bits | |||
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 | ||||
745 | static int | |||
746 | lo0bits | |||
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 | ||||
807 | static Bigint * | |||
808 | i2b | |||
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 | ||||
823 | static Bigint * | |||
824 | mult | |||
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
| |||
845 | c = a; | |||
846 | a = b; | |||
847 | b = c; | |||
848 | } | |||
849 | k = a->k; | |||
850 | wa = a->wds; | |||
851 | wb = b->wds; | |||
852 | wc = wa + wb; | |||
853 | if (wc > a->maxwds) { | |||
854 | k++; | |||
855 | } | |||
856 | c = Balloc(k); | |||
857 | for(x = c->x, xa = x + wc; x < xa; x++) { | |||
858 | *x = 0; | |||
859 | } | |||
860 | xa = a->x; | |||
861 | xae = xa + wa; | |||
862 | xb = b->x; | |||
863 | xbe = xb + wb; | |||
864 | xc0 = c->x; | |||
865 | #ifdef ULLong | |||
866 | for(; xb < xbe; xc0++) { | |||
867 | if (y = *xb++) { | |||
868 | x = xa; | |||
869 | xc = xc0; | |||
870 | carry = 0; | |||
871 | do { | |||
872 | z = *x++ * (ULLong)y + *xc + carry; | |||
873 | carry = z >> 32; | |||
874 | *xc++ = z & FFFFFFFF0xffffffffUL; | |||
875 | } | |||
876 | while(x < xae); | |||
877 | *xc = carry; | |||
878 | } | |||
879 | } | |||
880 | #else | |||
881 | #ifdef Pack_32 | |||
882 | for(; xb < xbe; xb++, xc0++) { | |||
883 | if (y = *xb & 0xffff) { | |||
884 | x = xa; | |||
885 | xc = xc0; | |||
886 | carry = 0; | |||
887 | do { | |||
888 | z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; | |||
| ||||
889 | carry = z >> 16; | |||
890 | z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; | |||
891 | carry = z2 >> 16; | |||
892 | Storeinc(xc, z2, z)(((unsigned short *)xc)[1] = (unsigned short)z2, ((unsigned short *)xc)[0] = (unsigned short)z, xc++); | |||
893 | } | |||
894 | while(x < xae); | |||
895 | *xc = carry; | |||
896 | } | |||
897 | if (y = *xb >> 16) { | |||
898 | x = xa; | |||
899 | xc = xc0; | |||
900 | carry = 0; | |||
901 | z2 = *xc; | |||
902 | do { | |||
903 | z = (*x & 0xffff) * y + (*xc >> 16) + carry; | |||
904 | carry = z >> 16; | |||
905 | Storeinc(xc, z, z2)(((unsigned short *)xc)[1] = (unsigned short)z, ((unsigned short *)xc)[0] = (unsigned short)z2, xc++); | |||
906 | z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; | |||
907 | carry = z2 >> 16; | |||
908 | } | |||
909 | while(x < xae); | |||
910 | *xc = z2; | |||
911 | } | |||
912 | } | |||
913 | #else | |||
914 | for(; xb < xbe; xc0++) { | |||
915 | if (y = *xb++) { | |||
916 | x = xa; | |||
917 | xc = xc0; | |||
918 | carry = 0; | |||
919 | do { | |||
920 | z = *x++ * y + *xc + carry; | |||
921 | carry = z >> 16; | |||
922 | *xc++ = z & 0xffff; | |||
923 | } | |||
924 | while(x < xae); | |||
925 | *xc = carry; | |||
926 | } | |||
927 | } | |||
928 | #endif | |||
929 | #endif | |||
930 | for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; | |||
931 | c->wds = wc; | |||
932 | return c; | |||
933 | } | |||
934 | ||||
935 | static Bigint *p5s; | |||
936 | ||||
937 | static Bigint * | |||
938 | pow5mult | |||
939 | #ifdef KR_headers | |||
940 | (b, k) Bigint *b; int k; | |||
941 | #else | |||
942 | (Bigint *b, int k) | |||
943 | #endif | |||
944 | { | |||
945 | Bigint *b1, *p5, *p51; | |||
946 | int i; | |||
947 | static int p05[3] = { 5, 25, 125 }; | |||
948 | ||||
949 | if (i = k & 3) { | |||
950 | b = multadd(b, p05[i-1], 0); | |||
951 | } | |||
952 | ||||
953 | if (!(k >>= 2)) { | |||
954 | return b; | |||
955 | } | |||
956 | if (!(p5 = p5s)) { | |||
957 | /* first time */ | |||
958 | #ifdef MULTIPLE_THREADS | |||
959 | ACQUIRE_DTOA_LOCK(1)PR_Lock(dtoa_lock[1]); | |||
960 | if (!(p5 = p5s)) { | |||
961 | p5 = p5s = i2b(625); | |||
962 | p5->next = 0; | |||
963 | } | |||
964 | FREE_DTOA_LOCK(1)PR_Unlock(dtoa_lock[1]); | |||
965 | #else | |||
966 | p5 = p5s = i2b(625); | |||
967 | p5->next = 0; | |||
968 | #endif | |||
969 | } | |||
970 | for(;;) { | |||
971 | if (k & 1) { | |||
972 | b1 = mult(b, p5); | |||
973 | Bfree(b); | |||
974 | b = b1; | |||
975 | } | |||
976 | if (!(k >>= 1)) { | |||
977 | break; | |||
978 | } | |||
979 | if (!(p51 = p5->next)) { | |||
980 | #ifdef MULTIPLE_THREADS | |||
981 | ACQUIRE_DTOA_LOCK(1)PR_Lock(dtoa_lock[1]); | |||
982 | if (!(p51 = p5->next)) { | |||
983 | p51 = p5->next = mult(p5,p5); | |||
984 | p51->next = 0; | |||
985 | } | |||
986 | FREE_DTOA_LOCK(1)PR_Unlock(dtoa_lock[1]); | |||
987 | #else | |||
988 | p51 = p5->next = mult(p5,p5); | |||
989 | p51->next = 0; | |||
990 | #endif | |||
991 | } | |||
992 | p5 = p51; | |||
993 | } | |||
994 | return b; | |||
995 | } | |||
996 | ||||
997 | static Bigint * | |||
998 | lshift | |||
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 | ||||
1062 | static int | |||
1063 | cmp | |||
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 | ||||
1101 | static Bigint * | |||
1102 | diff | |||
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 | ||||
1197 | static double | |||
1198 | ulp | |||
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 | ||||
1240 | static double | |||
1241 | b2d | |||
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 | |||
1301 | ret_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 | ||||
1312 | static Bigint * | |||
1313 | d2b | |||
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 | ||||
1455 | static double | |||
1456 | ratio | |||
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 | ||||
1499 | static CONSTconst double | |||
1500 | tens[] = { | |||
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 | ||||
1509 | static CONSTconst double | |||
1510 | #ifdef IEEE_Arith | |||
1511 | bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; | |||
1512 | static 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 | |||
1526 | bigtens[] = { 1e16, 1e32, 1e64 }; | |||
1527 | static CONSTconst double tinytens[] = { 1e-16, 1e-32, 1e-64 }; | |||
1528 | #define n_bigtens5 3 | |||
1529 | #else | |||
1530 | bigtens[] = { 1e16, 1e32 }; | |||
1531 | static 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 | ||||
1550 | static int | |||
1551 | match | |||
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 | |||
1574 | static void | |||
1575 | hexnan | |||
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 | ||||
1633 | PR_IMPLEMENT(double)__attribute__((visibility("default"))) double | |||
1634 | PR_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 | } | |||
1689 | break2: | |||
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) { | |||
1741 | have_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 | } | |||
1762 | dig_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 */ | |||
1842 | ret0: | |||
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 | */ | |||
1913 | vax_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) { | |||
1973 | ovfl: | |||
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) { | |||
2077 | undfl: | |||
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 | ||||
2102 | bd0 = s2b(s0, nd0, nd, y); | |||
2103 | ||||
2104 | for(;;) { | |||
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 | } | |||
2227 | apply_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]) { | |||
2351 | drop_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 | |||
2565 | cont: | |||
2566 | Bfree(bb); | |||
2567 | Bfree(bd); | |||
2568 | Bfree(bs); | |||
2569 | Bfree(delta); | |||
2570 | } | |||
2571 | #ifdef SET_INEXACT | |||
2572 | if (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 | } | |||
2579 | else if (!oldinexact) { | |||
2580 | clear_inexact(); | |||
2581 | } | |||
2582 | #endif | |||
2583 | #ifdef Avoid_Underflow | |||
2584 | if (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 | |||
2597 | if (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 | |||
2603 | retfree: | |||
2604 | Bfree(bb); | |||
2605 | Bfree(bd); | |||
2606 | Bfree(bs); | |||
2607 | Bfree(bd0); | |||
2608 | Bfree(delta); | |||
2609 | ret: | |||
2610 | if (se) { | |||
2611 | *se = (char *)s; | |||
2612 | } | |||
2613 | return sign ? -dval(rv)(rv).d : dval(rv)(rv).d; | |||
2614 | } | |||
2615 | ||||
2616 | static int | |||
2617 | quorem | |||
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 | |||
2742 | static char *dtoa_result; | |||
2743 | #endif | |||
2744 | ||||
2745 | static char * | |||
2746 | #ifdef KR_headers | |||
2747 | rv_alloc(i) int i; | |||
2748 | #else | |||
2749 | rv_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 | ||||
2769 | static char * | |||
2770 | #ifdef KR_headers | |||
2771 | nrv_alloc(s, rve, n) char *s, **rve; int n; | |||
2772 | #else | |||
2773 | nrv_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 | ||||
2794 | static void | |||
2795 | #ifdef KR_headers | |||
2796 | freedtoa(s) char *s; | |||
2797 | #else | |||
2798 | freedtoa(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 | ||||
2845 | static char * | |||
2846 | dtoa | |||
2847 | #ifdef KR_headers | |||
2848 | (dd, mode, ndigits, decpt, sign, rve) | |||
2849 | double 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) { | |||
| ||||
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) | |||
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) { | |||
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))) { | |||
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 | } | |||
3009 | else { | |||
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 | |||
3014 | : word1(d)(d).L[0] << 32 - i; | |||
3015 | dval(d2)(d2).d = x; | |||
3016 | word0(d2)(d2).L[1] -= 31*Exp_msk10x100000; /* adjust exponent */ | |||
3017 | i -= (Bias1023 + (P53-1) - 1) + 1; | |||
3018 | denorm = 1; | |||
3019 | } | |||
3020 | #endif | |||
3021 | ds = (dval(d2)(d2).d-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981; | |||
3022 | k = (int)ds; | |||
3023 | if (ds < 0. && ds != k) { | |||
3024 | k--; /* want k = floor(ds) */ | |||
3025 | } | |||
3026 | k_check = 1; | |||
3027 | if (k >= 0 && k <= Ten_pmax22) { | |||
3028 | if (dval(d)(d).d < tens[k]) { | |||
3029 | k--; | |||
3030 | } | |||
3031 | k_check = 0; | |||
3032 | } | |||
3033 | j = bbits - i - 1; | |||
3034 | if (j >= 0) { | |||
3035 | b2 = 0; | |||
3036 | s2 = j; | |||
3037 | } | |||
3038 | else { | |||
3039 | b2 = -j; | |||
3040 | s2 = 0; | |||
3041 | } | |||
3042 | if (k
| |||
3043 | b5 = 0; | |||
3044 | s5 = k; | |||
3045 | s2 += k; | |||
3046 | } | |||
3047 | else { | |||
3048 | b2 -= k; | |||
3049 | b5 = -k; | |||
3050 | s5 = 0; | |||
3051 | } | |||
3052 | if (mode < 0 || mode > 9) { | |||
3053 | mode = 0; | |||
3054 | } | |||
3055 | ||||
3056 | #ifndef SET_INEXACT | |||
3057 | #ifdef Check_FLT_ROUNDS | |||
3058 | try_quick = Rounding(__builtin_flt_rounds()) == 1; | |||
3059 | #else | |||
3060 | try_quick = 1; | |||
3061 | #endif | |||
3062 | #endif /*SET_INEXACT*/ | |||
3063 | ||||
3064 | if (mode > 5) { | |||
3065 | mode -= 4; | |||
3066 | try_quick = 0; | |||
3067 | } | |||
3068 | leftright = 1; | |||
3069 | switch(mode) { | |||
3070 | case 0: | |||
3071 | case 1: | |||
3072 | ilim = ilim1 = -1; | |||
3073 | i = 18; | |||
3074 | ndigits = 0; | |||
3075 | break; | |||
3076 | case 2: | |||
3077 | leftright = 0; | |||
3078 | /* no break */ | |||
3079 | case 4: | |||
3080 | if (ndigits <= 0) { | |||
3081 | ndigits = 1; | |||
3082 | } | |||
3083 | ilim = ilim1 = i = ndigits; | |||
3084 | break; | |||
3085 | case 3: | |||
3086 | leftright = 0; | |||
3087 | /* no break */ | |||
3088 | case 5: | |||
3089 | i = ndigits + k + 1; | |||
3090 | ilim = i; | |||
3091 | ilim1 = i - 1; | |||
3092 | if (i <= 0) { | |||
3093 | i = 1; | |||
3094 | } | |||
3095 | } | |||
3096 | s = s0 = rv_alloc(i); | |||
3097 | ||||
3098 | #ifdef Honor_FLT_ROUNDS | |||
3099 | if (mode > 1 && rounding != 1) { | |||
3100 | leftright = 0; | |||
3101 | } | |||
3102 | #endif | |||
3103 | ||||
3104 | if (ilim
| |||
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 | |||
3207 | fast_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 | ||||
3216 | if (be
| |||
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) { | |||
3253 | bump_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 | ||||
3268 | m2 = b2; | |||
3269 | m5 = b5; | |||
3270 | mhi = mlo = 0; | |||
3271 | if (leftright
| |||
3272 | i = | |||
3273 | #ifndef Sudden_Underflow | |||
3274 | denorm
| |||
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 | } | |||
3285 | if (m2 > 0 && s2 > 0) { | |||
3286 | i = m2 < s2 ? m2 : s2; | |||
3287 | b2 -= i; | |||
3288 | m2 -= i; | |||
3289 | s2 -= i; | |||
3290 | } | |||
3291 | if (b5 > 0) { | |||
3292 | if (leftright
| |||
3293 | if (m5
| |||
3294 | mhi = pow5mult(mhi, m5); | |||
3295 | b1 = mult(mhi, b); | |||
3296 | Bfree(b); | |||
3297 | b = b1; | |||
3298 | } | |||
3299 | if (j = b5 - m5) { | |||
3300 | b = pow5mult(b, j); | |||
3301 | } | |||
3302 | } | |||
3303 | else { | |||
3304 | b = pow5mult(b, b5); | |||
3305 | } | |||
3306 | } | |||
3307 | S = i2b(1); | |||
3308 | if (s5 > 0) { | |||
3309 | S = pow5mult(S, s5); | |||
3310 | } | |||
3311 | ||||
3312 | /* Check for special case that d is a normalized power of 2. */ | |||
3313 | ||||
3314 | spec_case = 0; | |||
3315 | if ((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 | |||
3340 | if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) { | |||
3341 | i = 32 - i; | |||
3342 | } | |||
3343 | #else | |||
3344 | if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) { | |||
3345 | i = 16 - i; | |||
3346 | } | |||
3347 | #endif | |||
3348 | if (i > 4) { | |||
3349 | i -= 4; | |||
3350 | b2 += i; | |||
3351 | m2 += i; | |||
3352 | s2 += i; | |||
3353 | } | |||
3354 | else if (i < 4) { | |||
3355 | i += 28; | |||
3356 | b2 += i; | |||
3357 | m2 += i; | |||
3358 | s2 += i; | |||
3359 | } | |||
3360 | if (b2 > 0) { | |||
3361 | b = lshift(b, b2); | |||
3362 | } | |||
3363 | if (s2 > 0) { | |||
3364 | S = lshift(S, s2); | |||
3365 | } | |||
3366 | if (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 | } | |||
3376 | if (ilim <= 0 && (mode == 3 || mode == 5)) { | |||
3377 | if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { | |||
3378 | /* no digits, fcvt style */ | |||
3379 | no_digits: | |||
3380 | k = -1 - ndigits; | |||
3381 | goto ret; | |||
3382 | } | |||
3383 | one_digit: | |||
3384 | *s++ = '1'; | |||
3385 | k++; | |||
3386 | goto ret; | |||
3387 | } | |||
3388 | if (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 | } | |||
3460 | accept_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 */ | |||
3471 | round_9_up: | |||
3472 | *s++ = '9'; | |||
3473 | goto roundoff; | |||
3474 | } | |||
3475 | *s++ = dig + 1; | |||
3476 | goto ret; | |||
3477 | } | |||
3478 | #ifdef Honor_FLT_ROUNDS | |||
3479 | keep_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 | } | |||
3495 | else | |||
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 | |||
3513 | switch(rounding) { | |||
3514 | case 0: goto trimzeros; | |||
3515 | case 2: goto roundoff; | |||
3516 | } | |||
3517 | #endif | |||
3518 | b = lshift(b, 1); | |||
3519 | j = cmp(b, S); | |||
3520 | if (j > 0 || j == 0 && dig & 1) { | |||
3521 | roundoff: | |||
3522 | while(*--s == '9') | |||
3523 | if (s == s0) { | |||
3524 | k++; | |||
3525 | *s++ = '1'; | |||
3526 | goto ret; | |||
3527 | } | |||
3528 | ++*s++; | |||
3529 | } | |||
3530 | else { | |||
3531 | #ifdef Honor_FLT_ROUNDS | |||
3532 | trimzeros: | |||
3533 | #endif | |||
3534 | while(*--s == '0'); | |||
3535 | s++; | |||
3536 | } | |||
3537 | ret: | |||
3538 | Bfree(S); | |||
3539 | if (mhi) { | |||
3540 | if (mlo && mlo != mhi) { | |||
3541 | Bfree(mlo); | |||
3542 | } | |||
3543 | Bfree(mhi); | |||
3544 | } | |||
3545 | ret1: | |||
3546 | #ifdef SET_INEXACT | |||
3547 | if (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 | } | |||
3554 | else if (!oldinexact) { | |||
3555 | clear_inexact(); | |||
3556 | } | |||
3557 | #endif | |||
3558 | Bfree(b); | |||
3559 | *s = 0; | |||
3560 | *decpt = k + 1; | |||
3561 | if (rve) { | |||
3562 | *rve = s; | |||
3563 | } | |||
3564 | return s0; | |||
3565 | } | |||
3566 | #ifdef __cplusplus | |||
3567 | } | |||
3568 | #endif | |||
3569 | ||||
3570 | PR_IMPLEMENT(PRStatus)__attribute__((visibility("default"))) PRStatus | |||
3571 | PR_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 | */ | |||
3616 | PR_IMPLEMENT(void)__attribute__((visibility("default"))) void | |||
3617 | PR_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 | } | |||
3695 | done: | |||
3696 | PR_DELETE(num){ PR_Free(num); (num) = ((void*)0); }; | |||
3697 | } |