d/freebsd-dev
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Import gdtoa sources dated 2011-03-04. This version includes a number

Browse Source 
of bugfixes, although I believe we already have local patches for the
ones people are likely to notice.

Per a request from arundel@, I also added the vendor's change log,
which is available separately from ftp://ftp.netlib.org/fp/.
This commit is contained in:
David Schultz 2011-03-09 06:14:33 +00:00
parent 9ea4d2a874
commit 21a2b1c905
Notes: svn2git 2020-12-20 02:59:44 +00:00 
svn path=/vendor/gdtoa/dist/; revision=219422
svn path=/vendor/gdtoa/20110304/; revision=219423; tag=vendor/gdtoa/20110304
52 changed files with 3629 additions and 758 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
README
View File

@ -353,5 +353,12 @@ you also compile with -DNO_LOCALE_CACHE, the details about the
current "decimal point" character string are cached and assumed not
to change during the program's execution.
On machines with a 64-bit long double and perhaps a 113-bit "quad"
type, you can invoke "make Printf" to add Printf (and variants, such
as Fprintf) to gdtoa.a. These are analogs, declared in stdio1.h, of
printf and fprintf, etc. in which %La, %Le, %Lf, and %Lg are for long
double and (if appropriate) %Lqa, %Lqe, %Lqf, and %Lqg are for quad
precision printing.
Please send comments to David M. Gay (dmg at acm dot org, with " at "
changed at "@" and " dot " changed to ".").

672
changes Normal file
View File

