summaryrefslogtreecommitdiff
path: root/src/lib/evil
diff options
context:
space:
mode:
authorVincent Torri <vincent dot torri at gmail dot com>2015-12-03 09:15:12 +0100
committerTom Hacohen <tom@stosb.com>2015-12-04 13:17:55 +0000
commitc66608dce06340b0dc779ea5101c8ea3774a2a7c (patch)
tree6721763c8ff1f4f55da03ea67fce8a6355048c27 /src/lib/evil
parent5411e23a0f0f17172cef3ad32658799a9174d4a3 (diff)
Evil: remove useless code
This code was used for GNU printf and al. but it is useless now
Diffstat (limited to 'src/lib/evil')
-rw-r--r--src/lib/evil/gdtoa/README357
-rw-r--r--src/lib/evil/gdtoa/README.mingw20
-rw-r--r--src/lib/evil/gdtoa/arithchk.c192
-rw-r--r--src/lib/evil/gdtoa/dmisc.c196
-rw-r--r--src/lib/evil/gdtoa/dtoa.c750
-rw-r--r--src/lib/evil/gdtoa/g__fmt.c142
-rw-r--r--src/lib/evil/gdtoa/g_dfmt.c90
-rw-r--r--src/lib/evil/gdtoa/g_ffmt.c88
-rw-r--r--src/lib/evil/gdtoa/g_xfmt.c143
-rw-r--r--src/lib/evil/gdtoa/gd_arith.h6
-rw-r--r--src/lib/evil/gdtoa/gd_qnan.h12
-rw-r--r--src/lib/evil/gdtoa/gdtoa.c733
-rw-r--r--src/lib/evil/gdtoa/gdtoa.h121
-rw-r--r--src/lib/evil/gdtoa/gdtoa_fltrnds.h18
-rw-r--r--src/lib/evil/gdtoa/gdtoaimp.h645
-rw-r--r--src/lib/evil/gdtoa/gethex.c340
-rw-r--r--src/lib/evil/gdtoa/gmisc.c76
-rw-r--r--src/lib/evil/gdtoa/hd_init.c49
-rw-r--r--src/lib/evil/gdtoa/hexnan.c139
-rw-r--r--src/lib/evil/gdtoa/misc.c860
-rw-r--r--src/lib/evil/gdtoa/qnan.c116
-rw-r--r--src/lib/evil/gdtoa/smisc.c149
-rw-r--r--src/lib/evil/gdtoa/strtodg.c979
-rw-r--r--src/lib/evil/gdtoa/strtof.c77
-rw-r--r--src/lib/evil/gdtoa/strtopx.c119
-rw-r--r--src/lib/evil/gdtoa/sum.c91
-rw-r--r--src/lib/evil/gdtoa/ulp.c61
27 files changed, 0 insertions, 6569 deletions
diff --git a/src/lib/evil/gdtoa/README b/src/lib/evil/gdtoa/README
deleted file mode 100644
index ce8be55888..0000000000
--- a/src/lib/evil/gdtoa/README
+++ /dev/null
@@ -1,357 +0,0 @@
1This directory contains source for a library of binary -> decimal
2and decimal -> binary conversion routines, for single-, double-,
3and extended-precision IEEE binary floating-point arithmetic, and
4other IEEE-like binary floating-point, including "double double",
5as in
6
7 T. J. Dekker, "A Floating-Point Technique for Extending the
8 Available Precision", Numer. Math. 18 (1971), pp. 224-242
9
10and
11
12 "Inside Macintosh: PowerPC Numerics", Addison-Wesley, 1994
13
14The conversion routines use double-precision floating-point arithmetic
15and, where necessary, high precision integer arithmetic. The routines
16are generalizations of the strtod and dtoa routines described in
17
18 David M. Gay, "Correctly Rounded Binary-Decimal and
19 Decimal-Binary Conversions", Numerical Analysis Manuscript
20 No. 90-10, Bell Labs, Murray Hill, 1990;
21 http://cm.bell-labs.com/cm/cs/what/ampl/REFS/rounding.ps.gz
22
23(based in part on papers by Clinger and Steele & White: see the
24references in the above paper).
25
26The present conversion routines should be able to use any of IEEE binary,
27VAX, or IBM-mainframe double-precision arithmetic internally, but I (dmg)
28have so far only had a chance to test them with IEEE double precision
29arithmetic.
30
31The core conversion routines are strtodg for decimal -> binary conversions
32and gdtoa for binary -> decimal conversions. These routines operate
33on arrays of unsigned 32-bit integers of type ULong, a signed 32-bit
34exponent of type Long, and arithmetic characteristics described in
35struct FPI; FPI, Long, and ULong are defined in gdtoa.h. File arith.h
36is supposed to provide #defines that cause gdtoa.h to define its
37types correctly. File arithchk.c is source for a program that
38generates a suitable arith.h on all systems where I've been able to
39test it.
40
41The core conversion routines are meant to be called by helper routines
42that know details of the particular binary arithmetic of interest and
43convert. The present directory provides helper routines for 5 variants
44of IEEE binary floating-point arithmetic, each indicated by one or
45two letters:
46
47 f IEEE single precision
48 d IEEE double precision
49 x IEEE extended precision, as on Intel 80x87
50 and software emulations of Motorola 68xxx chips
51 that do not pad the way the 68xxx does, but
52 only store 80 bits
53 xL IEEE extended precision, as on Motorola 68xxx chips
54 Q quad precision, as on Sun Sparc chips
55 dd double double, pairs of IEEE double numbers
56 whose sum is the desired value
57
58For decimal -> binary conversions, there are three families of
59helper routines: one for round-nearest (or the current rounding
60mode on IEEE-arithmetic systems that provide the C99 fegetround()
61function, if compiled with -DHonor_FLT_ROUNDS):
62
63 strtof
64 strtod
65 strtodd
66 strtopd
67 strtopf
68 strtopx
69 strtopxL
70 strtopQ
71
72one with rounding direction specified:
73
74 strtorf
75 strtord
76 strtordd
77 strtorx
78 strtorxL
79 strtorQ
80
81and one for computing an interval (at most one bit wide) that contains
82the decimal number:
83
84 strtoIf
85 strtoId
86 strtoIdd
87 strtoIx
88 strtoIxL
89 strtoIQ
90
91The latter call strtoIg, which makes one call on strtodg and adjusts
92the result to provide the desired interval. On systems where native
93arithmetic can easily make one-ulp adjustments on values in the
94desired floating-point format, it might be more efficient to use the
95native arithmetic. Routine strtodI is a variant of strtoId that
96illustrates one way to do this for IEEE binary double-precision
97arithmetic -- but whether this is more efficient remains to be seen.
98
99Functions strtod and strtof have "natural" return types, float and
100double -- strtod is specified by the C standard, and strtof appears
101in the stdlib.h of some systems, such as (at least some) Linux systems.
102The other functions write their results to their final argument(s):
103to the final two argument for the strtoI... (interval) functions,
104and to the final argument for the others (strtop... and strtor...).
105Where possible, these arguments have "natural" return types (double*
106or float*), to permit at least some type checking. In reality, they
107are viewed as arrays of ULong (or, for the "x" functions, UShort)
108values. On systems where long double is the appropriate type, one can
109pass long double* final argument(s) to these routines. The int value
110that these routines return is the return value from the call they make
111on strtodg; see the enum of possible return values in gdtoa.h.
112
113Source files g_ddfmt.c, misc.c, smisc.c, strtod.c, strtodg.c, and ulp.c
114should use true IEEE double arithmetic (not, e.g., double extended),
115at least for storing (and viewing the bits of) the variables declared
116"double" within them.
117
118One detail indicated in struct FPI is whether the target binary
119arithmetic departs from the IEEE standard by flushing denormalized
120numbers to 0. On systems that do this, the helper routines for
121conversion to double-double format (when compiled with
122Sudden_Underflow #defined) penalize the bottom of the exponent
123range so that they return a nonzero result only when the least
124significant bit of the less significant member of the pair of
125double values returned can be expressed as a normalized double
126value. An alternative would be to drop to 53-bit precision near
127the bottom of the exponent range. To get correct rounding, this
128would (in general) require two calls on strtodg (one specifying
129126-bit arithmetic, then, if necessary, one specifying 53-bit
130arithmetic).
131
132By default, the core routine strtodg and strtod set errno to ERANGE
133if the result overflows to +Infinity or underflows to 0. Compile
134these routines with NO_ERRNO #defined to inhibit errno assignments.
135
136Routine strtod is based on netlib's "dtoa.c from fp", and
137(f = strtod(s,se)) is more efficient for some conversions than, say,
138strtord(s,se,1,&f). Parts of strtod require true IEEE double
139arithmetic with the default rounding mode (round-to-nearest) and, on
140systems with IEEE extended-precision registers, double-precision
141(53-bit) rounding precision. If the machine uses (the equivalent of)
142Intel 80x87 arithmetic, the call
143 _control87(PC_53, MCW_PC);
144does this with many compilers. Whether this or another call is
145appropriate depends on the compiler; for this to work, it may be
146necessary to #include "float.h" or another system-dependent header
147file.
148
149Source file strtodnrp.c gives a strtod that does not require 53-bit
150rounding precision on systems (such as Intel IA32 systems) that may
151suffer double rounding due to use of extended-precision registers.
152For some conversions this variant of strtod is less efficient than the
153one in strtod.c when the latter is run with 53-bit rounding precision.
154
155The values that the strto* routines return for NaNs are determined by
156gd_qnan.h, which the makefile generates by running the program whose
157source is qnan.c. Note that the rules for distinguishing signaling
158from quiet NaNs are system-dependent. For cross-compilation, you need
159to determine arith.h and gd_qnan.h suitably, e.g., using the
160arithmetic of the target machine.
161
162C99's hexadecimal floating-point constants are recognized by the
163strto* routines (but this feature has not yet been heavily tested).
164Compiling with NO_HEX_FP #defined disables this feature.
165
166When compiled with -DINFNAN_CHECK, the strto* routines recognize C99's
167NaN and Infinity syntax. Moreover, unless No_Hex_NaN is #defined, the
168strto* routines also recognize C99's NaN(...) syntax: they accept
169(case insensitively) strings of the form NaN(x), where x is a string
170of hexadecimal digits and spaces; if there is only one string of
171hexadecimal digits, it is taken for the fraction bits of the resulting
172NaN; if there are two or more strings of hexadecimal digits, each
173string is assigned to the next available sequence of 32-bit words of
174fractions bits (starting with the most significant), right-aligned in
175each sequence.
176
177For binary -> decimal conversions, I've provided just one family
178of helper routines:
179
180 g_ffmt
181 g_dfmt
182 g_ddfmt
183 g_xfmt
184 g_xLfmt
185 g_Qfmt
186
187which do a "%g" style conversion either to a specified number of decimal
188places (if their ndig argument is positive), or to the shortest
189decimal string that rounds to the given binary floating-point value
190(if ndig <= 0). They write into a buffer supplied as an argument
191and return either a pointer to the end of the string (a null character)
192in the buffer, if the buffer was long enough, or 0. Other forms of
193conversion are easily done with the help of gdtoa(), such as %e or %f
194style and conversions with direction of rounding specified (so that, if
195desired, the decimal value is either >= or <= the binary value).
196On IEEE-arithmetic systems that provide the C99 fegetround() function,
197if compiled with -DHonor_FLT_ROUNDS, these routines honor the current
198rounding mode.
199
200For an example of more general conversions based on dtoa(), see
201netlib's "printf.c from ampl/solvers".
202
203For double-double -> decimal, g_ddfmt() assumes IEEE-like arithmetic
204of precision max(126, #bits(input)) bits, where #bits(input) is the
205number of mantissa bits needed to represent the sum of the two double
206values in the input.
207
208The makefile creates a library, gdtoa.a. To use the helper
209routines, a program only needs to include gdtoa.h. All the
210source files for gdtoa.a include a more extensive gdtoaimp.h;
211among other things, gdtoaimp.h has #defines that make "internal"
212names end in _D2A. To make a "system" library, one could modify
213these #defines to make the names start with __.
214
215Various comments about possible #defines appear in gdtoaimp.h,
216but for most purposes, arith.h should set suitable #defines.
217
218Systems with preemptive scheduling of multiple threads require some
219manual intervention. On such systems, it's necessary to compile
220dmisc.c, dtoa.c gdota.c, and misc.c with MULTIPLE_THREADS #defined,
221and to provide (or suitably #define) two locks, acquired by
222ACQUIRE_DTOA_LOCK(n) and freed by FREE_DTOA_LOCK(n) for n = 0 or 1.
223(The second lock, accessed in pow5mult, ensures lazy evaluation of
224only one copy of high powers of 5; omitting this lock would introduce
225a small probability of wasting memory, but would otherwise be harmless.)
226Routines that call dtoa or gdtoa directly must also invoke freedtoa(s)
227to free the value s returned by dtoa or gdtoa. It's OK to do so whether
228or not MULTIPLE_THREADS is #defined, and the helper g_*fmt routines
229listed above all do this indirectly (in gfmt_D2A(), which they all call).
230
231By default, there is a private pool of memory of length 2000 bytes
232for intermediate quantities, and MALLOC (see gdtoaimp.h) is called only
233if the private pool does not suffice. 2000 is large enough that MALLOC
234is called only under very unusual circumstances (decimal -> binary
235conversion of very long strings) for conversions to and from double
236precision. For systems with preemptively scheduled multiple threads
237or for conversions to extended or quad, it may be appropriate to
238#define PRIVATE_MEM nnnn, where nnnn is a suitable value > 2000.
239For extended and quad precisions, -DPRIVATE_MEM=20000 is probably
240plenty even for many digits at the ends of the exponent range.
241Use of the private pool avoids some overhead.
242
243Directory test provides some test routines. See its README.
244I've also tested this stuff (except double double conversions)
245with Vern Paxson's testbase program: see
246
247 V. Paxson and W. Kahan, "A Program for Testing IEEE Binary-Decimal
248 Conversion", manuscript, May 1991,
249 ftp://ftp.ee.lbl.gov/testbase-report.ps.Z .
250
251(The same ftp directory has source for testbase.)
252
253Some system-dependent additions to CFLAGS in the makefile:
254
255 HU-UX: -Aa -Ae
256 OSF (DEC Unix): -ieee_with_no_inexact
257 SunOS 4.1x: -DKR_headers -DBad_float_h
258
259If you want to put this stuff into a shared library and your
260operating system requires export lists for shared libraries,
261the following would be an appropriate export list:
262
263 dtoa
264 freedtoa
265 g_Qfmt
266 g_ddfmt
267 g_dfmt
268 g_ffmt
269 g_xLfmt
270 g_xfmt
271 gdtoa
272 strtoIQ
273 strtoId
274 strtoIdd
275 strtoIf
276 strtoIx
277 strtoIxL
278 strtod
279 strtodI
280 strtodg
281 strtof
282 strtopQ
283 strtopd
284 strtopdd
285 strtopf
286 strtopx
287 strtopxL
288 strtorQ
289 strtord
290 strtordd
291 strtorf
292 strtorx
293 strtorxL
294
295When time permits, I (dmg) hope to write in more detail about the
296present conversion routines; for now, this README file must suffice.
297Meanwhile, if you wish to write helper functions for other kinds of
298IEEE-like arithmetic, some explanation of struct FPI and the bits
299array may be helpful. Both gdtoa and strtodg operate on a bits array
300described by FPI *fpi. The bits array is of type ULong, a 32-bit
301unsigned integer type. Floating-point numbers have fpi->nbits bits,
302with the least significant 32 bits in bits[0], the next 32 bits in
303bits[1], etc. These numbers are regarded as integers multiplied by
3042^e (i.e., 2 to the power of the exponent e), where e is the second
305argument (be) to gdtoa and is stored in *exp by strtodg. The minimum
306and maximum exponent values fpi->emin and fpi->emax for normalized
307floating-point numbers reflect this arrangement. For example, the
308P754 standard for binary IEEE arithmetic specifies doubles as having
30953 bits, with normalized values of the form 1.xxxxx... times 2^(b-1023),
310with 52 bits (the x's) and the biased exponent b represented explicitly;
311b is an unsigned integer in the range 1 <= b <= 2046 for normalized
312finite doubles, b = 0 for denormals, and b = 2047 for Infinities and NaNs.
313To turn an IEEE double into the representation used by strtodg and gdtoa,
314we multiply 1.xxxx... by 2^52 (to make it an integer) and reduce the
315exponent e = (b-1023) by 52:
316
317 fpi->emin = 1 - 1023 - 52
318 fpi->emax = 1046 - 1023 - 52
319
320In various wrappers for IEEE double, we actually write -53 + 1 rather
321than -52, to emphasize that there are 53 bits including one implicit bit.
322Field fpi->rounding indicates the desired rounding direction, with
323possible values
324 FPI_Round_zero = toward 0,
325 FPI_Round_near = unbiased rounding -- the IEEE default,
326 FPI_Round_up = toward +Infinity, and
327 FPI_Round_down = toward -Infinity
328given in gdtoa.h.
329
330Field fpi->sudden_underflow indicates whether strtodg should return
331denormals or flush them to zero. Normal floating-point numbers have
332bit fpi->nbits in the bits array on. Denormals have it off, with
333exponent = fpi->emin. Strtodg provides distinct return values for normals
334and denormals; see gdtoa.h.
335
336Compiling g__fmt.c, strtod.c, and strtodg.c with -DUSE_LOCALE causes
337the decimal-point character to be taken from the current locale; otherwise
338it is '.'.
339
340Source files dtoa.c and strtod.c in this directory are derived from
341netlib's "dtoa.c from fp" and are meant to function equivalently.
342When compiled with Honor_FLT_ROUNDS #defined (on systems that provide
343FLT_ROUNDS and fegetround() as specified in the C99 standard), they
344honor the current rounding mode. Because FLT_ROUNDS is buggy on some
345(Linux) systems -- not reflecting calls on fesetround(), as the C99
346standard says it should -- when Honor_FLT_ROUNDS is #defined, the
347current rounding mode is obtained from fegetround() rather than from
348FLT_ROUNDS, unless Trust_FLT_ROUNDS is also #defined.
349
350Compile with -DUSE_LOCALE to use the current locale; otherwise
351decimal points are assumed to be '.'. With -DUSE_LOCALE, unless
352you also compile with -DNO_LOCALE_CACHE, the details about the
353current "decimal point" character string are cached and assumed not
354to change during the program's execution.
355
356Please send comments to David M. Gay (dmg at acm dot org, with " at "
357changed at "@" and " dot " changed to ".").
diff --git a/src/lib/evil/gdtoa/README.mingw b/src/lib/evil/gdtoa/README.mingw
deleted file mode 100644
index b41d1fdae0..0000000000
--- a/src/lib/evil/gdtoa/README.mingw
+++ /dev/null
@@ -1,20 +0,0 @@
1The gdtoa code here is based on David M. Gay's original
2gdtoa source at http://www.netlib.org/fp/ from Sep. 27,
32010. The major changes between the original source and
4the mingw port here include:
5
6* IBM, CRAY and VAX code removed.
7* KR_headers, ANSI, Void and Char ifdefs are removed.
8* gdtoa symbols are prepended with "__".
9* g_xfmt() uses __fpclassifyl() instead of interpreting
10 the flags bit-wise.
11* lo0bits() and hi0bits() of misc.c replaced by wrappers
12 to gcc's __builtin_clz()
13* The double/ulong union renamed from U to dbl_union
14 (grep'ped better..)
15* A few compiler warning fixes here and there.
16* A few other insignificant changes (if any..)
17
18MinGW specific compile-time definitions are at the top of
19gdtoaimp.h and gdtoa.h headers.
20
diff --git a/src/lib/evil/gdtoa/arithchk.c b/src/lib/evil/gdtoa/arithchk.c
deleted file mode 100644
index 901bab1fb3..0000000000
--- a/src/lib/evil/gdtoa/arithchk.c
+++ /dev/null
@@ -1,192 +0,0 @@
1/****************************************************************
2Copyright (C) 1997, 1998 Lucent Technologies
3All Rights Reserved
4
5Permission to use, copy, modify, and distribute this software and
6its documentation for any purpose and without fee is hereby
7granted, provided that the above copyright notice appear in all
8copies and that both that the copyright notice and this
9permission notice and warranty disclaimer appear in supporting
10documentation, and that the name of Lucent or any of its entities
11not be used in advertising or publicity pertaining to
12distribution of the software without specific, written prior
13permission.
14
15LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
16INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
17IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
18SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
19WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
20IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
21ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
22THIS SOFTWARE.
23****************************************************************/
24
25/* Try to deduce arith.h from arithmetic properties. */
26#ifdef MINGW_BUILD_GEN
27
28#include <stdio.h>
29
30 static int dalign;
31 typedef struct
32Akind {
33 char *name;
34 int kind;
35 } Akind;
36
37 static Akind
38IEEE_8087 = { "IEEE_8087", 1 },
39IEEE_MC68k = { "IEEE_MC68k", 2 },
40IBM = { "IBM", 3 },
41VAX = { "VAX", 4 },
42CRAY = { "CRAY", 5};
43
44 static Akind *
45Lcheck()
46{
47 union {
48 double d;
49 long L[2];
50 } u;
51 struct {
52 double d;
53 long L;
54 } x[2];
55
56 if (sizeof(x) > 2*(sizeof(double) + sizeof(long)))
57 dalign = 1;
58 u.L[0] = u.L[1] = 0;
59 u.d = 1e13;
60 if (u.L[0] == 1117925532 && u.L[1] == -448790528)
61 return &IEEE_MC68k;
62 if (u.L[1] == 1117925532 && u.L[0] == -448790528)
63 return &IEEE_8087;
64 if (u.L[0] == -2065213935 && u.L[1] == 10752)
65 return &VAX;
66 if (u.L[0] == 1267827943 && u.L[1] == 704643072)
67 return &IBM;
68 return 0;
69 }
70
71 static Akind *
72icheck()
73{
74 union {
75 double d;
76 int L[2];
77 } u;
78 struct {
79 double d;
80 int L;
81 } x[2];
82
83 if (sizeof(x) > 2*(sizeof(double) + sizeof(int)))
84 dalign = 1;
85 u.L[0] = u.L[1] = 0;
86 u.d = 1e13;
87 if (u.L[0] == 1117925532 && u.L[1] == -448790528)
88 return &IEEE_MC68k;
89 if (u.L[1] == 1117925532 && u.L[0] == -448790528)
90 return &IEEE_8087;
91 if (u.L[0] == -2065213935 && u.L[1] == 10752)
92 return &VAX;
93 if (u.L[0] == 1267827943 && u.L[1] == 704643072)
94 return &IBM;
95 return 0;
96 }
97
98char *emptyfmt = ""; /* avoid possible warning message with printf("") */
99
100 static Akind *
101ccheck()
102{
103 union {
104 double d;
105 long L;
106 } u;
107 long Cray1;
108
109 /* Cray1 = 4617762693716115456 -- without overflow on non-Crays */
110 Cray1 = printf(emptyfmt) < 0 ? 0 : 4617762;
111 if (printf(emptyfmt, Cray1) >= 0)
112 Cray1 = 1000000*Cray1 + 693716;
113 if (printf(emptyfmt, Cray1) >= 0)
114 Cray1 = 1000000*Cray1 + 115456;
115 u.d = 1e13;
116 if (u.L == Cray1)
117 return &CRAY;
118 return 0;
119 }
120
121 static int
122fzcheck()
123{
124 double a, b;
125 int i;
126
127 a = 1.;
128 b = .1;
129 for(i = 155;; b *= b, i >>= 1) {
130 if (i & 1) {
131 a *= b;
132 if (i == 1)
133 break;
134 }
135 }
136 b = a * a;
137 return b == 0.;
138 }
139
140 int
141main()
142{
143 Akind *a = 0;
144 int Ldef = 0;
145 FILE *f;
146
147#ifdef WRITE_ARITH_H /* for Symantec's buggy "make" */
148 f = fopen("arith.h", "w");
149 if (!f) {
150 printf("Cannot open arith.h\n");
151 return 1;
152 }
153#else
154 f = stdout;
155#endif
156
157 if (sizeof(double) == 2*sizeof(long))
158 a = Lcheck();
159 else if (sizeof(double) == 2*sizeof(int)) {
160 Ldef = 1;
161 a = icheck();
162 }
163 else if (sizeof(double) == sizeof(long))
164 a = ccheck();
165 if (a) {
166 fprintf(f, "#define %s\n#define Arith_Kind_ASL %d\n",
167 a->name, a->kind);
168 if (Ldef)
169 fprintf(f, "#define Long int\n#define Intcast (int)(long)\n");
170 if (dalign)
171 fprintf(f, "#define Double_Align\n");
172 if (sizeof(char*) == 8)
173 fprintf(f, "#define X64_bit_pointers\n");
174#ifndef NO_LONG_LONG
175 if (sizeof(long long) < 8)
176#endif
177 fprintf(f, "#define NO_LONG_LONG\n");
178 if (a->kind <= 2 && fzcheck())
179 fprintf(f, "#define Sudden_Underflow\n");
180#ifdef WRITE_ARITH_H /* for Symantec's buggy "make" */
181 fclose(f);
182#endif
183 return 0;
184 }
185 fprintf(f, "/* Unknown arithmetic */\n");
186#ifdef WRITE_ARITH_H /* for Symantec's buggy "make" */
187 fclose(f);
188#endif
189 return 1;
190 }
191#endif /* MINGW_BUILD_GEN */
192
diff --git a/src/lib/evil/gdtoa/dmisc.c b/src/lib/evil/gdtoa/dmisc.c
deleted file mode 100644
index c61e9fa8a5..0000000000
--- a/src/lib/evil/gdtoa/dmisc.c
+++ /dev/null
@@ -1,196 +0,0 @@
1/****************************************************************
2
3The author of this software is David M. Gay.
4
5Copyright (C) 1998 by Lucent Technologies
6All Rights Reserved
7
8Permission to use, copy, modify, and distribute this software and
9its documentation for any purpose and without fee is hereby
10granted, provided that the above copyright notice appear in all
11copies and that both that the copyright notice and this
12permission notice and warranty disclaimer appear in supporting
13documentation, and that the name of Lucent or any of its entities
14not be used in advertising or publicity pertaining to
15distribution of the software without specific, written prior
16permission.
17
18LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
25THIS SOFTWARE.
26
27****************************************************************/
28
29/* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
31
32#include "gdtoaimp.h"
33
34#ifndef MULTIPLE_THREADS
35char *dtoa_result;
36#endif
37
38char *rv_alloc (int i)
39{
40 int j, k, *r;
41
42 j = sizeof(ULong);
43 for(k = 0;
44 (int) (sizeof(Bigint) - sizeof(ULong) - sizeof(int)) + j <= i;
45 j <<= 1)
46 k++;
47 r = (int*)Balloc(k);
48 *r = k;
49 return
50#ifndef MULTIPLE_THREADS
51 dtoa_result =
52#endif
53 (char *)(r+1);
54}
55
56char *nrv_alloc (char *s, char **rve, int n)
57{
58 char *rv, *t;
59
60 t = rv = rv_alloc(n);
61 while((*t = *s++) !=0)
62 t++;
63 if (rve)
64 *rve = t;
65 return rv;
66}
67
68/* freedtoa(s) must be used to free values s returned by dtoa
69 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
70 * but for consistency with earlier versions of dtoa, it is optional
71 * when MULTIPLE_THREADS is not defined.
72 */
73
74void __freedtoa (char *s)
75{
76 Bigint *b = (Bigint *)((int *)s - 1);
77 b->maxwds = 1 << (b->k = *(int*)b);
78 Bfree(b);
79#ifndef MULTIPLE_THREADS
80 if (s == dtoa_result)
81 dtoa_result = 0;
82#endif
83}
84
85int quorem (Bigint *b, Bigint *S)
86{
87 int n;
88 ULong *bx, *bxe, q, *sx, *sxe;
89#ifdef ULLong
90 ULLong borrow, carry, y, ys;
91#else
92 ULong borrow, carry, y, ys;
93#ifdef Pack_32
94 ULong si, z, zs;
95#endif
96#endif
97
98 n = S->wds;
99#ifdef DEBUG
100 /*debug*/ if (b->wds > n)
101 /*debug*/ Bug("oversize b in quorem");
102#endif
103 if (b->wds < n)
104 return 0;
105 sx = S->x;
106 sxe = sx + --n;
107 bx = b->x;
108 bxe = bx + n;
109 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
110#ifdef DEBUG
111 /*debug*/ if (q > 9)
112 /*debug*/ Bug("oversized quotient in quorem");
113#endif
114 if (q) {
115 borrow = 0;
116 carry = 0;
117 do {
118#ifdef ULLong
119 ys = *sx++ * (ULLong)q + carry;
120 carry = ys >> 32;
121 y = *bx - (ys & 0xffffffffUL) - borrow;
122 borrow = y >> 32 & 1UL;
123 *bx++ = y & 0xffffffffUL;
124#else
125#ifdef Pack_32
126 si = *sx++;
127 ys = (si & 0xffff) * q + carry;
128 zs = (si >> 16) * q + (ys >> 16);
129 carry = zs >> 16;
130 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
131 borrow = (y & 0x10000) >> 16;
132 z = (*bx >> 16) - (zs & 0xffff) - borrow;
133 borrow = (z & 0x10000) >> 16;
134 Storeinc(bx, z, y);
135#else
136 ys = *sx++ * q + carry;
137 carry = ys >> 16;
138 y = *bx - (ys & 0xffff) - borrow;
139 borrow = (y & 0x10000) >> 16;
140 *bx++ = y & 0xffff;
141#endif
142#endif
143 } while(sx <= sxe);
144
145 if (!*bxe) {
146 bx = b->x;
147 while(--bxe > bx && !*bxe)
148 --n;
149 b->wds = n;
150 }
151 }
152
153 if (cmp(b, S) >= 0) {
154 q++;
155 borrow = 0;
156 carry = 0;
157 bx = b->x;
158 sx = S->x;
159 do {
160#ifdef ULLong
161 ys = *sx++ + carry;
162 carry = ys >> 32;
163 y = *bx - (ys & 0xffffffffUL) - borrow;
164 borrow = y >> 32 & 1UL;
165 *bx++ = y & 0xffffffffUL;
166#else
167#ifdef Pack_32
168 si = *sx++;
169 ys = (si & 0xffff) + carry;
170 zs = (si >> 16) + (ys >> 16);
171 carry = zs >> 16;
172 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
173 borrow = (y & 0x10000) >> 16;
174 z = (*bx >> 16) - (zs & 0xffff) - borrow;
175 borrow = (z & 0x10000) >> 16;
176 Storeinc(bx, z, y);
177#else
178 ys = *sx++ + carry;
179 carry = ys >> 16;
180 y = *bx - (ys & 0xffff) - borrow;
181 borrow = (y & 0x10000) >> 16;
182 *bx++ = y & 0xffff;
183#endif
184#endif
185 } while(sx <= sxe);
186
187 bx = b->x;
188 bxe = bx + n;
189 if (!*bxe) {
190 while(--bxe > bx && !*bxe)
191 --n;
192 b->wds = n;
193 }
194 }
195 return q;
196}
diff --git a/src/lib/evil/gdtoa/dtoa.c b/src/lib/evil/gdtoa/dtoa.c
deleted file mode 100644
index 2906bd99f6..0000000000
--- a/src/lib/evil/gdtoa/dtoa.c
+++ /dev/null
@@ -1,750 +0,0 @@
1/****************************************************************
2
3The author of this software is David M. Gay.
4
5Copyright (C) 1998, 1999 by Lucent Technologies
6All Rights Reserved
7
8Permission to use, copy, modify, and distribute this software and
9its documentation for any purpose and without fee is hereby
10granted, provided that the above copyright notice appear in all
11copies and that both that the copyright notice and this
12permission notice and warranty disclaimer appear in supporting
13documentation, and that the name of Lucent or any of its entities
14not be used in advertising or publicity pertaining to
15distribution of the software without specific, written prior
16permission.
17
18LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
25THIS SOFTWARE.
26
27****************************************************************/
28
29/* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
31
32#include "gdtoaimp.h"
33
34/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
35 *
36 * Inspired by "How to Print Floating-Point Numbers Accurately" by
37 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
38 *
39 * Modifications:
40 * 1. Rather than iterating, we use a simple numeric overestimate
41 * to determine k = floor(log10(d)). We scale relevant
42 * quantities using O(log2(k)) rather than O(k) multiplications.
43 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
44 * try to generate digits strictly left to right. Instead, we
45 * compute with fewer bits and propagate the carry if necessary
46 * when rounding the final digit up. This is often faster.
47 * 3. Under the assumption that input will be rounded nearest,
48 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
49 * That is, we allow equality in stopping tests when the
50 * round-nearest rule will give the same floating-point value
51 * as would satisfaction of the stopping test with strict
52 * inequality.
53 * 4. We remove common factors of powers of 2 from relevant
54 * quantities.
55 * 5. When converting floating-point integers less than 1e16,
56 * we use floating-point arithmetic rather than resorting
57 * to multiple-precision integers.
58 * 6. When asked to produce fewer than 15 digits, we first try
59 * to get by with floating-point arithmetic; we resort to
60 * multiple-precision integer arithmetic only if we cannot
61 * guarantee that the floating-point calculation has given
62 * the correctly rounded result. For k requested digits and
63 * "uniformly" distributed input, the probability is
64 * something like 10^(k-15) that we must resort to the Long
65 * calculation.
66 */
67
68#ifdef Honor_FLT_ROUNDS
69#undef Check_FLT_ROUNDS
70#define Check_FLT_ROUNDS
71#else
72#define Rounding Flt_Rounds
73#endif
74
75char *__dtoa (double d0, int mode, int ndigits, int *decpt, int *sign, char **rve)
76{
77 /* Arguments ndigits, decpt, sign are similar to those
78 of ecvt and fcvt; trailing zeros are suppressed from
79 the returned string. If not null, *rve is set to point
80 to the end of the return value. If d is +-Infinity or NaN,
81 then *decpt is set to 9999.
82
83 mode:
84 0 ==> shortest string that yields d when read in
85 and rounded to nearest.
86 1 ==> like 0, but with Steele & White stopping rule;
87 e.g. with IEEE P754 arithmetic , mode 0 gives
88 1e23 whereas mode 1 gives 9.999999999999999e22.
89 2 ==> max(1,ndigits) significant digits. This gives a
90 return value similar to that of ecvt, except
91 that trailing zeros are suppressed.
92 3 ==> through ndigits past the decimal point. This
93 gives a return value similar to that from fcvt,
94 except that trailing zeros are suppressed, and
95 ndigits can be negative.
96 4,5 ==> similar to 2 and 3, respectively, but (in
97 round-nearest mode) with the tests of mode 0 to
98 possibly return a shorter string that rounds to d.
99 With IEEE arithmetic and compilation with
100 -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
101 as modes 2 and 3 when FLT_ROUNDS != 1.
102 6-9 ==> Debugging modes similar to mode - 4: don't try
103 fast floating-point estimate (if applicable).
104
105 Values of mode other than 0-9 are treated as mode 0.
106
107 Sufficient space is allocated to the return value
108 to hold the suppressed trailing zeros.
109 */
110
111 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
112 j, j2, k, k0, k_check, leftright, m2, m5, s2, s5,
113 spec_case, try_quick;
114 Long L;
115#ifndef Sudden_Underflow
116 int denorm;
117 ULong x;
118#endif
119 Bigint *b, *b1, *delta, *mlo, *mhi, *S;
120 union _dbl_union d, d2, eps;
121 double ds;
122 char *s, *s0;
123#ifdef SET_INEXACT
124 int inexact, oldinexact;
125#endif
126#ifdef Honor_FLT_ROUNDS /*{*/
127 int Rounding;
128#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
129 Rounding = Flt_Rounds;
130#else /*}{*/
131 Rounding = 1;
132 switch(fegetround()) {
133 case FE_TOWARDZERO: Rounding = 0; break;
134 case FE_UPWARD: Rounding = 2; break;
135 case FE_DOWNWARD: Rounding = 3;
136 }
137#endif /*}}*/
138#endif /*}*/
139
140#ifndef MULTIPLE_THREADS
141 if (dtoa_result) {
142 __freedtoa(dtoa_result);
143 dtoa_result = 0;
144 }
145#endif
146 d.d = d0;
147 if (word0(&d) & Sign_bit) {
148 /* set sign for everything, including 0's and NaNs */
149 *sign = 1;
150 word0(&d) &= ~Sign_bit; /* clear sign bit */
151 }
152 else
153 *sign = 0;
154
155 if ((word0(&d) & Exp_mask) == Exp_mask)
156 {
157 /* Infinity or NaN */
158 *decpt = 9999;
159 if (!word1(&d) && !(word0(&d) & 0xfffff))
160 return nrv_alloc("Infinity", rve, 8);
161 return nrv_alloc("NaN", rve, 3);
162 }
163 if (!dval(&d)) {
164 *decpt = 1;
165 return nrv_alloc("0", rve, 1);
166 }
167
168#ifdef SET_INEXACT
169 try_quick = oldinexact = get_inexact();
170 inexact = 1;
171#endif
172#ifdef Honor_FLT_ROUNDS
173 if (Rounding >= 2) {
174 if (*sign)
175 Rounding = Rounding == 2 ? 0 : 2;
176 else
177 if (Rounding != 2)
178 Rounding = 0;
179 }
180#endif
181
182 b = d2b(dval(&d), &be, &bbits);
183#ifdef Sudden_Underflow
184 i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
185#else
186 if (( i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
187#endif
188 dval(&d2) = dval(&d);
189 word0(&d2) &= Frac_mask1;
190 word0(&d2) |= Exp_11;
191
192 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
193 * log10(x) = log(x) / log(10)
194 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
195 * log10(&d) = (i-Bias)*log(2)/log(10) + log10(&d2)
196 *
197 * This suggests computing an approximation k to log10(&d) by
198 *
199 * k = (i - Bias)*0.301029995663981
200 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
201 *
202 * We want k to be too large rather than too small.
203 * The error in the first-order Taylor series approximation
204 * is in our favor, so we just round up the constant enough
205 * to compensate for any error in the multiplication of
206 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
207 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
208 * adding 1e-13 to the constant term more than suffices.
209 * Hence we adjust the constant term to 0.1760912590558.
210 * (We could get a more accurate k by invoking log10,
211 * but this is probably not worthwhile.)
212 */
213
214 i -= Bias;
215#ifndef Sudden_Underflow
216 denorm = 0;
217 }
218 else {
219 /* d is denormalized */
220
221 i = bbits + be + (Bias + (P-1) - 1);
222 x = i > 32 ? word0(&d) << (64 - i) | word1(&d) >> (i - 32)
223 : word1(&d) << (32 - i);
224 dval(&d2) = x;
225 word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
226 i -= (Bias + (P-1) - 1) + 1;
227 denorm = 1;
228 }
229#endif
230 ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
231 k = (int)ds;
232 if (ds < 0. && ds != k)
233 k--; /* want k = floor(ds) */
234 k_check = 1;
235 if (k >= 0 && k <= Ten_pmax) {
236 if (dval(&d) < tens[k])
237 k--;
238 k_check = 0;
239 }
240 j = bbits - i - 1;
241 if (j >= 0) {
242 b2 = 0;
243 s2 = j;
244 }
245 else {
246 b2 = -j;
247 s2 = 0;
248 }
249 if (k >= 0) {
250 b5 = 0;
251 s5 = k;
252 s2 += k;
253 }
254 else {
255 b2 -= k;
256 b5 = -k;
257 s5 = 0;
258 }
259 if (mode < 0 || mode > 9)
260 mode = 0;
261
262#ifndef SET_INEXACT
263#ifdef Check_FLT_ROUNDS
264 try_quick = Rounding == 1;
265#else
266 try_quick = 1;
267#endif
268#endif /*SET_INEXACT*/
269
270 if (mode > 5) {
271 mode -= 4;
272 try_quick = 0;
273 }
274 leftright = 1;
275 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
276 /* silence erroneous "gcc -Wall" warning. */
277 switch(mode) {
278 case 0:
279 case 1:
280 i = 18;
281 ndigits = 0;
282 break;
283 case 2:
284 leftright = 0;
285 /* no break */
286 case 4:
287 if (ndigits <= 0)
288 ndigits = 1;
289 ilim = ilim1 = i = ndigits;
290 break;
291 case 3:
292 leftright = 0;
293 /* no break */
294 case 5:
295 i = ndigits + k + 1;
296 ilim = i;
297 ilim1 = i - 1;
298 if (i <= 0)
299 i = 1;
300 }
301 s = s0 = rv_alloc(i);
302
303#ifdef Honor_FLT_ROUNDS
304 if (mode > 1 && Rounding != 1)
305 leftright = 0;
306#endif
307
308 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
309
310 /* Try to get by with floating-point arithmetic. */
311
312 i = 0;
313 dval(&d2) = dval(&d);
314 k0 = k;
315 ilim0 = ilim;
316 ieps = 2; /* conservative */
317 if (k > 0) {
318 ds = tens[k&0xf];
319 j = k >> 4;
320 if (j & Bletch) {
321 /* prevent overflows */
322 j &= Bletch - 1;
323 dval(&d) /= bigtens[n_bigtens-1];
324 ieps++;
325 }
326 for(; j; j >>= 1, i++)
327 if (j & 1) {
328 ieps++;
329 ds *= bigtens[i];
330 }
331 dval(&d) /= ds;
332 }
333 else if (( j2 = -k )!=0) {
334 dval(&d) *= tens[j2 & 0xf];
335 for(j = j2 >> 4; j; j >>= 1, i++)
336 if (j & 1) {
337 ieps++;
338 dval(&d) *= bigtens[i];
339 }
340 }
341 if (k_check && dval(&d) < 1. && ilim > 0) {
342 if (ilim1 <= 0)
343 goto fast_failed;
344 ilim = ilim1;
345 k--;
346 dval(&d) *= 10.;
347 ieps++;
348 }
349 dval(&eps) = ieps*dval(&d) + 7.;
350 word0(&eps) -= (P-1)*Exp_msk1;
351 if (ilim == 0) {
352 S = mhi = 0;
353 dval(&d) -= 5.;
354 if (dval(&d) > dval(&eps))
355 goto one_digit;
356 if (dval(&d) < -dval(&eps))
357 goto no_digits;
358 goto fast_failed;
359 }
360#ifndef No_leftright
361 if (leftright) {
362 /* Use Steele & White method of only
363 * generating digits needed.
364 */
365 dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
366 for(i = 0;;) {
367 L = dval(&d);
368 dval(&d) -= L;
369 *s++ = '0' + (int)L;
370 if (dval(&d) < dval(&eps))
371 goto ret1;
372 if (1. - dval(&d) < dval(&eps))
373 goto bump_up;
374 if (++i >= ilim)
375 break;
376 dval(&eps) *= 10.;
377 dval(&d) *= 10.;
378 }
379 }
380 else {
381#endif
382 /* Generate ilim digits, then fix them up. */
383 dval(&eps) *= tens[ilim-1];
384 for(i = 1;; i++, dval(&d) *= 10.) {
385 L = (Long)(dval(&d));
386 if (!(dval(&d) -= L))
387 ilim = i;
388 *s++ = '0' + (int)L;
389 if (i == ilim) {
390 if (dval(&d) > 0.5 + dval(&eps))
391 goto bump_up;
392 else if (dval(&d) < 0.5 - dval(&eps)) {
393 while(*--s == '0');
394 s++;
395 goto ret1;
396 }
397 break;
398 }
399 }
400#ifndef No_leftright
401 }
402#endif
403 fast_failed:
404 s = s0;
405 dval(&d) = dval(&d2);
406 k = k0;
407 ilim = ilim0;
408 }
409
410 /* Do we have a "small" integer? */
411
412 if (be >= 0 && k <= Int_max) {
413 /* Yes. */
414 ds = tens[k];
415 if (ndigits < 0 && ilim <= 0) {
416 S = mhi = 0;
417 if (ilim < 0 || dval(&d) <= 5*ds)
418 goto no_digits;
419 goto one_digit;
420 }
421 for(i = 1;; i++, dval(&d) *= 10.) {
422 L = (Long)(dval(&d) / ds);
423 dval(&d) -= L*ds;
424#ifdef Check_FLT_ROUNDS
425 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
426 if (dval(&d) < 0) {
427 L--;
428 dval(&d) += ds;
429 }
430#endif
431 *s++ = '0' + (int)L;
432 if (!dval(&d)) {
433#ifdef SET_INEXACT
434 inexact = 0;
435#endif
436 break;
437 }
438 if (i == ilim) {
439#ifdef Honor_FLT_ROUNDS
440 if (mode > 1)
441 switch(Rounding) {
442 case 0: goto ret1;
443 case 2: goto bump_up;
444 }
445#endif
446 dval(&d) += dval(&d);
447#ifdef ROUND_BIASED
448 if (dval(&d) >= ds)
449#else
450 if (dval(&d) > ds || (dval(&d) == ds && L & 1))
451#endif
452 {
453 bump_up:
454 while(*--s == '9')
455 if (s == s0) {
456 k++;
457 *s = '0';
458 break;
459 }
460 ++*s++;
461 }
462 break;
463 }
464 }
465 goto ret1;
466 }
467
468 m2 = b2;
469 m5 = b5;
470 mhi = mlo = 0;
471 if (leftright) {
472 i =
473#ifndef Sudden_Underflow
474 denorm ? be + (Bias + (P-1) - 1 + 1) :
475#endif
476 1 + P - bbits;
477 b2 += i;
478 s2 += i;
479 mhi = i2b(1);
480 }
481 if (m2 > 0 && s2 > 0) {
482 i = m2 < s2 ? m2 : s2;
483 b2 -= i;
484 m2 -= i;
485 s2 -= i;
486 }
487 if (b5 > 0) {
488 if (leftright) {
489 if (m5 > 0) {
490 mhi = pow5mult(mhi, m5);
491 b1 = mult(mhi, b);
492 Bfree(b);
493 b = b1;
494 }
495 if (( j = b5 - m5 )!=0)
496 b = pow5mult(b, j);
497 }
498 else
499 b = pow5mult(b, b5);
500 }
501 S = i2b(1);
502 if (s5 > 0)
503 S = pow5mult(S, s5);
504
505 /* Check for special case that d is a normalized power of 2. */
506
507 spec_case = 0;
508 if ((mode < 2 || leftright)
509#ifdef Honor_FLT_ROUNDS
510 && Rounding == 1
511#endif
512 ) {
513 if (!word1(&d) && !(word0(&d) & Bndry_mask)
514#ifndef Sudden_Underflow
515 && word0(&d) & (Exp_mask & ~Exp_msk1)
516#endif
517 ) {
518 /* The special case */
519 b2 += Log2P;
520 s2 += Log2P;
521 spec_case = 1;
522 }
523 }
524
525 /* Arrange for convenient computation of quotients:
526 * shift left if necessary so divisor has 4 leading 0 bits.
527 *
528 * Perhaps we should just compute leading 28 bits of S once
529 * and for all and pass them and a shift to quorem, so it
530 * can do shifts and ors to compute the numerator for q.
531 */
532#ifdef Pack_32
533 if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f )!=0)
534 i = 32 - i;
535#else
536 if (( i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf )!=0)
537 i = 16 - i;
538#endif
539 if (i > 4) {
540 i -= 4;
541 b2 += i;
542 m2 += i;
543 s2 += i;
544 }
545 else if (i < 4) {
546 i += 28;
547 b2 += i;
548 m2 += i;
549 s2 += i;
550 }
551 if (b2 > 0)
552 b = lshift(b, b2);
553 if (s2 > 0)
554 S = lshift(S, s2);
555 if (k_check) {
556 if (cmp(b,S) < 0) {
557 k--;
558 b = multadd(b, 10, 0); /* we botched the k estimate */
559 if (leftright)
560 mhi = multadd(mhi, 10, 0);
561 ilim = ilim1;
562 }
563 }
564 if (ilim <= 0 && (mode == 3 || mode == 5)) {
565 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
566 /* no digits, fcvt style */
567 no_digits:
568 k = -1 - ndigits;
569 goto ret;
570 }
571 one_digit:
572 *s++ = '1';
573 k++;
574 goto ret;
575 }
576 if (leftright) {
577 if (m2 > 0)
578 mhi = lshift(mhi, m2);
579
580 /* Compute mlo -- check for special case
581 * that d is a normalized power of 2.
582 */
583
584 mlo = mhi;
585 if (spec_case) {
586 mhi = Balloc(mhi->k);
587 Bcopy(mhi, mlo);
588 mhi = lshift(mhi, Log2P);
589 }
590
591 for(i = 1;;i++) {
592 dig = quorem(b,S) + '0';
593 /* Do we yet have the shortest decimal string
594 * that will round to d?
595 */
596 j = cmp(b, mlo);
597 delta = diff(S, mhi);
598 j2 = delta->sign ? 1 : cmp(b, delta);
599 Bfree(delta);
600#ifndef ROUND_BIASED
601 if (j2 == 0 && mode != 1 && !(word1(&d) & 1)
602#ifdef Honor_FLT_ROUNDS
603 && Rounding >= 1
604#endif
605 ) {
606 if (dig == '9')
607 goto round_9_up;
608 if (j > 0)
609 dig++;
610#ifdef SET_INEXACT
611 else if (!b->x[0] && b->wds <= 1)
612 inexact = 0;
613#endif
614 *s++ = dig;
615 goto ret;
616 }
617#endif
618 if (j < 0 || (j == 0 && mode != 1
619#ifndef ROUND_BIASED
620 && !(word1(&d) & 1)
621#endif
622 )) {
623 if (!b->x[0] && b->wds <= 1) {
624#ifdef SET_INEXACT
625 inexact = 0;
626#endif
627 goto accept_dig;
628 }
629#ifdef Honor_FLT_ROUNDS
630 if (mode > 1)
631 switch(Rounding) {
632 case 0: goto accept_dig;
633 case 2: goto keep_dig;
634 }
635#endif /*Honor_FLT_ROUNDS*/
636 if (j2 > 0) {
637 b = lshift(b, 1);
638 j2 = cmp(b, S);
639#ifdef ROUND_BIASED
640 if (j2 >= 0 /*)*/
641#else
642 if ((j2 > 0 || (j2 == 0 && dig & 1))
643#endif
644 && dig++ == '9')
645 goto round_9_up;
646 }
647 accept_dig:
648 *s++ = dig;
649 goto ret;
650 }
651 if (j2 > 0) {
652#ifdef Honor_FLT_ROUNDS
653 if (!Rounding)
654 goto accept_dig;
655#endif
656 if (dig == '9') { /* possible if i == 1 */
657 round_9_up:
658 *s++ = '9';
659 goto roundoff;
660 }
661 *s++ = dig + 1;
662 goto ret;
663 }
664#ifdef Honor_FLT_ROUNDS
665 keep_dig:
666#endif
667 *s++ = dig;
668 if (i == ilim)
669 break;
670 b = multadd(b, 10, 0);
671 if (mlo == mhi)
672 mlo = mhi = multadd(mhi, 10, 0);
673 else {
674 mlo = multadd(mlo, 10, 0);
675 mhi = multadd(mhi, 10, 0);
676 }
677 }
678 }
679 else
680 for(i = 1;; i++) {
681 *s++ = dig = quorem(b,S) + '0';
682 if (!b->x[0] && b->wds <= 1) {
683#ifdef SET_INEXACT
684 inexact = 0;
685#endif
686 goto ret;
687 }
688 if (i >= ilim)
689 break;
690 b = multadd(b, 10, 0);
691 }
692
693 /* Round off last digit */
694
695#ifdef Honor_FLT_ROUNDS
696 switch(Rounding) {
697 case 0: goto trimzeros;
698 case 2: goto roundoff;
699 }
700#endif
701 b = lshift(b, 1);
702 j = cmp(b, S);
703#ifdef ROUND_BIASED
704 if (j >= 0)
705#else
706 if (j > 0 || (j == 0 && dig & 1))
707#endif
708 {
709 roundoff:
710 while(*--s == '9')
711 if (s == s0) {
712 k++;
713 *s++ = '1';
714 goto ret;
715 }
716 ++*s++;
717 }
718 else {
719#ifdef Honor_FLT_ROUNDS
720 trimzeros:
721#endif
722 while(*--s == '0');
723 s++;
724 }
725 ret:
726 Bfree(S);
727 if (mhi) {
728 if (mlo && mlo != mhi)
729 Bfree(mlo);
730 Bfree(mhi);
731 }
732 ret1:
733#ifdef SET_INEXACT
734 if (inexact) {
735 if (!oldinexact) {
736 word0(&d) = Exp_1 + (70 << Exp_shift);
737 word1(&d) = 0;
738 dval(&d) += 1.;
739 }
740 }
741 else if (!oldinexact)
742 clear_inexact();
743#endif
744 Bfree(b);
745 *s = 0;
746 *decpt = k + 1;
747 if (rve)
748 *rve = s;
749 return s0;
750}
diff --git a/src/lib/evil/gdtoa/g__fmt.c b/src/lib/evil/gdtoa/g__fmt.c
deleted file mode 100644
index 49bd95a845..0000000000
--- a/src/lib/evil/gdtoa/g__fmt.c
+++ /dev/null
@@ -1,142 +0,0 @@
1/****************************************************************
2
3The author of this software is David M. Gay.
4
5Copyright (C) 1998 by Lucent Technologies
6All Rights Reserved
7
8Permission to use, copy, modify, and distribute this software and
9its documentation for any purpose and without fee is hereby
10granted, provided that the above copyright notice appear in all
11copies and that both that the copyright notice and this
12permission notice and warranty disclaimer appear in supporting
13documentation, and that the name of Lucent or any of its entities
14not be used in advertising or publicity pertaining to
15distribution of the software without specific, written prior
16permission.
17
18LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
25THIS SOFTWARE.
26
27****************************************************************/
28
29/* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
31
32#include "gdtoaimp.h"
33
34#ifdef USE_LOCALE
35#include "locale.h"
36#endif
37
38char *__g__fmt (char *b, char *s, char *se, int decpt, ULong sign, size_t blen)
39{
40 int i, j, k;
41 char *be, *s0;
42 size_t len;
43#ifdef USE_LOCALE
44#ifdef NO_LOCALE_CACHE
45 char *decimalpoint = localeconv()->decimal_point;
46 size_t dlen = strlen(decimalpoint);
47#else
48 char *decimalpoint;
49 static char *decimalpoint_cache;
50 static size_t dlen;
51 if (!(s0 = decimalpoint_cache)) {
52 s0 = localeconv()->decimal_point;
53 dlen = strlen(s0);
54 if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
55 strcpy(decimalpoint_cache, s0);
56 s0 = decimalpoint_cache;
57 }
58 }
59 decimalpoint = s0;
60#endif
61#else
62#define dlen 0
63#endif
64 s0 = s;
65 len = (se-s) + dlen + 6; /* 6 = sign + e+dd + trailing null */
66 if (blen < len)
67 goto ret0;
68 be = b + blen - 1;
69 if (sign)
70 *b++ = '-';
71 if (decpt <= -4 || decpt > se - s + 5) {
72 *b++ = *s++;
73 if (*s) {
74#ifdef USE_LOCALE
75 while((*b = *decimalpoint++))
76 ++b;
77#else
78 *b++ = '.';
79#endif
80 while((*b = *s++) !=0)
81 b++;
82 }
83 *b++ = 'e';
84 /* sprintf(b, "%+.2d", decpt - 1); */
85 if (--decpt < 0) {
86 *b++ = '-';
87 decpt = -decpt;
88 }
89 else
90 *b++ = '+';
91 for(j = 2, k = 10; 10*k <= decpt; j++, k *= 10){}
92 for(;;) {
93 i = decpt / k;
94 if (b >= be)
95 goto ret0;
96 *b++ = i + '0';
97 if (--j <= 0)
98 break;
99 decpt -= i*k;
100 decpt *= 10;
101 }
102 *b = 0;
103 }
104 else if (decpt <= 0) {
105#ifdef USE_LOCALE
106 while((*b = *decimalpoint++))
107 ++b;
108#else
109 *b++ = '.';
110#endif
111 if (be < b - decpt + (se - s))
112 goto ret0;
113 for(; decpt < 0; decpt++)
114 *b++ = '0';
115 while((*b = *s++) != 0)
116 b++;
117 }
118 else {
119 while((*b = *s++) != 0) {
120 b++;
121 if (--decpt == 0 && *s) {
122#ifdef USE_LOCALE
123 while((*b = *decimalpoint++))
124 ++b;
125#else
126 *b++ = '.';
127#endif
128 }
129 }
130 if (b + decpt > be) {
131 ret0:
132 b = 0;
133 goto ret;
134 }
135 for(; decpt > 0; decpt--)
136 *b++ = '0';
137 *b = 0;
138 }
139 ret:
140 __freedtoa(s0);
141 return b;
142}
diff --git a/src/lib/evil/gdtoa/g_dfmt.c b/src/lib/evil/gdtoa/g_dfmt.c
deleted file mode 100644
index 50ed708a6c..0000000000
--- a/src/lib/evil/gdtoa/g_dfmt.c
+++ /dev/null
@@ -1,90 +0,0 @@
1/****************************************************************
2
3The author of this software is David M. Gay.
4
5Copyright (C) 1998 by Lucent Technologies
6All Rights Reserved
7
8Permission to use, copy, modify, and distribute this software and
9its documentation for any purpose and without fee is hereby
10granted, provided that the above copyright notice appear in all
11copies and that both that the copyright notice and this
12permission notice and warranty disclaimer appear in supporting
13documentation, and that the name of Lucent or any of its entities
14not be used in advertising or publicity pertaining to
15distribution of the software without specific, written prior
16permission.
17
18LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
25THIS SOFTWARE.
26
27****************************************************************/
28
29/* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
31
32#include "gdtoaimp.h"
33
34char *__g_dfmt (char *buf, double *d, int ndig, size_t bufsize)
35{
36 static FPI fpi0 = { 53, 1-1023-53+1, 2046-1023-53+1, 1, 0, Int_max };
37 char *b, *s, *se;
38 ULong bits[2], *L, sign;
39 int decpt, ex, i, mode;
40#ifdef Honor_FLT_ROUNDS
41#include "gdtoa_fltrnds.h"
42#else
43#define fpi &fpi0
44#endif
45
46 if (ndig < 0)
47 ndig = 0;
48 if ((int) bufsize < ndig + 10)
49 return 0;
50
51 L = (ULong*)d;
52 sign = L[_0] & 0x80000000L;
53 if ((L[_0] & 0x7ff00000) == 0x7ff00000) {
54 /* Infinity or NaN */
55 if (bufsize < 10)
56 return 0;
57 if (L[_0] & 0xfffff || L[_1]) {
58 return strcp(buf, "NaN");
59 }
60 b = buf;
61 if (sign)
62 *b++ = '-';
63 return strcp(b, "Infinity");
64 }
65 if (L[_1] == 0 && (L[_0] ^ sign) == 0 /*d == 0.*/) {
66 b = buf;
67#ifndef IGNORE_ZERO_SIGN
68 if (L[_0] & 0x80000000L)
69 *b++ = '-';
70#endif
71 *b++ = '0';
72 *b = 0;
73 return b;
74 }
75 bits[0] = L[_1];
76 bits[1] = L[_0] & 0xfffff;
77 if ( (ex = (L[_0] >> 20) & 0x7ff) !=0)
78 bits[1] |= 0x100000;
79 else
80 ex = 1;
81 ex -= 0x3ff + 52;
82 mode = 2;
83 if (ndig <= 0)
84 mode = 0;
85 i = STRTOG_Normal;
86 if (sign)
87 i = STRTOG_Normal | STRTOG_Neg;
88 s = __gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
89 return __g__fmt(buf, s, se, decpt, sign, bufsize);
90}
diff --git a/src/lib/evil/gdtoa/g_ffmt.c b/src/lib/evil/gdtoa/g_ffmt.c
deleted file mode 100644
index f3f7c2419c..0000000000
--- a/src/lib/evil/gdtoa/g_ffmt.c
+++ /dev/null
@@ -1,88 +0,0 @@
1/****************************************************************
2
3The author of this software is David M. Gay.
4
5Copyright (C) 1998 by Lucent Technologies
6All Rights Reserved
7
8Permission to use, copy, modify, and distribute this software and
9its documentation for any purpose and without fee is hereby
10granted, provided that the above copyright notice appear in all
11copies and that both that the copyright notice and this
12permission notice and warranty disclaimer appear in supporting
13documentation, and that the name of Lucent or any of its entities
14not be used in advertising or publicity pertaining to
15distribution of the software without specific, written prior
16permission.
17
18LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
25THIS SOFTWARE.
26
27****************************************************************/
28
29/* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
31
32#include "gdtoaimp.h"
33
34char *__g_ffmt (char *buf, float *f, int ndig, size_t bufsize)
35{
36 static FPI fpi0 = { 24, 1-127-24+1, 254-127-24+1, 1, 0, 6 };
37 char *b, *s, *se;
38 ULong bits[1], *L, sign;
39 int decpt, ex, i, mode;
40#ifdef Honor_FLT_ROUNDS
41#include "gdtoa_fltrnds.h"
42#else
43#define fpi &fpi0
44#endif
45
46 if (ndig < 0)
47 ndig = 0;
48 if ((int) bufsize < ndig + 10)
49 return 0;
50
51 L = (ULong*)f;
52 sign = L[0] & 0x80000000L;
53 if ((L[0] & 0x7f800000) == 0x7f800000) {
54 /* Infinity or NaN */
55 if (L[0] & 0x7fffff) {
56 return strcp(buf, "NaN");
57 }
58 b = buf;
59 if (sign)
60 *b++ = '-';
61 return strcp(b, "Infinity");
62 }
63 if (*f == 0.) {
64 b = buf;
65#ifndef IGNORE_ZERO_SIGN
66 if (L[0] & 0x80000000L)
67 *b++ = '-';
68#endif
69 *b++ = '0';
70 *b = 0;
71 return b;
72 }
73 bits[0] = L[0] & 0x7fffff;
74 if ( (ex = (L[0] >> 23) & 0xff) !=0)
75 bits[0] |= 0x800000;
76 else
77 ex = 1;
78 ex -= 0x7f + 23;
79 mode = 2;
80 if (ndig <= 0) {
81 if (bufsize < 16)
82 return 0;
83 mode = 0;
84 }
85 i = STRTOG_Normal;
86 s = __gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
87 return __g__fmt(buf, s, se, decpt, sign, bufsize);
88}
diff --git a/src/lib/evil/gdtoa/g_xfmt.c b/src/lib/evil/gdtoa/g_xfmt.c
deleted file mode 100644
index da11a5881b..0000000000
--- a/src/lib/evil/gdtoa/g_xfmt.c
+++ /dev/null
@@ -1,143 +0,0 @@
1/****************************************************************
2
3The author of this software is David M. Gay.
4
5Copyright (C) 1998 by Lucent Technologies
6All Rights Reserved
7
8Permission to use, copy, modify, and distribute this software and
9its documentation for any purpose and without fee is hereby
10granted, provided that the above copyright notice appear in all
11copies and that both that the copyright notice and this
12permission notice and warranty disclaimer appear in supporting
13documentation, and that the name of Lucent or any of its entities
14not be used in advertising or publicity pertaining to
15distribution of the software without specific, written prior
16permission.
17
18LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
25THIS SOFTWARE.
26
27****************************************************************/
28
29/* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
31
32#include "gdtoaimp.h"
33
34#undef _0
35#undef _1
36
37/* one or the other of IEEE_MC68k or IEEE_8087 should be #defined */
38
39#ifdef IEEE_MC68k
40#define _0 0
41#define _1 1
42#define _2 2
43#define _3 3
44#define _4 4
45#endif
46#ifdef IEEE_8087
47#define _0 4
48#define _1 3
49#define _2 2
50#define _3 1
51#define _4 0
52#endif
53
54char *__g_xfmt (char *buf, void *V, int ndig, size_t bufsize)
55{
56 static FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, 0, Int_max };
57 char *b, *s, *se;
58 ULong bits[2], sign;
59 UShort *L;
60 int decpt, ex, i, mode;
61#if defined(__MINGW32__) || defined(__MINGW64__)
62 int fptype = __fpclassifyl (*(long double*) V);
63#endif /* MinGW */
64#ifdef Honor_FLT_ROUNDS
65#include "gdtoa_fltrnds.h"
66#else
67#define fpi &fpi0
68#endif
69
70 if (ndig < 0)
71 ndig = 0;
72 if ((int) bufsize < ndig + 10)
73 return 0;
74
75 L = (UShort *)V;
76 sign = L[_0] & 0x8000;
77 ex = L[_0] & 0x7fff;
78 bits[1] = (L[_1] << 16) | L[_2];
79 bits[0] = (L[_3] << 16) | L[_4];
80
81#if defined(__MINGW32__) || defined(__MINGW64__)
82 if (fptype & FP_NAN) {
83 /* NaN or Inf */
84 if (fptype & FP_NORMAL) {
85 b = buf;
86 *b++ = sign ? '-': '+';
87 strncpy (b, "Infinity", ndig ? ndig : 8);
88 return (buf + strlen (buf));
89 }
90 strncpy (buf, "NaN", ndig ? ndig : 3);
91 return (buf + strlen (buf));
92 }
93 else if (fptype & FP_NORMAL) {
94 /* Normal or subnormal */
95 if (fptype & FP_ZERO) {
96 i = STRTOG_Denormal;
97 ex = 1;
98 }
99 else
100 i = STRTOG_Normal;
101 }
102#else
103 if (ex != 0) {
104 if (ex == 0x7fff) {
105 /* Infinity or NaN */
106 if (bits[0] | bits[1])
107 b = strcp(buf, "NaN");
108 else {
109 b = buf;
110 if (sign)
111 *b++ = '-';
112 b = strcp(b, "Infinity");
113 }
114 return b;
115 }
116 i = STRTOG_Normal;
117 }
118 else if (bits[0] | bits[1]) {
119 i = STRTOG_Denormal;
120 ex = 1;
121 }
122#endif
123 else {
124 /* i = STRTOG_Zero; */
125 b = buf;
126#ifndef IGNORE_ZERO_SIGN
127 if (sign)
128 *b++ = '-';
129#endif
130 *b++ = '0';
131 *b = 0;
132 return b;
133 }
134 ex -= 0x3fff + 63;
135 mode = 2;
136 if (ndig <= 0) {
137 if (bufsize < 32)
138 return 0;
139 mode = 0;
140 }
141 s = __gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
142 return __g__fmt(buf, s, se, decpt, sign, bufsize);
143}
diff --git a/src/lib/evil/gdtoa/gd_arith.h b/src/lib/evil/gdtoa/gd_arith.h
deleted file mode 100644
index f86e633b37..0000000000
--- a/src/lib/evil/gdtoa/gd_arith.h
+++ /dev/null
@@ -1,6 +0,0 @@
1#define IEEE_8087
2#define Arith_Kind_ASL 1
3#define Double_Align
4#ifdef _WIN64
5#define X64_bit_pointers
6#endif /* w64 */
diff --git a/src/lib/evil/gdtoa/gd_qnan.h b/src/lib/evil/gdtoa/gd_qnan.h
deleted file mode 100644
index 87eba8fb31..0000000000
--- a/src/lib/evil/gdtoa/gd_qnan.h
+++ /dev/null
@@ -1,12 +0,0 @@
1#define f_QNAN 0xffc00000
2#define d_QNAN0 0x0
3#define d_QNAN1 0xfff80000
4#define ld_QNAN0 0x0
5#define ld_QNAN1 0xc0000000
6#define ld_QNAN2 0xffff
7#define ld_QNAN3 0x0
8#define ldus_QNAN0 0x0
9#define ldus_QNAN1 0x0
10#define ldus_QNAN2 0x0
11#define ldus_QNAN3 0xc000
12#define ldus_QNAN4 0xffff
diff --git a/src/lib/evil/gdtoa/gdtoa.c b/src/lib/evil/gdtoa/gdtoa.c
deleted file mode 100644
index cf9c290fb4..0000000000
--- a/src/lib/evil/gdtoa/gdtoa.c
+++ /dev/null
@@ -1,733 +0,0 @@
1/****************************************************************
2
3The author of this software is David M. Gay.
4
5Copyright (C) 1998, 1999 by Lucent Technologies
6All Rights Reserved
7
8Permission to use, copy, modify, and distribute this software and
9its documentation for any purpose and without fee is hereby
10granted, provided that the above copyright notice appear in all
11copies and that both that the copyright notice and this
12permission notice and warranty disclaimer appear in supporting
13documentation, and that the name of Lucent or any of its entities
14not be used in advertising or publicity pertaining to
15distribution of the software without specific, written prior
16permission.
17
18LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
25THIS SOFTWARE.
26
27****************************************************************/
28
29/* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
31
32#include "gdtoaimp.h"
33
34static Bigint *bitstob (ULong *bits, int nbits, int *bbits)
35{
36 int i, k;
37 Bigint *b;
38 ULong *be, *x, *x0;
39
40 i = ULbits;
41 k = 0;
42 while(i < nbits) {
43 i <<= 1;
44 k++;
45 }
46#ifndef Pack_32
47 if (!k)
48 k = 1;
49#endif
50 b = Balloc(k);
51 be = bits + ((nbits - 1) >> kshift);
52 x = x0 = b->x;
53 do {
54 *x++ = *bits & ALL_ON;
55#ifdef Pack_16
56 *x++ = (*bits >> 16) & ALL_ON;
57#endif
58 } while(++bits <= be);
59 i = x - x0;
60 while(!x0[--i])
61 if (!i) {
62 b->wds = 0;
63 *bbits = 0;
64 goto ret;
65 }
66 b->wds = i + 1;
67 *bbits = i*ULbits + 32 - hi0bits(b->x[i]);
68 ret:
69 return b;
70}
71
72/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
73 *
74 * Inspired by "How to Print Floating-Point Numbers Accurately" by
75 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
76 *
77 * Modifications:
78 * 1. Rather than iterating, we use a simple numeric overestimate
79 * to determine k = floor(log10(d)). We scale relevant
80 * quantities using O(log2(k)) rather than O(k) multiplications.
81 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
82 * try to generate digits strictly left to right. Instead, we
83 * compute with fewer bits and propagate the carry if necessary
84 * when rounding the final digit up. This is often faster.
85 * 3. Under the assumption that input will be rounded nearest,
86 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
87 * That is, we allow equality in stopping tests when the
88 * round-nearest rule will give the same floating-point value
89 * as would satisfaction of the stopping test with strict
90 * inequality.
91 * 4. We remove common factors of powers of 2 from relevant
92 * quantities.
93 * 5. When converting floating-point integers less than 1e16,
94 * we use floating-point arithmetic rather than resorting
95 * to multiple-precision integers.
96 * 6. When asked to produce fewer than 15 digits, we first try
97 * to get by with floating-point arithmetic; we resort to
98 * multiple-precision integer arithmetic only if we cannot
99 * guarantee that the floating-point calculation has given
100 * the correctly rounded result. For k requested digits and
101 * "uniformly" distributed input, the probability is
102 * something like 10^(k-15) that we must resort to the Long
103 * calculation.
104 */
105
106char *__gdtoa (FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits,
107 int *decpt, char **rve)
108{
109 /* Arguments ndigits and decpt are similar to the second and third
110 arguments of ecvt and fcvt; trailing zeros are suppressed from
111 the returned string. If not null, *rve is set to point
112 to the end of the return value. If d is +-Infinity or NaN,
113 then *decpt is set to 9999.
114 be = exponent: value = (integer represented by bits) * (2 to the power of be).
115
116 mode:
117 0 ==> shortest string that yields d when read in
118 and rounded to nearest.
119 1 ==> like 0, but with Steele & White stopping rule;
120 e.g. with IEEE P754 arithmetic , mode 0 gives
121 1e23 whereas mode 1 gives 9.999999999999999e22.
122 2 ==> max(1,ndigits) significant digits. This gives a
123 return value similar to that of ecvt, except
124 that trailing zeros are suppressed.
125 3 ==> through ndigits past the decimal point. This
126 gives a return value similar to that from fcvt,
127 except that trailing zeros are suppressed, and
128 ndigits can be negative.
129 4-9 should give the same return values as 2-3, i.e.,
130 4 <= mode <= 9 ==> same return as mode
131 2 + (mode & 1). These modes are mainly for
132 debugging; often they run slower but sometimes
133 faster than modes 2-3.
134 4,5,8,9 ==> left-to-right digit generation.
135 6-9 ==> don't try fast floating-point estimate
136 (if applicable).
137
138 Values of mode other than 0-9 are treated as mode 0.
139
140 Sufficient space is allocated to the return value
141 to hold the suppressed trailing zeros.
142 */
143
144 int bbits, b2, b5, be0, dig, i, ieps, ilim, ilim0, ilim1, inex;
145 int j, j2, k, k0, k_check, kind, leftright, m2, m5, nbits;
146 int rdir, s2, s5, spec_case, try_quick;
147 Long L;
148 Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S;
149 double d2, ds;
150 char *s, *s0;
151 union _dbl_union d, eps;
152
153#ifndef MULTIPLE_THREADS
154 if (dtoa_result) {
155 __freedtoa(dtoa_result);
156 dtoa_result = 0;
157 }
158#endif
159 inex = 0;
160 kind = *kindp &= ~STRTOG_Inexact;
161 switch(kind & STRTOG_Retmask) {
162 case STRTOG_Zero:
163 goto ret_zero;
164 case STRTOG_Normal:
165 case STRTOG_Denormal:
166 break;
167 case STRTOG_Infinite:
168 *decpt = -32768;
169 return nrv_alloc("Infinity", rve, 8);
170 case STRTOG_NaN:
171 *decpt = -32768;
172 return nrv_alloc("NaN", rve, 3);
173 default:
174 return 0;
175 }
176 b = bitstob(bits, nbits = fpi->nbits, &bbits);
177 be0 = be;
178 if ( (i = trailz(b)) !=0) {
179 rshift(b, i);
180 be += i;
181 bbits -= i;
182 }
183 if (!b->wds) {
184 Bfree(b);
185 ret_zero:
186 *decpt = 1;
187 return nrv_alloc("0", rve, 1);
188 }
189
190 dval(&d) = b2d(b, &i);
191 i = be + bbits - 1;
192 word0(&d) &= Frac_mask1;
193 word0(&d) |= Exp_11;
194
195 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
196 * log10(x) = log(x) / log(10)
197 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
198 * log10(&d) = (i-Bias)*log(2)/log(10) + log10(d2)
199 *
200 * This suggests computing an approximation k to log10(&d) by
201 *
202 * k = (i - Bias)*0.301029995663981
203 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
204 *
205 * We want k to be too large rather than too small.
206 * The error in the first-order Taylor series approximation
207 * is in our favor, so we just round up the constant enough
208 * to compensate for any error in the multiplication of
209 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
210 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
211 * adding 1e-13 to the constant term more than suffices.
212 * Hence we adjust the constant term to 0.1760912590558.
213 * (We could get a more accurate k by invoking log10,
214 * but this is probably not worthwhile.)
215 */
216 ds = (dval(&d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
217
218 /* correct assumption about exponent range */
219 if ((j = i) < 0)
220 j = -j;
221 if ((j -= 1077) > 0)
222 ds += j * 7e-17;
223
224 k = (int)ds;
225 if (ds < 0. && ds != k)
226 k--; /* want k = floor(ds) */
227 k_check = 1;
228 word0(&d) += (be + bbits - 1) << Exp_shift;
229 if (k >= 0 && k <= Ten_pmax) {
230 if (dval(&d) < tens[k])
231 k--;
232 k_check = 0;
233 }
234 j = bbits - i - 1;
235 if (j >= 0) {
236 b2 = 0;
237 s2 = j;
238 }
239 else {
240 b2 = -j;
241 s2 = 0;
242 }
243 if (k >= 0) {
244 b5 = 0;
245 s5 = k;
246 s2 += k;
247 }
248 else {
249 b2 -= k;
250 b5 = -k;
251 s5 = 0;
252 }
253 if (mode < 0 || mode > 9)
254 mode = 0;
255 try_quick = 1;
256 if (mode > 5) {
257 mode -= 4;
258 try_quick = 0;
259 }
260 else if (i >= -4 - Emin || i < Emin)
261 try_quick = 0;
262 leftright = 1;
263 ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
264 /* silence erroneous "gcc -Wall" warning. */
265 switch(mode) {
266 case 0:
267 case 1:
268 i = (int)(nbits * .30103) + 3;
269 ndigits = 0;
270 break;
271 case 2:
272 leftright = 0;
273 /* no break */
274 case 4:
275 if (ndigits <= 0)
276 ndigits = 1;
277 ilim = ilim1 = i = ndigits;
278 break;
279 case 3:
280 leftright = 0;
281 /* no break */
282 case 5:
283 i = ndigits + k + 1;
284 ilim = i;
285 ilim1 = i - 1;
286 if (i <= 0)
287 i = 1;
288 }
289 s = s0 = rv_alloc(i);
290
291 if ( (rdir = fpi->rounding - 1) !=0) {
292 if (rdir < 0)
293 rdir = 2;
294 if (kind & STRTOG_Neg)
295 rdir = 3 - rdir;
296 }
297
298 /* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */
299
300 if (ilim >= 0 && ilim <= Quick_max && try_quick && !rdir
301#ifndef IMPRECISE_INEXACT
302 && k == 0
303#endif
304 ) {
305
306 /* Try to get by with floating-point arithmetic. */
307
308 i = 0;
309 d2 = dval(&d);
310 k0 = k;
311 ilim0 = ilim;
312 ieps = 2; /* conservative */
313 if (k > 0) {
314 ds = tens[k&0xf];
315 j = k >> 4;
316 if (j & Bletch) {
317 /* prevent overflows */
318 j &= Bletch - 1;
319 dval(&d) /= bigtens[n_bigtens-1];
320 ieps++;
321 }
322 for(; j; j >>= 1, i++)
323 if (j & 1) {
324 ieps++;
325 ds *= bigtens[i];
326 }
327 }
328 else {
329 ds = 1.;
330 if ( (j2 = -k) !=0) {
331 dval(&d) *= tens[j2 & 0xf];
332 for(j = j2 >> 4; j; j >>= 1, i++)
333 if (j & 1) {
334 ieps++;
335 dval(&d) *= bigtens[i];
336 }
337 }
338 }
339 if (k_check && dval(&d) < 1. && ilim > 0) {
340 if (ilim1 <= 0)
341 goto fast_failed;
342 ilim = ilim1;
343 k--;
344 dval(&d) *= 10.;
345 ieps++;
346 }
347 dval(&eps) = ieps*dval(&d) + 7.;
348 word0(&eps) -= (P-1)*Exp_msk1;
349 if (ilim == 0) {
350 S = mhi = 0;
351 dval(&d) -= 5.;
352 if (dval(&d) > dval(&eps))
353 goto one_digit;
354 if (dval(&d) < -dval(&eps))
355 goto no_digits;
356 goto fast_failed;
357 }
358#ifndef No_leftright
359 if (leftright) {
360 /* Use Steele & White method of only
361 * generating digits needed.
362 */
363 dval(&eps) = ds*0.5/tens[ilim-1] - dval(&eps);
364 for(i = 0;;) {
365 L = (Long)(dval(&d)/ds);
366 dval(&d) -= L*ds;
367 *s++ = '0' + (int)L;
368 if (dval(&d) < dval(&eps)) {
369 if (dval(&d))
370 inex = STRTOG_Inexlo;
371 goto ret1;
372 }
373 if (ds - dval(&d) < dval(&eps))
374 goto bump_up;
375 if (++i >= ilim)
376 break;
377 dval(&eps) *= 10.;
378 dval(&d) *= 10.;
379 }
380 }
381 else {
382#endif
383 /* Generate ilim digits, then fix them up. */
384 dval(&eps) *= tens[ilim-1];
385 for(i = 1;; i++, dval(&d) *= 10.) {
386 if ( (L = (Long)(dval(&d)/ds)) !=0)
387 dval(&d) -= L*ds;
388 *s++ = '0' + (int)L;
389 if (i == ilim) {
390 ds *= 0.5;
391 if (dval(&d) > ds + dval(&eps))
392 goto bump_up;
393 else if (dval(&d) < ds - dval(&eps)) {
394 if (dval(&d))
395 inex = STRTOG_Inexlo;
396 goto clear_trailing0;
397 }
398 break;
399 }
400 }
401#ifndef No_leftright
402 }
403#endif
404 fast_failed:
405 s = s0;
406 dval(&d) = d2;
407 k = k0;
408 ilim = ilim0;
409 }
410
411 /* Do we have a "small" integer? */
412
413 if (be >= 0 && k <= fpi->int_max) {
414 /* Yes. */
415 ds = tens[k];
416 if (ndigits < 0 && ilim <= 0) {
417 S = mhi = 0;
418 if (ilim < 0 || dval(&d) <= 5*ds)
419 goto no_digits;
420 goto one_digit;
421 }
422 for(i = 1;; i++, dval(&d) *= 10.) {
423 L = dval(&d) / ds;
424 dval(&d) -= L*ds;
425#ifdef Check_FLT_ROUNDS
426 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
427 if (dval(&d) < 0) {
428 L--;
429 dval(&d) += ds;
430 }
431#endif
432 *s++ = '0' + (int)L;
433 if (dval(&d) == 0.)
434 break;
435 if (i == ilim) {
436 if (rdir) {
437 if (rdir == 1)
438 goto bump_up;
439 inex = STRTOG_Inexlo;
440 goto ret1;
441 }
442 dval(&d) += dval(&d);
443#ifdef ROUND_BIASED
444 if (dval(&d) >= ds)
445#else
446 if (dval(&d) > ds || (dval(&d) == ds && L & 1))
447#endif
448 {
449 bump_up:
450 inex = STRTOG_Inexhi;
451 while(*--s == '9')
452 if (s == s0) {
453 k++;
454 *s = '0';
455 break;
456 }
457 ++*s++;
458 }
459 else {
460 inex = STRTOG_Inexlo;
461 clear_trailing0:
462 while(*--s == '0'){}
463 ++s;
464 }
465 break;
466 }
467 }
468 goto ret1;
469 }
470
471 m2 = b2;
472 m5 = b5;
473 mhi = mlo = 0;
474 if (leftright) {
475 i = nbits - bbits;
476 if (be - i++ < fpi->emin && mode != 3 && mode != 5) {
477 /* denormal */
478 i = be - fpi->emin + 1;
479 if (mode >= 2 && ilim > 0 && ilim < i)
480 goto small_ilim;
481 }
482 else if (mode >= 2) {
483 small_ilim:
484 j = ilim - 1;
485 if (m5 >= j)
486 m5 -= j;
487 else {
488 s5 += j -= m5;
489 b5 += j;
490 m5 = 0;
491 }
492 if ((i = ilim) < 0) {
493 m2 -= i;
494 i = 0;
495 }
496 }
497 b2 += i;
498 s2 += i;
499 mhi = i2b(1);
500 }
501 if (m2 > 0 && s2 > 0) {
502 i = m2 < s2 ? m2 : s2;
503 b2 -= i;
504 m2 -= i;
505 s2 -= i;
506 }
507 if (b5 > 0) {
508 if (leftright) {
509 if (m5 > 0) {
510 mhi = pow5mult(mhi, m5);
511 b1 = mult(mhi, b);
512 Bfree(b);
513 b = b1;
514 }
515 if ( (j = b5 - m5) !=0)
516 b = pow5mult(b, j);
517 }
518 else
519 b = pow5mult(b, b5);
520 }
521 S = i2b(1);
522 if (s5 > 0)
523 S = pow5mult(S, s5);
524
525 /* Check for special case that d is a normalized power of 2. */
526
527 spec_case = 0;
528 if (mode < 2) {
529 if (bbits == 1 && be0 > fpi->emin + 1) {
530 /* The special case */
531 b2++;
532 s2++;
533 spec_case = 1;
534 }
535 }
536
537 /* Arrange for convenient computation of quotients:
538 * shift left if necessary so divisor has 4 leading 0 bits.
539 *
540 * Perhaps we should just compute leading 28 bits of S once
541 * and for all and pass them and a shift to quorem, so it
542 * can do shifts and ors to compute the numerator for q.
543 */
544 i = ((s5 ? hi0bits(S->x[S->wds-1]) : ULbits - 1) - s2 - 4) & kmask;
545 m2 += i;
546 if ((b2 += i) > 0)
547 b = lshift(b, b2);
548 if ((s2 += i) > 0)
549 S = lshift(S, s2);
550 if (k_check) {
551 if (cmp(b,S) < 0) {
552 k--;
553 b = multadd(b, 10, 0); /* we botched the k estimate */
554 if (leftright)
555 mhi = multadd(mhi, 10, 0);
556 ilim = ilim1;
557 }
558 }
559 if (ilim <= 0 && mode > 2) {
560 if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
561 /* no digits, fcvt style */
562 no_digits:
563 k = -1 - ndigits;
564 inex = STRTOG_Inexlo;
565 goto ret;
566 }
567 one_digit:
568 inex = STRTOG_Inexhi;
569 *s++ = '1';
570 k++;
571 goto ret;
572 }
573 if (leftright) {
574 if (m2 > 0)
575 mhi = lshift(mhi, m2);
576
577 /* Compute mlo -- check for special case
578 * that d is a normalized power of 2.
579 */
580
581 mlo = mhi;
582 if (spec_case) {
583 mhi = Balloc(mhi->k);
584 Bcopy(mhi, mlo);
585 mhi = lshift(mhi, 1);
586 }
587
588 for(i = 1;;i++) {
589 dig = quorem(b,S) + '0';
590 /* Do we yet have the shortest decimal string
591 * that will round to d?
592 */
593 j = cmp(b, mlo);
594 delta = diff(S, mhi);
595 j2 = delta->sign ? 1 : cmp(b, delta);
596 Bfree(delta);
597#ifndef ROUND_BIASED
598 if (j2 == 0 && !mode && !(bits[0] & 1) && !rdir) {
599 if (dig == '9')
600 goto round_9_up;
601 if (j <= 0) {
602 if (b->wds > 1 || b->x[0])
603 inex = STRTOG_Inexlo;
604 }
605 else {
606 dig++;
607 inex = STRTOG_Inexhi;
608 }
609 *s++ = dig;
610 goto ret;
611 }
612#endif
613 if (j < 0 || (j == 0 && !mode
614#ifndef ROUND_BIASED
615 && !(bits[0] & 1)
616#endif
617 )) {
618 if (rdir && (b->wds > 1 || b->x[0])) {
619 if (rdir == 2) {
620 inex = STRTOG_Inexlo;
621 goto accept;
622 }
623 while (cmp(S,mhi) > 0) {
624 *s++ = dig;
625 mhi1 = multadd(mhi, 10, 0);
626 if (mlo == mhi)
627 mlo = mhi1;
628 mhi = mhi1;
629 b = multadd(b, 10, 0);
630 dig = quorem(b,S) + '0';
631 }
632 if (dig++ == '9')
633 goto round_9_up;
634 inex = STRTOG_Inexhi;
635 goto accept;
636 }
637 if (j2 > 0) {
638 b = lshift(b, 1);
639 j2 = cmp(b, S);
640#ifdef ROUND_BIASED
641 if (j2 >= 0 /*)*/
642#else
643 if ((j2 > 0 || (j2 == 0 && dig & 1))
644#endif
645 && dig++ == '9')
646 goto round_9_up;
647 inex = STRTOG_Inexhi;
648 }
649 if (b->wds > 1 || b->x[0])
650 inex = STRTOG_Inexlo;
651 accept:
652 *s++ = dig;
653 goto ret;
654 }
655 if (j2 > 0 && rdir != 2) {
656 if (dig == '9') { /* possible if i == 1 */
657 round_9_up:
658 *s++ = '9';
659 inex = STRTOG_Inexhi;
660 goto roundoff;
661 }
662 inex = STRTOG_Inexhi;
663 *s++ = dig + 1;
664 goto ret;
665 }
666 *s++ = dig;
667 if (i == ilim)
668 break;
669 b = multadd(b, 10, 0);
670 if (mlo == mhi)
671 mlo = mhi = multadd(mhi, 10, 0);
672 else {
673 mlo = multadd(mlo, 10, 0);
674 mhi = multadd(mhi, 10, 0);
675 }
676 }
677 }
678 else
679 for(i = 1;; i++) {
680 *s++ = dig = quorem(b,S) + '0';
681 if (i >= ilim)
682 break;
683 b = multadd(b, 10, 0);
684 }
685
686 /* Round off last digit */
687
688 if (rdir) {
689 if (rdir == 2 || (b->wds <= 1 && !b->x[0]))
690 goto chopzeros;
691 goto roundoff;
692 }
693 b = lshift(b, 1);
694 j = cmp(b, S);
695#ifdef ROUND_BIASED
696 if (j >= 0)
697#else
698 if (j > 0 || (j == 0 && dig & 1))
699#endif
700 {
701 roundoff:
702 inex = STRTOG_Inexhi;
703 while(*--s == '9')
704 if (s == s0) {
705 k++;
706 *s++ = '1';
707 goto ret;
708 }
709 ++*s++;
710 }
711 else {
712 chopzeros:
713 if (b->wds > 1 || b->x[0])
714 inex = STRTOG_Inexlo;
715 while(*--s == '0'){}
716 ++s;
717 }
718 ret:
719 Bfree(S);
720 if (mhi) {
721 if (mlo && mlo != mhi)
722 Bfree(mlo);
723 Bfree(mhi);
724 }
725 ret1:
726 Bfree(b);
727 *s = 0;
728 *decpt = k + 1;
729 if (rve)
730 *rve = s;
731 *kindp |= inex;
732 return s0;
733}
diff --git a/src/lib/evil/gdtoa/gdtoa.h b/src/lib/evil/gdtoa/gdtoa.h
deleted file mode 100644
index ef72bf99ed..0000000000
--- a/src/lib/evil/gdtoa/gdtoa.h
+++ /dev/null
@@ -1,121 +0,0 @@
1/****************************************************************
2
3The author of this software is David M. Gay.
4
5Copyright (C) 1998 by Lucent Technologies
6All Rights Reserved
7
8Permission to use, copy, modify, and distribute this software and
9its documentation for any purpose and without fee is hereby
10granted, provided that the above copyright notice appear in all
11copies and that both that the copyright notice and this
12permission notice and warranty disclaimer appear in supporting
13documentation, and that the name of Lucent or any of its entities
14not be used in advertising or publicity pertaining to
15distribution of the software without specific, written prior
16permission.
17
18LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
25THIS SOFTWARE.
26
27****************************************************************/
28
29/* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
31
32/* Modified by Danny Smith for inclusion in libmingwex.a
33 Aug 2006 */
34
35/* Modified by Vincent Torri for inclusion in Evil
36 February 2012 */
37
38#ifndef GDTOA_H_INCLUDED
39#define GDTOA_H_INCLUDED
40
41#include "gd_arith.h"
42#include <stddef.h> /* for size_t */
43
44#if defined(__MINGW32__) || defined(__MINGW64__)
45/* keep the 'Long' definition as 'long' for compatibility
46 * with older/other software. long in w64 is 32 bits anyway..
47 */
48#define Long long /* Windows long is 32 bit */
49#undef NO_LONG_LONG /* we have long long type */
50#endif /* MinGW */
51
52#ifndef Long
53#define Long int
54#endif
55#ifndef ULong
56typedef unsigned Long ULong;
57#endif
58#ifndef UShort
59typedef unsigned short UShort;
60#endif
61
62enum { /* return values from strtodg */
63 STRTOG_Zero = 0,
64 STRTOG_Normal = 1,
65 STRTOG_Denormal = 2,
66 STRTOG_Infinite = 3,
67 STRTOG_NaN = 4,
68 STRTOG_NaNbits = 5,
69 STRTOG_NoNumber = 6,
70 STRTOG_Retmask = 7,
71
72 /* The following may be or-ed into one of the above values. */
73
74 STRTOG_Neg = 0x08, /* does not affect STRTOG_Inexlo or STRTOG_Inexhi */
75 STRTOG_Inexlo = 0x10, /* returned result rounded toward zero */
76 STRTOG_Inexhi = 0x20, /* returned result rounded away from zero */
77 STRTOG_Inexact = 0x30,
78 STRTOG_Underflow= 0x40,
79 STRTOG_Overflow = 0x80
80};
81
82typedef struct
83FPI {
84 int nbits;
85 int emin;
86 int emax;
87 int rounding;
88 int sudden_underflow;
89 int int_max;
90} FPI;
91
92enum { /* FPI.rounding values: same as FLT_ROUNDS */
93 FPI_Round_zero = 0,
94 FPI_Round_near = 1,
95 FPI_Round_up = 2,
96 FPI_Round_down = 3
97};
98
99#ifdef __cplusplus
100extern "C" {
101#endif
102
103extern char* __dtoa (double d, int mode, int ndigits, int *decpt,
104 int *sign, char **rve);
105extern char* __gdtoa (FPI *fpi, int be, ULong *bits, int *kindp,
106 int mode, int ndigits, int *decpt, char **rve);
107extern void __freedtoa (char *);
108
109extern float __evil_strtof (const char *, char **);
110extern long double __evil_strtold (const char *, char **);
111extern int __strtodg (const char *, char **, FPI *, Long *, ULong *);
112
113extern char* __g__fmt (char*, char*, char*, int, ULong, size_t);
114extern char* __g_dfmt (char*, double*, int, size_t);
115extern char* __g_ffmt (char*, float*, int, size_t);
116extern char* __g_xfmt (char*, void*, int, size_t);
117
118#ifdef __cplusplus
119}
120#endif
121#endif /* GDTOA_H_INCLUDED */
diff --git a/src/lib/evil/gdtoa/gdtoa_fltrnds.h b/src/lib/evil/gdtoa/gdtoa_fltrnds.h
deleted file mode 100644
index 28c474e2c2..0000000000
--- a/src/lib/evil/gdtoa/gdtoa_fltrnds.h
+++ /dev/null
@@ -1,18 +0,0 @@
1 FPI *fpi, fpi1;
2 int Rounding;
3#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
4 Rounding = Flt_Rounds;
5#else /*}{*/
6 Rounding = 1;
7 switch(fegetround()) {
8 case FE_TOWARDZERO: Rounding = 0; break;
9 case FE_UPWARD: Rounding = 2; break;
10 case FE_DOWNWARD: Rounding = 3;
11 }
12#endif /*}}*/
13 fpi = &fpi0;
14 if (Rounding != 1) {
15 fpi1 = fpi0;
16 fpi = &fpi1;
17 fpi1.rounding = Rounding;
18 }
diff --git a/src/lib/evil/gdtoa/gdtoaimp.h b/src/lib/evil/gdtoa/gdtoaimp.h
deleted file mode 100644
index 3a06c2072f..0000000000
--- a/src/lib/evil/gdtoa/gdtoaimp.h
+++ /dev/null
@@ -1,645 +0,0 @@
1/****************************************************************
2
3The author of this software is David M. Gay.
4
5Copyright (C) 1998-2000 by Lucent Technologies
6All Rights Reserved
7
8Permission to use, copy, modify, and distribute this software and
9its documentation for any purpose and without fee is hereby
10granted, provided that the above copyright notice appear in all
11copies and that both that the copyright notice and this
12permission notice and warranty disclaimer appear in supporting
13documentation, and that the name of Lucent or any of its entities
14not be used in advertising or publicity pertaining to
15distribution of the software without specific, written prior
16permission.
17
18LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
25THIS SOFTWARE.
26
27****************************************************************/
28
29/* This is a variation on dtoa.c that converts arbitary binary
30 floating-point formats to and from decimal notation. It uses
31 double-precision arithmetic internally, so there are still
32 various #ifdefs that adapt the calculations to the native
33 double-precision arithmetic (any of IEEE, VAX D_floating,
34 or IBM mainframe arithmetic).
35
36 Please send bug reports to David M. Gay (dmg at acm dot org,
37 with " at " changed at "@" and " dot " changed to ".").
38 */
39
40/* On a machine with IEEE extended-precision registers, it is
41 * necessary to specify double-precision (53-bit) rounding precision
42 * before invoking strtod or dtoa. If the machine uses (the equivalent
43 * of) Intel 80x87 arithmetic, the call
44 * _control87(PC_53, MCW_PC);
45 * does this with many compilers. Whether this or another call is
46 * appropriate depends on the compiler; for this to work, it may be
47 * necessary to #include "float.h" or another system-dependent header
48 * file.
49 */
50
51/* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
52 *
53 * This strtod returns a nearest machine number to the input decimal
54 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
55 * broken by the IEEE round-even rule. Otherwise ties are broken by
56 * biased rounding (add half and chop).
57 *
58 * Inspired loosely by William D. Clinger's paper "How to Read Floating
59 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 112-126].
60 *
61 * Modifications:
62 *
63 * 1. We only require IEEE, IBM, or VAX double-precision
64 * arithmetic (not IEEE double-extended).
65 * 2. We get by with floating-point arithmetic in a case that
66 * Clinger missed -- when we're computing d * 10^n
67 * for a small integer d and the integer n is not too
68 * much larger than 22 (the maximum integer k for which
69 * we can represent 10^k exactly), we may be able to
70 * compute (d*10^k) * 10^(e-k) with just one roundoff.
71 * 3. Rather than a bit-at-a-time adjustment of the binary
72 * result in the hard case, we use floating-point
73 * arithmetic to determine the adjustment to within
74 * one bit; only in really hard cases do we need to
75 * compute a second residual.
76 * 4. Because of 3., we don't need a large table of powers of 10
77 * for ten-to-e (just some small tables, e.g. of 10^k
78 * for 0 <= k <= 22).
79 */
80
81/*
82 * #define IEEE_8087 for IEEE-arithmetic machines where the least
83 * significant byte has the lowest address.
84 * #define IEEE_MC68k for IEEE-arithmetic machines where the most
85 * significant byte has the lowest address.
86 * #define Long int on machines with 32-bit ints and 64-bit longs.
87 * #define Sudden_Underflow for IEEE-format machines without gradual
88 * underflow (i.e., that flush to zero on underflow).
89 * #define IBM for IBM mainframe-style floating-point arithmetic.
90 * #define VAX for VAX-style floating-point arithmetic (D_floating).
91 * #define No_leftright to omit left-right logic in fast floating-point
92 * computation of dtoa and gdtoa. This will cause modes 4 and 5 to be
93 * treated the same as modes 2 and 3 for some inputs.
94 * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
95 * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
96 * that use extended-precision instructions to compute rounded
97 * products and quotients) with IBM.
98 * #define ROUND_BIASED for IEEE-format with biased rounding and arithmetic
99 * that rounds toward +Infinity.
100 * #define ROUND_BIASED_without_Round_Up for IEEE-format with biased
101 * rounding when the underlying floating-point arithmetic uses
102 * unbiased rounding. This prevent using ordinary floating-point
103 * arithmetic when the result could be computed with one rounding error.
104 * #define Inaccurate_Divide for IEEE-format with correctly rounded
105 * products but inaccurate quotients, e.g., for Intel i860.
106 * #define NO_LONG_LONG on machines that do not have a "long long"
107 * integer type (of >= 64 bits). On such machines, you can
108 * #define Just_16 to store 16 bits per 32-bit Long when doing
109 * high-precision integer arithmetic. Whether this speeds things
110 * up or slows things down depends on the machine and the number
111 * being converted. If long long is available and the name is
112 * something other than "long long", #define Llong to be the name,
113 * and if "unsigned Llong" does not work as an unsigned version of
114 * Llong, #define #ULLong to be the corresponding unsigned type.
115 * #define KR_headers for old-style C function headers.
116 * #define Bad_float_h if your system lacks a float.h or if it does not
117 * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
118 * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
119 * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
120 * if memory is available and otherwise does something you deem
121 * appropriate. If MALLOC is undefined, malloc will be invoked
122 * directly -- and assumed always to succeed.
123 * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
124 * memory allocations from a private pool of memory when possible.
125 * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
126 * unless #defined to be a different length. This default length
127 * suffices to get rid of MALLOC calls except for unusual cases,
128 * such as decimal-to-binary conversion of a very long string of
129 * digits. When converting IEEE double precision values, the
130 * longest string gdtoa can return is about 751 bytes long. For
131 * conversions by strtod of strings of 800 digits and all gdtoa
132 * conversions of IEEE doubles in single-threaded executions with
133 * 8-byte pointers, PRIVATE_MEM >= 7400 appears to suffice; with
134 * 4-byte pointers, PRIVATE_MEM >= 7112 appears adequate.
135 * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
136 * #defined automatically on IEEE systems. On such systems,
137 * when INFNAN_CHECK is #defined, strtod checks
138 * for Infinity and NaN (case insensitively).
139 * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
140 * strtodg also accepts (case insensitively) strings of the form
141 * NaN(x), where x is a string of hexadecimal digits (optionally
142 * preceded by 0x or 0X) and spaces; if there is only one string
143 * of hexadecimal digits, it is taken for the fraction bits of the
144 * resulting NaN; if there are two or more strings of hexadecimal
145 * digits, each string is assigned to the next available sequence
146 * of 32-bit words of fractions bits (starting with the most
147 * significant), right-aligned in each sequence.
148 * Unless GDTOA_NON_PEDANTIC_NANCHECK is #defined, input "NaN(...)"
149 * is consumed even when ... has the wrong form (in which case the
150 * "(...)" is consumed but ignored).
151 * #define MULTIPLE_THREADS if the system offers preemptively scheduled
152 * multiple threads. In this case, you must provide (or suitably
153 * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
154 * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
155 * in pow5mult, ensures lazy evaluation of only one copy of high
156 * powers of 5; omitting this lock would introduce a small
157 * probability of wasting memory, but would otherwise be harmless.)
158 * You must also invoke freedtoa(s) to free the value s returned by
159 * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
160 * #define IMPRECISE_INEXACT if you do not care about the setting of
161 * the STRTOG_Inexact bits in the special case of doing IEEE double
162 * precision conversions (which could also be done by the strtog in
163 * dtoa.c).
164 * #define NO_HEX_FP to disable recognition of C9x's hexadecimal
165 * floating-point constants.
166 * #define -DNO_ERRNO to suppress setting errno (in strtod.c and
167 * strtodg.c).
168 * #define NO_STRING_H to use private versions of memcpy.
169 * On some K&R systems, it may also be necessary to
170 * #define DECLARE_SIZE_T in this case.
171 * #define USE_LOCALE to use the current locale's decimal_point value.
172 */
173
174#ifndef GDTOAIMP_H_INCLUDED
175#define GDTOAIMP_H_INCLUDED
176#include "gdtoa.h"
177#include "gd_qnan.h"
178
179#if defined(__MINGW32__) || defined(__MINGW64__)
180#define MULTIPLE_THREADS 1
181#define USE_LOCALE 1
182#define NO_LOCALE_CACHE 1
183#endif /* MinGW */
184
185#ifdef Honor_FLT_ROUNDS
186#include <fenv.h>
187#endif
188
189#ifdef DEBUG
190#include <stdio.h>
191#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
192#endif
193
194#include <stdlib.h>
195#include <string.h>
196
197#ifdef MALLOC
198extern void *MALLOC (size_t);
199#else
200#define MALLOC malloc
201#endif
202
203#undef IEEE_Arith
204#undef Avoid_Underflow
205#ifdef IEEE_MC68k
206#define IEEE_Arith
207#endif
208#ifdef IEEE_8087
209#define IEEE_Arith
210#endif
211
212#include <errno.h>
213
214#ifdef NO_ERRNO
215#define SET_ERRNO(x)
216#else
217#define SET_ERRNO(x) \
218 errno = (x)
219#endif
220
221#ifdef Bad_float_h
222
223#ifdef IEEE_Arith
224#define DBL_DIG 15
225#define DBL_MAX_10_EXP 308
226#define DBL_MAX_EXP 1024
227#define FLT_RADIX 2
228#define DBL_MAX 1.7976931348623157e+308
229#endif
230
231#ifdef IBM
232#define DBL_DIG 16
233#define DBL_MAX_10_EXP 75
234#define DBL_MAX_EXP 63
235#define FLT_RADIX 16
236#define DBL_MAX 7.2370055773322621e+75
237#endif
238
239#ifdef VAX
240#define DBL_DIG 16
241#define DBL_MAX_10_EXP 38
242#define DBL_MAX_EXP 127
243#define FLT_RADIX 2
244#define DBL_MAX 1.7014118346046923e+38
245#define n_bigtens 2
246#endif
247
248#ifndef LONG_MAX
249#define LONG_MAX 2147483647
250#endif
251
252#else /* ifndef Bad_float_h */
253#include <float.h>
254#endif /* Bad_float_h */
255
256#ifdef IEEE_Arith
257#define Scale_Bit 0x10
258#define n_bigtens 5
259#endif
260
261#ifdef IBM
262#define n_bigtens 3
263#endif
264
265#ifdef VAX
266#define n_bigtens 2
267#endif
268
269#ifndef __MATH_H__
270#include <math.h>
271#endif
272
273#ifdef __cplusplus
274extern "C" {
275#endif
276
277#if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
278Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
279#endif
280
281typedef union _dbl_union { double d; ULong L[2]; } dbl_union;
282
283#ifdef IEEE_8087
284#define word0(x) (x)->L[1]
285#define word1(x) (x)->L[0]
286#else
287#define word0(x) (x)->L[0]
288#define word1(x) (x)->L[1]
289#endif
290#define dval(x) (x)->d
291
292/* The following definition of Storeinc is appropriate for MIPS processors.
293 * An alternative that might be better on some machines is
294 * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
295 */
296#if defined(IEEE_8087) + defined(VAX)
297#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
298((unsigned short *)a)[0] = (unsigned short)c, a++)
299#else
300#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
301((unsigned short *)a)[1] = (unsigned short)c, a++)
302#endif
303
304/* #define P DBL_MANT_DIG */
305/* Ten_pmax = floor(P*log(2)/log(5)) */
306/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
307/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
308/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
309
310#ifdef IEEE_Arith
311#define Exp_shift 20
312#define Exp_shift1 20
313#define Exp_msk1 0x100000
314#define Exp_msk11 0x100000
315#define Exp_mask 0x7ff00000
316#define P 53
317#define Bias 1023
318#define Emin (-1022)
319#define Exp_1 0x3ff00000
320#define Exp_11 0x3ff00000
321#define Ebits 11
322#define Frac_mask 0xfffff
323#define Frac_mask1 0xfffff
324#define Ten_pmax 22
325#define Bletch 0x10
326#define Bndry_mask 0xfffff
327#define Bndry_mask1 0xfffff
328#define LSB 1
329#define Sign_bit 0x80000000
330#define Log2P 1
331#define Tiny0 0
332#define Tiny1 1
333#define Quick_max 14
334#define Int_max 14
335
336#ifndef Flt_Rounds
337#ifdef FLT_ROUNDS
338#define Flt_Rounds FLT_ROUNDS
339#else
340#define Flt_Rounds 1
341#endif
342#endif /*Flt_Rounds*/
343
344#else /* ifndef IEEE_Arith */
345#undef Sudden_Underflow
346#define Sudden_Underflow
347#ifdef IBM
348#undef Flt_Rounds
349#define Flt_Rounds 0
350#define Exp_shift 24
351#define Exp_shift1 24
352#define Exp_msk1 0x1000000
353#define Exp_msk11 0x1000000
354#define Exp_mask 0x7f000000
355#define P 14
356#define Bias 65
357#define Exp_1 0x41000000
358#define Exp_11 0x41000000
359#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
360#define Frac_mask 0xffffff
361#define Frac_mask1 0xffffff
362#define Bletch 4
363#define Ten_pmax 22
364#define Bndry_mask 0xefffff
365#define Bndry_mask1 0xffffff
366#define LSB 1
367#define Sign_bit 0x80000000
368#define Log2P 4
369#define Tiny0 0x100000
370#define Tiny1 0
371#define Quick_max 14
372#define Int_max 15
373#else /* VAX */
374#undef Flt_Rounds
375#define Flt_Rounds 1
376#define Exp_shift 23
377#define Exp_shift1 7
378#define Exp_msk1 0x80
379#define Exp_msk11 0x800000
380#define Exp_mask 0x7f80
381#define P 56
382#define Bias 129
383#define Exp_1 0x40800000
384#define Exp_11 0x4080
385#define Ebits 8
386#define Frac_mask 0x7fffff
387#define Frac_mask1 0xffff007f
388#define Ten_pmax 24
389#define Bletch 2
390#define Bndry_mask 0xffff007f
391#define Bndry_mask1 0xffff007f
392#define LSB 0x10000
393#define Sign_bit 0x8000
394#define Log2P 1
395#define Tiny0 0x80
396#define Tiny1 0
397#define Quick_max 15
398#define Int_max 15
399#endif /* IBM, VAX */
400#endif /* IEEE_Arith */
401
402#ifndef IEEE_Arith
403#define ROUND_BIASED
404#else
405#ifdef ROUND_BIASED_without_Round_Up
406#undef ROUND_BIASED
407#define ROUND_BIASED
408#endif
409#endif
410
411#ifdef RND_PRODQUOT
412#define rounded_product(a,b) a = rnd_prod(a, b)
413#define rounded_quotient(a,b) a = rnd_quot(a, b)
414extern double rnd_prod(double, double), rnd_quot(double, double);
415#else
416#define rounded_product(a,b) a *= b
417#define rounded_quotient(a,b) a /= b
418#endif
419
420#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
421#define Big1 0xffffffff
422
423#undef Pack_16
424#ifndef Pack_32
425#define Pack_32
426#endif
427
428#ifdef NO_LONG_LONG
429#undef ULLong
430#ifdef Just_16
431#undef Pack_32
432#define Pack_16
433/* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
434 * This makes some inner loops simpler and sometimes saves work
435 * during multiplications, but it often seems to make things slightly
436 * slower. Hence the default is now to store 32 bits per Long.
437 */
438#endif
439#else /* long long available */
440#ifndef Llong
441#define Llong long long
442#endif
443#ifndef ULLong
444#define ULLong unsigned Llong
445#endif
446#endif /* NO_LONG_LONG */
447
448#ifdef Pack_32
449#define ULbits 32
450#define kshift 5
451#define kmask 31
452#define ALL_ON 0xffffffff
453#else
454#define ULbits 16
455#define kshift 4
456#define kmask 15
457#define ALL_ON 0xffff
458#endif
459
460#ifndef MULTIPLE_THREADS
461#define ACQUIRE_DTOA_LOCK(n) /*nothing*/
462#define FREE_DTOA_LOCK(n) /*nothing*/
463#endif
464
465#define Kmax 9
466
467#define Bigint __Bigint
468struct
469Bigint {
470 struct Bigint *next;
471 int k, maxwds, sign, wds;
472 ULong x[1];
473};
474typedef struct Bigint Bigint;
475
476#ifdef NO_STRING_H
477#ifdef DECLARE_SIZE_T
478typedef unsigned int size_t;
479#endif
480extern void memcpy_D2A (void*, const void*, size_t);
481#define Bcopy(x,y) memcpy_D2A(&x->sign,&y->sign,y->wds*sizeof(ULong) + 2*sizeof(int))
482#else /* !NO_STRING_H */
483#define Bcopy(x,y) memcpy(&x->sign,&y->sign,y->wds*sizeof(ULong) + 2*sizeof(int))
484#endif /* NO_STRING_H */
485
486#ifdef __GNUC__
487static inline int
488__lo0bits_D2A (ULong *y)
489{
490 int ret = __builtin_ctz(*y);
491 *y = *y >> ret;
492 return ret;
493}
494
495static inline int
496__hi0bits_D2A (ULong y)
497{
498 return __builtin_clz(y);
499}
500#endif
501
502#define Balloc __Balloc_D2A
503#define Bfree __Bfree_D2A
504#define ULtoQ __ULtoQ_D2A
505#define ULtof __ULtof_D2A
506#define ULtod __ULtod_D2A
507#define ULtodd __ULtodd_D2A
508#define ULtox __ULtox_D2A
509#define ULtoxL __ULtoxL_D2A
510#define any_on __any_on_D2A
511#define b2d __b2d_D2A
512#define bigtens __bigtens_D2A
513#define cmp __cmp_D2A
514#define copybits __copybits_D2A
515#define d2b __d2b_D2A
516#define decrement __decrement_D2A
517#define diff __diff_D2A
518#define dtoa_result __dtoa_result_D2A
519#define gethex __gethex_D2A
520#define hexdig __hexdig_D2A
521#define hexnan __hexnan_D2A
522#define hi0bits_D2A __hi0bits_D2A
523#define hi0bits(x) __hi0bits_D2A((ULong)(x))
524#define i2b __i2b_D2A
525#define increment __increment_D2A
526#define lo0bits __lo0bits_D2A
527#define lshift __lshift_D2A
528#define match __match_D2A
529#define mult __mult_D2A
530#define multadd __multadd_D2A
531#define nrv_alloc __nrv_alloc_D2A
532#define pow5mult __pow5mult_D2A
533#define quorem __quorem_D2A
534#define ratio __ratio_D2A
535#define rshift __rshift_D2A
536#define rv_alloc __rv_alloc_D2A
537#define s2b __s2b_D2A
538#define set_ones __set_ones_D2A
539#define strcp_D2A __strcp_D2A
540#define strcp __strcp_D2A
541#define strtoIg __strtoIg_D2A
542#define sum __sum_D2A
543#define tens __tens_D2A
544#define tinytens __tinytens_D2A
545#define tinytens __tinytens_D2A
546#define trailz __trailz_D2A
547#define ulp __ulp_D2A
548
549#define hexdig_init_D2A __mingw_hexdig_init_D2A
550
551extern char *dtoa_result;
552extern const double bigtens[], tens[], tinytens[];
553extern unsigned char hexdig[];
554
555extern Bigint *Balloc (int);
556extern void Bfree (Bigint*);
557extern void ULtof (ULong*, ULong*, Long, int);
558extern void ULtod (ULong*, ULong*, Long, int);
559extern void ULtodd (ULong*, ULong*, Long, int);
560extern void ULtoQ (ULong*, ULong*, Long, int);
561extern void ULtox (UShort*, ULong*, Long, int);
562extern void ULtoxL (ULong*, ULong*, Long, int);
563extern ULong any_on (Bigint*, int);
564extern double b2d (Bigint*, int*);
565extern int cmp (Bigint*, Bigint*);
566extern void copybits (ULong*, int, Bigint*);
567extern Bigint *d2b (double, int*, int*);
568extern void decrement (Bigint*);
569extern Bigint *diff (Bigint*, Bigint*);
570extern int gethex (const char**, FPI*, Long*, Bigint**, int);
571extern void hexdig_init_D2A(void);
572extern int hexnan (const char**, FPI*, ULong*);
573#ifndef __GNUC__
574extern int hi0bits_D2A (ULong);
575#endif
576extern Bigint *i2b (int);
577extern Bigint *increment (Bigint*);
578#ifndef __GNUC__
579extern int lo0bits (ULong*);
580#endif
581extern Bigint *lshift (Bigint*, int);
582extern int match (const char**, char*);
583extern Bigint *mult (Bigint*, Bigint*);
584extern Bigint *multadd (Bigint*, int, int);
585extern char *nrv_alloc (char*, char **, int);
586extern Bigint *pow5mult (Bigint*, int);
587extern int quorem (Bigint*, Bigint*);
588extern double ratio (Bigint*, Bigint*);
589extern void rshift (Bigint*, int);
590extern char *rv_alloc (int);
591extern Bigint *s2b (const char*, int, int, ULong, int);
592extern Bigint *set_ones (Bigint*, int);
593extern char *strcp (char*, const char*);
594extern Bigint *sum (Bigint*, Bigint*);
595extern int trailz (Bigint*);
596extern double ulp (dbl_union *);
597
598#ifdef __cplusplus
599}
600#endif
601/*
602 * NAN_WORD0 and NAN_WORD1 are only referenced in strtod.c. Prior to
603 * 20050115, they used to be hard-wired here (to 0x7ff80000 and 0,
604 * respectively), but now are determined by compiling and running
605 * qnan.c to generate gd_qnan.h, which specifies d_QNAN0 and d_QNAN1.
606 * Formerly gdtoaimp.h recommended supplying suitable -DNAN_WORD0=...
607 * and -DNAN_WORD1=... values if necessary. This should still work.
608 * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
609 */
610#ifdef IEEE_Arith
611#ifndef NO_INFNAN_CHECK
612#undef INFNAN_CHECK
613#define INFNAN_CHECK
614#endif
615#ifdef IEEE_MC68k
616#define _0 0
617#define _1 1
618#ifndef NAN_WORD0
619#define NAN_WORD0 d_QNAN0
620#endif
621#ifndef NAN_WORD1
622#define NAN_WORD1 d_QNAN1
623#endif
624#else
625#define _0 1
626#define _1 0
627#ifndef NAN_WORD0
628#define NAN_WORD0 d_QNAN1
629#endif
630#ifndef NAN_WORD1
631#define NAN_WORD1 d_QNAN0
632#endif
633#endif
634#else
635#undef INFNAN_CHECK
636#endif
637
638#undef SI
639#ifdef Sudden_Underflow
640#define SI 1
641#else
642#define SI 0
643#endif
644
645#endif /* GDTOAIMP_H_INCLUDED */
diff --git a/src/lib/evil/gdtoa/gethex.c b/src/lib/evil/gdtoa/gethex.c
deleted file mode 100644
index 276175453e..0000000000
--- a/src/lib/evil/gdtoa/gethex.c
+++ /dev/null
@@ -1,340 +0,0 @@
1/****************************************************************
2
3The author of this software is David M. Gay.
4
5Copyright (C) 1998 by Lucent Technologies
6All Rights Reserved
7
8Permission to use, copy, modify, and distribute this software and
9its documentation for any purpose and without fee is hereby
10granted, provided that the above copyright notice appear in all
11copies and that both that the copyright notice and this
12permission notice and warranty disclaimer appear in supporting
13documentation, and that the name of Lucent or any of its entities
14not be used in advertising or publicity pertaining to
15distribution of the software without specific, written prior
16permission.
17
18LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
25THIS SOFTWARE.
26
27****************************************************************/
28
29/* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
31
32#include "gdtoaimp.h"
33
34#ifdef USE_LOCALE
35#include "locale.h"
36#endif
37
38int gethex (const char **sp, FPI *fpi, Long *expo, Bigint **bp, int sign)
39{
40 Bigint *b;
41 const unsigned char *decpt, *s0, *s, *s1;
42 int big, esign, havedig, irv, j, k, n, n0, nbits, up, zret;
43 ULong L, lostbits, *x;
44 Long e, e1;
45#ifdef USE_LOCALE
46 int i;
47 const unsigned char *decimalpoint;
48#ifdef NO_LOCALE_CACHE
49 decimalpoint = (unsigned char *)localeconv()->decimal_point;
50#else
51 static unsigned char *decimalpoint_cache;
52 if (!(s0 = decimalpoint_cache)) {
53 s0 = (unsigned char *)localeconv()->decimal_point;
54 decimalpoint_cache = (unsigned char *)
55 MALLOC(strlen((char *)s0) + 1);
56 if (decimalpoint_cache) {
57 strcpy((char *)decimalpoint_cache, (char *)s0);
58 s0 = decimalpoint_cache;
59 }
60 }
61 decimalpoint = s0;
62#endif
63#endif
64
65 if (!hexdig['0'])
66 hexdig_init_D2A();
67 *bp = 0;
68 havedig = 0;
69 s0 = *(const unsigned char **)sp + 2;
70 while(s0[havedig] == '0')
71 havedig++;
72 s0 += havedig;
73 s = s0;
74 decpt = 0;
75 zret = 0;
76 e = 0;
77 if (hexdig[*s])
78 havedig++;
79 else {
80 zret = 1;
81#ifdef USE_LOCALE
82 for(i = 0; decimalpoint[i]; ++i) {
83 if (s[i] != decimalpoint[i])
84 goto pcheck;
85 }
86 decpt = s += i;
87#else
88 if (*s != '.')
89 goto pcheck;
90 decpt = ++s;
91#endif
92 if (!hexdig[*s])
93 goto pcheck;
94 while(*s == '0')
95 s++;
96 if (hexdig[*s])
97 zret = 0;
98 havedig = 1;
99 s0 = s;
100 }
101 while(hexdig[*s])
102 s++;
103#ifdef USE_LOCALE
104 if (*s == *decimalpoint && !decpt) {
105 for(i = 1; decimalpoint[i]; ++i) {
106 if (s[i] != decimalpoint[i])
107 goto pcheck;
108 }
109 decpt = s += i;
110#else
111 if (*s == '.' && !decpt) {
112 decpt = ++s;
113#endif
114 while(hexdig[*s])
115 s++;
116 }/*}*/
117 if (decpt)
118 e = -(((Long)(s-decpt)) << 2);
119 pcheck:
120 s1 = s;
121 big = esign = 0;
122 switch(*s) {
123 case 'p':
124 case 'P':
125 switch(*++s) {
126 case '-':
127 esign = 1;
128 /* no break */
129 case '+':
130 s++;
131 }
132 if ((n = hexdig[*s]) == 0 || n > 0x19) {
133 s = s1;
134 break;
135 }
136 e1 = n - 0x10;
137 while((n = hexdig[*++s]) !=0 && n <= 0x19) {
138 if (e1 & 0xf8000000)
139 big = 1;
140 e1 = 10*e1 + n - 0x10;
141 }
142 if (esign)
143 e1 = -e1;
144 e += e1;
145 }
146 *sp = (char*)s;
147 if (!havedig)
148 *sp = (char*)s0 - 1;
149 if (zret)
150 return STRTOG_Zero;
151 if (big) {
152 if (esign) {
153 switch(fpi->rounding) {
154 case FPI_Round_up:
155 if (sign)
156 break;
157 goto ret_tiny;
158 case FPI_Round_down:
159 if (!sign)
160 break;
161 goto ret_tiny;
162 }
163 goto retz;
164 ret_tiny:
165 b = Balloc(0);
166 b->wds = 1;
167 b->x[0] = 1;
168 goto dret;
169 }
170 switch(fpi->rounding) {
171 case FPI_Round_near:
172 goto ovfl1;
173 case FPI_Round_up:
174 if (!sign)
175 goto ovfl1;
176 goto ret_big;
177 case FPI_Round_down:
178 if (sign)
179 goto ovfl1;
180 goto ret_big;
181 }
182 ret_big:
183 nbits = fpi->nbits;
184 n0 = n = nbits >> kshift;
185 if (nbits & kmask)
186 ++n;
187 for(j = n, k = 0; j >>= 1; ++k);
188 *bp = b = Balloc(k);
189 b->wds = n;
190 for(j = 0; j < n0; ++j)
191 b->x[j] = ALL_ON;
192 if (n > n0)
193 b->x[j] = ULbits >> (ULbits - (nbits & kmask));
194 *expo = fpi->emin;
195 return STRTOG_Normal | STRTOG_Inexlo;
196 }
197 n = s1 - s0 - 1;
198 for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
199 k++;
200 b = Balloc(k);
201 x = b->x;
202 n = 0;
203 L = 0;
204#ifdef USE_LOCALE
205 for(i = 0; decimalpoint[i+1]; ++i);
206#endif
207 while(s1 > s0) {
208#ifdef USE_LOCALE
209 if (*--s1 == decimalpoint[i]) {
210 s1 -= i;
211 continue;
212 }
213#else
214 if (*--s1 == '.')
215 continue;
216#endif
217 if (n == ULbits) {
218 *x++ = L;
219 L = 0;
220 n = 0;
221 }
222 L |= (hexdig[*s1] & 0x0f) << n;
223 n += 4;
224 }
225 *x++ = L;
226 b->wds = n = x - b->x;
227 n = ULbits*n - hi0bits(L);
228 nbits = fpi->nbits;
229 lostbits = 0;
230 x = b->x;
231 if (n > nbits) {
232 n -= nbits;
233 if (any_on(b,n)) {
234 lostbits = 1;
235 k = n - 1;
236 if (x[k>>kshift] & 1 << (k & kmask)) {
237 lostbits = 2;
238 if (k > 0 && any_on(b,k))
239 lostbits = 3;
240 }
241 }
242 rshift(b, n);
243 e += n;
244 }
245 else if (n < nbits) {
246 n = nbits - n;
247 b = lshift(b, n);
248 e -= n;
249 x = b->x;
250 }
251 if (e > fpi->emax) {
252 ovfl:
253 Bfree(b);