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