@ -0,0 +1,672 @@
Sun Jun 30 13:48:26 EDT 1991:
dtoa.c: adjust dtoa to allow negative ndigits for modes 3,5,7,9
(fixed-point mode); fix rounding bug in these modes when the input
d (to be converted) satisfies 10^-(ndigits+1) <= |d| < 10^-ndigits ,
i.e., when the result, before rounding, would be empty but might
round to one digit. Adjust the decpt returned in these modes when
the result is empty (i.e., when |d| <= 5 * 10^-ndigits).
Tue Jul 2 21:44:00 EDT 1991
Correct an inefficiency introduced 2 days ago in dtoa's handling of
integers in modes 0, 1.
Mon Sep 9 23:29:38 EDT 1991
dtoa.c: remove superfluous declaration of size_t.
Sun Oct 6 15:34:15 EDT 1991
dtoa.c: fix another bug in modes 3,5,7,9 when the result, before
rounding, would be empty, but rounds to one digit: *decpt was low by
one.
Sat Jan 18 12:30:04 EST 1992
dtoa.c: add some #ifdef KR_headers lines relevant only if IBM is
defined; for input decimal strings representing numbers too large, have
strtod return HUGE_VAL only if __STDC__ is defined; otherwise have it
return +-Infinity for IEEE arithmetic, +- the largest machine number
for IBM and VAX arithmetic. (If __STDC__ is not defined, HUGE_VAL may
not be defined either, or it may be wrong.)
Mon Apr 27 23:13:43 EDT 1992
dtoa.c: tweak strtod (one-line addition) so the end-pointer = start
pointer when the input has, e.g., only white space.
Thu May 7 18:04:46 EDT 1992
dtoa.c: adjust treatment of exponent field (in strtod) to behave
reasonably with huge numbers and 16-bit ints.
Fri Jun 19 08:29:02 EDT 1992
dtoa.c: fix a botch in placement of #ifdef __cplusplus (which only
matters if you're using a C++ compiler).
Wed Oct 21 11:23:07 EDT 1992
dtoa.c: add #ifdef Bad_float_h lines for systems with missing or
inferior float.h .
Thu Apr 22 07:54:48 EDT 1993
dtoa.c: change < to <= in line 2059:
< for(result_k = 0; sizeof(Bigint) - sizeof(unsigned long) + j < i;
---
> for(result_k = 0; sizeof(Bigint) - sizeof(unsigned long) + j <= i;
With 32-bit ints, the former could give too small a block for the return
value when, e.g., mode = 2 or 4 and ndigits = 24 (16 for 16-bit ints).
Mon Jun 21 12:56:42 EDT 1993
dtoa.c: tweak to work with 32-bit ints and 64-bit longs
when compiled with -DLong=int .
Wed Jan 26 11:09:16 EST 1994
dtoa.c: fix bug in strtod's handling of numbers with very
negative exponents (e.g. 1.8826e-512), which should underflow to 0;
fix storage leak in strtod with underflows and overflows near
the underflow and overflow thresholds.
Mon Feb 28 11:37:30 EST 1994
dtoa.c:
85a86,89
> * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
> * if memory is available and otherwise does something you deem
> * appropriate. If MALLOC is undefined, malloc will be invoked
> * directly -- and assumed always to succeed.
87a92,95
> #ifndef MALLOC
> #define MALLOC malloc
> #endif
>
352c360
< rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(Long));
---
> rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(Long));
Thu Mar 3 16:56:39 EST 1994
dtoa.c: if MALLOC is #defined, declare it.
Wed Jan 4 15:45:34 EST 1995
dtoa.c: add CONST qualification to tens, bigtens, tinytens (for use
on embedded systems with little spare RAM).
Fri Mar 1 08:55:39 EST 1996
g_fmt.c: honor the sign of 0 and return the first argument (buf).
Sat Jul 6 07:59:28 EDT 1996
dtoa.c: cosmetic changes: "ULong" rather than "unsigned Long";
update comments to reflect AT&T breakup.
Mon Aug 5 23:31:24 EDT 1996
dtoa.c: add comment about invoking _control87(PC_53, MCW_PC)
(or the equivalent) on 80x87 machines before calling strtod or dtoa.
Tue Dec 17 15:01:56 EST 1996
dtoa.c: new #define possibilities: #define INFNAN_CHECK to have
strtod check (case insensitively) for "Infinity" and "NaN" on machines
with IEEE arithmetic; #define MULTIPLE_THREADS if the system offers
preemptively scheduled multiple threads, in which case you must supply
routines ACQUIRE_DTOA_LOCK(n) and FREE_DTOA_LOCK(n) (n = 0 or 1).
New void freedtoa(char*) for freeing values returned by dtoa; use of
freedtoa() is required if MULTIPLE_THREADS is #defined, and is merely
recommended otherwise.
g_fmt.c: adjusted to invoke freedtoa().
Wed Feb 12 00:40:01 EST 1997
dtoa.c: strtod: on IEEE systems, scale to avoid intermediate
underflows when the result does not underflow; compiling with
-DNO_IEEE_Scale restores the old logic. Fix a bug, revealed by
input string 2.2250738585072012e-308, in treating input just less
than the smallest normalized number. (The bug introduced an extra
ULP of error in this special case.)
Tue May 12 11:13:04 EDT 1998
dtoa.c: strtod: fix a glitch introduced with the scaling of 19970212
that caused one-bit rounding errors in certain denormal numbers, such
as 8.44291197326099e-309, which was read as 8.442911973260987e-309.
Remove #ifdef Unsigned_Shifts logic in favor of unsigned arithmetic.
Unless compiled with -DNO_LONG_LONG, use 64-bit arithmetic where
possible.
Fri May 15 07:49:07 EDT 1998
dtoa.c: strtod: fix another glitch with scaling to avoid underflow
with IEEE arithmetic, again revealed by the input string
2.2250738585072012e-308, which was rounded to the largest denormal
rather than the smallest normal double precision number.
Wed Aug 5 23:27:26 EDT 1998
gdtoa.tar.gz: tweaks in response to comments from Shawn C. Sheridan
(with no effect on the resulting .o files except when strtod.c is
compiled with -DNO_ERRNO); bigtens --> bigtens_D2A (a symbol meant
to be private to gdtoa.a).
Sat Sep 12 17:05:15 EDT 1998
gdtoa.tar.gz: more changes in response to comments from Shawn C.
Sheridan (including repair of a glitch in g_ffmt.c). For consistency
and possible convenience, there are some new functions and some name
changes to existing ones:
Old New
--- g_xLfmt
strtoQ strtopQ
--- strtopd
strtodd strtopdd
--- strtopf
strtox strtopx
--- strtopxL
--- strtorxL
--- strtoIxL
Functions strtopd and strtopf are variations of strtod and strtof,
respectively, which write their results to their final (pointer)
arguments. Functions strtorf and strtord are now analogous to the
other strtor* functions in that they now have a final pointer
argument to which they write their results, and they return the
int value they get from strtodg.
The xL family (g_xLfmt, strto[Irp]xL) is a variation of the old x
family (for 80-bit IEEE double-extended precision) that assumes the
storage layout of the Motorola 68881's double-extended format: 80
interesting bits stored in 3 unsigned 32-bit ints (with a "hole", 16
zero bits, in the word that holds the sign and exponent). The x
family now deals with 80-bit (5 unsigned 16-bit ints) rather than
96-bit arrays (3 unsigned 32-bit ints) that hold its 80-bit
double-extended values. (This relaxes the alignment requirements of
the x family and results in strto[Ipr]x writing 80 rather than 96 bits
to their final arguments.)
Each g_*fmt routine now returns a pointer to the null character
that terminates the strings it writes, rather than a pointer to
the beginning of that string (the first argument). These routines
still return 0 (NULL) if the first argument is too short.
The second argument to g_dfmt is now pointer (to a double) rather
than a double value.
Thu Oct 29 21:54:00 EST 1998
dtoa.c: Fix bug in strtod under -DSudden_Underflow and (the default)
-DAvoid_Underflow: some numbers that should have suffered sudden
underflow were scaled inappropriately (giving nonzero return values).
Example: "1e-320" gave -3.6304123742133376e+280 rather than 0.
Mon Nov 2 15:41:16 EST 1998
dtoa.c: tweak to remove LL suffixes from numeric constants (for
compilers that offer a 64-bit long long type but do not recognize the
LL constants prescribed by C9x, the proposed update to the ANSI/ISO C
standard). Thanks to Earl Chew for pointing out the existence of such
compilers.
gdtoa.tar.gz: renamed gdtoa.tgz and updated to incorporate the above
changes (of 29 Oct. and 2 Nov. 1998) to dtoa.c.
Thu Mar 25 17:56:44 EST 1999
dtoa.c, gdtoa.tgz: fix a bug in strtod's reading of 4.9e-324:
it returned 0 rather than the smallest denormal.
Mon Apr 12 10:39:25 EDT 1999
gdtoa.tgz: test/ftest.c: change %.7g to %.8g throughout.
Fri Aug 20 19:17:52 EDT 1999
gdtoa.tgz: gdtoa.c: fix two bugs reported by David Chase (thanks!):
1. An adjustment for denormalized numbers around 503 was off by one.
2. A check for "The special case" around line 551 omitted the condition
that we not be at the bottom of the exponent range.
Mon Sep 13 10:53:34 EDT 1999
dtoa.c: computationally invisible tweak for the benefit of people
who actually read the code:
2671c2671
< && word0(d) & Exp_mask
---
> && word0(d) & (Exp_mask & Exp_mask << 1)
I.e., in dtoa(), the "special case" does not arise for the smallest
normalized IEEE double. Thanks to Waldemar Horwat for pointing this
out and suggesting the modified test above. Also, some tweaks for
compilation with -DKR_headers.
gdtoa.tgz: gdtoa.c: analogous change:
552c552
< if (bbits == 1 && be0 > fpi->emin) {
---
> if (bbits == 1 && be0 > fpi->emin + 1) {
This has no effect on the g*fmt.c routines, but might affect the
computation of the shortest decimal string that rounds to the
smallest normalized floating-point number of other precisions.
gdota.tgz: test/d.out test/dI.out test/dd.out: updated to reflect
previous changes (of 19990820); test/*.c: most test programs modified
to permit #hex input. See the comments.
Fri Sep 17 01:39:25 EDT 1999
Try again to update dtoa.c: somehow dtoa.c got put back to a version
from 3 years ago after this "changes" file was updated on 13 Sept. 1999.
One more tweak to omit a warning on some systems:
2671c2671
< && word0(d) & (Exp_mask & Exp_mask << 1)
---
> && word0(d) & (Exp_mask & ~Exp_msk1)
Plus changes to avoid trouble with aggressively optimizing compilers
(e.g., gcc 2.95.1 under -O2). On some systems, these changes do not
affect the resulting machine code; on others, the old way of viewing
a double as a pair of ULongs is available with -DYES_ALIAS.
Tue Sep 21 09:21:25 EDT 1999
gdtoa.tgz: changes analogous to those of 17 Sept. 1999 to dtoa.c to
avoid trouble with aggressively optimizing compilers.
Wed Dec 15 13:14:38 EST 1999
dtoa.c: tweak to bypass a bug with HUGE_VAL on HP systems.
Mon Jan 17 18:32:52 EST 2000
dtoa.c and gdtoa.tgz: strtod: set errno = ERANGE on all inputs that
underflow to zero (not just those sufficiently less than the smallest
positive denormalized number).
gdtoa.tgz: README: point out that compiling with -DNO_ERRNO inhibits
errno assignments (by strtod and the core converter, strtodg).
Tue Jan 18 16:35:31 EST 2000
dtoa.c and gdtoa.tgz: strtod: modify the test inserted yesterday so
it may work correctly with buggy 80x87 compilers. (The change matters,
e.g., to Microsoft Visual C++ 4.2 and 6.0.)
Thu Nov 2 21:00:45 EST 2000
dtoa.c and gdtoa.tgz:
1. Fix bug in test for exact half-way cases of denormalized numbers
(without -DNO_IEEE_Scale).
2. Compilation with -DNO_ERRNO prevents strtod from assigning
errno = ERANGE when the result overflows or underflows to 0.
3. With IEEE arithmetic and no -DNO_IEEE_Scale, adjust scaling so
ulp(d) never returns a denormalized number. This and other tweaks
permit strtod and dtoa to work correctly on machines that flush
underflows to zero but otherwise use IEEE arithmetic without
Sudden_Underflow being #defined (and with strtod simply returning 0
instead of denormalized numbers).
4. Compilations with -DUSE_LOCALE causes strtod to use the current
locale's decimal_point value.
5. Under compilations with -DINFNAN_CHECK, strtod and strtodg (case
insensitively) treat "inf" the same as "infinity" and, unless
compiled with -DNo_Hex_NaN, accept "nan(x)", where x is a string of
hexadecimal digits and spaces, as a NaN whose value is constructed
from x (as explained in comments near the top of dtoa.c and in
gdtoaimp.h).
6. The default PRIVATE_MEM is increased slightly (to 2304), and comments
near the top of dtoa.c provide more discussion of PRIVATE_MEM.
7. Meanings of dtoa modes 4,5,8,9 changed. See comments in dtoa.c and
gdtoa.c; modes 4 and 5 may now provide shorter strings that round
(in round-nearest mode) to the given double value. (Paxson's
testbase program is unhappy with this new rounding, as it can
introduce an error of more than one base-10 ulp when 17 or more
decimal digits are requested.)
8. With IEEE arithmetic, compilation with -DHonor_FLT_ROUNDS causes
strtod and dtoa to round according to FLT_ROUNDS:
0 ==> towards 0,
1 ==> nearest,
2 ==> towards +Infinity,
3 ==> towards -Infinity.
9. With IEEE arithmetic, compilation with -DSET_INEXACT causes extra
computation (and sometimes slower conversions in dtoa and strtod,
particularly for dtoa in cases where otherwise some simple floating-
point computations would suffice) to set the IEEE inexact flag
correctly. As comments in dtoa.c explain in more detail, this
requires compilation in an environment (such as #include "dtoa.c"
in suitable source) that provides
int get_inexact(void);
void clear_inexact(void);
10. On input "-x", return 0 rather than -0.
gdtoa.tgz: gethex.c: adjust logic for reading hex constants to accord
with current wording of the C99 standard. Previously, I thought hex
constants read by strtod and friends had to have either a decimal point
or an exponent field; p. 307 of the C99 standard states that both are
optional. Because of the complexity of this reading, it is available
only in the variant of strtod that appears in gdtoa.tgz.
strtodg (gdtoa.tgz): New return value STRTOG_NaNbits (with
STRTOG_NoNumber renumbered). Allow STRTOG_Neg bit in strtodg returns
for STRTOG_NaN and STRTOG_NaNbits.
gdtoa.tgz: Fix uninitialized variable bug in g_Qfmt.c's handling of NaNs.
Mon Nov 13 14:00:05 EST 2000
gdtoa.tgz: strtodg: fix a storage leak and an apparently rare infinite
loop with a boundary case of directed rounding. Example input to
gdtoa/test/Qtest where the loop bug bit:
r 3
35184372088831.999999999999999999999999999999999999
This was revealed by testbase for quad precision Solaris arithmetic;
it did not show up in several other testbase configurations.
Wed Feb 7 12:56:11 EST 2001
dtoa.c: fix bug (possible infinite loop, e.g., with
2.47032822920623272e-324) introduced 20001113 in handling the special
case of a power of 2 to be rounded down one ulp. Add test (required
by changes of 20001113) for the extra special case of 2^-1075 (half
the smallest denormal).
gdtoa.tgz: corresponding adjustments to strtod.c and strtodg.c.
Tue Mar 13 00:46:09 EST 2001
gdtoa.tgz: gdtoa/strtodg.c: fix bug in handling values exactly half
an ulp less than the smallest normal floating-point number;
gdtoa/*test.c: fix glitch in handling "n ..." lines (intended to
change "ndig").
Wed Mar 6 10:13:52 EST 2002
gdtoa.tgz: add gdtoa/test/strtodt.c and gdtoa/test/testnos3 to test
strtod on hard cases posted by Fred Tydeman to comp.arch.arithmetic on
26 Feb. 1996. Add comment to gdtoa/README about strtod requiring true
IEEE arithmetic (with 53-bit rounding precision on 80x87 chips).
In gdtoa/test, automate selection of expected output files [xQ]*.out.
Wed Mar 5 10:35:41 EST 2003
gdtoa.tgz: fix a bug in strtod's handling of 0-valued 0x... "decimal"
strings. A fault was possible. Thanks to David Shultz for reporting
this bug.
Tue Mar 18 09:38:28 EST 2003
gdtoa.tgz: fix a glitch in strtodg.c with -DUSE_LOCALE; add #ifdef
USE_LOCALE lines to g__fmt.c (to affect binary --> decimal conversions
via the g*fmt routines), and add comments about -DUSE_LOCALE to README.
In short, compiling strtod.c, strtodg.c, and g__fmt.c with -DUSE_LOCALE
causes them to determine the decimal-point character from the current
locale. (Otherwise it is '.'.)
Fri Mar 21 16:36:27 EST 2003
gdtoa.tgz: gethex.c: add #ifdef USE_LOCAL logic; strtod.c: fix a
glitch in handling 0x... input (the return from gethex was ignored).
Wed Mar 26 15:35:10 EST 2003
gdtoa.tgz: gethex.c: pedantic (and normally invisible) change:
use unsigned char for decimalpoint variable (under -DUSE_LOCALE).
Sat Jan 17 23:58:52 MST 2004
gdtoa.tgz: gethex.c: supply missing parens in test for whether a
denormal result should be zero, correct logic for rounding up when the
result is denormal, and when returning zero or Infinity, set *bp = 0;
strtod.c: switch on gethex(...) & STRTOG_Retmask rather than just on
gethex(), and only copybits(..., bb) when bb is nonzero. This
mattered for underflows and overflows in 0x notation.
Thu Mar 25 22:34:56 MST 2004
dtoa.c and gdtoa.c/misc.c: change "(!x & 1)" to "(!x)" to avoid
confusion by human readers -- the object code is unaffected (with
reasonable compilers).
Mon Apr 12 00:44:22 MDT 2004
dtoa.c and gdtoa.tar.gz: update contact info. for dmg and correct
page numbers in comment on Steele & White (1990).
gdtoa.tgz: add strtodnrp.c for a variant of strtod that is slower
but does not require 53-bit rounding precision on Intel IA32 systems.
Tue Apr 13 00:28:14 MDT 2004
gdtoa.tgz: strtod.c: fix glitch when both INFNAN_CHECK and No_Hex_NaN
are #defined. Thanks to David Mendenhall for pointing this bug out.
Wed Jan 5 22:39:17 MST 2005
gdtoa.tgz:
gethex.c: fix the bug reported by Stefan Farfeleder of ignoring a
binary-exponent-part if the converted number is zero.
strto[pr]x.c: fix bug reported by Stefan Farfeleder in setting the
exponent of denormals (which should be 0, not 1).
g_xfmt.c: fix a corresponding bug with denormals.
strtodg.c: fix a bug under IBM (base 16) arithemtic reported
by Greg Alexander: a correction to the binary exponent for changes to
the exponent of a native double value for avoiding overflow had to be
multiplied by 4 ("e2 <<= 2;").
Various files: minor tweaks for portability.
Sat Jan 15 15:36:03 MST 2005
gdtoa.tgz: gethex.c: fix a glitch introduced last week (and reported
by Stefan Farfelder) with 0x forms with no nonzero digits before the "."
character, e.g., 0x.1 (which was rendered as 0 rather than .0625).
gdtoa.tgz: many files: add automatic computation of gd_qnan.h for
giving the system-dependent format of a quiet NaN (the one generated
for Infinity - Infinity). Tweak test/makefile so differences in the
spelling of Infinity ("INF" or "Inf" on some systems) do not matter.
Fix bug in strtod.c and strtodg.c under which, e.g., -.nan was read
as NaN rather than unacceptable input (causing return 0). Adjust
comments in README about nan(...). Fix glitch in test/dt.c.
Sun Jan 16 18:22:13 MST 2005
gdtoa.tgz: strtodg.c: fix long-standing bug in handling input
that rounds up to 2^nbits, e.g., strtof("16777215.5"). Thanks to
Edward Moy for reporting this problem.
gdtoa.tgz: Fix some bugs with -DJust_16.
Thu Sep 22 22:40:16 MDT 2005
gdtoa.tgz:
strtod.c: unless prevented by -DNO_FENV_H, include C99's fenv.h
and with hex input, get the current rounding mode from fegetround().
With decimal input, strtod honors the rounding mode automatically.
Thanks to David Schultz (das at FreeBSD dot ORG) for pointing
strtodg.c: fix a bug with handling numbers very near the largest
possible one, which were sometimes incorrectly converted to Infinity.
Thanks to Edward Moy (emoy at apple dot com) for pointing this out.
g_Qfmt.c: change strcpy to strcp. Thanks to J. T. Conklin
(jtc at acorntoolworks dot com) for pointing this out.
test/xtest.c: fix some subscript bugs.
test/x.ou0, test/x.ou1, test/xL.: update in response to the above fix to
test/xtest.c.
test/makefile: add -lm to some link lines (needed for fegetround).
Sun Jan 21 20:26:44 MST 2007
gdtoa.tgz:
strtodg.c: fix a botch in the test of whether to increase rvbits
before terminating the big for(;;) loop with dsign true: change
if (hi0bits(rvb->x[(rvb->wds - 1) >> kshift])
!= j)
rvbits++;
to
if (hi0bits(rvb->x[rvb->wds - 1]) != j)
rvbits++;
Example of input where this bug bit: 1.9e27. Thanks to Edward Moy
<emoy@apple.com> for providing this example. Also, simplify the
preceding computation of j.
test/README: add comment that strtodt needs to operate with 53-bit
rounding precision on Intel x86 systems, and add a pointer to Paxson's
paper.
Sat Mar 15 11:44:31 MDT 2008
dtoa.c and gdtoa.tgz: with -DINFNAN_CHECK and without
-DGDOTA_NON_PEDANTIC_NANCHECK, conform to the ill-advised prescription
in the C99 standard of consuming (...) in "nan(...)" even when ...
is not of the expected form. Allow an optional 0x or 0X to precede
the string of hex digits in the expected form of ... .
gdtoa.tgz: gethex.c: have, e.g., strtod("0xyz",&se) set se to "xyz".
Previously it was incorrectly set to "0xyz".
Thu Aug 28 22:37:35 MDT 2008
dtoa.c and gdtoa.tgz: Fix a bug in strtod when compiled with
-DHonor_FLT_ROUNDS: in rounding modes other than "to nearest",
strtod looped on input larger than and within a factor of 2 of
the largest finite floating-point number. Since FLT_ROUNDS is buggy
on some (Linux) systems in that it does not reflect calls on
fesetround(), when Honor_FLT_ROUNDS is #defined, get the curren
rounding mode from fegetround() rather than FLT_ROUNDS, unless
Trust_FLT_ROUNDS is also #defined.
gdtoa/test/getround.c in gdtoa.tgz: simply report the current
rounding mode when the input line is "r" by itself. (Previously it
did so, but also complained of invalid input.)
gdtoa/gethex.c: fix an off-by-one bug in a rounding test; detect and
deal with huge exponents (positive or negative). This affected the
reading of hexadecimal floating-point values (0x...). Also set errno
to ERANGE on out-of-range values (unless compiled with -DNO_ERRNO).
gdtoa/strtod.c: adjust scaling of tinytens[4] (as in dtoa.c) to
avoid double rounding when dealing with numbers near the bottom of
the exponent range.
Sat Aug 30 23:37:07 MDT 2008
gdtoa/gethex.c: ensure *bp is set to something (NULL if nothing else).
Bring gdtoa/xsum0.out and gdtoa/test/xsum0.out up to date.
Tue Sep 9 22:08:30 MDT 2008
gdtoa/strto*.c and gdtoa/*fmt.c: if compiled with -DUSE_LOCALE, use
the current locale's decimal point character string.
gdtoa/gdtoa.c: trim trailing zeros in a missed case (e.g., ndigits = 6
on 1020302).
dtoa.c and gdtoa/strtod.c: on systems with IEEE arithmetic (and without
NO_ERRNO being defined) set ERANGE for denormal values as well as real
underflows.
gdtoa/strtodg.c: fix an off-by-one bug in rounding to the largest
representable magnitude when nbits is a multiple of 32.
gdtoa/*fmt.c and gdtoa/gdtoa.h: bufsize changed from unsigned to size_t.
gdtoaimp.h, *fmt.c: change in calling sequence to internal g__fmt(),
which now explicitly checks bufsize.
Relevant routines (see README) honor the current rounding mode if
compiled with -DHonor_FLT_ROUNDS on IEEE-arithmetic systems that provide
the C99 fegetround() function.
gdtoa/test/getround.c can optionally be compiled with
-DHonor_FLT_ROUNDS and/or -DUSE_MY_LOCALE for manual testing of gdtoa.a
compiled with -DHonor_FLT_ROUNDS or -DUSE_LOCALE.
Fri Oct 10 20:07:15 MDT 2008
gdtoa/gethex.c: fix a bug reading hexadecimal floating-point values
starting with "0xd" for a nonzero digit d (such as "0x1.0002p3"). The
bug caused the values to be read as zero with endptr set incorrectly.
Tue Oct 28 00:14:08 MDT 2008
gdtoa/strtod.c: fix a comment glitch (with commented {}).
Tue Nov 11 23:05:25 MST 2008
gdtoa: fix a glitch in the strto* routines when compiled with
-DUSE_LOCALE and the locale's decimal-point string is two or more
characters long. Wrong conversions were then possible.
Fri Dec 5 18:20:36 MST 2008
gdtoa.tgz: fix bugs with reading C99-style hexadecimal floating-point
values when compiled with -DPack_16; on IEEE-arithmetic systems, make
INFNAN_CHECK the default unless NO_INFNAN_CHECK is #defined. (This is
consistent with dtoa.c, which has worked this way for a while.)
dtoa.c: add recognition of C99-style hexadecimal floating-point
values (unless compiled with NO_HEX_FP is #defined).
Thu Dec 11 23:10:23 MST 2008
dtoa.c: omit an unused variable.
Fri Jan 2 22:45:33 MST 2009
dtoa.c: tweak to banish some compiler warnings.
Sun Mar 1 20:57:22 MST 2009
dtoa.c, gdtoa/{g__fmt.c, gethex.c, strtod.c, strtodg.c}: change malloc
to MALLOC.
dtoa.c and gdtoa/gdtoaimp.h and gdtoa/misc.c: reduce Kmax, and use
MALLOC and FREE or free for huge blocks, which are possible only in
pathological cases, such as dtoa calls in mode 3 with thousands of
digits requested, or strtod() calls with thousand of digits. For the
latter case, I have an alternate approach that runs much faster
and uses less memory, but finding time to get it ready for distribution
may take a while.
Mon Mar 16 00:32:43 MDT 2009
dtoa.c: Fix a bug under -DUSE_LOCALE in handling "decimal point"
strings more than one character long.
dtoa.c and gdtoa/misc.c: Remove a buggy test activated with
-DDEBUG.
dtoa.c and gdtoa/gdtoa.c: simplify logic for "4 leading 0 bits".
dtoa.c: Add logic (that can be disabled with -DNO_STRTOD_BIGCOMP
and that) to strtod for more efficiently handling a very long input
string. It proceeds by initially truncating the input string, then if
necessary comparing the whole string with a decimal expansion to
decide close cases. This logic is only used for input more than
STRTOD_DIGLIM digits long (default 40), and for now only applies to
IEEE arithmetic (for want of other kinds of platforms on which to run
tests). This only appears worthwhile for absurdly long input strings,
so a corresponding update to gdtoa does not seem warranted.
dtoa.c, gdtoa.tgz: tweaks (mostly adding unnecessary parens) to
silence "gcc -Wall" warnings. Aside from a couple of minor changes
to banish erroneous warnings about uninitialized variables, the tweaks
do not affect the generated object code.
Sat Apr 11 23:25:58 MDT 2009
dtoa.c: fix glitch in compiling with -DNo_Hex_NaN and the bug of
accepting (e.g.) ".nan" or ".inf" as NaN or Infinity.
gdtoa.tgz: tweaks to silence warnings from "gcc -Wstrict-aliasing=2";
update xsum0.out files.
Sun Apr 19 23:40:24 MDT 2009
dtoa.c, gdtoa/misc.c: do not attempt to allocate large memory blocks
from the private memory pool (which was an unlikely event, but a bug).
gdtoa/strtopx.c, gdtoa/strtopxL.c, gdtoa/strtorx.c, gdtoa/strtorxL.c:
supply explicit bit for Infinity. Note that the Q routines (which do
not supply this bit) are appropriate for Sparc quad precision (probably
known as long double with most current compilers).
Wed Dec 9 08:14:52 MST 2009
gdtoa.tgz: add gdtoa/printf.c* and modify makefile so "make Printf"
adds a printf to gdtoa.a (to be accessed with #include "stdio1.h" to
get gdtoa/stdio1.h, which you might install in some standard place).
On Intel/AMD i386, x86_64, and Sparc systems, this adds formats %La,
%Le, %Lf and %Lg to handle long double. On x86_64 systems, it also
adds %Lqa, %Lqe, %Lqf and %Lqg to handle 128-bit bit types (called
__float128 by gcc and _Quad by the Intel compiler). In gdtoa/test,
"make pf_test" tests this printf (provided the system is an i386,
x86_64, or Sparc system). On x86_64 systems, "make pf_testLq" tests
the %Lq... formats (briefly).
Mon Jan 11 22:25:17 MST 2010
dtoa.c: fix a minor performance bug and, under compilation with -DDEBUG,
an erroneous error message "oversize b in quorem" in strtod's processing
of some input that underflows to zero. Also fix a bug in bigcomp()'s
handling of numbers that will scale to denormal values. The increments
that bigcomp applied were ignoring the effects of denormalization.
Sat Jan 23 00:25:54 MST 2010
dtoa.c: Fix some glitches in recently introduced changes meant to
speed up returns in pedantic cases. In quorem, adjust #ifdef DEBUG
stuff so it does not complain when bigcomp() calls quorem on input
near the smallest representable number and rounding up by a bit causes
a quorem return > 9 (which in this case is not a bug). Fix a memory
leak in the unlikely case of overflow only being detected after some
high-precision integer computations. Fix an off-by-one bug in
handling a large number of digits with a few nonzero digits, followed
by many zeros, and then some nonzero digits. (This does not happen
with sensible input.) Fix an off-by-one bug in a recently introduced
quick test for underflow (i.e., zero result) on input at the bottom of
the exponent range. Thanks to Mark Dickinson for pointing these bugs
out.
dtoa.c and gdtoa/strtod.c: Fix an obscure bug in strtod's handling
of some inputs of many digits at the bottom of the exponent range:
results were sometimes off by a bit when gdtoa/strtod.c or dtoa.c was
compiled without -DNO_IEEE_SCALE and, for dtoa.c, when compiled with
-DNO_STRTOD_BIGCOMP.
gdtoa/test/testnos3: add some examples that went wrong before
the present changes.
Sat Jan 23 23:29:02 MST 2010
dtoa.c: more tweaks relevant only to absurd input.
Tue Feb 2 23:05:34 MST 2010
dtoa.c: add test for setting errno = ERANGE when input of many digits
is rounded to Infinity or underflows to zero. Fix a memory leak in
such instances.
gdtoa/strtod.c: make some corresponding changes.
Wed Jul 7 09:25:46 MDT 2010
dtoa.c: adjust to use bigcomp when necessary when compiled with
-DHonor_FLT_ROUNDS (and without -DNO_STRTOD_BIGCOMP), and the rounding
mode is torwards +Infinity. An input (supplied by Rick Regan
<exploringbinary@gmail.com>) where this matters is
1.100000000000000088817841970012523233890533447265626
gdtoa/strtod.c: fix errors (introduced 20090411) when compiled
with -DHonor_FLT_ROUNDS.
Wed Sep 15 09:00:26 MDT 2010
dtoa.c, gdtoa/dtoa.c, gdtoa/gdtoa.c: fix bugs with -DROUND_BIASED
pointed out by Jay Foad.
Mon Sep 27 13:43:30 MDT 2010
gdtoa/gdtoa.c: fix a glitch (not revealed by compilation) in the
changes of 15 Sept. 2010.
Fri Nov 5 13:02:41 MDT 2010
dtoa.c: fix a bug related to bigcomp: decimal strings with all
zeros before the decimal point more than 40 significant digits that
required use of bigcomp might be converted very incorrectly.
Example: .010000000000000000057612911342378542997169 .
Thanks to Rick Regan <exploringbinary@gmail.com> for reporting the
symptoms and providing an example.
20110403:
dtoa.c, gdtoa/gdtoaimp.h, gdtoa/strtod.c: if
ROUND_BIASED_without_Round_Up is #defined, assume ROUND_BIASED and
omit the quick computation that would use ordinary arithmetic to
compute the correctly rounded result with one rounding error. If you
want biased rounding with IEEE-style format "double" and will operate
with rounding toward +Infinity, it suffices to #define ROUND_BIASED
(and thus retain the quick computation when it is appropriate).
gdtoa/gdtoa.h: change default Long from long to int (with the goal
of portability when compiling without -DLong=... specified). On some
64-bit systems, long is a 64-bit type; we need a 32-bit type here.
dtoa.c, gdtoa/gdtoa.c: fix a glith with ndigits with mode = 4 at
the bottom of the exponent range, e.g., 1e-323.

168
dtoa.c
View File

@ -75,10 +75,10 @@ THIS SOFTWARE.
char *
dtoa
#ifdef KR_headers
(d, mode, ndigits, decpt, sign, rve)
double d; int mode, ndigits, *decpt, *sign; char **rve;
(d0, mode, ndigits, decpt, sign, rve)
double d0; int mode, ndigits, *decpt, *sign; char **rve;
#else
(double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
(double d0, int mode, int ndigits, int *decpt, int *sign, char **rve)
#endif
{
/* Arguments ndigits, decpt, sign are similar to those
@ -124,7 +124,8 @@ dtoa
ULong x;
#endif
Bigint *b, *b1, *delta, *mlo, *mhi, *S;
double d2, ds, eps;
U d, d2, eps;
double ds;
char *s, *s0;
#ifdef SET_INEXACT
int inexact, oldinexact;
@ -149,35 +150,35 @@ dtoa
dtoa_result = 0;
}
#endif
if (word0(d) & Sign_bit) {
d.d = d0;
if (word0(&d) & Sign_bit) {
/* set sign for everything, including 0's and NaNs */
*sign = 1;
word0(d) &= ~Sign_bit; /* clear sign bit */
word0(&d) &= ~Sign_bit; /* clear sign bit */
}
else
*sign = 0;
#if defined(IEEE_Arith) + defined(VAX)
#ifdef IEEE_Arith
if ((word0(d) & Exp_mask) == Exp_mask)
if ((word0(&d) & Exp_mask) == Exp_mask)
#else
if (word0(d) == 0x8000)
if (word0(&d) == 0x8000)
#endif
{
/* Infinity or NaN */
*decpt = 9999;
#ifdef IEEE_Arith
if (!word1(d) && !(word0(d) & 0xfffff))
if (!word1(&d) && !(word0(&d) & 0xfffff))
return nrv_alloc("Infinity", rve, 8);
#endif
return nrv_alloc("NaN", rve, 3);
}
#endif
#ifdef IBM
dval(d) += 0; /* normalize */
dval(&d) += 0; /* normalize */
#endif
if (!dval(d)) {
if (!dval(&d)) {
*decpt = 1;
return nrv_alloc("0", rve, 1);
}
@ -196,26 +197,26 @@ dtoa
}
#endif
b = d2b(dval(d), &be, &bbits);
b = d2b(dval(&d), &be, &bbits);
#ifdef Sudden_Underflow
i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
#else
if (( i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
if (( i = (int)(word0(&d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)) )!=0) {
#endif
dval(d2) = dval(d);
word0(d2) &= Frac_mask1;
word0(d2) |= Exp_11;
dval(&d2) = dval(&d);
word0(&d2) &= Frac_mask1;
word0(&d2) |= Exp_11;
#ifdef IBM
if (( j = 11 - hi0bits(word0(d2) & Frac_mask) )!=0)
dval(d2) /= 1 << j;
if (( j = 11 - hi0bits(word0(&d2) & Frac_mask) )!=0)
dval(&d2) /= 1 << j;
#endif
/* log(x) ~=~ log(1.5) + (x-1.5)/1.5
* log10(x) = log(x) / log(10)
* ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
* log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
* log10(&d) = (i-Bias)*log(2)/log(10) + log10(&d2)
*
* This suggests computing an approximation k to log10(d) by
* This suggests computing an approximation k to log10(&d) by
*
* k = (i - Bias)*0.301029995663981
* + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
@ -244,21 +245,21 @@ dtoa
/* d is denormalized */
i = bbits + be + (Bias + (P-1) - 1);
x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32
: word1(d) << 32 - i;
dval(d2) = x;
word0(d2) -= 31*Exp_msk1; /* adjust exponent */
x = i > 32 ? word0(&d) << (64 - i) | word1(&d) >> (i - 32)
: word1(&d) << (32 - i);
dval(&d2) = x;
word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
i -= (Bias + (P-1) - 1) + 1;
denorm = 1;
}
#endif
ds = (dval(d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
k = (int)ds;
if (ds < 0. && ds != k)
k--; /* want k = floor(ds) */
k_check = 1;
if (k >= 0 && k <= Ten_pmax) {
if (dval(d) < tens[k])
if (dval(&d) < tens[k])
k--;
k_check = 0;
}
@ -297,10 +298,11 @@ dtoa
try_quick = 0;
}
leftright = 1;
ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
/* silence erroneous "gcc -Wall" warning. */
switch(mode) {
case 0:
case 1:
ilim = ilim1 = -1;
i = 18;
ndigits = 0;
break;
@ -334,7 +336,7 @@ dtoa
/* Try to get by with floating-point arithmetic. */
i = 0;
dval(d2) = dval(d);
dval(&d2) = dval(&d);
k0 = k;
ilim0 = ilim;
ieps = 2; /* conservative */
@ -344,7 +346,7 @@ dtoa
if (j & Bletch) {
/* prevent overflows */
j &= Bletch - 1;
dval(d) /= bigtens[n_bigtens-1];
dval(&d) /= bigtens[n_bigtens-1];
ieps++;
}
for(; j; j >>= 1, i++)
@ -352,32 +354,32 @@ dtoa
ieps++;
ds *= bigtens[i];
}
dval(d) /= ds;
dval(&d) /= ds;
}
else if (( j1 = -k )!=0) {
dval(d) *= tens[j1 & 0xf];
dval(&d) *= tens[j1 & 0xf];
for(j = j1 >> 4; j; j >>= 1, i++)
if (j & 1) {
ieps++;
dval(d) *= bigtens[i];
dval(&d) *= bigtens[i];
}
}
if (k_check && dval(d) < 1. && ilim > 0) {
if (k_check && dval(&d) < 1. && ilim > 0) {
if (ilim1 <= 0)
goto fast_failed;
ilim = ilim1;
k--;
dval(d) *= 10.;
dval(&d) *= 10.;
ieps++;
}
dval(eps) = ieps*dval(d) + 7.;
word0(eps) -= (P-1)*Exp_msk1;
dval(&eps) = ieps*dval(&d) + 7.;
word0(&eps) -= (P-1)*Exp_msk1;
if (ilim == 0) {
S = mhi = 0;
dval(d) -= 5.;
if (dval(d) > dval(eps))
dval(&d) -= 5.;
if (dval(&d) > dval(&eps))
goto one_digit;
if (dval(d) < -dval(eps))
if (dval(&d) < -dval(&eps))
goto no_digits;
goto fast_failed;
}
@ -386,34 +388,34 @@ dtoa
/* Use Steele & White method of only
* generating digits needed.
*/
dval(eps) = 0.5/tens[ilim-1] - dval(eps);
dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
for(i = 0;;) {
L = dval(d);
dval(d) -= L;
L = dval(&d);
dval(&d) -= L;
*s++ = '0' + (int)L;
if (dval(d) < dval(eps))
if (dval(&d) < dval(&eps))
goto ret1;
if (1. - dval(d) < dval(eps))
if (1. - dval(&d) < dval(&eps))
goto bump_up;
if (++i >= ilim)
break;
dval(eps) *= 10.;
dval(d) *= 10.;
dval(&eps) *= 10.;
dval(&d) *= 10.;
}
}
else {
#endif
/* Generate ilim digits, then fix them up. */
dval(eps) *= tens[ilim-1];
for(i = 1;; i++, dval(d) *= 10.) {
L = (Long)(dval(d));
if (!(dval(d) -= L))
dval(&eps) *= tens[ilim-1];
for(i = 1;; i++, dval(&d) *= 10.) {
L = (Long)(dval(&d));
if (!(dval(&d) -= L))
ilim = i;
*s++ = '0' + (int)L;
if (i == ilim) {
if (dval(d) > 0.5 + dval(eps))
if (dval(&d) > 0.5 + dval(&eps))
goto bump_up;
else if (dval(d) < 0.5 - dval(eps)) {
else if (dval(&d) < 0.5 - dval(&eps)) {
while(*--s == '0');
s++;
goto ret1;
@ -426,7 +428,7 @@ dtoa
#endif
fast_failed:
s = s0;
dval(d) = dval(d2);
dval(&d) = dval(&d2);
k = k0;
ilim = ilim0;
}
@ -438,22 +440,22 @@ dtoa
ds = tens[k];
if (ndigits < 0 && ilim <= 0) {
S = mhi = 0;
if (ilim < 0 || dval(d) <= 5*ds)
if (ilim < 0 || dval(&d) <= 5*ds)
goto no_digits;
goto one_digit;
}
for(i = 1;; i++, dval(d) *= 10.) {
L = (Long)(dval(d) / ds);
dval(d) -= L*ds;
for(i = 1;; i++, dval(&d) *= 10.) {
L = (Long)(dval(&d) / ds);
dval(&d) -= L*ds;
#ifdef Check_FLT_ROUNDS
/* If FLT_ROUNDS == 2, L will usually be high by 1 */
if (dval(d) < 0) {
if (dval(&d) < 0) {
L--;
dval(d) += ds;
dval(&d) += ds;
}
#endif
*s++ = '0' + (int)L;
if (!dval(d)) {
if (!dval(&d)) {
#ifdef SET_INEXACT
inexact = 0;
#endif
@ -467,8 +469,13 @@ dtoa
case 2: goto bump_up;
}
#endif
dval(d) += dval(d);
if (dval(d) > ds || dval(d) == ds && L & 1) {
dval(&d) += dval(&d);
#ifdef ROUND_BIASED
if (dval(&d) >= ds)
#else
if (dval(&d) > ds || (dval(&d) == ds && L & 1))
#endif
{
bump_up:
while(*--s == '9')
if (s == s0) {
@ -533,9 +540,9 @@ dtoa
&& Rounding == 1
#endif
) {
if (!word1(d) && !(word0(d) & Bndry_mask)
if (!word1(&d) && !(word0(&d) & Bndry_mask)
#ifndef Sudden_Underflow
&& word0(d) & (Exp_mask & ~Exp_msk1)
&& word0(&d) & (Exp_mask & ~Exp_msk1)
#endif
) {
/* The special case */
@ -621,7 +628,7 @@ dtoa
j1 = delta->sign ? 1 : cmp(b, delta);
Bfree(delta);
#ifndef ROUND_BIASED
if (j1 == 0 && mode != 1 && !(word1(d) & 1)
if (j1 == 0 && mode != 1 && !(word1(&d) & 1)
#ifdef Honor_FLT_ROUNDS
&& Rounding >= 1
#endif
@ -638,11 +645,11 @@ dtoa
goto ret;
}
#endif
if (j < 0 || j == 0 && mode != 1
if (j < 0 || (j == 0 && mode != 1
#ifndef ROUND_BIASED
&& !(word1(d) & 1)
&& !(word1(&d) & 1)
#endif
) {
)) {
if (!b->x[0] && b->wds <= 1) {
#ifdef SET_INEXACT
inexact = 0;
@ -659,7 +666,11 @@ dtoa
if (j1 > 0) {
b = lshift(b, 1);
j1 = cmp(b, S);
if ((j1 > 0 || j1 == 0 && dig & 1)
#ifdef ROUND_BIASED
if (j1 >= 0 /*)*/
#else
if ((j1 > 0 || (j1 == 0 && dig & 1))
#endif
&& dig++ == '9')
goto round_9_up;
}
@ -719,7 +730,12 @@ dtoa
#endif
b = lshift(b, 1);
j = cmp(b, S);
if (j > 0 || j == 0 && dig & 1) {
#ifdef ROUND_BIASED
if (j >= 0)
#else
if (j > 0 || (j == 0 && dig & 1))
#endif
{
roundoff:
while(*--s == '9')
if (s == s0) {
@ -730,7 +746,9 @@ dtoa
++*s++;
}
else {
#ifdef Honor_FLT_ROUNDS
trimzeros:
#endif
while(*--s == '0');
s++;
}
@ -745,9 +763,9 @@ dtoa
#ifdef SET_INEXACT
if (inexact) {
if (!oldinexact) {
word0(d) = Exp_1 + (70 << Exp_shift);
word1(d) = 0;
dval(d) += 1.;
word0(&d) = Exp_1 + (70 << Exp_shift);
word1(&d) = 0;
dval(&d) += 1.;
}
}
else if (!oldinexact)

2
g__fmt.c
View File

@ -56,7 +56,7 @@ g__fmt(char *b, char *s, char *se, int decpt, ULong sign, size_t blen)
if (!(s0 = decimalpoint_cache)) {
s0 = localeconv()->decimal_point;
dlen = strlen(s0);
if ((decimalpoint_cache = (char*)malloc(strlen(s0) + 1))) {
if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
strcpy(decimalpoint_cache, s0);
s0 = decimalpoint_cache;
}

23
g_ddfmt.c
View File

@ -33,9 +33,9 @@ THIS SOFTWARE.
char *
#ifdef KR_headers
g_ddfmt(buf, dd, ndig, bufsize) char *buf; double *dd; int ndig; size_t bufsize;
g_ddfmt(buf, dd0, ndig, bufsize) char *buf; double *dd0; int ndig; size_t bufsize;
#else
g_ddfmt(char *buf, double *dd, int ndig, size_t bufsize)
g_ddfmt(char *buf, double *dd0, int ndig, size_t bufsize)
#endif
{
FPI fpi;
@ -43,7 +43,7 @@ g_ddfmt(char *buf, double *dd, int ndig, size_t bufsize)
ULong *L, bits0[4], *bits, *zx;
int bx, by, decpt, ex, ey, i, j, mode;
Bigint *x, *y, *z;
double ddx[2];
U *dd, ddx[2];
#ifdef Honor_FLT_ROUNDS /*{{*/
int Rounding;
#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
@ -63,7 +63,8 @@ g_ddfmt(char *buf, double *dd, int ndig, size_t bufsize)
if (bufsize < 10 || bufsize < ndig + 8)
return 0;
L = (ULong*)dd;
dd = (U*)dd0;
L = dd->L;
if ((L[_0] & 0x7ff00000L) == 0x7ff00000L) {
/* Infinity or NaN */
if (L[_0] & 0xfffff || L[_1]) {
@ -88,7 +89,7 @@ g_ddfmt(char *buf, double *dd, int ndig, size_t bufsize)
goto nanret;
goto infret;
}
if (dd[0] + dd[1] == 0.) {
if (dval(&dd[0]) + dval(&dd[1]) == 0.) {
b = buf;
#ifndef IGNORE_ZERO_SIGN
if (L[_0] & L[2+_0] & 0x80000000L)
@ -99,16 +100,16 @@ g_ddfmt(char *buf, double *dd, int ndig, size_t bufsize)
return b;
}
if ((L[_0] & 0x7ff00000L) < (L[2+_0] & 0x7ff00000L)) {
ddx[1] = dd[0];
ddx[0] = dd[1];
dval(&ddx[1]) = dval(&dd[0]);
dval(&ddx[0]) = dval(&dd[1]);
dd = ddx;
L = (ULong*)dd;
L = dd->L;
}
z = d2b(dd[0], &ex, &bx);
if (dd[1] == 0.)
z = d2b(dval(&dd[0]), &ex, &bx);
if (dval(&dd[1]) == 0.)
goto no_y;
x = z;
y = d2b(dd[1], &ey, &by);
y = d2b(dval(&dd[1]), &ey, &by);
if ( (i = ex - ey) !=0) {
if (i > 0) {
x = lshift(x, i);

2
g_dfmt.c
View File

@ -88,6 +88,8 @@ g_dfmt(char *buf, double *d, int ndig, size_t bufsize)
if (ndig <= 0)
mode = 0;
i = STRTOG_Normal;
if (sign)
i = STRTOG_Normal | STRTOG_Neg;
s = gdtoa(fpi, ex, bits, &i, mode, ndig, &decpt, &se);
return g__fmt(buf, s, se, decpt, sign, bufsize);
}

172
gdtoa.c
View File

@ -123,6 +123,7 @@ gdtoa
the returned string. If not null, *rve is set to point
to the end of the return value. If d is +-Infinity or NaN,
then *decpt is set to 9999.
be = exponent: value = (integer represented by bits) * (2 to the power of be).
mode:
0 ==> shortest string that yields d when read in
@ -157,8 +158,9 @@ gdtoa
int rdir, s2, s5, spec_case, try_quick;
Long L;
Bigint *b, *b1, *delta, *mlo, *mhi, *mhi1, *S;
double d, d2, ds, eps;
double d2, ds;
char *s, *s0;
U d, eps;
#ifndef MULTIPLE_THREADS
if (dtoa_result) {
@ -197,21 +199,21 @@ gdtoa
return nrv_alloc("0", rve, 1);
}
dval(d) = b2d(b, &i);
dval(&d) = b2d(b, &i);
i = be + bbits - 1;
word0(d) &= Frac_mask1;
word0(d) |= Exp_11;
word0(&d) &= Frac_mask1;
word0(&d) |= Exp_11;
#ifdef IBM
if ( (j = 11 - hi0bits(word0(d) & Frac_mask)) !=0)
dval(d) /= 1 << j;
if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0)
dval(&d) /= 1 << j;
#endif
/* log(x) ~=~ log(1.5) + (x-1.5)/1.5
* log10(x) = log(x) / log(10)
* ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
* log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
* log10(&d) = (i-Bias)*log(2)/log(10) + log10(d2)
*
* This suggests computing an approximation k to log10(d) by
* This suggests computing an approximation k to log10(&d) by
*
* k = (i - Bias)*0.301029995663981
* + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
@ -231,7 +233,7 @@ gdtoa
i <<= 2;
i += j;
#endif
ds = (dval(d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
ds = (dval(&d)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
/* correct assumption about exponent range */
if ((j = i) < 0)
@ -246,13 +248,13 @@ gdtoa
#ifdef IBM
j = be + bbits - 1;
if ( (j1 = j & 3) !=0)
dval(d) *= 1 << j1;
word0(d) += j << Exp_shift - 2 & Exp_mask;
dval(&d) *= 1 << j1;
word0(&d) += j << Exp_shift - 2 & Exp_mask;
#else
word0(d) += (be + bbits - 1) << Exp_shift;
word0(&d) += (be + bbits - 1) << Exp_shift;
#endif
if (k >= 0 && k <= Ten_pmax) {
if (dval(d) < tens[k])
if (dval(&d) < tens[k])
k--;
k_check = 0;
}
@ -282,11 +284,14 @@ gdtoa
mode -= 4;
try_quick = 0;
}
else if (i >= -4 - Emin || i < Emin)
try_quick = 0;
leftright = 1;
ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
/* silence erroneous "gcc -Wall" warning. */
switch(mode) {
case 0:
case 1:
ilim = ilim1 = -1;
i = (int)(nbits * .30103) + 3;
ndigits = 0;
break;
@ -328,10 +333,10 @@ gdtoa
/* Try to get by with floating-point arithmetic. */
i = 0;
d2 = dval(d);
d2 = dval(&d);
#ifdef IBM
if ( (j = 11 - hi0bits(word0(d) & Frac_mask)) !=0)
dval(d) /= 1 << j;
if ( (j = 11 - hi0bits(word0(&d) & Frac_mask)) !=0)
dval(&d) /= 1 << j;
#endif
k0 = k;
ilim0 = ilim;
@ -342,7 +347,7 @@ gdtoa
if (j & Bletch) {
/* prevent overflows */
j &= Bletch - 1;
dval(d) /= bigtens[n_bigtens-1];
dval(&d) /= bigtens[n_bigtens-1];
ieps++;
}
for(; j; j >>= 1, i++)
@ -354,30 +359,30 @@ gdtoa
else {
ds = 1.;
if ( (j1 = -k) !=0) {
dval(d) *= tens[j1 & 0xf];
dval(&d) *= tens[j1 & 0xf];
for(j = j1 >> 4; j; j >>= 1, i++)
if (j & 1) {
ieps++;
dval(d) *= bigtens[i];
dval(&d) *= bigtens[i];
}
}
}
if (k_check && dval(d) < 1. && ilim > 0) {
if (k_check && dval(&d) < 1. && ilim > 0) {
if (ilim1 <= 0)
goto fast_failed;
ilim = ilim1;
k--;
dval(d) *= 10.;
dval(&d) *= 10.;
ieps++;
}
dval(eps) = ieps*dval(d) + 7.;
word0(eps) -= (P-1)*Exp_msk1;
dval(&eps) = ieps*dval(&d) + 7.;
word0(&eps) -= (P-1)*Exp_msk1;
if (ilim == 0) {
S = mhi = 0;
dval(d) -= 5.;
if (dval(d) > dval(eps))
dval(&d) -= 5.;
if (dval(&d) > dval(&eps))
goto one_digit;
if (dval(d) < -dval(eps))
if (dval(&d) < -dval(&eps))
goto no_digits;
goto fast_failed;
}
@ -386,38 +391,38 @@ gdtoa
/* Use Steele & White method of only
* generating digits needed.
*/
dval(eps) = ds*0.5/tens[ilim-1] - dval(eps);
dval(&eps) = ds*0.5/tens[ilim-1] - dval(&eps);
for(i = 0;;) {
L = (Long)(dval(d)/ds);
dval(d) -= L*ds;
L = (Long)(dval(&d)/ds);
dval(&d) -= L*ds;
*s++ = '0' + (int)L;
if (dval(d) < dval(eps)) {
if (dval(d))
if (dval(&d) < dval(&eps)) {
if (dval(&d))
inex = STRTOG_Inexlo;
goto ret1;
}
if (ds - dval(d) < dval(eps))
if (ds - dval(&d) < dval(&eps))
goto bump_up;
if (++i >= ilim)
break;
dval(eps) *= 10.;
dval(d) *= 10.;
dval(&eps) *= 10.;
dval(&d) *= 10.;
}
}
else {
#endif
/* Generate ilim digits, then fix them up. */
dval(eps) *= tens[ilim-1];
for(i = 1;; i++, dval(d) *= 10.) {
if ( (L = (Long)(dval(d)/ds)) !=0)
dval(d) -= L*ds;
dval(&eps) *= tens[ilim-1];
for(i = 1;; i++, dval(&d) *= 10.) {
if ( (L = (Long)(dval(&d)/ds)) !=0)
dval(&d) -= L*ds;
*s++ = '0' + (int)L;
if (i == ilim) {
ds *= 0.5;
if (dval(d) > ds + dval(eps))
if (dval(&d) > ds + dval(&eps))
goto bump_up;
else if (dval(d) < ds - dval(eps)) {
if (dval(d))
else if (dval(&d) < ds - dval(&eps)) {
if (dval(&d))
inex = STRTOG_Inexlo;
goto clear_trailing0;
}
@ -429,7 +434,7 @@ gdtoa
#endif
fast_failed:
s = s0;
dval(d) = d2;
dval(&d) = d2;
k = k0;
ilim = ilim0;
}
@ -441,22 +446,22 @@ gdtoa
ds = tens[k];
if (ndigits < 0 && ilim <= 0) {
S = mhi = 0;
if (ilim < 0 || dval(d) <= 5*ds)
if (ilim < 0 || dval(&d) <= 5*ds)
goto no_digits;
goto one_digit;
}
for(i = 1;; i++, dval(d) *= 10.) {
L = dval(d) / ds;
dval(d) -= L*ds;
for(i = 1;; i++, dval(&d) *= 10.) {
L = dval(&d) / ds;
dval(&d) -= L*ds;
#ifdef Check_FLT_ROUNDS
/* If FLT_ROUNDS == 2, L will usually be high by 1 */
if (dval(d) < 0) {
if (dval(&d) < 0) {
L--;
dval(d) += ds;
dval(&d) += ds;
}
#endif
*s++ = '0' + (int)L;
if (dval(d) == 0.)
if (dval(&d) == 0.)
break;
if (i == ilim) {
if (rdir) {
@ -465,8 +470,13 @@ gdtoa
inex = STRTOG_Inexlo;
goto ret1;
}
dval(d) += dval(d);
if (dval(d) > ds || dval(d) == ds && L & 1) {
dval(&d) += dval(&d);
#ifdef ROUND_BIASED
if (dval(&d) >= ds)
#else
if (dval(&d) > ds || (dval(&d) == ds && L & 1))
#endif
{
bump_up:
inex = STRTOG_Inexhi;
while(*--s == '9')
@ -493,13 +503,15 @@ gdtoa
m5 = b5;
mhi = mlo = 0;
if (leftright) {
if (mode < 2) {
i = nbits - bbits;
if (be - i++ < fpi->emin)
/* denormal */
i = be - fpi->emin + 1;
i = nbits - bbits;
if (be - i++ < fpi->emin && mode != 3 && mode != 5) {
/* denormal */
i = be - fpi->emin + 1;
if (mode >= 2 && ilim > 0 && ilim < i)
goto small_ilim;
}
else {
else if (mode >= 2) {
small_ilim:
j = ilim - 1;
if (m5 >= j)
m5 -= j;
@ -560,28 +572,11 @@ gdtoa
* and for all and pass them and a shift to quorem, so it
* can do shifts and ors to compute the numerator for q.
*/
#ifdef Pack_32
if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) !=0)
i = 32 - i;
#else
if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) !=0)
i = 16 - i;
#endif
if (i > 4) {
i -= 4;
b2 += i;
m2 += i;
s2 += i;
}
else if (i < 4) {
i += 28;
b2 += i;
m2 += i;
s2 += i;
}
if (b2 > 0)
i = ((s5 ? hi0bits(S->x[S->wds-1]) : ULbits - 1) - s2 - 4) & kmask;
m2 += i;
if ((b2 += i) > 0)
b = lshift(b, b2);
if (s2 > 0)
if ((s2 += i) > 0)
S = lshift(S, s2);
if (k_check) {
if (cmp(b,S) < 0) {
@ -646,11 +641,11 @@ gdtoa
goto ret;
}
#endif
if (j < 0 || j == 0 && !mode
if (j < 0 || (j == 0 && !mode
#ifndef ROUND_BIASED
&& !(bits[0] & 1)
#endif
) {
)) {
if (rdir && (b->wds > 1 || b->x[0])) {
if (rdir == 2) {
inex = STRTOG_Inexlo;
@ -673,7 +668,11 @@ gdtoa
if (j1 > 0) {
b = lshift(b, 1);
j1 = cmp(b, S);
if ((j1 > 0 || j1 == 0 && dig & 1)
#ifdef ROUND_BIASED
if (j1 >= 0 /*)*/
#else
if ((j1 > 0 || (j1 == 0 && dig & 1))
#endif
&& dig++ == '9')
goto round_9_up;
inex = STRTOG_Inexhi;
@ -718,13 +717,18 @@ gdtoa
/* Round off last digit */
if (rdir) {
if (rdir == 2 || b->wds <= 1 && !b->x[0])
if (rdir == 2 || (b->wds <= 1 && !b->x[0]))
goto chopzeros;
goto roundoff;
}
b = lshift(b, 1);
j = cmp(b, S);
if (j > 0 || j == 0 && dig & 1) {
#ifdef ROUND_BIASED
if (j >= 0)
#else
if (j > 0 || (j == 0 && dig & 1))
#endif
{
roundoff:
inex = STRTOG_Inexhi;
while(*--s == '9')

2
gdtoa.h
View File

@ -36,7 +36,7 @@ THIS SOFTWARE.
#include <stddef.h> /* for size_t */
#ifndef Long
#define Long long
#define Long int
#endif
#ifndef ULong
typedef unsigned Long ULong;

49
gdtoaimp.h
View File

@ -94,7 +94,12 @@ THIS SOFTWARE.
* #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
* that use extended-precision instructions to compute rounded
* products and quotients) with IBM.
* #define ROUND_BIASED for IEEE-format with biased rounding.
* #define ROUND_BIASED for IEEE-format with biased rounding and arithmetic
* that rounds toward +Infinity.
* #define ROUND_BIASED_without_Round_Up for IEEE-format with biased
* rounding when the underlying floating-point arithmetic uses
* unbiased rounding. This prevent using ordinary floating-point
* arithmetic when the result could be computed with one rounding error.
* #define Inaccurate_Divide for IEEE-format with correctly rounded
* products but inaccurate quotients, e.g., for Intel i860.
* #define NO_LONG_LONG on machines that do not have a "long long"
@ -113,7 +118,12 @@ THIS SOFTWARE.
* #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
* if memory is available and otherwise does something you deem
* appropriate. If MALLOC is undefined, malloc will be invoked
* directly -- and assumed always to succeed.
* directly -- and assumed always to succeed. Similarly, if you
* want something other than the system's free() to be called to
* recycle memory acquired from MALLOC, #define FREE to be the
* name of the alternate routine. (FREE or free is only called in
* pathological cases, e.g., in a gdtoa call after a gdtoa return in
* mode 3 with thousands of digits requested.)
* #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
* memory allocations from a private pool of memory when possible.
* When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
@ -162,11 +172,6 @@ THIS SOFTWARE.
* #define NO_STRING_H to use private versions of memcpy.
* On some K&R systems, it may also be necessary to
* #define DECLARE_SIZE_T in this case.
* #define YES_ALIAS to permit aliasing certain double values with
* arrays of ULongs. This leads to slightly better code with
* some compilers and was always used prior to 19990916, but it
* is not strictly legal and can cause trouble with aggressively
* optimizing compilers (e.g., gcc 2.95.1 under -O2).
* #define USE_LOCALE to use the current locale's decimal_point value.
*/
@ -270,25 +275,14 @@ Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
typedef union { double d; ULong L[2]; } U;
#ifdef YES_ALIAS
#define dval(x) x
#ifdef IEEE_8087
#define word0(x) ((ULong *)&x)[1]
#define word1(x) ((ULong *)&x)[0]
#define word0(x) (x)->L[1]
#define word1(x) (x)->L[0]
#else
#define word0(x) ((ULong *)&x)[0]
#define word1(x) ((ULong *)&x)[1]
#define word0(x) (x)->L[0]
#define word1(x) (x)->L[1]
#endif
#else /* !YES_ALIAS */
#ifdef IEEE_8087
#define word0(x) ((U*)&x)->L[1]
#define word1(x) ((U*)&x)->L[0]
#else
#define word0(x) ((U*)&x)->L[0]
#define word1(x) ((U*)&x)->L[1]
#endif
#define dval(x) ((U*)&x)->d
#endif /* YES_ALIAS */
#define dval(x) (x)->d
/* The following definition of Storeinc is appropriate for MIPS processors.
* An alternative that might be better on some machines is
@ -402,6 +396,11 @@ typedef union { double d; ULong L[2]; } U;
#ifndef IEEE_Arith
#define ROUND_BIASED
#else
#ifdef ROUND_BIASED_without_Round_Up
#undef ROUND_BIASED
#define ROUND_BIASED
#endif
#endif
#ifdef RND_PRODQUOT
@ -462,7 +461,7 @@ extern double rnd_prod(double, double), rnd_quot(double, double);
#define FREE_DTOA_LOCK(n) /*nothing*/
#endif
#define Kmax 15
#define Kmax 9
struct
Bigint {
@ -575,7 +574,7 @@ extern void memcpy_D2A ANSI((void*, const void*, size_t));
extern double strtod ANSI((const char *s00, char **se));
extern Bigint *sum ANSI((Bigint*, Bigint*));
extern int trailz ANSI((Bigint*));
extern double ulp ANSI((double));
extern double ulp ANSI((U*));
#ifdef __cplusplus
}

10
gethex.c
View File

@ -57,7 +57,7 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
static unsigned char *decimalpoint_cache;
if (!(s0 = decimalpoint_cache)) {
s0 = (unsigned char*)localeconv()->decimal_point;
if ((decimalpoint_cache = (char*)malloc(strlen(s0) + 1))) {
if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
strcpy(decimalpoint_cache, s0);
s0 = decimalpoint_cache;
}
@ -199,7 +199,7 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
return STRTOG_Normal | STRTOG_Inexlo;
}
n = s1 - s0 - 1;
for(k = 0; n > (1 << kshift-2) - 1; n >>= 1)
for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
k++;
b = Balloc(k);
x = b->x;
@ -314,7 +314,7 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
break;
case FPI_Round_near:
if (lostbits & 2
&& (lostbits & 1) | x[0] & 1)
&& (lostbits | x[0]) & 1)
up = 1;
break;
case FPI_Round_up:
@ -333,8 +333,8 @@ gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
irv = STRTOG_Normal;
}
else if (b->wds > k
|| (n = nbits & kmask) !=0
&& hi0bits(x[k-1]) < 32-n) {
|| ((n = nbits & kmask) !=0
&& hi0bits(x[k-1]) < 32-n)) {
rshift(b,1);
if (++e > fpi->emax)
goto ovfl;

6
hexnan.c
View File

@ -77,7 +77,7 @@ hexnan( CONST char **sp, FPI *fpi, ULong *x0)
if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')
&& *(CONST unsigned char*)(s+3) > ' ')
s += 2;
while(c = *(CONST unsigned char*)++s) {
while((c = *(CONST unsigned char*)++s)) {
if (!(h = hexdig[c])) {
if (c <= ' ') {
if (hd0 < havedig) {
@ -109,7 +109,7 @@ hexnan( CONST char **sp, FPI *fpi, ULong *x0)
*sp = s + 1;
break;
}
} while(c = *++s);
} while((c = *++s));
#endif
return STRTOG_NaN;
}
@ -120,7 +120,7 @@ hexnan( CONST char **sp, FPI *fpi, ULong *x0)
i = 1;
*--x = 0;
}
*x = (*x << 4) | h & 0xf;
*x = (*x << 4) | (h & 0xf);
}
if (!havedig)
return STRTOG_NaN;

41
makefile
View File

@ -30,6 +30,8 @@ CFLAGS = -g
.c.o:
$(CC) -c $(CFLAGS) $*.c
# invoke "make Printf" to add printf.o to gdtoa.a (if desired)
all: arith.h gd_qnan.h gdtoa.a
arith.h: arithchk.c
@ -42,27 +44,36 @@ gd_qnan.h: arith.h qnan.c
./a.out >gd_qnan.h
rm -f a.out qnan.o
gdtoa.a: dmisc.c dtoa.c g_Qfmt.c g__fmt.c g_ddfmt.c g_dfmt.c g_ffmt.c\
g_xLfmt.c g_xfmt.c gdtoa.c gethex.c gmisc.c hd_init.c hexnan.c\
misc.c smisc.c strtoIQ.c strtoId.c strtoIdd.c strtoIf.c strtoIg.c\
strtoIx.c strtoIxL.c strtod.c strtodI.c strtodg.c strtof.c strtopQ.c\
strtopd.c strtopdd.c strtopf.c strtopx.c strtopxL.c strtorQ.c\
strtord.c strtordd.c strtorf.c strtorx.c strtorxL.c sum.c ulp.c
gdtoa.a: dmisc.c dtoa.c g_Qfmt.c g__fmt.c g_ddfmt.c g_dfmt.c\
g_ffmt.c g_xLfmt.c g_xfmt.c gdtoa.c gethex.c gmisc.c hd_init.c\
hexnan.c misc.c smisc.c strtoIQ.c strtoId.c strtoIdd.c\
strtoIf.c strtoIg.c strtoIx.c strtoIxL.c strtod.c strtodI.c\
strtodg.c strtof.c strtopQ.c strtopd.c strtopdd.c strtopf.c\
strtopx.c strtopxL.c strtorQ.c strtord.c strtordd.c strtorf.c\
strtorx.c strtorxL.c sum.c ulp.c
$(CC) -c $(CFLAGS) $?
x=`echo $? | sed 's/\.c/.o/g'` && ar ruv gdtoa.a $$x && rm $$x
ranlib gdtoa.a || true
Printf: all printf.c
$(CC) -c $(CFLAGS) printf.c
ar ruv gdtoa.a printf.o
rm printf.o
touch Printf
# If your system lacks ranlib, you do not need it.
xs0 = README arithchk.c dmisc.c dtoa.c g_Qfmt.c g__fmt.c g_ddfmt.c g_dfmt.c\
g_ffmt.c g_xLfmt.c g_xfmt.c gdtoa.c gdtoa.h gdtoa_fltrnds.h gdtoaimp.h\
gethex.c gmisc.c hd_init.c hexnan.c makefile misc.c qnan.c smisc.c\
strtoIQ.c strtoId.c strtoIdd.c strtoIf.c strtoIg.c strtoIx.c strtoIxL.c\
strtod.c strtodI.c strtodg.c strtodnrp.c strtof.c strtopQ.c strtopd.c\
strtopdd.c strtopf.c strtopx.c strtopxL.c strtorQ.c strtord.c strtordd.c\
strtorf.c strtorx.c strtorxL.c sum.c ulp.c
xs0 = README arithchk.c dmisc.c dtoa.c g_Qfmt.c g__fmt.c g_ddfmt.c\
g_dfmt.c g_ffmt.c g_xLfmt.c g_xfmt.c gdtoa.c gdtoa.h\
gdtoa_fltrnds.h gdtoaimp.h gethex.c gmisc.c hd_init.c hexnan.c\
makefile misc.c printf.c printf.c0 qnan.c smisc.c stdio1.h\
strtoIQ.c strtoId.c strtoIdd.c strtoIf.c strtoIg.c strtoIx.c\
strtoIxL.c strtod.c strtodI.c strtodg.c strtodnrp.c strtof.c\
strtopQ.c strtopd.c strtopdd.c strtopf.c strtopx.c strtopxL.c\
strtorQ.c strtord.c strtordd.c strtorf.c strtorx.c strtorxL.c\
sum.c ulp.c
# "make xsum.out" to check for transmission errors; source for xsum is
# "make -r xsum.out" to check for transmission errors; source for xsum is
# netlib's "xsum.c from f2c", e.g.,
# ftp://netlib.bell-labs.com/netlib/f2c/xsum.c.gz
@ -71,4 +82,4 @@ xsum.out: xsum0.out $(xs0)
cmp xsum0.out xsum1.out && mv xsum1.out xsum.out || diff xsum[01].out
clean:
rm -f arith.h gd_qnan.h *.[ao] xsum.out xsum1.out
rm -f arith.h gd_qnan.h *.[ao] Printf xsum.out xsum1.out

82
misc.c
View File

@ -55,7 +55,9 @@ Balloc
#endif
ACQUIRE_DTOA_LOCK(0);
if ( (rv = freelist[k]) !=0) {
/* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
/* but this case seems very unlikely. */
if (k <= Kmax && (rv = freelist[k]) !=0) {
freelist[k] = rv->next;
}
else {
@ -65,7 +67,7 @@ Balloc
#else
len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
/sizeof(double);
if (pmem_next - private_mem + len <= PRIVATE_mem) {
if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
rv = (Bigint*)pmem_next;
pmem_next += len;
}
@ -89,10 +91,18 @@ Bfree
#endif
{
if (v) {
ACQUIRE_DTOA_LOCK(0);
v->next = freelist[v->k];
freelist[v->k] = v;
FREE_DTOA_LOCK(0);
if (v->k > Kmax)
#ifdef FREE
FREE((void*)v);
#else
free((void*)v);
#endif
else {
ACQUIRE_DTOA_LOCK(0);
v->next = freelist[v->k];
freelist[v->k] = v;
FREE_DTOA_LOCK(0);
}
}
}
@ -104,8 +114,8 @@ lo0bits
(ULong *y)
#endif
{
register int k;
register ULong x = *y;
int k;
ULong x = *y;
if (x & 7) {
if (x & 1)
@ -204,12 +214,12 @@ multadd
int
hi0bits_D2A
#ifdef KR_headers
(x) register ULong x;
(x) ULong x;
#else
(register ULong x)
(ULong x)
#endif
{
register int k = 0;
int k = 0;
if (!(x & 0xffff0000)) {
k = 16;
@ -612,12 +622,12 @@ b2d
{
ULong *xa, *xa0, w, y, z;
int k;
double d;
U d;
#ifdef VAX
ULong d0, d1;
#else
#define d0 word0(d)
#define d1 word1(d)
#define d0 word0(&d)
#define d1 word1(&d)
#endif
xa0 = a->x;
@ -630,16 +640,16 @@ b2d
*e = 32 - k;
#ifdef Pack_32
if (k < Ebits) {
d0 = Exp_1 | y >> Ebits - k;
d0 = Exp_1 | y >> (Ebits - k);
w = xa > xa0 ? *--xa : 0;
d1 = y << (32-Ebits) + k | w >> Ebits - k;
d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
goto ret_d;
}
z = xa > xa0 ? *--xa : 0;
if (k -= Ebits) {
d0 = Exp_1 | y << k | z >> 32 - k;
d0 = Exp_1 | y << k | z >> (32 - k);
y = xa > xa0 ? *--xa : 0;
d1 = z << k | y >> 32 - k;
d1 = z << k | y >> (32 - k);
}
else {
d0 = Exp_1 | y;
@ -663,10 +673,10 @@ b2d
#endif
ret_d:
#ifdef VAX
word0(d) = d0 >> 16 | d0 << 16;
word1(d) = d1 >> 16 | d1 << 16;
word0(&d) = d0 >> 16 | d0 << 16;
word1(&d) = d1 >> 16 | d1 << 16;
#endif
return dval(d);
return dval(&d);
}
#undef d0
#undef d1
@ -674,12 +684,13 @@ b2d
Bigint *
d2b
#ifdef KR_headers
(d, e, bits) double d; int *e, *bits;
(dd, e, bits) double dd; int *e, *bits;
#else
(double d, int *e, int *bits)
(double dd, int *e, int *bits)
#endif
{
Bigint *b;
U d;
#ifndef Sudden_Underflow
int i;
#endif
@ -687,11 +698,14 @@ d2b
ULong *x, y, z;
#ifdef VAX
ULong d0, d1;
d0 = word0(d) >> 16 | word0(d) << 16;
d1 = word1(d) >> 16 | word1(d) << 16;
#else
#define d0 word0(d)
#define d1 word1(d)
#define d0 word0(&d)
#define d1 word1(&d)
#endif
d.d = dd;
#ifdef VAX
d0 = word0(&d) >> 16 | word0(&d) << 16;
d1 = word1(&d) >> 16 | word1(&d) << 16;
#endif
#ifdef Pack_32
@ -715,7 +729,7 @@ d2b
#ifdef Pack_32
if ( (y = d1) !=0) {
if ( (k = lo0bits(&y)) !=0) {
x[0] = y | z << 32 - k;
x[0] = y | z << (32 - k);
z >>= k;
}
else
@ -726,10 +740,6 @@ d2b
b->wds = (x[1] = z) !=0 ? 2 : 1;
}
else {
#ifdef DEBUG
if (!z)
Bug("Zero passed to d2b");
#endif
k = lo0bits(&z);
x[0] = z;
#ifndef Sudden_Underflow
@ -788,7 +798,7 @@ d2b
#endif
#ifdef IBM
*e = (de - Bias - (P-1) << 2) + k;
*bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
*bits = 4*P + 8 - k - hi0bits(word0(&d) & Frac_mask);
#else
*e = de - Bias - (P-1) + k;
*bits = P - k;
@ -841,7 +851,7 @@ strcp_D2A(a, b) char *a; char *b;
strcp_D2A(char *a, CONST char *b)
#endif
{
while(*a = *b++)
while((*a = *b++))
a++;
return a;
}
@ -855,8 +865,8 @@ memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
memcpy_D2A(void *a1, void *b1, size_t len)
#endif
{
register char *a = (char*)a1, *ae = a + len;
register char *b = (char*)b1, *a0 = a;
char *a = (char*)a1, *ae = a + len;
char *b = (char*)b1, *a0 = a;
while(a < ae)
*a++ = *b++;
return a0;

8
printf.c Normal file
View File

@ -0,0 +1,8 @@
#ifdef __sun
#define Use_GDTOA_Qtype
#else
#if defined(__i386) || defined(__x86_64)
#define Use_GDTOA_for_i386_long_double
#endif
#endif
#include "printf.c0"

1669
printf.c0 Normal file
View File

File diff suppressed because it is too large Load Diff

20
smisc.c
View File

@ -77,33 +77,33 @@ ratio
(Bigint *a, Bigint *b)
#endif
{
double da, db;
U da, db;
int k, ka, kb;
dval(da) = b2d(a, &ka);
dval(db) = b2d(b, &kb);
dval(&da) = b2d(a, &ka);
dval(&db) = b2d(b, &kb);
k = ka - kb + ULbits*(a->wds - b->wds);
#ifdef IBM
if (k > 0) {
word0(da) += (k >> 2)*Exp_msk1;
word0(&da) += (k >> 2)*Exp_msk1;
if (k &= 3)
dval(da) *= 1 << k;
dval(&da) *= 1 << k;
}
else {
k = -k;
word0(db) += (k >> 2)*Exp_msk1;
word0(&db) += (k >> 2)*Exp_msk1;
if (k &= 3)
dval(db) *= 1 << k;
dval(&db) *= 1 << k;
}
#else
if (k > 0)
word0(da) += k*Exp_msk1;
word0(&da) += k*Exp_msk1;
else {
k = -k;
word0(db) += k*Exp_msk1;
word0(&db) += k*Exp_msk1;
}
#endif
return dval(da) / dval(db);
return dval(&da) / dval(&db);
}
#ifdef INFNAN_CHECK

103
stdio1.h Normal file
View File

@ -0,0 +1,103 @@
/****************************************************************
Copyright (C) 1997-1999 Lucent Technologies
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation, and that the name of Lucent or any of its entities
not be used in advertising or publicity pertaining to
distribution of the software without specific, written prior
permission.
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
THIS SOFTWARE.
****************************************************************/
/* stdio1.h -- for using Printf, Fprintf, Sprintf while
* retaining the system-supplied printf, fprintf, sprintf.
*/
#ifndef STDIO1_H_included
#define STDIO1_H_included
#ifndef STDIO_H_included /* allow suppressing stdio.h */
#include <stdio.h> /* in case it's already included, */
#endif /* e.g., by cplex.h */
#ifdef KR_headers
#ifndef _SIZE_T
#define _SIZE_T
typedef unsigned int size_t;
#endif
#define ANSI(x) ()
#include "varargs.h"
#ifndef Char
#define Char char
#endif
#else
#define ANSI(x) x
#include "stdarg.h"
#ifndef Char
#define Char void
#endif
#endif
#ifndef NO_STDIO1
#ifdef __cplusplus
extern "C" {
#endif
extern int Fprintf ANSI((FILE*, const char*, ...));
extern int Printf ANSI((const char*, ...));
extern int Sprintf ANSI((char*, const char*, ...));
extern int Snprintf ANSI((char*, size_t, const char*, ...));
extern void Perror ANSI((const char*));
extern int Vfprintf ANSI((FILE*, const char*, va_list));
extern int Vsprintf ANSI((char*, const char*, va_list));
extern int Vsnprintf ANSI((char*, size_t, const char*, va_list));
#ifdef PF_BUF
extern FILE *stderr_ASL;
extern void (*pfbuf_print_ASL) ANSI((char*));
extern char *pfbuf_ASL;
extern void fflush_ASL ANSI((FILE*));
#ifdef fflush
#define old_fflush_ASL fflush
#undef fflush
#endif
#define fflush fflush_ASL
#endif
#ifdef __cplusplus
}
#endif
#undef printf
#undef fprintf
#undef sprintf
#undef perror
#undef vfprintf
#undef vsprintf
#define printf Printf
#define fprintf Fprintf
#undef snprintf /* for MacOSX */
#undef vsnprintf /* for MacOSX */
#define snprintf Snprintf
#define sprintf Sprintf
#define perror Perror
#define vfprintf Vfprintf
#define vsnprintf Vsnprintf
#define vsprintf Vsprintf
#endif /* NO_STDIO1 */
#endif /* STDIO1_H_included */

2
strtoIg.c
View File

@ -74,7 +74,7 @@ strtoIg(CONST char *s00, char **se, FPI *fpi, Long *exp, Bigint **B, int *rvp)
goto swapcheck;
}
if (b1->wds > nw
|| nb1 && b1->x[nw1] & 1L << nb1) {
|| (nb1 && b1->x[nw1] & 1L << nb1)) {
if (++e1 > fpi->emax)
rv1 = STRTOG_Infinite | STRTOG_Inexhi;
rshift(b1, 1);

446
strtod.c
View File

@ -57,6 +57,28 @@ static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
#define Rounding Flt_Rounds
#endif
#ifdef Avoid_Underflow /*{*/
static double
sulp
#ifdef KR_headers
(x, scale) U *x; int scale;
#else
(U *x, int scale)
#endif
{
U u;
double rv;
int i;
rv = ulp(x);
if (!scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0)
return rv; /* Is there an example where i <= 0 ? */
word0(&u) = Exp_1 + (i << Exp_shift);
word1(&u) = 0;
return rv * u.d;
}
#endif /*}*/
double
strtod
#ifdef KR_headers
@ -71,10 +93,14 @@ strtod
int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
CONST char *s, *s0, *s1;
double aadj, aadj1, adj, rv, rv0;
double aadj;
Long L;
U adj, aadj1, rv, rv0;
ULong y, z;
Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
#ifdef Avoid_Underflow
ULong Lsb, Lsb1;
#endif
#ifdef SET_INEXACT
int inexact, oldinexact;
#endif
@ -88,7 +114,7 @@ strtod
static int dplen;
if (!(s0 = decimalpoint_cache)) {
s0 = localeconv()->decimal_point;
if ((decimalpoint_cache = (char*)malloc(strlen(s0) + 1))) {
if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
strcpy(decimalpoint_cache, s0);
s0 = decimalpoint_cache;
}
@ -115,7 +141,7 @@ strtod
#endif /*}*/
sign = nz0 = nz = decpt = 0;
dval(rv) = 0.;
dval(&rv) = 0.;
for(s = s00;;s++) switch(*s) {
case '-':
sign = 1;
@ -147,20 +173,12 @@ strtod
case 'x':
case 'X':
{
#if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD) /*{{*/
#ifdef Honor_FLT_ROUNDS
FPI fpi1 = fpi;
#ifdef Honor_FLT_ROUNDS /*{{*/
fpi1.rounding = Rounding;
#else /*}{*/
switch(fegetround()) {
case FE_TOWARDZERO: fpi1.rounding = 0; break;
case FE_UPWARD: fpi1.rounding = 2; break;
case FE_DOWNWARD: fpi1.rounding = 3;
}
#endif /*}}*/
#else /*}{*/
#else
#define fpi1 fpi
#endif /*}}*/
#endif
switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) {
case STRTOG_NoNumber:
s = s00;
@ -285,8 +303,8 @@ strtod
--s;
if (!match(&s,"inity"))
++s;
word0(rv) = 0x7ff00000;
word1(rv) = 0;
word0(&rv) = 0x7ff00000;
word1(&rv) = 0;
goto ret;
}
break;
@ -297,13 +315,13 @@ strtod
if (*s == '(' /*)*/
&& hexnan(&s, &fpinan, bits)
== STRTOG_NaNbits) {
word0(rv) = 0x7ff00000 | bits[1];
word1(rv) = bits[0];
word0(&rv) = 0x7ff00000 | bits[1];
word1(&rv) = bits[0];
}
else {
#endif
word0(rv) = NAN_WORD0;
word1(rv) = NAN_WORD1;
word0(&rv) = NAN_WORD0;
word1(&rv) = NAN_WORD1;
#ifndef No_Hex_NaN
}
#endif
@ -327,13 +345,13 @@ strtod
if (!nd0)
nd0 = nd;
k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
dval(rv) = y;
dval(&rv) = y;
if (k > 9) {
#ifdef SET_INEXACT
if (k > DBL_DIG)
oldinexact = get_inexact();
#endif
dval(rv) = tens[k - 9] * dval(rv) + z;
dval(&rv) = tens[k - 9] * dval(&rv) + z;
}
bd0 = 0;
if (nd <= DBL_DIG
@ -345,6 +363,7 @@ strtod
) {
if (!e)
goto ret;
#ifndef ROUND_BIASED_without_Round_Up
if (e > 0) {
if (e <= Ten_pmax) {
#ifdef VAX
@ -353,11 +372,11 @@ strtod
#ifdef Honor_FLT_ROUNDS
/* round correctly FLT_ROUNDS = 2 or 3 */
if (sign) {
rv = -rv;
rv.d = -rv.d;
sign = 0;
}
#endif
/* rv = */ rounded_product(dval(rv), tens[e]);
/* rv = */ rounded_product(dval(&rv), tens[e]);
goto ret;
#endif
}
@ -369,25 +388,25 @@ strtod
#ifdef Honor_FLT_ROUNDS
/* round correctly FLT_ROUNDS = 2 or 3 */
if (sign) {
rv = -rv;
rv.d = -rv.d;
sign = 0;
}
#endif
e -= i;
dval(rv) *= tens[i];
dval(&rv) *= tens[i];
#ifdef VAX
/* VAX exponent range is so narrow we must
* worry about overflow here...
*/
vax_ovfl_check:
word0(rv) -= P*Exp_msk1;
/* rv = */ rounded_product(dval(rv), tens[e]);
if ((word0(rv) & Exp_mask)
word0(&rv) -= P*Exp_msk1;
/* rv = */ rounded_product(dval(&rv), tens[e]);
if ((word0(&rv) & Exp_mask)
> Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
goto ovfl;
word0(rv) += P*Exp_msk1;
word0(&rv) += P*Exp_msk1;
#else
/* rv = */ rounded_product(dval(rv), tens[e]);
/* rv = */ rounded_product(dval(&rv), tens[e]);
#endif
goto ret;
}
@ -397,14 +416,15 @@ strtod
#ifdef Honor_FLT_ROUNDS
/* round correctly FLT_ROUNDS = 2 or 3 */
if (sign) {
rv = -rv;
rv.d = -rv.d;
sign = 0;
}
#endif
/* rv = */ rounded_quotient(dval(rv), tens[-e]);
/* rv = */ rounded_quotient(dval(&rv), tens[-e]);
goto ret;
}
#endif
#endif /* ROUND_BIASED_without_Round_Up */
}
e1 += nd - k;
@ -432,67 +452,73 @@ strtod
if (e1 > 0) {
if ( (i = e1 & 15) !=0)
dval(rv) *= tens[i];
dval(&rv) *= tens[i];
if (e1 &= ~15) {
if (e1 > DBL_MAX_10_EXP) {
ovfl:
#ifndef NO_ERRNO
errno = ERANGE;
#endif
/* Can't trust HUGE_VAL */
#ifdef IEEE_Arith
#ifdef Honor_FLT_ROUNDS
switch(Rounding) {
case 0: /* toward 0 */
case 3: /* toward -infinity */
word0(rv) = Big0;
word1(rv) = Big1;
word0(&rv) = Big0;
word1(&rv) = Big1;
break;
default:
word0(rv) = Exp_mask;
word1(rv) = 0;
word0(&rv) = Exp_mask;
word1(&rv) = 0;
}
#else /*Honor_FLT_ROUNDS*/
word0(rv) = Exp_mask;
word1(rv) = 0;
word0(&rv) = Exp_mask;
word1(&rv) = 0;
#endif /*Honor_FLT_ROUNDS*/
#ifdef SET_INEXACT
/* set overflow bit */
dval(rv0) = 1e300;
dval(rv0) *= dval(rv0);
dval(&rv0) = 1e300;
dval(&rv0) *= dval(&rv0);
#endif
#else /*IEEE_Arith*/
word0(rv) = Big0;
word1(rv) = Big1;
word0(&rv) = Big0;
word1(&rv) = Big1;
#endif /*IEEE_Arith*/
if (bd0)
goto retfree;
range_err:
if (bd0) {
Bfree(bb);
Bfree(bd);
Bfree(bs);
Bfree(bd0);
Bfree(delta);
}
#ifndef NO_ERRNO
errno = ERANGE;
#endif
goto ret;
}
e1 >>= 4;
for(j = 0; e1 > 1; j++, e1 >>= 1)
if (e1 & 1)
dval(rv) *= bigtens[j];
dval(&rv) *= bigtens[j];
/* The last multiplication could overflow. */
word0(rv) -= P*Exp_msk1;
dval(rv) *= bigtens[j];
if ((z = word0(rv) & Exp_mask)
word0(&rv) -= P*Exp_msk1;
dval(&rv) *= bigtens[j];
if ((z = word0(&rv) & Exp_mask)
> Exp_msk1*(DBL_MAX_EXP+Bias-P))
goto ovfl;
if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
/* set to largest number */
/* (Can't trust DBL_MAX) */
word0(rv) = Big0;
word1(rv) = Big1;
word0(&rv) = Big0;
word1(&rv) = Big1;
}
else
word0(rv) += P*Exp_msk1;
word0(&rv) += P*Exp_msk1;
}
}
else if (e1 < 0) {
e1 = -e1;
if ( (i = e1 & 15) !=0)
dval(rv) /= tens[i];
dval(&rv) /= tens[i];
if (e1 >>= 4) {
if (e1 >= 1 << n_bigtens)
goto undfl;
@ -501,44 +527,39 @@ strtod
scale = 2*P;
for(j = 0; e1 > 0; j++, e1 >>= 1)
if (e1 & 1)
dval(rv) *= tinytens[j];
if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
dval(&rv) *= tinytens[j];
if (scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
>> Exp_shift)) > 0) {
/* scaled rv is denormal; zap j low bits */
if (j >= 32) {
word1(rv) = 0;
word1(&rv) = 0;
if (j >= 53)
word0(rv) = (P+2)*Exp_msk1;
word0(&rv) = (P+2)*Exp_msk1;
else
word0(rv) &= 0xffffffff << j-32;
word0(&rv) &= 0xffffffff << (j-32);
}
else
word1(rv) &= 0xffffffff << j;
word1(&rv) &= 0xffffffff << j;
}
#else
for(j = 0; e1 > 1; j++, e1 >>= 1)
if (e1 & 1)
dval(rv) *= tinytens[j];
dval(&rv) *= tinytens[j];
/* The last multiplication could underflow. */
dval(rv0) = dval(rv);
dval(rv) *= tinytens[j];
if (!dval(rv)) {
dval(rv) = 2.*dval(rv0);
dval(rv) *= tinytens[j];
dval(&rv0) = dval(&rv);
dval(&rv) *= tinytens[j];
if (!dval(&rv)) {
dval(&rv) = 2.*dval(&rv0);
dval(&rv) *= tinytens[j];
#endif
if (!dval(rv)) {
if (!dval(&rv)) {
undfl:
dval(rv) = 0.;
#ifndef NO_ERRNO
errno = ERANGE;
#endif
if (bd0)
goto retfree;
goto ret;
dval(&rv) = 0.;
goto range_err;
}
#ifndef Avoid_Underflow
word0(rv) = Tiny0;
word1(rv) = Tiny1;
word0(&rv) = Tiny0;
word1(&rv) = Tiny1;
/* The refinement below will clean
* this approximation up.
*/
@ -556,7 +577,7 @@ strtod
for(;;) {
bd = Balloc(bd0->k);
Bcopy(bd, bd0);
bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
bb = d2b(dval(&rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
bs = i2b(1);
if (e >= 0) {
@ -577,12 +598,19 @@ strtod
bs2++;
#endif
#ifdef Avoid_Underflow
Lsb = LSB;
Lsb1 = 0;
j = bbe - scale;
i = j + bbbits - 1; /* logb(rv) */
if (i < Emin) /* denormal */
j += P - Emin;
else
j = P + 1 - bbbits;
j = P + 1 - bbbits;
if (i < Emin) { /* denormal */
i = Emin - i;
j -= i;
if (i < 32)
Lsb <<= i;
else
Lsb1 = Lsb << (i-32);
}
#else /*Avoid_Underflow*/
#ifdef Sudden_Underflow
#ifdef IBM
@ -592,7 +620,7 @@ strtod
#endif
#else /*Sudden_Underflow*/
j = bbe;
i = j + bbbits - 1; /* logb(rv) */
i = j + bbbits - 1; /* logb(&rv) */
if (i < Emin) /* denormal */
j += P - Emin;
else
@ -643,15 +671,15 @@ strtod
}
if (Rounding) {
if (dsign) {
adj = 1.;
dval(&adj) = 1.;
goto apply_adj;
}
}
else if (!dsign) {
adj = -1.;
if (!word1(rv)
&& !(word0(rv) & Frac_mask)) {
y = word0(rv) & Exp_mask;
dval(&adj) = -1.;
if (!word1(&rv)
&& !(word0(&rv) & Frac_mask)) {
y = word0(&rv) & Exp_mask;
#ifdef Avoid_Underflow
if (!scale || y > 2*P*Exp_msk1)
#else
@ -660,66 +688,66 @@ strtod
{
delta = lshift(delta,Log2P);
if (cmp(delta, bs) <= 0)
adj = -0.5;
dval(&adj) = -0.5;
}
}
apply_adj:
#ifdef Avoid_Underflow
if (scale && (y = word0(rv) & Exp_mask)
if (scale && (y = word0(&rv) & Exp_mask)
<= 2*P*Exp_msk1)
word0(adj) += (2*P+1)*Exp_msk1 - y;
word0(&adj) += (2*P+1)*Exp_msk1 - y;
#else
#ifdef Sudden_Underflow
if ((word0(rv) & Exp_mask) <=
if ((word0(&rv) & Exp_mask) <=
P*Exp_msk1) {
word0(rv) += P*Exp_msk1;
dval(rv) += adj*ulp(dval(rv));
word0(rv) -= P*Exp_msk1;
word0(&rv) += P*Exp_msk1;
dval(&rv) += adj*ulp(&rv);
word0(&rv) -= P*Exp_msk1;
}
else
#endif /*Sudden_Underflow*/
#endif /*Avoid_Underflow*/
dval(rv) += adj*ulp(dval(rv));
dval(&rv) += adj.d*ulp(&rv);
}
break;
}
adj = ratio(delta, bs);
if (adj < 1.)
adj = 1.;
if (adj <= 0x7ffffffe) {
/* adj = Rounding ? ceil(adj) : floor(adj); */
y = adj;
if (y != adj) {
dval(&adj) = ratio(delta, bs);
if (adj.d < 1.)
dval(&adj) = 1.;
if (adj.d <= 0x7ffffffe) {
/* dval(&adj) = Rounding ? ceil(&adj) : floor(&adj); */
y = adj.d;
if (y != adj.d) {
if (!((Rounding>>1) ^ dsign))
y++;
adj = y;
dval(&adj) = y;
}
}
#ifdef Avoid_Underflow
if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
word0(adj) += (2*P+1)*Exp_msk1 - y;
if (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
word0(&adj) += (2*P+1)*Exp_msk1 - y;
#else
#ifdef Sudden_Underflow
if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
word0(rv) += P*Exp_msk1;
adj *= ulp(dval(rv));
if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
word0(&rv) += P*Exp_msk1;
dval(&adj) *= ulp(&rv);
if (dsign)
dval(rv) += adj;
dval(&rv) += adj;
else
dval(rv) -= adj;
word0(rv) -= P*Exp_msk1;
dval(&rv) -= adj;
word0(&rv) -= P*Exp_msk1;
goto cont;
}
#endif /*Sudden_Underflow*/
#endif /*Avoid_Underflow*/
adj *= ulp(dval(rv));
dval(&adj) *= ulp(&rv);
if (dsign) {
if (word0(rv) == Big0 && word1(rv) == Big1)
if (word0(&rv) == Big0 && word1(&rv) == Big1)
goto ovfl;
dval(rv) += adj;
dval(&rv) += adj.d;
}
else
dval(rv) -= adj;
dval(&rv) -= adj.d;
goto cont;
}
#endif /*Honor_FLT_ROUNDS*/
@ -728,12 +756,12 @@ strtod
/* Error is less than half an ulp -- check for
* special case of mantissa a power of two.
*/
if (dsign || word1(rv) || word0(rv) & Bndry_mask
if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
#ifdef IEEE_Arith
#ifdef Avoid_Underflow
|| (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
|| (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
#else
|| (word0(rv) & Exp_mask) <= Exp_msk1
|| (word0(&rv) & Exp_mask) <= Exp_msk1
#endif
#endif
) {
@ -758,32 +786,34 @@ strtod
if (i == 0) {
/* exactly half-way between */
if (dsign) {
if ((word0(rv) & Bndry_mask1) == Bndry_mask1
&& word1(rv) == (
if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
&& word1(&rv) == (
#ifdef Avoid_Underflow
(scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
(scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
#endif
0xffffffff)) {
/*boundary case -- increment exponent*/
word0(rv) = (word0(rv) & Exp_mask)
if (word0(&rv) == Big0 && word1(&rv) == Big1)
goto ovfl;
word0(&rv) = (word0(&rv) & Exp_mask)
+ Exp_msk1
#ifdef IBM
| Exp_msk1 >> 4
#endif
;
word1(rv) = 0;
word1(&rv) = 0;
#ifdef Avoid_Underflow
dsign = 0;
#endif
break;
}
}
else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
drop_down:
/* boundary case -- decrement exponent */
#ifdef Sudden_Underflow /*{{*/
L = word0(rv) & Exp_mask;
L = word0(&rv) & Exp_mask;
#ifdef IBM
if (L < Exp_msk1)
#else
@ -798,7 +828,7 @@ strtod
#else /*Sudden_Underflow}{*/
#ifdef Avoid_Underflow
if (scale) {
L = word0(rv) & Exp_mask;
L = word0(&rv) & Exp_mask;
if (L <= (2*P+1)*Exp_msk1) {
if (L > (P+2)*Exp_msk1)
/* round even ==> */
@ -809,10 +839,10 @@ strtod
}
}
#endif /*Avoid_Underflow*/
L = (word0(rv) & Exp_mask) - Exp_msk1;
L = (word0(&rv) & Exp_mask) - Exp_msk1;
#endif /*Sudden_Underflow}}*/
word0(rv) = L | Bndry_mask1;
word1(rv) = 0xffffffff;
word0(&rv) = L | Bndry_mask1;
word1(&rv) = 0xffffffff;
#ifdef IBM
goto cont;
#else
@ -820,16 +850,33 @@ strtod
#endif
}
#ifndef ROUND_BIASED
if (!(word1(rv) & LSB))
#ifdef Avoid_Underflow
if (Lsb1) {
if (!(word0(&rv) & Lsb1))
break;
}
else if (!(word1(&rv) & Lsb))
break;
#else
if (!(word1(&rv) & LSB))
break;
#endif
#endif
if (dsign)
dval(rv) += ulp(dval(rv));
#ifdef Avoid_Underflow
dval(&rv) += sulp(&rv, scale);
#else
dval(&rv) += ulp(&rv);
#endif
#ifndef ROUND_BIASED
else {
dval(rv) -= ulp(dval(rv));
#ifdef Avoid_Underflow
dval(&rv) -= sulp(&rv, scale);
#else
dval(&rv) -= ulp(&rv);
#endif
#ifndef Sudden_Underflow
if (!dval(rv))
if (!dval(&rv))
goto undfl;
#endif
}
@ -841,14 +888,14 @@ strtod
}
if ((aadj = ratio(delta, bs)) <= 2.) {
if (dsign)
aadj = aadj1 = 1.;
else if (word1(rv) || word0(rv) & Bndry_mask) {
aadj = dval(&aadj1) = 1.;
else if (word1(&rv) || word0(&rv) & Bndry_mask) {
#ifndef Sudden_Underflow
if (word1(rv) == Tiny1 && !word0(rv))
if (word1(&rv) == Tiny1 && !word0(&rv))
goto undfl;
#endif
aadj = 1.;
aadj1 = -1.;
dval(&aadj1) = -1.;
}
else {
/* special case -- power of FLT_RADIX to be */
@ -858,45 +905,45 @@ strtod
aadj = 1./FLT_RADIX;
else
aadj *= 0.5;
aadj1 = -aadj;
dval(&aadj1) = -aadj;
}
}
else {
aadj *= 0.5;
aadj1 = dsign ? aadj : -aadj;
dval(&aadj1) = dsign ? aadj : -aadj;
#ifdef Check_FLT_ROUNDS
switch(Rounding) {
case 2: /* towards +infinity */
aadj1 -= 0.5;
dval(&aadj1) -= 0.5;
break;
case 0: /* towards 0 */
case 3: /* towards -infinity */
aadj1 += 0.5;
dval(&aadj1) += 0.5;
}
#else
if (Flt_Rounds == 0)
aadj1 += 0.5;
dval(&aadj1) += 0.5;
#endif /*Check_FLT_ROUNDS*/
}
y = word0(rv) & Exp_mask;
y = word0(&rv) & Exp_mask;
/* Check for overflow */
if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
dval(rv0) = dval(rv);
word0(rv) -= P*Exp_msk1;
adj = aadj1 * ulp(dval(rv));
dval(rv) += adj;
if ((word0(rv) & Exp_mask) >=
dval(&rv0) = dval(&rv);
word0(&rv) -= P*Exp_msk1;
dval(&adj) = dval(&aadj1) * ulp(&rv);
dval(&rv) += dval(&adj);
if ((word0(&rv) & Exp_mask) >=
Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
if (word0(rv0) == Big0 && word1(rv0) == Big1)
if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
goto ovfl;
word0(rv) = Big0;
word1(rv) = Big1;
word0(&rv) = Big0;
word1(&rv) = Big1;
goto cont;
}
else
word0(rv) += P*Exp_msk1;
word0(&rv) += P*Exp_msk1;
}
else {
#ifdef Avoid_Underflow
@ -905,58 +952,58 @@ strtod
if ((z = aadj) <= 0)
z = 1;
aadj = z;
aadj1 = dsign ? aadj : -aadj;
dval(&aadj1) = dsign ? aadj : -aadj;
}
word0(aadj1) += (2*P+1)*Exp_msk1 - y;
word0(&aadj1) += (2*P+1)*Exp_msk1 - y;
}
adj = aadj1 * ulp(dval(rv));
dval(rv) += adj;
dval(&adj) = dval(&aadj1) * ulp(&rv);
dval(&rv) += dval(&adj);
#else
#ifdef Sudden_Underflow
if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
dval(rv0) = dval(rv);
word0(rv) += P*Exp_msk1;
adj = aadj1 * ulp(dval(rv));
dval(rv) += adj;
if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
dval(&rv0) = dval(&rv);
word0(&rv) += P*Exp_msk1;
dval(&adj) = dval(&aadj1) * ulp(&rv);
dval(&rv) += adj;
#ifdef IBM
if ((word0(rv) & Exp_mask) < P*Exp_msk1)
if ((word0(&rv) & Exp_mask) < P*Exp_msk1)
#else
if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
#endif
{
if (word0(rv0) == Tiny0
&& word1(rv0) == Tiny1)
if (word0(&rv0) == Tiny0
&& word1(&rv0) == Tiny1)
goto undfl;
word0(rv) = Tiny0;
word1(rv) = Tiny1;
word0(&rv) = Tiny0;
word1(&rv) = Tiny1;
goto cont;
}
else
word0(rv) -= P*Exp_msk1;
word0(&rv) -= P*Exp_msk1;
}
else {
adj = aadj1 * ulp(dval(rv));
dval(rv) += adj;
dval(&adj) = dval(&aadj1) * ulp(&rv);
dval(&rv) += adj;
}
#else /*Sudden_Underflow*/
/* Compute adj so that the IEEE rounding rules will
* correctly round rv + adj in some half-way cases.
* If rv * ulp(rv) is denormalized (i.e.,
/* Compute dval(&adj) so that the IEEE rounding rules will
* correctly round rv + dval(&adj) in some half-way cases.
* If rv * ulp(&rv) is denormalized (i.e.,
* y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
* trouble from bits lost to denormalization;
* example: 1.2e-307 .
*/
if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
aadj1 = (double)(int)(aadj + 0.5);
dval(&aadj1) = (double)(int)(aadj + 0.5);
if (!dsign)
aadj1 = -aadj1;
dval(&aadj1) = -dval(&aadj1);
}
adj = aadj1 * ulp(dval(rv));
dval(rv) += adj;
dval(&adj) = dval(&aadj1) * ulp(&rv);
dval(&rv) += adj;
#endif /*Sudden_Underflow*/
#endif /*Avoid_Underflow*/
}
z = word0(rv) & Exp_mask;
z = word0(&rv) & Exp_mask;
#ifndef SET_INEXACT
#ifdef Avoid_Underflow
if (!scale)
@ -966,7 +1013,7 @@ strtod
L = (Long)aadj;
aadj -= L;
/* The tolerances below are conservative. */
if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
if (aadj < .4999999 || aadj > .5000001)
break;
}
@ -980,12 +1027,17 @@ strtod
Bfree(bs);
Bfree(delta);
}
Bfree(bb);
Bfree(bd);
Bfree(bs);
Bfree(bd0);
Bfree(delta);
#ifdef SET_INEXACT
if (inexact) {
if (!oldinexact) {
word0(rv0) = Exp_1 + (70 << Exp_shift);
word1(rv0) = 0;
dval(rv0) += 1.;
word0(&rv0) = Exp_1 + (70 << Exp_shift);
word1(&rv0) = 0;
dval(&rv0) += 1.;
}
}
else if (!oldinexact)
@ -993,36 +1045,30 @@ strtod
#endif
#ifdef Avoid_Underflow
if (scale) {
word0(rv0) = Exp_1 - 2*P*Exp_msk1;
word1(rv0) = 0;
dval(rv) *= dval(rv0);
word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
word1(&rv0) = 0;
dval(&rv) *= dval(&rv0);
#ifndef NO_ERRNO
/* try to avoid the bug of testing an 8087 register value */
#ifdef IEEE_Arith
if (!(word0(rv) & Exp_mask))
if (!(word0(&rv) & Exp_mask))
#else
if (word0(rv) == 0 && word1(rv) == 0)
if (word0(&rv) == 0 && word1(&rv) == 0)
#endif
errno = ERANGE;
#endif
}
#endif /* Avoid_Underflow */
#ifdef SET_INEXACT
if (inexact && !(word0(rv) & Exp_mask)) {
if (inexact && !(word0(&rv) & Exp_mask)) {
/* set underflow bit */
dval(rv0) = 1e-300;
dval(rv0) *= dval(rv0);
dval(&rv0) = 1e-300;
dval(&rv0) *= dval(&rv0);
}
#endif
retfree:
Bfree(bb);
Bfree(bd);
Bfree(bs);
Bfree(bd0);
Bfree(delta);
ret:
if (se)
*se = (char *)s;
return sign ? -dval(rv) : dval(rv);
return sign ? -dval(&rv) : dval(&rv);
}

88
strtodI.c
View File

@ -33,16 +33,16 @@ THIS SOFTWARE.
static double
#ifdef KR_headers
ulpdown(d) double *d;
ulpdown(d) U *d;
#else
ulpdown(double *d)
ulpdown(U *d)
#endif
{
double u;
ULong *L = (ULong*)d;
ULong *L = d->L;
u = ulp(*d);
if (!(L[_1] | L[_0] & 0xfffff)
u = ulp(d);
if (!(L[_1] | (L[_0] & 0xfffff))
&& (L[_0] & 0x7ff00000) > 0x00100000)
u *= 0.5;
return u;
@ -59,10 +59,6 @@ strtodI(CONST char *s, char **sp, double *dd)
ULong bits[2], sign;
Long exp;
int j, k;
typedef union {
double d[2];
ULong L[4];
} U;
U *u;
k = strtodg(s, sp, &fpi, &exp, bits);
@ -70,17 +66,17 @@ strtodI(CONST char *s, char **sp, double *dd)
sign = k & STRTOG_Neg ? 0x80000000L : 0;
switch(k & STRTOG_Retmask) {
case STRTOG_NoNumber:
u->d[0] = u->d[1] = 0.;
dval(&u[0]) = dval(&u[1]) = 0.;
break;
case STRTOG_Zero:
u->d[0] = u->d[1] = 0.;
dval(&u[0]) = dval(&u[1]) = 0.;
#ifdef Sudden_Underflow
if (k & STRTOG_Inexact) {
if (sign)
u->L[_0] = 0x80100000L;
word0(&u[0]) = 0x80100000L;
else
u->L[2+_0] = 0x100000L;
word0(&u[1]) = 0x100000L;
}
break;
#else
@ -88,80 +84,80 @@ strtodI(CONST char *s, char **sp, double *dd)
#endif
case STRTOG_Denormal:
u->L[_1] = bits[0];
u->L[_0] = bits[1];
word1(&u[0]) = bits[0];
word0(&u[0]) = bits[1];
goto contain;
case STRTOG_Normal:
u->L[_1] = bits[0];
u->L[_0] = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20);
word1(&u[0]) = bits[0];
word0(&u[0]) = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20);
contain:
j = k & STRTOG_Inexact;
if (sign) {
u->L[_0] |= sign;
word0(&u[0]) |= sign;
j = STRTOG_Inexact - j;
}
switch(j) {
case STRTOG_Inexlo:
#ifdef Sudden_Underflow
if ((u->L[_0] & 0x7ff00000) < 0x3500000) {
u->L[2+_0] = u->L[_0] + 0x3500000;
u->L[2+_1] = u->L[_1];
u->d[1] += ulp(u->d[1]);
u->L[2+_0] -= 0x3500000;
if (!(u->L[2+_0] & 0x7ff00000)) {
u->L[2+_0] = sign;
u->L[2+_1] = 0;
word0(&u[1]) = word0(&u[0]) + 0x3500000;
word1(&u[1]) = word1(&u[0]);
dval(&u[1]) += ulp(&u[1]);
word0(&u[1]) -= 0x3500000;
if (!(word0(&u[1]) & 0x7ff00000)) {
word0(&u[1]) = sign;
word1(&u[1]) = 0;
}
}
else
#endif
u->d[1] = u->d[0] + ulp(u->d[0]);
dval(&u[1]) = dval(&u[0]) + ulp(&u[0]);
break;
case STRTOG_Inexhi:
u->d[1] = u->d[0];
dval(&u[1]) = dval(&u[0]);
#ifdef Sudden_Underflow
if ((u->L[_0] & 0x7ff00000) < 0x3500000) {
u->L[_0] += 0x3500000;
u->d[0] -= ulpdown(u->d);
u->L[_0] -= 0x3500000;
if (!(u->L[_0] & 0x7ff00000)) {
u->L[_0] = sign;
u->L[_1] = 0;
if ((word0(&u[0]) & 0x7ff00000) < 0x3500000) {
word0(&u[0]) += 0x3500000;
dval(&u[0]) -= ulpdown(u);
word0(&u[0]) -= 0x3500000;
if (!(word0(&u[0]) & 0x7ff00000)) {
word0(&u[0]) = sign;
word1(&u[0]) = 0;
}
}
else
#endif
u->d[0] -= ulpdown(u->d);
dval(&u[0]) -= ulpdown(u);
break;
default:
u->d[1] = u->d[0];
dval(&u[1]) = dval(&u[0]);
}
break;
case STRTOG_Infinite:
u->L[_0] = u->L[2+_0] = sign | 0x7ff00000;
u->L[_1] = u->L[2+_1] = 0;
word0(&u[0]) = word0(&u[1]) = sign | 0x7ff00000;
word1(&u[0]) = word1(&u[1]) = 0;
if (k & STRTOG_Inexact) {
if (sign) {
u->L[2+_0] = 0xffefffffL;
u->L[2+_1] = 0xffffffffL;
word0(&u[1]) = 0xffefffffL;
word1(&u[1]) = 0xffffffffL;
}
else {
u->L[_0] = 0x7fefffffL;
u->L[_1] = 0xffffffffL;
word0(&u[0]) = 0x7fefffffL;
word1(&u[0]) = 0xffffffffL;
}
}
break;
case STRTOG_NaN:
u->L[0] = u->L[2] = d_QNAN0;
u->L[1] = u->L[3] = d_QNAN1;
u->L[0] = (u+1)->L[0] = d_QNAN0;
u->L[1] = (u+1)->L[1] = d_QNAN1;
break;
case STRTOG_NaNbits:
u->L[_0] = u->L[2+_0] = 0x7ff00000 | sign | bits[1];
u->L[_1] = u->L[2+_1] = bits[0];
word0(&u[0]) = word0(&u[1]) = 0x7ff00000 | sign | bits[1];
word1(&u[0]) = word1(&u[1]) = bits[0];
}
return k;
}

121
strtodg.c
View File

@ -172,9 +172,9 @@ set_ones(Bigint *b, int n)
rvOK
#ifdef KR_headers
(d, fpi, exp, bits, exact, rd, irv)
double d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv;
U *d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv;
#else
(double d, FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv)
(U *d, FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv)
#endif
{
Bigint *b;
@ -182,7 +182,7 @@ rvOK
int bdif, e, j, k, k1, nb, rv;
carry = rv = 0;
b = d2b(d, &e, &bdif);
b = d2b(dval(d), &e, &bdif);
bdif -= nb = fpi->nbits;
e += bdif;
if (bdif <= 0) {
@ -291,9 +291,9 @@ rvOK
static int
#ifdef KR_headers
mantbits(d) double d;
mantbits(d) U *d;
#else
mantbits(double d)
mantbits(U *d)
#endif
{
ULong L;
@ -327,8 +327,9 @@ strtodg
int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign;
int sudden_underflow;
CONST char *s, *s0, *s1;
double adj, adj0, rv, tol;
double adj0, tol;
Long L;
U adj, rv;
ULong *b, *be, y, z;
Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
#ifdef USE_LOCALE /*{{*/
@ -341,7 +342,7 @@ strtodg
static int dplen;
if (!(s0 = decimalpoint_cache)) {
s0 = localeconv()->decimal_point;
if ((decimalpoint_cache = (char*)malloc(strlen(s0) + 1))) {
if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
strcpy(decimalpoint_cache, s0);
s0 = decimalpoint_cache;
}
@ -355,7 +356,7 @@ strtodg
irv = STRTOG_Zero;
denorm = sign = nz0 = nz = 0;
dval(rv) = 0.;
dval(&rv) = 0.;
rvb = 0;
nbits = fpi->nbits;
for(s = s00;;s++) switch(*s) {
@ -547,13 +548,13 @@ strtodg
if (!nd0)
nd0 = nd;
k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
dval(rv) = y;
dval(&rv) = y;
if (k > 9)
dval(rv) = tens[k - 9] * dval(rv) + z;
dval(&rv) = tens[k - 9] * dval(&rv) + z;
bd0 = 0;
if (nbits <= P && nd <= DBL_DIG) {
if (!e) {
if (rvOK(dval(rv), fpi, exp, bits, 1, rd, &irv))
if (rvOK(&rv, fpi, exp, bits, 1, rd, &irv))
goto ret;
}
else if (e > 0) {
@ -561,9 +562,9 @@ strtodg
#ifdef VAX
goto vax_ovfl_check;
#else
i = fivesbits[e] + mantbits(dval(rv)) <= P;
/* rv = */ rounded_product(dval(rv), tens[e]);
if (rvOK(dval(rv), fpi, exp, bits, i, rd, &irv))
i = fivesbits[e] + mantbits(&rv) <= P;
/* rv = */ rounded_product(dval(&rv), tens[e]);
if (rvOK(&rv, fpi, exp, bits, i, rd, &irv))
goto ret;
e1 -= e;
goto rv_notOK;
@ -576,32 +577,32 @@ strtodg
*/
e2 = e - i;
e1 -= i;
dval(rv) *= tens[i];
dval(&rv) *= tens[i];
#ifdef VAX
/* VAX exponent range is so narrow we must
* worry about overflow here...
*/
vax_ovfl_check:
dval(adj) = dval(rv);
word0(adj) -= P*Exp_msk1;
/* adj = */ rounded_product(dval(adj), tens[e2]);
if ((word0(adj) & Exp_mask)
dval(&adj) = dval(&rv);
word0(&adj) -= P*Exp_msk1;
/* adj = */ rounded_product(dval(&adj), tens[e2]);
if ((word0(&adj) & Exp_mask)
> Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
goto rv_notOK;
word0(adj) += P*Exp_msk1;
dval(rv) = dval(adj);
word0(&adj) += P*Exp_msk1;
dval(&rv) = dval(&adj);
#else
/* rv = */ rounded_product(dval(rv), tens[e2]);
/* rv = */ rounded_product(dval(&rv), tens[e2]);
#endif
if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv))
goto ret;
e1 -= e2;
}
}
#ifndef Inaccurate_Divide
else if (e >= -Ten_pmax) {
/* rv = */ rounded_quotient(dval(rv), tens[-e]);
if (rvOK(dval(rv), fpi, exp, bits, 0, rd, &irv))
/* rv = */ rounded_quotient(dval(&rv), tens[-e]);
if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv))
goto ret;
e1 -= e;
}
@ -615,45 +616,45 @@ strtodg
e2 = 0;
if (e1 > 0) {
if ( (i = e1 & 15) !=0)
dval(rv) *= tens[i];
dval(&rv) *= tens[i];
if (e1 &= ~15) {
e1 >>= 4;
while(e1 >= (1 << n_bigtens-1)) {
e2 += ((word0(rv) & Exp_mask)
while(e1 >= (1 << (n_bigtens-1))) {
e2 += ((word0(&rv) & Exp_mask)
>> Exp_shift1) - Bias;
word0(rv) &= ~Exp_mask;
word0(rv) |= Bias << Exp_shift1;
dval(rv) *= bigtens[n_bigtens-1];
e1 -= 1 << n_bigtens-1;
word0(&rv) &= ~Exp_mask;
word0(&rv) |= Bias << Exp_shift1;
dval(&rv) *= bigtens[n_bigtens-1];
e1 -= 1 << (n_bigtens-1);
}
e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
word0(rv) &= ~Exp_mask;
word0(rv) |= Bias << Exp_shift1;
e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
word0(&rv) &= ~Exp_mask;
word0(&rv) |= Bias << Exp_shift1;
for(j = 0; e1 > 0; j++, e1 >>= 1)
if (e1 & 1)
dval(rv) *= bigtens[j];
dval(&rv) *= bigtens[j];
}
}
else if (e1 < 0) {
e1 = -e1;
if ( (i = e1 & 15) !=0)
dval(rv) /= tens[i];
dval(&rv) /= tens[i];
if (e1 &= ~15) {
e1 >>= 4;
while(e1 >= (1 << n_bigtens-1)) {
e2 += ((word0(rv) & Exp_mask)
while(e1 >= (1 << (n_bigtens-1))) {
e2 += ((word0(&rv) & Exp_mask)
>> Exp_shift1) - Bias;
word0(rv) &= ~Exp_mask;
word0(rv) |= Bias << Exp_shift1;
dval(rv) *= tinytens[n_bigtens-1];
e1 -= 1 << n_bigtens-1;
word0(&rv) &= ~Exp_mask;
word0(&rv) |= Bias << Exp_shift1;
dval(&rv) *= tinytens[n_bigtens-1];
e1 -= 1 << (n_bigtens-1);
}
e2 += ((word0(rv) & Exp_mask) >> Exp_shift1) - Bias;
word0(rv) &= ~Exp_mask;
word0(rv) |= Bias << Exp_shift1;
e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
word0(&rv) &= ~Exp_mask;
word0(&rv) |= Bias << Exp_shift1;
for(j = 0; e1 > 0; j++, e1 >>= 1)
if (e1 & 1)
dval(rv) *= tinytens[j];
dval(&rv) *= tinytens[j];
}
}
#ifdef IBM
@ -664,7 +665,7 @@ strtodg
*/
e2 <<= 2;
#endif
rvb = d2b(dval(rv), &rve, &rvbits); /* rv = rvb * 2^rve */
rvb = d2b(dval(&rv), &rve, &rvbits); /* rv = rvb * 2^rve */
rve += e2;
if ((j = rvbits - nbits) > 0) {
rshift(rvb, j);
@ -842,7 +843,7 @@ strtodg
}
else
irv = STRTOG_Normal | STRTOG_Inexhi;
if (bbbits < nbits && !denorm || !(rvb->x[0] & 1))
if ((bbbits < nbits && !denorm) || !(rvb->x[0] & 1))
break;
if (dsign) {
rvb = increment(rvb);
@ -859,7 +860,7 @@ strtodg
}
break;
}
if ((dval(adj) = ratio(delta, bs)) <= 2.) {
if ((dval(&adj) = ratio(delta, bs)) <= 2.) {
adj1:
inex = STRTOG_Inexlo;
if (dsign) {
@ -873,15 +874,15 @@ strtodg
irv = STRTOG_Underflow | STRTOG_Inexlo;
break;
}
adj0 = dval(adj) = 1.;
adj0 = dval(&adj) = 1.;
}
else {
adj0 = dval(adj) *= 0.5;
adj0 = dval(&adj) *= 0.5;
if (dsign) {
asub = 0;
inex = STRTOG_Inexlo;
}
if (dval(adj) < 2147483647.) {
if (dval(&adj) < 2147483647.) {
L = adj0;
adj0 -= L;
switch(rd) {
@ -900,12 +901,12 @@ strtodg
inex = STRTOG_Inexact - inex;
}
}
dval(adj) = L;
dval(&adj) = L;
}
}
y = rve + rvbits;
/* adj *= ulp(dval(rv)); */
/* adj *= ulp(dval(&rv)); */
/* if (asub) rv -= adj; else rv += adj; */
if (!denorm && rvbits < nbits) {
@ -913,7 +914,7 @@ strtodg
rve -= j;
rvbits = nbits;
}
ab = d2b(dval(adj), &abe, &abits);
ab = d2b(dval(&adj), &abe, &abits);
if (abe < 0)
rshift(ab, -abe);
else if (abe > 0)
@ -967,15 +968,15 @@ strtodg
z = rve + rvbits;
if (y == z && L) {
/* Can we stop now? */
tol = dval(adj) * 5e-16; /* > max rel error */
dval(adj) = adj0 - .5;
if (dval(adj) < -tol) {
tol = dval(&adj) * 5e-16; /* > max rel error */
dval(&adj) = adj0 - .5;
if (dval(&adj) < -tol) {
if (adj0 > tol) {
irv |= inex;
break;
}
}
else if (dval(adj) > tol && adj0 < 1. - tol) {
else if (dval(&adj) > tol && adj0 < 1. - tol) {
irv |= inex;
break;
}

2
strtof.c
View File

@ -58,7 +58,7 @@ strtof(CONST char *s, char **sp)
case STRTOG_Normal:
case STRTOG_NaNbits:
u.L[0] = bits[0] & 0x7fffff | exp + 0x7f + 23 << 23;
u.L[0] = (bits[0] & 0x7fffff) | ((exp + 0x7f + 23) << 23);
break;
case STRTOG_Denormal:

36
strtopdd.c
View File

@ -67,8 +67,8 @@ strtopdd(CONST char *s, char **sp, double *dd)
case STRTOG_Normal:
u->L[_1] = (bits[1] >> 21 | bits[2] << 11) & 0xffffffffL;
u->L[_0] = bits[2] >> 21 | bits[3] << 11 & 0xfffff
| exp + 0x3ff + 105 << 20;
u->L[_0] = (bits[2] >> 21) | ((bits[3] << 11) & 0xfffff)
| ((exp + 0x3ff + 105) << 20);
exp += 0x3ff + 52;
if (bits[1] &= 0x1fffff) {
i = hi0bits(bits[1]) - 11;
@ -79,7 +79,7 @@ strtopdd(CONST char *s, char **sp, double *dd)
else
exp -= i;
if (i > 0) {
bits[1] = bits[1] << i | bits[0] >> 32-i;
bits[1] = bits[1] << i | bits[0] >> (32-i);
bits[0] = bits[0] << i & 0xffffffffL;
}
}
@ -92,11 +92,11 @@ strtopdd(CONST char *s, char **sp, double *dd)
else
exp -= i;
if (i < 32) {
bits[1] = bits[0] >> 32 - i;
bits[1] = bits[0] >> (32 - i);
bits[0] = bits[0] << i & 0xffffffffL;
}
else {
bits[1] = bits[0] << i - 32;
bits[1] = bits[0] << (i - 32);
bits[0] = 0;
}
}
@ -105,7 +105,7 @@ strtopdd(CONST char *s, char **sp, double *dd)
break;
}
u->L[2+_1] = bits[0];
u->L[2+_0] = bits[1] & 0xfffff | exp << 20;
u->L[2+_0] = (bits[1] & 0xfffff) | (exp << 20);
break;
case STRTOG_Denormal:
@ -124,10 +124,10 @@ strtopdd(CONST char *s, char **sp, double *dd)
nearly_normal:
i = hi0bits(bits[3]) - 11; /* i >= 12 */
j = 32 - i;
u->L[_0] = (bits[3] << i | bits[2] >> j) & 0xfffff
| 65 - i << 20;
u->L[_0] = ((bits[3] << i | bits[2] >> j) & 0xfffff)
| ((65 - i) << 20);
u->L[_1] = (bits[2] << i | bits[1] >> j) & 0xffffffffL;
u->L[2+_0] = bits[1] & (1L << j) - 1;
u->L[2+_0] = bits[1] & ((1L << j) - 1);
u->L[2+_1] = bits[0];
break;
@ -136,34 +136,34 @@ strtopdd(CONST char *s, char **sp, double *dd)
if (i < 0) {
j = -i;
i += 32;
u->L[_0] = bits[2] >> j & 0xfffff | (33 + j) << 20;
u->L[_1] = (bits[2] << i | bits[1] >> j) & 0xffffffffL;
u->L[2+_0] = bits[1] & (1L << j) - 1;
u->L[_0] = (bits[2] >> j & 0xfffff) | (33 + j) << 20;
u->L[_1] = ((bits[2] << i) | (bits[1] >> j)) & 0xffffffffL;
u->L[2+_0] = bits[1] & ((1L << j) - 1);
u->L[2+_1] = bits[0];
break;
}
if (i == 0) {
u->L[_0] = bits[2] & 0xfffff | 33 << 20;
u->L[_0] = (bits[2] & 0xfffff) | (33 << 20);
u->L[_1] = bits[1];
u->L[2+_0] = 0;
u->L[2+_1] = bits[0];
break;
}
j = 32 - i;
u->L[_0] = (bits[2] << i | bits[1] >> j) & 0xfffff
| j + 1 << 20;
u->L[_0] = (((bits[2] << i) | (bits[1] >> j)) & 0xfffff)
| ((j + 1) << 20);
u->L[_1] = (bits[1] << i | bits[0] >> j) & 0xffffffffL;
u->L[2+_0] = 0;
u->L[2+_1] = bits[0] & (1L << j) - 1;
u->L[2+_1] = bits[0] & ((1L << j) - 1);
break;
hardly_normal:
j = 11 - hi0bits(bits[1]);
i = 32 - j;
u->L[_0] = bits[1] >> j & 0xfffff | j + 1 << 20;
u->L[_0] = (bits[1] >> j & 0xfffff) | ((j + 1) << 20);
u->L[_1] = (bits[1] << i | bits[0] >> j) & 0xffffffffL;
u->L[2+_0] = 0;
u->L[2+_1] = bits[0] & (1L << j) - 1;
u->L[2+_1] = bits[0] & ((1L << j) - 1);
break;
case STRTOG_Infinite:

2
strtopf.c
View File

@ -58,7 +58,7 @@ strtopf(CONST char *s, char **sp, float *f)
case STRTOG_Normal:
case STRTOG_NaNbits:
L[0] = bits[0] & 0x7fffff | exp + 0x7f + 23 << 23;
L[0] = (bits[0] & 0x7fffff) | ((exp + 0x7f + 23) << 23);
break;
case STRTOG_Denormal:

3
strtopx.c
View File

@ -92,7 +92,8 @@ strtopx(CONST char *s, char **sp, void *V)
case STRTOG_Infinite:
L[_0] = 0x7fff;
L[_1] = L[_2] = L[_3] = L[_4] = 0;
L[_1] = 0x8000;
L[_2] = L[_3] = L[_4] = 0;
break;
case STRTOG_NaN:

3
strtopxL.c
View File

@ -82,7 +82,8 @@ strtopxL(CONST char *s, char **sp, void *V)
case STRTOG_Infinite:
L[_0] = 0x7fff << 16;
L[_1] = L[_2] = 0;
L[_1] = 0x80000000;
L[_2] = 0;
break;
case STRTOG_NaN:

34
strtordd.c
View File

@ -48,8 +48,8 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
case STRTOG_Normal:
L[_1] = (bits[1] >> 21 | bits[2] << 11) & (ULong)0xffffffffL;
L[_0] = bits[2] >> 21 | bits[3] << 11 & 0xfffff
| exp + 0x3ff + 105 << 20;
L[_0] = (bits[2] >> 21) | (bits[3] << 11 & 0xfffff)
| ((exp + 0x3ff + 105) << 20);
exp += 0x3ff + 52;
if (bits[1] &= 0x1fffff) {
i = hi0bits(bits[1]) - 11;
@ -60,7 +60,7 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
else
exp -= i;
if (i > 0) {
bits[1] = bits[1] << i | bits[0] >> 32-i;
bits[1] = bits[1] << i | bits[0] >> (32-i);
bits[0] = bits[0] << i & (ULong)0xffffffffL;
}
}
@ -73,11 +73,11 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
else
exp -= i;
if (i < 32) {
bits[1] = bits[0] >> 32 - i;
bits[1] = bits[0] >> (32 - i);
bits[0] = bits[0] << i & (ULong)0xffffffffL;
}
else {
bits[1] = bits[0] << i - 32;
bits[1] = bits[0] << (i - 32);
bits[0] = 0;
}
}
@ -86,7 +86,7 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
break;
}
L[2+_1] = bits[0];
L[2+_0] = bits[1] & 0xfffff | exp << 20;
L[2+_0] = (bits[1] & 0xfffff) | (exp << 20);
break;
case STRTOG_Denormal:
@ -105,10 +105,10 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
nearly_normal:
i = hi0bits(bits[3]) - 11; /* i >= 12 */
j = 32 - i;
L[_0] = (bits[3] << i | bits[2] >> j) & 0xfffff
| 65 - i << 20;
L[_0] = ((bits[3] << i | bits[2] >> j) & 0xfffff)
| ((65 - i) << 20);
L[_1] = (bits[2] << i | bits[1] >> j) & 0xffffffffL;
L[2+_0] = bits[1] & ((ULong)1L << j) - 1;
L[2+_0] = bits[1] & (((ULong)1L << j) - 1);
L[2+_1] = bits[0];
break;
@ -117,34 +117,34 @@ ULtodd(ULong *L, ULong *bits, Long exp, int k)
if (i < 0) {
j = -i;
i += 32;
L[_0] = bits[2] >> j & 0xfffff | (33 + j) << 20;
L[_0] = (bits[2] >> j & 0xfffff) | ((33 + j) << 20);
L[_1] = (bits[2] << i | bits[1] >> j) & 0xffffffffL;
L[2+_0] = bits[1] & ((ULong)1L << j) - 1;
L[2+_0] = bits[1] & (((ULong)1L << j) - 1);
L[2+_1] = bits[0];
break;
}
if (i == 0) {
L[_0] = bits[2] & 0xfffff | 33 << 20;
L[_0] = (bits[2] & 0xfffff) | (33 << 20);
L[_1] = bits[1];
L[2+_0] = 0;
L[2+_1] = bits[0];
break;
}
j = 32 - i;
L[_0] = (bits[2] << i | bits[1] >> j) & 0xfffff
| j + 1 << 20;
L[_0] = (((bits[2] << i) | (bits[1] >> j)) & 0xfffff)
| ((j + 1) << 20);
L[_1] = (bits[1] << i | bits[0] >> j) & 0xffffffffL;
L[2+_0] = 0;
L[2+_1] = bits[0] & (1L << j) - 1;
L[2+_1] = bits[0] & ((1L << j) - 1);
break;
hardly_normal:
j = 11 - hi0bits(bits[1]);
i = 32 - j;
L[_0] = bits[1] >> j & 0xfffff | j + 1 << 20;
L[_0] = (bits[1] >> j & 0xfffff) | ((j + 1) << 20);
L[_1] = (bits[1] << i | bits[0] >> j) & 0xffffffffL;
L[2+_0] = 0;
L[2+_1] = bits[0] & ((ULong)1L << j) - 1;
L[2+_1] = bits[0] & (((ULong)1L << j) - 1);
break;
case STRTOG_Infinite:

2
strtorf.c
View File

@ -46,7 +46,7 @@ ULtof(ULong *L, ULong *bits, Long exp, int k)
case STRTOG_Normal:
case STRTOG_NaNbits:
L[0] = bits[0] & 0x7fffff | exp + 0x7f + 23 << 23;
L[0] = (bits[0] & 0x7fffff) | ((exp + 0x7f + 23) << 23);
break;
case STRTOG_Denormal:

3
strtorx.c
View File

@ -80,7 +80,8 @@ ULtox(UShort *L, ULong *bits, Long exp, int k)
case STRTOG_Infinite:
L[_0] = 0x7fff;
L[_1] = L[_2] = L[_3] = L[_4] = 0;
L[_1] = 0x8000;
L[_2] = L[_3] = L[_4] = 0;
break;
case STRTOG_NaN:

3
strtorxL.c
View File

@ -70,7 +70,8 @@ ULtoxL(ULong *L, ULong *bits, Long exp, int k)
case STRTOG_Infinite:
L[_0] = 0x7fff << 16;
L[_1] = L[_2] = 0;
L[_1] = 0x80000000;
L[_2] = 0;
break;
case STRTOG_NaN:

23
test/Qtest.c
View File

@ -81,9 +81,14 @@ THIS SOFTWARE.
main(Void)
{
char *s, *s1, *se, *se1;
int i, dItry, ndig = 0, r = 1;
union { long double d; ULong bits[4]; } u, v[2];
int Ltest, i, dItry, ndig = 0, r = 1;
union { long double d; ULong bits[4]; } u, v[2], w;
w.bits[0] = w.bits[3] = 0;
w.d = 1.;
u.d = 3.;
w.d = w.d / u.d;
Ltest = sizeof(long double) == 16 && w.bits[0] && w.bits[3];
while( (s = fgets(ibuf, sizeof(ibuf), stdin)) !=0) {
while(*s <= ' ')
if (!*s++)
@ -95,7 +100,7 @@ main(Void)
continue;
case 'n':
i = s[1];
if (i <= ' ' || i >= '0' && i <= '9') {
if (i <= ' ' || (i >= '0' && i <= '9')) {
ndig = atoi(s+1);
continue;
}
@ -113,8 +118,8 @@ main(Void)
}
}
printf("\nInput: %s", ibuf);
printf(" --> f = #%lx %lx %lx %lx\n", u.bits[_0],
u.bits[_1], u.bits[_2], u.bits[_3]);
printf(" --> f = #%lx %lx %lx %lx\n", U u.bits[_0],
U u.bits[_1], U u.bits[_2], U u.bits[_3]);
goto fmt_test;
}
dItry = 1;
@ -128,7 +133,7 @@ main(Void)
printf("with bits = #%lx %lx %lx %lx\n",
U u.bits[_0], U u.bits[_1], U u.bits[_2], U u.bits[_3]);
fmt_test:
if (sizeof(long double) == 16)
if (Ltest)
printf("printf(\"%%.35Lg\") gives %.35Lg\n", u.d);
se = g_Qfmt(obuf, u.bits, ndig, sizeof(obuf));
printf("g_Qfmt(%d) gives %d bytes: \"%s\"\n\n", ndig,
@ -145,7 +150,7 @@ main(Void)
printf("fI[0] == fI[1] = #%lx %lx %lx %lx\n",
U v[0].bits[_0], U v[0].bits[_1],
U v[0].bits[_2], U v[0].bits[_3]);
if (sizeof(long double) == 16)
if (Ltest)
printf("= %.35Lg\n", v[0].d);
}
}
@ -153,12 +158,12 @@ main(Void)
printf("fI[0] = #%lx %lx %lx %lx\n",
U v[0].bits[_0], U v[0].bits[_1],
U v[0].bits[_2], U v[0].bits[_3]);
if (sizeof(long double) == 16)
if (Ltest)
printf("= %.35Lg\n", v[0].d);
printf("fI[1] = #%lx %lx %lx %lx\n",
U v[1].bits[_0], U v[1].bits[_1],
U v[1].bits[_2], U v[1].bits[_3]);
if (sizeof(long double) == 16)
if (Ltest)
printf("= %.35Lg\n", v[1].d);
if (!memcmp(v[0].bits, u.bits, 16))
printf("fI[0] == strtod\n");

6
test/README
View File

@ -69,6 +69,12 @@ You can also or alternatively compile getround.c with
-DUSE_MY_LOCALE (when ../gdtoa.a is compiled with -DUSE_LOCALE)
to test multi-byte decimal points.
If in the parent directory, you have sucessfully invoked "make Printf"
to add a "printf" (called Printf and accessed via ../stdio1.h), then
here you can use "make pf_test" and (if you have both a 64-bit long
double and a 113-bit "quad" double type) "make pf_testLq" for a brief
test of %g and %a variants in Printf.
These are simple test programs, not meant for exhaustive testing,
but for manually testing "interesting" cases. Paxson's testbase
is good for more exhaustive testing, in part with random inputs.

2
test/ddtest.c
View File

@ -94,7 +94,7 @@ main(Void)
continue;
case 'n':
i = s[1];
if (i <= ' ' || i >= '0' && i <= '9') {
if (i <= ' ' || (i >= '0' && i <= '9')) {
ndig = atoi(s+1);
continue;
}

112
test/dt.c
View File

@ -44,6 +44,7 @@ THIS SOFTWARE.
#include <stdio.h>
#include "gdtoa.h"
int STRTOD_DIGLIM = 24;
#ifdef KR_headers
#define Void /*void*/
#else
@ -59,13 +60,18 @@ extern "C" double atof(const char*);
extern double atof ANSI((char*));
#endif
#endif
typedef union { double d; ULong L[2]; } U;
#ifdef IEEE_8087
#define word0(x) ((ULong *)&x)[1]
#define word1(x) ((ULong *)&x)[0]
#define word0(x) (x)->L[1]
#define word1(x) (x)->L[0]
#else
#define word0(x) ((ULong *)&x)[0]
#define word1(x) ((ULong *)&x)[1]
#define word0(x) (x)->L[0]
#define word1(x) (x)->L[1]
#endif
#define dval(x) (x)->d
#include "errno.h"
#ifdef __cplusplus
@ -93,14 +99,14 @@ g_fmt(char *b, double x)
if (sign)
*b++ = '-';
if (decpt == 9999) /* Infinity or Nan */ {
while(*b++ = *s++);
while((*b++ = *s++));
return;
}
if (decpt <= -4 || decpt > se - s + 5) {
*b++ = *s++;
if (*s) {
*b++ = '.';
while(*b = *s++)
while((*b = *s++))
b++;
}
*b++ = 'e';
@ -126,10 +132,10 @@ g_fmt(char *b, double x)
*b++ = '.';
for(; decpt < 0; decpt++)
*b++ = '0';
while(*b++ = *s++);
while((*b++ = *s++));
}
else {
while(*b = *s++) {
while((*b = *s++)) {
b++;
if (--decpt == 0 && *s)
*b++ = '.';
@ -148,40 +154,41 @@ baderrno(Void)
fflush(stderr);
}
#define U (unsigned long)
#define UL (unsigned long)
static void
#ifdef KR_headers
check(d) double d;
check(d) U *d;
#else
check(double d)
check(U *d)
#endif
{
char buf[64];
int decpt, sign;
char *s, *se;
double d1;
U d1;
s = dtoa(d, 0, 0, &decpt, &sign, &se);
s = dtoa(dval(d), 0, 0, &decpt, &sign, &se);
sprintf(buf, "%s%s%se%d", sign ? "-" : "",
decpt == 9999 ? "" : ".", s, decpt);
errno = 0;
d1 = strtod(buf, (char **)0);
dval(&d1) = strtod(buf, (char **)0);
if (errno)
baderrno();
if (d != d1) {
if (dval(d) != dval(&d1)) {
printf("sent d = %.17g = 0x%lx %lx, buf = %s\n",
d, U word0(d), U word1(d), buf);
dval(d), UL word0(d), UL word1(d), buf);
printf("got d1 = %.17g = 0x%lx %lx\n",
d1, U word0(d1), U word1(d1));
dval(&d1), UL word0(&d1), UL word1(&d1));
}
}
int
main(Void){
main(Void)
{
U d, d1;
char buf[2048], buf1[32];
char *fmt, *s, *s1, *se;
double d, d1;
int decpt, sign;
int mode = 0, ndigits = 17;
ULong x, y;
@ -196,8 +203,8 @@ main(Void){
}
printf("Input: %s", buf);
if (*buf == '#') {
x = word0(d);
y = word1(d);
x = word0(&d);
y = word1(&d);
/* sscanf(buf+1, "%lx %lx:%d %d", &x, &y, &mode, &ndigits); */
x = (ULong)strtoul(s1 = buf+1, &se, 16);
if (se > s1) {
@ -205,72 +212,79 @@ main(Void){
if (se > s1)
sscanf(se, ":%d %d", &mode, &ndigits);
}
word0(d) = x;
word1(d) = y;
word0(&d) = x;
word1(&d) = y;
fmt = "Output: d =\n%.17g = 0x%lx %lx\n";
}
else if (*buf == '*') {
x = strtoul(buf,&s,10);
if (!*s && x > 18)
STRTOD_DIGLIM = (int)x;
printf("STRTOD_DIGLIM = %lu\n", UL x);
continue;
}
else {
errno = 0;
d = strtod(buf,&se);
dval(&d) = strtod(buf,&se);
if (*se == ':')
sscanf(se+1,"%d %d", &mode, &ndigits);
d1 = atof(buf);
dval(&d1) = atof(buf);
fmt = "Output: d =\n%.17g = 0x%lx %lx, se = %s";
if (errno)
baderrno();
}
printf(fmt, d, U word0(d), U word1(d), se);
g_fmt(buf1, d);
printf(fmt, dval(&d), UL word0(&d), UL word1(&d), se);
g_fmt(buf1, dval(&d));
printf("\tg_fmt gives \"%s\"\n", buf1);
if (*buf != '#' && d != d1)
if (*buf != '#' && dval(&d) != dval(&d1))
printf("atof gives\n\
d1 = %.17g = 0x%lx %lx\nversus\n\
d = %.17g = 0x%lx %lx\n", d1, U word0(d1), U word1(d1),
d, U word0(d), U word1(d));
check(d);
s = dtoa(d, mode, ndigits, &decpt, &sign, &se);
d = %.17g = 0x%lx %lx\n", dval(&d1), UL word0(&d1), UL word1(&d1),
dval(&d), UL word0(&d), UL word1(&d));
check(&d);
s = dtoa(dval(&d), mode, ndigits, &decpt, &sign, &se);
printf("\tdtoa(mode = %d, ndigits = %d):\n", mode, ndigits);
printf("\tdtoa returns sign = %d, decpt = %d, %d digits:\n%s\n",
sign, decpt, se-s, s);
x = word1(d);
sign, decpt, (int)(se-s), s);
x = word1(&d);
if (x != 0xffffffff
&& (word0(d) & 0x7ff00000) != 0x7ff00000) {
&& (word0(&d) & 0x7ff00000) != 0x7ff00000) {
#ifdef VAX
z = x << 16 | x >> 16;
z++;
z = z << 16 | z >> 16;
word1(d) = z;
word1(&d) = z;
#else
word1(d) = x + 1;
word1(&d) = x + 1;
#endif
printf("\tnextafter(d,+Inf) = %.17g = 0x%lx %lx:\n",
d, U word0(d), U word1(d));
g_fmt(buf1, d);
dval(&d), UL word0(&d), UL word1(&d));
g_fmt(buf1, dval(&d));
printf("\tg_fmt gives \"%s\"\n", buf1);
s = dtoa(d, mode, ndigits, &decpt, &sign, &se);
s = dtoa(dval(&d), mode, ndigits, &decpt, &sign, &se);
printf(
"\tdtoa returns sign = %d, decpt = %d, %d digits:\n%s\n",
sign, decpt, se-s, s);
check(d);
sign, decpt, (int)(se-s), s);
check(&d);
}
if (x) {
#ifdef VAX
z = x << 16 | x >> 16;
z--;
z = z << 16 | z >> 16;
word1(d) = z;
word1(&d) = z;
#else
word1(d) = x - 1;
word1(&d) = x - 1;
#endif
printf("\tnextafter(d,-Inf) = %.17g = 0x%lx %lx:\n",
d, U word0(d), U word1(d));
g_fmt(buf1, d);
dval(&d), UL word0(&d), UL word1(&d));
g_fmt(buf1, dval(&d));
printf("\tg_fmt gives \"%s\"\n", buf1);
s = dtoa(d, mode, ndigits, &decpt, &sign, &se);
s = dtoa(dval(&d), mode, ndigits, &decpt, &sign, &se);
printf(
"\tdtoa returns sign = %d, decpt = %d, %d digits:\n%s\n",
sign, decpt, se-s, s);
check(d);
sign, decpt, (int)(se-s), s);
check(&d);
}
}
return 0;

2
test/dtest.c
View File

@ -78,7 +78,7 @@ main(Void)
continue;
case 'n':
i = s[1];
if (i <= ' ' || i >= '0' && i <= '9') {
if (i <= ' ' || (i >= '0' && i <= '9')) {
ndig = atoi(s+1);
continue;
}

2
test/ftest.c
View File

@ -77,7 +77,7 @@ main(Void)
continue;
case 'n':
i = s[1];
if (i <= ' ' || i >= '0' && i <= '9') {
if (i <= ' ' || (i >= '0' && i <= '9')) {
ndig = atoi(s+1);
continue;
}

2
test/getround.c
View File

@ -55,7 +55,7 @@ getround(int r, char *s)
r, dir[r]);
return r;
}
}
}
i = atoi(s);
if (i >= 0 && i < 4) {
printf("Rounding mode for strtor... ");

31
test/makefile
View File

@ -91,6 +91,10 @@ strtodt = strtodt.o $A
strtodt: $(strtodt)
$(CC) -o strtodt $(strtodt) $L
pftest = pftest.o $A
pftest: ../Printf $(pftest)
$(CC) -o pftest $(pftest) $L
## On Intel (and Intel-like) systems using extended-precision registers
## for double-precision (C type double) computations that sometimes suffer
## double rounding errors, the test below involving strtodt generally shows
@ -138,16 +142,31 @@ tests: Q.out x.out xL.out dt dItest ddtest dtest ftest Qtest xLtest xtest ddtest
(cd bad; for i in *; do cmp -s $$i ../obad/$$i && rm $$i;done; cd ..; rmdir bad)
touch tests
xs0 = README Qtest.c dItest.c ddtest.c dtest.c dt.c ftest.c getround.c \
strtoIdSI.c strtoIddSI.c strtodISI.c strtodt.c strtopddSI.c \
strtorddSI.c xLtest.c xQtest.c xtest.c rtestnos testnos testnos1 \
testnos3 dI.out dIsi.out ddsi.out dd.out dtst.out d.out f.out \
x.ou0 xL.ou0 x.ou1 xL.ou1 Q.ou0 Q.ou1 makefile
# To test Printf in ../gdtoa.a, "make pf_test" and perhaps "make pf_testLq"
# (if both long double and quad are desired and available).
pf_test: pftest
./pftest <pftestnos >zap
cmp pftest.out zap && rm zap
pf_testLq: pftest
./pftest <pfLqtestnos >zap
cmp pftestLq.out zap && rm zap
xs0 = README Q.ou0 Q.ou1 Qtest.c d.out dI.out dIsi.out dItest.c dd.out\
ddsi.out ddtest.c dt.c dtest.c dtst.out f.out ftest.c\
getround.c makefile pfLqtestnos pftest.c pftestQ.out\
pftestx.out pftestLq.out pftestnos rtestnos strtoIdSI.c\
strtoIddSI.c strtodISI.c strtodt.c strtopddSI.c strtorddSI.c\
testnos testnos1 testnos3 x.ou0 x.ou1 xL.ou0 xL.ou1 xLtest.c\
xQtest.c xtest.c
# invoke "make -r xsum.out"
xsum.out: xsum0.out $(xs0)
xsum $(xs0) >xsum1.out
cmp xsum0.out xsum1.out && mv xsum1.out xsum.out || diff xsum[01].out
clean:
rm -f *.[ao] dt *test *testsi strtodt strtodtnrp xsum.out xsum1.out tests zap x.out xL.out Q.out
rm -f *.[ao] dt *test *testsi pftest.out strtodt strtodtnrp xsum.out\
xsum1.out tests zap x.out xL.out Q.out
rm -rf bad

13
test/pfLqtestnos Normal file
View File

@ -0,0 +1,13 @@
%.8a 1.23
%.7a
%.2a
%.La 1.23
%.20Lqa 1.23
%La 3e27
%La 1.7e27
%Lqa 3e48
%Lqa 1e49
1.5
2.5
4.5
8.5

158
test/pftest.c Normal file
View File

@ -0,0 +1,158 @@
/****************************************************************
The author of this software is David M. Gay.
Copyright (C) 2009 by David M. Gay
All Rights Reserved
Permission to use, copy, modify, and distribute this software and
its documentation for any purpose and without fee is hereby
granted, provided that the above copyright notice appear in all
source-code copies and that both that the copyright notice and this
permission notice and warranty disclaimer appear in supporting
documentation.
THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, INDIRECT OR
CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF
USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR
OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
PERFORMANCE OF THIS SOFTWARE.
****************************************************************/
#include "stdio1.h"
#include "gdtoa.h"
#include <string.h>
#undef allow_Quad
#undef want_Quad
#undef want_Ux
#define want_LD
typedef union Ud {double x; unsigned int u[2]; } Ud;
#ifdef __x86_64 /*{{*/
#define want_Ux
#ifndef NO_GDTOA_i386_Quad /*{*/
typedef union UQ {__float128 x; unsigned int u[4]; } UQ;
#define allow_Quad(x) x
#define want_Quad
#endif /*}*/
#else /*}{*/
#ifdef __i386 /*{{*/
#define want_Ux
#else /*}{*/
#ifdef __sparc /*{{*/
typedef union UQ {long double x; unsigned int u[4]; } Ux;
#else /*}{*/
#ifdef __INTEL_COMPILER /*{*/
#undef want_Quad
#undef want_Ux
#undef want_LD
#endif /*}*/
#endif /*}}*/
#endif /*}}*/
#endif /*}}*/
#ifndef allow_Quad
#define allow_Quad(x) /*nothing*/
#endif
#ifdef want_Ux /*{{*/
typedef union Ux {long double x; unsigned short u[5]; } Ux;
#else /*}{*/
#ifdef __sparc
#define want_Ux
#endif
#endif /*}}*/
int
main(void)
{
Ud d;
allow_Quad(UQ q;)
char *b, buf[256], fmt[32], *s;
#ifdef want_Ux
Ux x;
x.x = 0.;
#endif
int k;
k = 0;
strcpy(fmt, "%.g");
d.x = 0.;
allow_Quad(q.x = 0.;)
while(fgets(buf, sizeof(buf), stdin)) {
for(b = buf; *b && *b != '\n'; ++b);
*b = 0;
if (b == buf)
continue;
b = buf;
if (*b == '%') {
for(k = 0; *b > ' '; ++b)
#ifdef want_LD /*{{*/
switch(*b) {
case 'L':
k = 1;
#ifdef want_Quad
break;
case 'q':
if (k >= 1)
k = 2;
#endif
}
#else /*}{*/
;
#endif /*}}*/
if (*b)
*b++ = 0;
if (b - buf < sizeof(fmt)) {
strcpy(fmt, buf);
}
}
if (*b) {
switch(k) {
case 0:
d.x = strtod(b,&s);
break;
case 1:
#ifdef want_Ux
#ifdef __sparc
strtopQ(b,&s,&x.x);
#else
strtopx(b,&s,&x.x);
#endif
#else
strtopQ(b,&s,&q.x);
#endif
break;
allow_Quad(case 2: strtopQ(b,&s,&q.x);)
}
if (*s)
printf("Ignoring \"%s\"\n", s);
}
switch(k) {
case 0:
printf("d.x = %.g = #%x %x; %s ==> ", d.x, d.u[1], d.u[0], fmt);
printf(fmt, d.x);
break;
case 1:
#ifdef __sparc
printf("x.x = %.Lg = #%x %x %x %x; %s ==> ", x.x,
x.u[0], x.u[1], x.u[2], x.u[3], fmt);
#else
printf("x.x = %.Lg = #%x %x %x %x %x; %s ==> ", x.x,
x.u[4], x.u[3], x.u[2], x.u[1], x.u[0], fmt);
#endif
printf(fmt, x.x);
#ifdef want_Quad
break;
case 2:
printf("q.x = %.Lqg = #%x %x %x %x; %s ==> ", q.x,
q.u[3], q.u[2], q.u[1], q.u[0], fmt);
printf(fmt, q.x);
#endif
}
putchar('\n');
}
return 0;
}

13
test/pftestLq.out Normal file
View File

@ -0,0 +1,13 @@
d.x = 1.23 = #3ff3ae14 7ae147ae; %.8a ==> 0x1.3ae147aep+0
d.x = 1.23 = #3ff3ae14 7ae147ae; %.7a ==> 0x1.3ae147bp+0
d.x = 1.23 = #3ff3ae14 7ae147ae; %.2a ==> 0x1.3bp+0
x.x = 1.23 = #3fff 9d70 a3d7 a3d 70a4; %.La ==> 0xap-3
q.x = 1.23 = #3fff3ae1 47ae147a e147ae14 7ae147ae; %.20Lqa ==> 0x1.3ae147ae147ae147ae14p+0
x.x = 3e+27 = #405a 9b18 ab5d f718 b6c; %La ==> 0x9.b18ab5df7180b6cp+88
x.x = 1.7e+27 = #4059 afc6 a015 291b 4024; %La ==> 0xa.fc6a015291b4024p+87
q.x = 3e+48 = #40a006be 53879565 60c1e1a9 9c13ee2; %Lqa ==> 0x1.06be5387956560c1e1a909c13ee2p+161
q.x = 1e+49 = #40a1b5e7 e08ca3a8 f6987819 baecbe22; %Lqa ==> 0x1.b5e7e08ca3a8f6987819baecbe22p+162
q.x = 1.5 = #3fff8000 0 0 0; %Lqa ==> 0x1.8p+0
q.x = 2.5 = #40004000 0 0 0; %Lqa ==> 0x1.4p+1
q.x = 4.5 = #40012000 0 0 0; %Lqa ==> 0x1.2p+2
q.x = 8.5 = #40021000 0 0 0; %Lqa ==> 0x1.1p+3

15
test/pftestQ.out Normal file
View File

@ -0,0 +1,15 @@
d.x = 1.23 = #7ae147ae 3ff3ae14; %.8a ==> 0x1.3ae147aep+0
d.x = 1.23 = #7ae147ae 3ff3ae14; %.7a ==> 0x1.3ae147bp+0
d.x = 1.23 = #7ae147ae 3ff3ae14; %.2a ==> 0x1.3bp+0
x.x = 1.23 = #3fff3ae1 47ae147a e147ae14 7ae147ae; %.La ==> 0x1p+0
x.x = 1.23 = #3fff3ae1 47ae147a e147ae14 7ae147ae; %.20La ==> 0x1.3ae147ae147ae147ae14p+0
x.x = 3e+27 = #405a3631 56bbee30 16d70000 0; %La ==> 0x1.363156bbee3016d7p+91
x.x = 1.7e+27 = #40595f8d 402a5236 80490000 0; %La ==> 0x1.5f8d402a52368049p+90
d.x = 1.5 = #0 3ff80000; %a ==> 0x1.8p+0
d.x = 2.5 = #0 40040000; %a ==> 0x1.4p+1
d.x = 4.5 = #0 40120000; %a ==> 0x1.2p+2
d.x = 8.5 = #0 40210000; %a ==> 0x1.1p+3
x.x = 1.5 = #3fff8000 0 0 0; %La ==> 0x1.8p+0
x.x = 2.5 = #40004000 0 0 0; %La ==> 0x1.4p+1
x.x = 4.5 = #40012000 0 0 0; %La ==> 0x1.2p+2
x.x = 8.5 = #40021000 0 0 0; %La ==> 0x1.1p+3

15
test/pftestnos Normal file
View File

@ -0,0 +1,15 @@
%.8a 1.23
%.7a
%.2a
%.La 1.23
%.20La 1.23
%La 3e27
%La 1.7e27
%a 1.5
2.5
4.5
8.5
%La 1.5
2.5
4.5
8.5

15
test/pftestx.out Normal file
View File

@ -0,0 +1,15 @@
d.x = 1.23 = #3ff3ae14 7ae147ae; %.8a ==> 0x1.3ae147aep+0
d.x = 1.23 = #3ff3ae14 7ae147ae; %.7a ==> 0x1.3ae147bp+0
d.x = 1.23 = #3ff3ae14 7ae147ae; %.2a ==> 0x1.3bp+0
x.x = 1.23 = #3fff 9d70 a3d7 a3d 70a4; %.La ==> 0xap-3
x.x = 1.23 = #3fff 9d70 a3d7 a3d 70a4; %.20La ==> 0x9.d70a3d70a3d70a4p-3
x.x = 3e+27 = #405a 9b18 ab5d f718 b6c; %La ==> 0x9.b18ab5df7180b6cp+88
x.x = 1.7e+27 = #4059 afc6 a015 291b 4024; %La ==> 0xa.fc6a015291b4024p+87
d.x = 1.5 = #3ff80000 0; %a ==> 0x1.8p+0
d.x = 2.5 = #40040000 0; %a ==> 0x1.4p+1
d.x = 4.5 = #40120000 0; %a ==> 0x1.2p+2
d.x = 8.5 = #40210000 0; %a ==> 0x1.1p+3
x.x = 1.5 = #3fff c000 0 0 0; %La ==> 0xcp-3
x.x = 2.5 = #4000 a000 0 0 0; %La ==> 0xap-2
x.x = 4.5 = #4001 9000 0 0 0; %La ==> 0x9p-1
x.x = 8.5 = #4002 8800 0 0 0; %La ==> 0x8.8p+0

8
test/strtodt.c
View File

@ -49,6 +49,8 @@ THIS SOFTWARE.
ULong L[2];
} U;
#define UL (unsigned long)
static int
process(char *fname, FILE *f)
{
@ -81,7 +83,7 @@ process(char *fname, FILE *f)
if (b.L[W0] != a.L[0] || b.L[W1] != a.L[1]) {
n++;
printf("Line %d of %s: got %lx %lx; expected %lx %lx\n",
line, fname, b.L[W0], b.L[W1], a.L[0], a.L[1]);
line, fname, UL b.L[W0], UL b.L[W1], UL a.L[0], UL a.L[1]);
}
}
return n;
@ -120,8 +122,8 @@ main(int argc, char **argv)
if (argc <= 1)
n = process("<stdin>", stdin);
else
while(s = *++argv)
if (f = fopen(s,"r")) {
while((s = *++argv))
if ((f = fopen(s,"r"))) {
n += process(s, f);
fclose(f);
}

24
test/testnos3
View File

@ -326,3 +326,27 @@
# first is 2^-1075 (half the smallest denormal)
2.4703282292062327208828439643411068618252990130716238221279284125033775363510437593264991818081799618989828234772285886546332835517796989819938739800539093906315035659515570226392290858392449105184435931802849936536152500319370457678249219365623669863658480757001585769269903706311928279558551332927834338409351978015531246597263579574622766465272827220056374006485499977096599470454020828166226237857393450736339007967761930577506740176324673600968951340535537458516661134223766678604162159680461914467291840300530057530849048765391711386591646239524912623653881879636239373280423891018672348497668235089863388587925628302755995657524455507255189313690836254779186948667994968324049705821028513185451396213837722826145437693412532098591327667236328125e-324 0 0
2.47032822920623272e-324 0 0
# examples reported by Mark Dickinson of bugs in the bigcomp() logic introduced
# 20090316 in dtoa.c to speed handling of absurdly long input:
12579816049008305546974391768996369464963024663104e-357 90bbd 7412d19f
17489628565202117263145367596028389348922981857013e-357 c938e 9000492f
18487398785991994634182916638542680759613590482273e-357 d4b3a ee198863
32002864200581033134358724675198044527469366773928e-358 24d1e ed8448e3
99999999999999994487665465554760717039532578546e-47 3ff00000 0
1.0000000000000000100000000000000000000001e44 4911efc6 59cf7d4c
1000000000000000000000000000000000000000e-16 44b52d02 c7e14af6
10000000000000000000000000000000000000000e-17 44b52d02 c7e14af6
10.900000000000000012345678912345678912345 4025cccc cccccccd
104308485241983990666713401708072175773165034278685682646111762292409330928739751702404658197872319129036519947435319418387839758990478549477777586673075945844895981012024387992135617064532141489278815239849108105951619997829153633535314849999674266169258928940692239684771590065027025835804863585454872499320500023126142553932654370362024104462255244034053203998964360882487378334860197725139151265590832887433736189468858614521708567646743455601905935595381852723723645799866672558576993978025033590728687206296379801363024094048327273913079612469982585674824156000783167963081616214710691759864332339239688734656548790656486646106983450809073750535624894296242072010195710276073042036425579852459556183541199012652571123898996574563824424330960027873516082763671875e-1075 78026 65fd9600
99037485700245683102805043437346965248029601286431e-373 0 2
99617639833743863161109961162881027406769510558457e-373 0 2
98852915025769345295749278351563179840130565591462e-372 0 14
99059944827693569659153042769690930905148015876788e-373 0 2
98914979205069368270421829889078356254059760327101e-372 0 14
0.999999999999999999999999999999999999999999999e23 44b52d02 c7e14af6
991633793189150720000000000000000000000000000000000000000e-33 44ea3f92 6bad90c6
37652435753827922121470370984740152789920e234 78f1667a c9e75d61
999999999999999996790597280027956716285163e-42 3ff00000 0
9483973038658180570348795755328802873667739881500874740826641664593613312413122937394311083577538394191754403820631172036846773125424639263833553383990195662207006139342261292777056851379062046720e0 68a03d69 82f2f936
20209005503919489280000000000000000000000000000000000000000e-40 43bc0bae 57e880e6

8
test/xLtest.c
View File

@ -83,7 +83,7 @@ main(Void)
int dItry, i, ndig = 0, r = 1;
union { long double d; ULong bits[3]; } u, v[2];
while(s = fgets(ibuf, sizeof(ibuf), stdin)) {
while((s = fgets(ibuf, sizeof(ibuf), stdin))) {
while(*s <= ' ')
if (!*s++)
continue;
@ -94,7 +94,7 @@ main(Void)
continue;
case 'n':
i = s[1];
if (i <= ' ' || i >= '0' && i <= '9') {
if (i <= ' ' || (i >= '0' && i <= '9')) {
ndig = atoi(s+1);
continue;
}
@ -109,8 +109,8 @@ main(Void)
u.bits[_2] = (ULong)strtoul(s1=se, &se, 16);
}
printf("\nInput: %s", ibuf);
printf(" --> f = #%lx %lx %lx\n", u.bits[_0],
u.bits[_1], u.bits[_2]);
printf(" --> f = #%lx %lx %lx\n", U u.bits[_0],
U u.bits[_1], U u.bits[_2]);
goto fmt_test;
}
dItry = 1;

22
test/xQtest.c
View File

@ -31,16 +31,28 @@ THIS SOFTWARE.
int
main(void)
{
union { long double d; unsigned int bits[4]; } u, w;
switch(sizeof(long double)) {
case 16:
w.bits[0] = w.bits[3] = 0;
w.d = 1.;
u.d = 3.;
w.d = w.d / u.d;
if (w.bits[0] && w.bits[3])
printf("cp x.ou0 x.out; cp xL.ou0 xL.out;"
" cp Q.ou1 Q.out; cp pftestQ.out pftest.out\n");
else
printf("cp x.ou0 x.out; cp xL.ou0 xL.out;"
" cp Q.ou0 Q.out; cp pftestx.out pftest.out\n");
break;
case 10:
case 12:
printf("cp x.ou1 x.out; cp xL.ou1 xL.out; cp Q.ou0 Q.out\n");
break;
case 16:
printf("cp x.ou0 x.out; cp xL.ou0 xL.out; cp Q.ou1 Q.out\n");
printf("cp x.ou1 x.out; cp xL.ou1 xL.out; cp Q.ou0 Q.out;"
" cp pftestx.out pftest.out\n");
break;
default:
printf("cp x.ou0 x.out; cp xL.ou0 xL.out; cp Q.ou0 Q.out\n");
printf("cp x.ou0 x.out; cp xL.ou0 xL.out; cp Q.ou0 Q.out;"
" cp pftestx.out pftest.out\n");
}
return 0;
}

54
test/xsum0.out
View File

@ -1,35 +1,41 @@
README e3adb571 3095
Qtest.c e8353ffc 5046
README efba0d5d 3412
Q.ou0 e4592b85 28742
Q.ou1 ea0b344d 39572
Qtest.c efe2b3f4 5116
d.out f271efc9 28131
dI.out d522eef 4369
dIsi.out 1dd6d02f 4350
dItest.c e33800ce 2371
ddtest.c f9d06e7b 4984
dtest.c ee533ac3 4078
dt.c 7eeda57 6384
ftest.c ec8a6654 3999
getround.c fe659fe7 2503
dd.out e262456e 40923
ddsi.out 1f94bbe2 10251
ddtest.c ef71cbf3 4986
dt.c e562c302 6808
dtest.c 9a5d01 4080
dtst.out e284ac98 23711
f.out 9013e91 21537
ftest.c 1c824a88 4001
getround.c ee95ed1 2502
makefile f714a641 5634
pfLqtestnos ffb9723 99
pftest.c ea314d7f 3452
pftestQ.out 198434fd 830
pftestx.out 1ccea5dd 788
pftestLq.out 1691dbfc 845
pftestnos ecbc9be6 101
rtestnos f94bcdf6 336
strtoIdSI.c 7bfb88b 49
strtoIddSI.c 72e8852 50
strtodISI.c ed08b740 49
strtodt.c aaf94bc 3330
strtodt.c f1aa53af 3374
strtopddSI.c 13e7138d 50
strtorddSI.c f7e4b1d5 50
xLtest.c f3f96ad1 4833
xQtest.c efdea3a2 1549
xtest.c ee81e661 4830
rtestnos f94bcdf6 336
testnos e89999d6 485
testnos1 7e16229 294
testnos3 fa5c8aca 11998
dI.out d522eef 4369
dIsi.out 1dd6d02f 4350
ddsi.out 1f94bbe2 10251
dd.out e262456e 40923
dtst.out e284ac98 23711
d.out f271efc9 28131
f.out 9013e91 21537
testnos3 f5ae7ef3 14403
x.ou0 1402f834 25372
xL.ou0 faa3a741 26363
x.ou1 f1af5a00 34581
xL.ou0 faa3a741 26363
xL.ou1 e349e5c 37165
Q.ou0 e4592b85 28742
Q.ou1 ea0b344d 39572
makefile 13d36c85 5148
xLtest.c ee11f673 4843
xQtest.c efbe29be 1912
xtest.c ec9f5deb 4834

4
test/xtest.c
View File

@ -85,7 +85,7 @@ main(Void)
int i, dItry, ndig = 0, r = 1;
union { long double d; UShort bits[5]; } u, v[2];
while(s = fgets(ibuf, sizeof(ibuf), stdin)) {
while((s = fgets(ibuf, sizeof(ibuf), stdin))) {
while(*s <= ' ')
if (!*s++)
continue;
@ -96,7 +96,7 @@ main(Void)
continue;
case 'n':
i = s[1];
if (i <= ' ' || i >= '0' && i <= '9') {
if (i <= ' ' || (i >= '0' && i <= '9')) {
ndig = atoi(s+1);
continue;
}

20
ulp.c
View File

@ -34,13 +34,13 @@ THIS SOFTWARE.
double
ulp
#ifdef KR_headers
(x) double x;
(x) U *x;
#else
(double x)
(U *x)
#endif
{
Long L;
double a;
U a;
L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
#ifndef Sudden_Underflow
@ -49,22 +49,22 @@ ulp
#ifdef IBM
L |= Exp_msk1 >> 4;
#endif
word0(a) = L;
word1(a) = 0;
word0(&a) = L;
word1(&a) = 0;
#ifndef Sudden_Underflow
}
else {
L = -L >> Exp_shift;
if (L < Exp_shift) {
word0(a) = 0x80000 >> L;
word1(a) = 0;
word0(&a) = 0x80000 >> L;
word1(&a) = 0;
}
else {
word0(a) = 0;
word0(&a) = 0;
L -= Exp_shift;
word1(a) = L >= 31 ? 1 : 1 << 31 - L;
word1(&a) = L >= 31 ? 1 : 1 << (31 - L);
}
}
#endif
return a;
return dval(&a);
}

57
xsum0.out
View File

@ -1,49 +1,52 @@
README ee32ce1c 15386
README 664d441 15750
arithchk.c ebbe5bc7 4075
dmisc.c c8daa18 4682
dtoa.c 14f5b9e1 17138
dtoa.c f6ed7306 17542
g_Qfmt.c f6aad16c 2926
g__fmt.c 1b8e8cd8 3401
g_ddfmt.c eb4d0889 4090
g_dfmt.c 10d694eb 2584
g__fmt.c e961947e 3401
g_ddfmt.c 1f75ef12 4161
g_dfmt.c 24be3a8 2629
g_ffmt.c 197652f1 2516
g_xLfmt.c 3fd29c5 2773
g_xfmt.c f6a580a 2862
gdtoa.c eabc1bc5 17024
gdtoa.h f3c171cc 4945
gdtoa.c e4651ba8 17355
gdtoa.h f966845d 4944
gdtoa_fltrnds.h 1aaf5112 421
gdtoaimp.h 490a372 19893
gethex.c e7327648 7106
gdtoaimp.h 395c83e 20044
gethex.c 144ecbb4 7107
gmisc.c 1859d016 2084
hd_init.c efdbe921 1797
hexnan.c 13f362b7 3462
makefile 3cd7ae 2933
misc.c 1757f7fc 14252
hexnan.c 891b04e 3468
makefile ffae86d0 3148
misc.c 1d84f517 14416
printf.c 15c97ab 158
printf.c0 efece124 29696
qnan.c efd33d64 3417
smisc.c 1064c213 3694
smisc.c feb7290 3699
stdio1.h e6e5dc11 2909
strtoIQ.c 1809dfcf 1939
strtoId.c f41ddac2 1931
strtoIdd.c f13e3bc3 2105
strtoIf.c f12c6af4 1875
strtoIg.c fd8c1bcf 3589
strtoIg.c fe47bc72 3591
strtoIx.c e50f716d 1960
strtoIxL.c ea0b821b 1931
strtod.c e7e4addc 21731
strtodI.c 1c2440ce 3915
strtodg.c f74828cc 21164
strtod.c 1dabf6 22744
strtodI.c f1dbb0af 4004
strtodg.c eacc7424 21188
strtodnrp.c af895e9 2538
strtof.c 1db524d 2155
strtof.c e1cb24dd 2161
strtopQ.c f244b012 2645
strtopd.c 1c3c9ce 1751
strtopdd.c 125d6c19 4580
strtopf.c 1939f590 2149
strtopx.c f78a816e 2700
strtopxL.c b77b701 2455
strtopdd.c 177251c9 4640
strtopf.c f305c6ea 2155
strtopx.c 1529e1f5 2710
strtopxL.c 69418e3 2469
strtorQ.c 9360a0b 2885
strtord.c af5c50e 2491
strtordd.c 1b266865 4936
strtorf.c f0d86e2b 2396
strtorx.c f19a56af 2947
strtorxL.c 167fe87c 2704
strtordd.c e3d20ab 4992
strtorf.c fbe2ea18 2402
strtorx.c 168faee5 2957
strtorxL.c 42c6cec 2718
sum.c f525bad9 2494
ulp.c 1e2e148f 1864
ulp.c edb56c0b 1866