Merge gdtoa-20110304.

31 files changed, 3183 insertions, 646 deletions

diff --git a/contrib/gdtoa/README b/contrib/gdtoa/README
index ce8be55..0c034f1 100644
--- a/contrib/gdtoa/README
+++ b/contrib/gdtoa/README
@@ -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 ".").
diff --git a/contrib/gdtoa/changes b/contrib/gdtoa/changes
new file mode 100644
index 0000000..3eb8138
--- /dev/null
+++ b/contrib/gdtoa/changes
@@ -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.
diff --git a/contrib/gdtoa/dtoa.c b/contrib/gdtoa/dtoa.c
index 48fdf5e..6e5de30 100644
--- a/contrib/gdtoa/dtoa.c
+++ b/contrib/gdtoa/dtoa.c
@@ -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)
diff --git a/contrib/gdtoa/g__fmt.c b/contrib/gdtoa/g__fmt.c
index fbccb7d..6c4b3d2 100644
--- a/contrib/gdtoa/g__fmt.c
+++ b/contrib/gdtoa/g__fmt.c
@@ -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;
 			}
diff --git a/contrib/gdtoa/g_ddfmt.c b/contrib/gdtoa/g_ddfmt.c
index b65d39d..1256655 100644
--- a/contrib/gdtoa/g_ddfmt.c
+++ b/contrib/gdtoa/g_ddfmt.c
@@ -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);
diff --git a/contrib/gdtoa/g_dfmt.c b/contrib/gdtoa/g_dfmt.c
index 23d8b24..8367868 100644
--- a/contrib/gdtoa/g_dfmt.c
+++ b/contrib/gdtoa/g_dfmt.c
@@ -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);
 	}
diff --git a/contrib/gdtoa/gdtoa.c b/contrib/gdtoa/gdtoa.c
index 3b64250..e15bb2a 100644
--- a/contrib/gdtoa/gdtoa.c
+++ b/contrib/gdtoa/gdtoa.c
@@ -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')
diff --git a/contrib/gdtoa/gdtoa.h b/contrib/gdtoa/gdtoa.h
index e59ebf6..f004ec3 100644
--- a/contrib/gdtoa/gdtoa.h
+++ b/contrib/gdtoa/gdtoa.h
@@ -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;
diff --git a/contrib/gdtoa/gdtoaimp.h b/contrib/gdtoa/gdtoaimp.h
index 9991ffa..75cdc4e 100644
--- a/contrib/gdtoa/gdtoaimp.h
+++ b/contrib/gdtoa/gdtoaimp.h
@@ -96,7 +96,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"
@@ -115,7 +120,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,
@@ -164,11 +174,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.
  */
 
@@ -287,25 +292,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
@@ -419,6 +413,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
@@ -659,7 +658,7 @@ extern void memcpy_D2A ANSI((void*, const void*, size_t));
  extern int strtorxL ANSI((CONST char *, char **, int, void *));
  extern Bigint *sum ANSI((Bigint*, Bigint*));
  extern int trailz ANSI((Bigint*));
- extern double ulp ANSI((double));
+ extern double ulp ANSI((U*));
 
 #ifdef __cplusplus
 }
diff --git a/contrib/gdtoa/gethex.c b/contrib/gdtoa/gethex.c
index c8ac8ba..a9982c9 100644
--- a/contrib/gdtoa/gethex.c
+++ b/contrib/gdtoa/gethex.c
@@ -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;
diff --git a/contrib/gdtoa/hexnan.c b/contrib/gdtoa/hexnan.c
index 58952db..a443721 100644
--- a/contrib/gdtoa/hexnan.c
+++ b/contrib/gdtoa/hexnan.c
@@ -29,10 +29,6 @@ THIS SOFTWARE.
 /* Please send bug reports to David M. Gay (dmg at acm dot org,
  * with " at " changed at "@" and " dot " changed to ".").	*/
 
-/* $FreeBSD$ */
-
-#include <ctype.h>
-
 #include "gdtoaimp.h"
 
  static void
@@ -75,14 +71,14 @@ hexnan( CONST char **sp, FPI *fpi, ULong *x0)
 	x1 = xe = x;
 	havedig = hd0 = i = 0;
 	s = *sp;
-
-	/* FreeBSD local: Accept (but ignore) the '0x' prefix. */
-	if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X'))
+	/* allow optional initial 0x or 0X */
+	while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
+		++s;
+	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 0
 			if (c <= ' ') {
 				if (hd0 < havedig) {
 					if (x < x1 && i < 8)
@@ -96,14 +92,26 @@ hexnan( CONST char **sp, FPI *fpi, ULong *x0)
 					x1 = x;
 					i = 0;
 					}
+				while(*(CONST unsigned char*)(s+1) <= ' ')
+					++s;
+				if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')
+				 && *(CONST unsigned char*)(s+3) > ' ')
+					s += 2;
 				continue;
 				}
 			if (/*(*/ c == ')' && havedig) {
 				*sp = s + 1;
 				break;
 				}
+#ifndef GDTOA_NON_PEDANTIC_NANCHECK
+			do {
+				if (/*(*/ c == ')') {
+					*sp = s + 1;
+					break;
+					}
+				} while((c = *++s));
 #endif
-			break;
+			return STRTOG_NaN;
 			}
 		havedig++;
 		if (++i > 8) {
@@ -112,9 +120,11 @@ hexnan( CONST char **sp, FPI *fpi, ULong *x0)
 			i = 1;
 			*--x = 0;
 			}
-		*x = (*x << 4) | h & 0xf;
+		*x = (*x << 4) | (h & 0xf);
 		}
-	if (havedig && x < x1 && i < 8)
+	if (!havedig)
+		return STRTOG_NaN;
+	if (x < x1 && i < 8)
 		L_shift(x, x1, i);
 	if (x > x0) {
 		x1 = x0;
@@ -128,7 +138,6 @@ hexnan( CONST char **sp, FPI *fpi, ULong *x0)
 		if ( (i = nbits & (ULbits-1)) !=0)
 			*xe &= ((ULong)0xffffffff) >> (ULbits - i);
 		}
-	if (havedig) {
 	for(x1 = xe;; --x1) {
 		if (*x1 != 0)
 			break;
@@ -137,22 +146,5 @@ hexnan( CONST char **sp, FPI *fpi, ULong *x0)
 			break;
 			}
 		}
-	}
-
-	/*
-	 * FreeBSD local: Accept all the sequences allowed by C99 and update
-	 * the tail pointer correctly. Don't accept any invalid sequences.
-	 */
-	if (c == '\0')	/* nan() calls this, too; tolerate a missing ')' */
-		return STRTOG_NaNbits;
-	if (c != ')') {
-		while(c = *(CONST unsigned char*)++s) {
-			if (c == ')')
-				break;
-			if (!isalnum(c) && c != '_')
-				return STRTOG_NaNbits;
-		}
-	}
-	*sp = s + 1;
 	return STRTOG_NaNbits;
 	}
diff --git a/contrib/gdtoa/makefile b/contrib/gdtoa/makefile
index 367eb9c..b1f18cd 100644
--- a/contrib/gdtoa/makefile
+++ b/contrib/gdtoa/makefile
@@ -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
diff --git a/contrib/gdtoa/misc.c b/contrib/gdtoa/misc.c
index 8d2888e..e5f7b04 100644
--- a/contrib/gdtoa/misc.c
+++ b/contrib/gdtoa/misc.c
@@ -92,7 +92,11 @@ Bfree
 {
 	if (v) {
 		if (v->k > Kmax)
+#ifdef FREE
+			FREE((void*)v);
+#else
 			free((void*)v);
+#endif
 		else {
 			ACQUIRE_DTOA_LOCK(0);
 			v->next = freelist[v->k];
@@ -110,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)
@@ -210,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;
@@ -618,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;
@@ -636,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;
@@ -669,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
@@ -680,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
@@ -693,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
@@ -721,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
@@ -732,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
@@ -794,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;
@@ -847,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;
 	}
@@ -861,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;
diff --git a/contrib/gdtoa/printf.c b/contrib/gdtoa/printf.c
new file mode 100644
index 0000000..02f1f3c
--- /dev/null
+++ b/contrib/gdtoa/printf.c
@@ -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"
diff --git a/contrib/gdtoa/printf.c0 b/contrib/gdtoa/printf.c0
new file mode 100644
index 0000000..c33bcc1
--- /dev/null
+++ b/contrib/gdtoa/printf.c0
@@ -0,0 +1,1669 @@
+/****************************************************************
+Copyright (C) 1997, 1999, 2001 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.
+****************************************************************/
+
+/* This implements most of ANSI C's printf, fprintf, and sprintf,
+ * omitting L, with %.0g and %.0G giving the shortest decimal string
+ * that rounds to the number being converted, and with negative
+ * precisions allowed for %f.
+ */
+
+#ifdef KR_headers
+#include "varargs.h"
+#else
+#include "stddef.h"
+#include "stdarg.h"
+#include "stdlib.h"
+#endif
+
+#ifdef Use_GDTOA_for_i386_long_double /*{{*/
+#include "gdtoa.h"
+#else /*}{*/
+#ifndef NO_PRINTF_A_FMT /*{*/
+#include "gdtoa.h"
+#endif /*}*/
+#endif /*}}*/
+
+#ifdef __i386
+#define NO_GDTOA_i386_Quad
+#endif
+
+#ifdef Use_GDTOA_for_i386_long_double /*{*/
+#ifndef NO_GDTOA_i386_Quad /*{*/
+#define GDTOA_both
+#define Use_GDTOA_Qtype
+#ifdef __ICC__ /* or __INTEL_COMPILER__ or __INTEL_COMPILER ?? */
+#define GDTOA_Qtype _Quad
+#else
+#define GDTOA_Qtype __float128
+#endif
+#endif /*} NO_GDTOA_i386_Quad */
+#endif /*} Use_GDTOA_for_i386_long_double */
+
+#ifdef Use_GDTOA_Qtype /*{*/
+#ifndef GDTOA_H_INCLUDED
+#include "gdtoa.h"
+#endif
+#ifndef GDTOA_Qtype
+#define GDTOA_Qtype long double
+#endif
+#endif /*}*/
+
+#ifndef GDTOA_H_INCLUDED /*{*/
+
+ enum {	/* return values from strtodg */
+	STRTOG_Zero	= 0,
+	STRTOG_Normal	= 1,
+	STRTOG_Denormal	= 2,
+	STRTOG_Infinite	= 3,
+	STRTOG_NaN	= 4,
+	STRTOG_NaNbits	= 5,
+	STRTOG_NoNumber	= 6,
+	STRTOG_Retmask	= 7};
+
+ typedef struct
+FPI {
+	int nbits;
+	int emin;
+	int emax;
+	int rounding;
+	int sudden_underflow;
+	} FPI;
+
+enum {	/* FPI.rounding values: same as FLT_ROUNDS */
+	FPI_Round_zero = 0,
+	FPI_Round_near = 1,
+	FPI_Round_up = 2,
+	FPI_Round_down = 3
+	};
+#endif /*}*/
+
+#ifdef NO_PRINTF_A_FMT /*{{*/
+#define WANT_A_FMT(x) /*nothing*/
+#else /*}{*/
+#define WANT_A_FMT(x) x
+#endif /*}}*/
+
+#include "stdio1.h"
+#include "string.h"
+#include "errno.h"
+
+ typedef struct
+Finfo {
+	union { FILE *cf; char *sf; } u;
+	char *ob0, *obe1;
+	size_t lastlen;
+	} Finfo;
+
+ typedef char *(*pgdtoa) ANSI((FPI*, int be, ULong *bits, int *kind, int mode, int ndigits, int *decpt, char **rve));
+
+ typedef struct
+FPBits {
+	ULong bits[4];	/* sufficient for quad; modify if considering wider types */
+	FPI *fpi;
+	pgdtoa gdtoa;
+	int sign;
+	int ex;	/* exponent */
+	int kind;
+	} FPBits;
+
+ typedef union U
+{
+	double d;
+	long double ld;
+#ifdef GDTOA_Qtype
+	GDTOA_Qtype Qd;
+#endif
+	unsigned int ui[4];
+	unsigned short us[5];
+	} U;
+
+ typedef char *(*Putfunc) ANSI((Finfo*, int*));
+ typedef void (*Fpbits) ANSI((U*, FPBits*));
+
+/* Would have preferred typedef void (*Fpbits)(va_list*, FPBits*)
+ * but gcc is buggy in this regard.
+ */
+
+#ifdef Use_GDTOA_for_i386_long_double /*{*/
+
+#ifdef IEEE_MC68k
+#define _0 0
+#define _1 1
+#define _2 2
+#define _3 3
+#define _4 4
+#endif
+#ifdef IEEE_8087
+#define _0 4
+#define _1 3
+#define _2 2
+#define _3 1
+#define _4 0
+#endif
+
+ static void
+xfpbits(U *u, FPBits *b)
+{
+	ULong *bits;
+	int ex, i;
+	static FPI fpi0 = { 64, 1-16383-64+1, 32766 - 16383 - 64 + 1, 1, 0 };
+
+	b->fpi = &fpi0;
+	b->gdtoa = gdtoa;
+	b->sign = u->us[_0] & 0x8000;
+	bits = b->bits;
+	bits[1] = (u->us[_1] << 16) | u->us[_2];
+	bits[0] = (u->us[_3] << 16) | u->us[_4];
+	if ( (ex = u->us[_0] & 0x7fff) !=0) {
+		i = STRTOG_Normal;
+		if (ex == 0x7fff)
+			/* Infinity or NaN */
+			i = bits[0] | bits[1] ? STRTOG_NaN : STRTOG_Infinite;
+		}
+	else if (bits[0] | bits[1]) {
+		i = STRTOG_Denormal;
+		ex = 1;
+		}
+	else
+		i = STRTOG_Zero;
+	b->kind = i;
+	b->ex = ex - (0x3fff + 63);
+	}
+
+#undef _0
+#undef _1
+#undef _2
+#undef _3
+#undef _4
+#define GDTOA_LD_fpbits xfpbits
+#endif /*} Use_GDTOA_for_i386_long_double */
+
+#ifdef Use_GDTOA_Qtype /*{*/
+#include "gdtoa.h"
+#ifndef GDTOA_Qtype
+#define GDTOA_Qtype long double
+#endif
+#ifdef GDTOA_LD_fpbits
+#define GDTOA_Q_fpbits Qfpbits
+#else
+#define GDTOA_LD_fpbits Qfpbits
+#endif
+
+#ifdef IEEE_MC68k
+#define _0 0
+#define _1 1
+#define _2 2
+#define _3 3
+#endif
+#ifdef IEEE_8087
+#define _0 3
+#define _1 2
+#define _2 1
+#define _3 0
+#endif
+
+ static void
+Qfpbits(U *u, FPBits *b)
+{
+	ULong *bits;
+	int ex, i;
+	static FPI fpi0 = { 113, 1-16383-113+1, 32766 - 16383 - 113 + 1, 1, 0 };
+
+	b->fpi = &fpi0;
+	b->gdtoa = gdtoa;
+	b->sign = u->ui[_0] & 0x80000000L;
+	bits = b->bits;
+	bits[3] = u->ui[_0] & 0xffff;
+	bits[2] = u->ui[_1];
+	bits[1] = u->ui[_2];
+	bits[0] = u->ui[_3];
+	if ( (ex = (u->ui[_0] & 0x7fff0000L) >> 16) !=0) {
+		if (ex == 0x7fff) {
+			/* Infinity or NaN */
+			i = bits[0] | bits[1] | bits[2] | bits[3]
+				? STRTOG_NaN : STRTOG_Infinite;
+			}
+		else {
+			i = STRTOG_Normal;
+			bits[3] |= 0x10000;
+			}
+		}
+	else if (bits[0] | bits[1] | bits[2] | bits[3]) {
+		i = STRTOG_Denormal;
+		ex = 1;
+		}
+	else
+		i = STRTOG_Zero;
+	b->kind = i;
+	b->ex = ex - (0x3fff + 112);
+	}
+
+#undef _0
+#undef _1
+#undef _2
+#undef _3
+#endif /*} GDTOA_Qtype */
+
+#ifdef KR_headers
+#define Const /* const */
+#define Voidptr char*
+#ifndef size_t__
+#define size_t int
+#define size_t__
+#endif
+
+#else
+
+#define Const const
+#define Voidptr void*
+
+#endif
+
+#undef MESS
+#ifndef Stderr
+#define Stderr stderr
+#endif
+
+#ifdef _windows_
+#undef PF_BUF
+#define MESS
+#include "mux0.h"
+#define stdout_or_err(f) (f == stdout)
+#else
+#define stdout_or_err(f) (f == Stderr || f == stdout)
+#endif
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+ extern char *dtoa ANSI((double, int, int, int*, int*, char **));
+ extern void freedtoa ANSI((char*));
+
+
+
+#ifdef USE_ULDIV
+/* This is for avoiding 64-bit divisions on the DEC Alpha, since */
+/* they are not portable among variants of OSF1 (DEC's Unix). */
+
+#define ULDIV(a,b) uldiv_ASL(a,(unsigned long)(b))
+
+#ifndef LLBITS
+#define LLBITS 6
+#endif
+#ifndef ULONG
+#define ULONG unsigned long
+#endif
+
+ static int
+klog(ULONG x)
+{
+	int k, rv = 0;
+
+	if (x > 1L)
+	    for(k = 1 << LLBITS-1;;) {
+		if (x >= (1L << k)) {
+			rv |= k;
+			x >>= k;
+			}
+		if (!(k >>= 1))
+			break;
+		}
+	return rv;
+	}
+
+ ULONG
+uldiv_ASL(ULONG a, ULONG b)
+{
+	int ka;
+	ULONG c, k;
+	static ULONG b0;
+	static int kb;
+
+	if (a < b)
+		return 0;
+	if (b != b0) {
+		b0 = b;
+		kb = klog(b);
+		}
+	k = 1;
+	if ((ka = klog(a) - kb) > 0) {
+		k <<= ka;
+		b <<= ka;
+		}
+	c = 0;
+	for(;;) {
+		if (a >= b) {
+			a -= b;
+			c |= k;
+			}
+		if (!(k >>= 1))
+			break;
+		a <<= 1;
+		}
+	return c;
+	}
+
+#else
+#define ULDIV(a,b) a / b
+#endif /* USE_ULDIV */
+
+#ifdef PF_BUF
+FILE *stderr_ASL = (FILE*)&stderr_ASL;
+void (*pfbuf_print_ASL) ANSI((char*));
+char *pfbuf_ASL;
+static char *pfbuf_next;
+static size_t pfbuf_len;
+extern Char *mymalloc_ASL ANSI((size_t));
+extern Char *myralloc_ASL ANSI((void *, size_t));
+
+#undef fflush
+#ifdef old_fflush_ASL
+#define fflush old_fflush_ASL
+#endif
+
+ void
+fflush_ASL(FILE *f)
+{
+	if (f == stderr_ASL) {
+		if (pfbuf_ASL && pfbuf_print_ASL) {
+			(*pfbuf_print_ASL)(pfbuf_ASL);
+			free(pfbuf_ASL);
+			pfbuf_ASL = 0;
+			}
+		}
+	else
+		fflush(f);
+	}
+
+ static void
+pf_put(char *buf, int len)
+{
+	size_t x, y;
+	if (!pfbuf_ASL) {
+		x = len + 256;
+		if (x < 512)
+			x = 512;
+		pfbuf_ASL = pfbuf_next = (char*)mymalloc_ASL(pfbuf_len = x);
+		}
+	else if ((y = (pfbuf_next - pfbuf_ASL) + len) >= pfbuf_len) {
+		x = pfbuf_len;
+		while((x <<= 1) <= y);
+		y = pfbuf_next - pfbuf_ASL;
+		pfbuf_ASL = (char*)myralloc_ASL(pfbuf_ASL, x);
+		pfbuf_next = pfbuf_ASL + y;
+		pfbuf_len = x;
+		}
+	memcpy(pfbuf_next, buf, len);
+	pfbuf_next += len;
+	*pfbuf_next = 0;
+	}
+
+ static char *
+pfput(Finfo *f, int *rvp)
+{
+	int n;
+	char *ob0 = f->ob0;
+	*rvp += n = (int)(f->obe1 - ob0);
+	pf_put(ob0, n);
+	return ob0;
+	}
+#endif /* PF_BUF */
+
+ static char *
+Fput
+#ifdef KR_headers
+	(f, rvp) Finfo *f; int *rvp;
+#else
+	(Finfo *f, int *rvp)
+#endif
+{
+	char *ob0 = f->ob0;
+
+	*rvp += f->obe1 - ob0;
+	*f->obe1 = 0;
+	fputs(ob0, f->u.cf);
+	return ob0;
+	}
+
+
+#ifdef _windows_
+int stdout_fileno_ASL = 1;
+
+ static char *
+Wput
+#ifdef KR_headers
+	(f, rvp) Finfo *f; int *rvp;
+#else
+	(Finfo *f, int *rvp)
+#endif
+{
+	char *ob0 = f->ob0;
+
+	*rvp += f->obe1 - ob0;
+	*f->obe1 = 0;
+	mwrite(ob0, f->obe1 - ob0);
+	return ob0;
+	}
+#endif /*_windows_*/
+
+
+#ifdef IEEE_MC68k
+#define _0 0
+#define _1 1
+#endif
+#ifdef IEEE_8087
+#define _0 1
+#define _1 0
+#endif
+
+ static void
+dfpbits(U *u, FPBits *b)
+{
+	ULong *bits;
+	int ex, i;
+	static FPI fpi0 = { 53, 1-1023-53+1, 2046-1023-53+1, 1, 0 };
+
+	b->fpi = &fpi0;
+	b->gdtoa = gdtoa;
+	b->sign = u->ui[_0] & 0x80000000L;
+	bits = b->bits;
+	bits[1] = u->ui[_0] & 0xfffff;
+	bits[0] = u->ui[_1];
+	if ( (ex = (u->ui[_0] & 0x7ff00000L) >> 20) !=0) {
+		if (ex == 0x7ff) {
+			/* Infinity or NaN */
+			i = bits[0] | bits[1] ? STRTOG_NaN : STRTOG_Infinite;
+			}
+		else {
+			i = STRTOG_Normal;
+			bits[1] |= 0x100000;
+			}
+		}
+	else if (bits[0] | bits[1]) {
+		i = STRTOG_Denormal;
+		ex = 1;
+		}
+	else
+		i = STRTOG_Zero;
+	b->kind = i;
+	b->ex = ex - (0x3ff + 52);
+	}
+
+#undef _0
+#undef _1
+
+#ifdef Honor_FLT_ROUNDS /*{{*/
+#ifdef Trust_FLT_ROUNDS /*{{*/
+#define RoundCheck if (Rounding == -1) Rounding = Flt_Rounds; if (Rounding != 1){\
+	fpi1 = *fpb.fpi; fpi1.rounding = Rounding; fpb.fpi = &fpi1;}
+#else /*}{*/
+#define RoundCheck if (Rounding == -1) { Rounding = 1; switch((fegetround()) {\
+	  case FE_TOWARDZERO:	Rounding = 0; break;\
+	  case FE_UPWARD:	Rounding = 2; break;\
+	  case FE_DOWNWARD:	Rounding = 3; }}\
+	if (Rounding != 1){\
+	fpi1 = *fpb.fpi; fpi1.rounding = Rounding; fpb.fpi = &fpi1;}
+#endif /*}}*/
+#else /*}{*/
+#define RoundCheck /*nothing*/
+#endif /*}}*/
+
+#ifndef NO_PRINTF_A_FMT /*{*/
+ static int
+fpiprec(FPBits *b)	/* return number of hex digits minus 1, or 0 for zero */
+{
+	FPI *fpi;
+	ULong *bits;
+	int i, j, k, m;
+
+	if (b->kind == STRTOG_Zero)
+		return b->ex = 0;
+	fpi = b->fpi;
+	bits = b->bits;
+	for(k = (fpi->nbits - 1) >> 2; k > 0; --k)
+		if ((bits[k >> 3] >> 4*(k & 7)) & 0xf) {
+			m = k >> 3;
+			for(i = 0; i <= m; ++i)
+				if (bits[i]) {
+					if (i > 0) {
+						k -= 8*i;
+						b->ex += 32*i;
+						for(j = i; j <= m; ++j)
+							bits[j-i] = bits[j];
+						}
+					break;
+					}
+			for(i = 0; i < 28 && !((bits[0] >> i) & 0xf); i += 4);
+			if (i) {
+				b->ex += i;
+				m = k >> 3;
+				k -= (i >> 2);
+				for(j = 0;;++j) {
+					bits[j] >>= i;
+					if (j == m)
+						break;
+					bits[j] |= bits[j+1] << (32 - i);
+					}
+				}
+			break;
+			}
+	return k;
+	}
+
+ static int
+bround(FPBits *b, int prec, int prec1)	/* round to prec hex digits after the "." */
+{					/* prec1 = incoming precision (after ".") */
+	FPI *fpi = b->fpi;
+	ULong *bits, t;
+	int i, inc, j, k, m, n;
+#ifdef Honor_FLT_ROUNDS
+	int rounding = fpi->rounding;
+
+	if (rounding > FPI_Round_near && b->sign)
+		rounding = FPI_Round_up + FPI_Round_down - rounding;
+	if (rounding == FPI_Round_down)
+		rounding = FPI_Round_zero;
+#endif
+	m = prec1 - prec;
+	bits = b->bits;
+	inc = 0;
+#ifdef Honor_FLT_ROUNDS
+	switch(rounding) {
+		case FPI_Round_up:
+			for(i = 0; i < m; i += 8)
+				if (bits[i>>3])
+					goto inc1;
+			if ((j = i - m) > 0 && bits[(i-8)>>3] << j*4)
+				goto inc1;
+			break;
+		case FPI_Round_near:
+#endif
+			k = m - 1;
+			if ((t = bits[k >> 3] >> (j = (k&7)*4)) & 8) {
+				if (t & 7)
+					goto inc1;
+				if (j && bits[k >> 3] << (32 - j))
+					goto inc1;
+				while(k >= 8) {
+					k -= 8;
+					if (bits[k>>3]) {
+ inc1:
+						inc = 1;
+						goto haveinc;
+						}
+					}
+				}
+#ifdef Honor_FLT_ROUNDS
+		}
+#endif
+ haveinc:
+	b->ex += m*4;
+	i = m >> 3;
+	k = prec1 >> 3;
+	j = i;
+	if ((n = 4*(m & 7)))
+		for(;; ++j) {
+			bits[j-i] = bits[j] >> n;
+			if (j == k)
+				break;
+			bits[j-i] |= bits[j+1] << (32-n);
+			}
+	else
+		for(;; ++j) {
+			bits[j-i] = bits[j];
+			if (j == k)
+				break;
+			}
+	k = prec >> 3;
+	if (inc) {
+		for(j = 0; !(++bits[j] & 0xffffffff); ++j);
+		if (j > k) {
+ onebit:
+			bits[0] = 1;
+			b->ex += 4*prec;
+			return 1;
+			}
+		if ((j = prec & 7) < 7 && bits[k] >> (j+1)*4)
+			goto onebit;
+		}
+	for(i = 0; !(bits[i >> 3] & (0xf << 4*(i&7))); ++i);
+	if (i) {
+		b->ex += 4*i;
+		prec -= i;
+		j = i >> 3;
+		i &= 7;
+		i *= 4;
+		for(m = j; ; ++m) {
+			bits[m-j] = bits[m] >> i;
+			if (m == k)
+				break;
+			bits[m-j] |= bits[m+1] << (32 - i);
+			}
+		}
+	return prec;
+	}
+#endif /*}NO_PRINTF_A_FMT*/
+
+#define put(x) { *outbuf++ = x; if (outbuf == obe) outbuf = (*fput)(f,&rv); }
+
+ static int
+x_sprintf
+#ifdef KR_headers
+	(obe, fput, f, fmt, ap)
+	char *obe, *fmt; Finfo *f; Putfunc fput; va_list ap;
+#else
+	(char *obe, Putfunc fput, Finfo *f, const char *fmt, va_list ap)
+#endif
+{
+	FPBits fpb;
+	Fpbits fpbits;
+	U u;
+	char *digits, *ob0, *outbuf, *s, *s0, *se;
+	Const char *fmt0;
+	char buf[32];
+	long i;
+	unsigned long j, ul;
+	double x;
+	int alt, base, c, decpt, dot, conv, i1, k, lead0, left,
+		len, prec, prec1, psign, rv, sign, width;
+	long Ltmp, *ip;
+	short sh;
+	unsigned short us;
+	unsigned int ui;
+#ifdef Honor_FLT_ROUNDS
+	FPI fpi1;
+	int Rounding = -1;
+#endif
+#ifndef NO_PRINTF_A_FMT /*{*/
+	int bex, bw;
+#endif /*} NO_PRINTF_A_FMT */
+	static char hex[] = "0123456789abcdefpx";
+	static char Hex[] = "0123456789ABCDEFPX";
+
+	ob0 = outbuf = f->ob0;
+	rv = 0;
+	for(;;) {
+		for(;;) {
+			switch(c = *fmt++) {
+				case 0:
+					goto done;
+				case '%':
+					break;
+				default:
+					put(c)
+					continue;
+				}
+			break;
+			}
+		alt=dot=lead0=left=len=prec=psign=sign=width=0;
+		fpbits = dfpbits;
+		fmt0 = fmt;
+ fmtloop:
+		switch(conv = *fmt++) {
+			case ' ':
+			case '+':
+				sign = conv;
+				goto fmtloop;
+			case '-':
+				if (dot)
+					psign = 1;
+				else
+					left = 1;
+				goto fmtloop;
+			case '#':
+				alt = 1;
+				goto fmtloop;
+			case '0':
+				if (!lead0 && !dot) {
+					lead0 = 1;
+					goto fmtloop;
+					}
+			case '1':
+			case '2':
+			case '3':
+			case '4':
+			case '5':
+			case '6':
+			case '7':
+			case '8':
+			case '9':
+				k = conv - '0';
+				while((c = *fmt) >= '0' && c <= '9') {
+					k = 10*k + c - '0';
+					fmt++;
+					}
+				if (dot)
+					prec = psign ? -k : k;
+				else
+					width = k;
+				goto fmtloop;
+			case 'h':
+				len = 2;
+				goto fmtloop;
+			case 'L':
+#ifdef GDTOA_LD_fpbits /*{*/
+				fpbits = GDTOA_LD_fpbits;
+#ifdef GDTOA_Q_fpbits
+				if (*fmt == 'q') {
+					++fmt;
+					fpbits = Qfpbits;
+					}
+#endif
+#endif /*}*/
+				goto fmtloop;
+			case 'l':
+				len = 1;
+				goto fmtloop;
+			case '.':
+				dot = 1;
+				goto fmtloop;
+			case '*':
+				k = va_arg(ap, int);
+				if (dot)
+					prec = k;
+				else {
+					if (k < 0) {
+						sign = '-';
+						k = -k;
+						}
+					width = k;
+					}
+				goto fmtloop;
+			case 'c':
+				c = va_arg(ap, int);
+				put(c)
+				continue;
+			case '%':
+				put(conv)
+				continue;
+			case 'u':
+				switch(len) {
+				  case 0:
+					ui = va_arg(ap, int);
+					i = ui;
+					break;
+				  case 1:
+					i = va_arg(ap, long);
+					break;
+				  case 2:
+					us = va_arg(ap, int);
+					i = us;
+				  }
+				sign = 0;
+				goto have_i;
+			case 'i':
+			case 'd':
+				switch(len) {
+				  case 0:
+					k = va_arg(ap, int);
+					i = k;
+					break;
+				  case 1:
+					i = va_arg(ap, long);
+					break;
+				  case 2:
+					sh = va_arg(ap, int);
+					i = sh;
+				  }
+				if (i < 0) {
+					sign = '-';
+					i = -i;
+					}
+ have_i:
+				base = 10;
+				ul = i;
+				digits = hex;
+ baseloop:
+				if (dot)
+					lead0 = 0;
+				s = buf;
+				if (!ul)
+					alt = 0;
+				do {
+					j = ULDIV(ul, base);
+					*s++ = digits[ul - base*j];
+					}
+					while((ul = j));
+				prec -= c = s - buf;
+				if (alt && conv == 'o' && prec <= 0)
+					prec = 1;
+				if ((width -= c) > 0) {
+					if (prec > 0)
+						width -= prec;
+					if (sign)
+						width--;
+					if (alt == 2)
+						width--;
+					}
+				if (left) {
+					if (alt == 2)
+						put('0') /* for 0x */
+					if (sign)
+						put(sign)
+					while(--prec >= 0)
+						put('0')
+					do put(*--s)
+						while(s > buf);
+					while(--width >= 0)
+						put(' ')
+					continue;
+					}
+				if (width > 0) {
+					if (lead0) {
+						if (alt == 2)
+							put('0')
+						if (sign)
+							put(sign)
+						while(--width >= 0)
+							put('0')
+						goto s_loop;
+						}
+					else
+						while(--width >= 0)
+							put(' ')
+					}
+				if (alt == 2)
+					put('0')
+				if (sign)
+					put(sign)
+ s_loop:
+				while(--prec >= 0)
+					put('0')
+				do put(*--s)
+					while(s > buf);
+				continue;
+			case 'n':
+				ip = va_arg(ap, long*);
+				if (!ip)
+					ip = &Ltmp;
+				c = outbuf - ob0 + rv;
+				switch(len) {
+				  case 0:
+					*(int*)ip = c;
+					break;
+				  case 1:
+					*ip = c;
+					break;
+				  case 2:
+					*(short*)ip = c;
+				  }
+				break;
+			case 'p':
+				len = alt = 1;
+				/* no break */
+			case 'x':
+				digits = hex;
+				goto more_x;
+			case 'X':
+				digits = Hex;
+ more_x:
+				if (alt) {
+					alt = 2;
+					sign = conv;
+					}
+				else
+					sign = 0;
+				base = 16;
+ get_u:
+				switch(len) {
+				  case 0:
+					ui = va_arg(ap, int);
+					ul = ui;
+					break;
+				  case 1:
+					ul = va_arg(ap, long);
+					break;
+				  case 2:
+					us = va_arg(ap, int);
+					ul = us;
+				  }
+				if (!ul)
+					sign = alt = 0;
+				goto baseloop;
+			case 'o':
+				base = 8;
+				digits = hex;
+				goto get_u;
+			case 's':
+				s0 = 0;
+				s = va_arg(ap, char*);
+				if (!s)
+					s = "<NULL>";
+				if (prec < 0)
+					prec = 0;
+ have_s:
+				if (dot) {
+					for(c = 0; c < prec; c++)
+						if (!s[c])
+							break;
+					prec = c;
+					}
+				else
+					prec = strlen(s);
+				width -= prec;
+				if (!left)
+					while(--width >= 0)
+						put(' ')
+				while(--prec >= 0)
+					put(*s++)
+				while(--width >= 0)
+					put(' ')
+				if (s0)
+					freedtoa(s0);
+				continue;
+			case 'f':
+				if (!dot)
+					prec = 6;
+#ifdef GDTOA_H_INCLUDED
+				if (fpbits == dfpbits) {
+#endif
+				x = va_arg(ap, double);
+				s = s0 = dtoa(x, 3, prec, &decpt, &fpb.sign, &se);
+#ifdef GDTOA_H_INCLUDED
+					}
+				else {
+#ifdef GDTOA_both
+					if (fpbits == GDTOA_LD_fpbits)
+						u.ld = va_arg(ap, long double);
+					else
+						u.Qd = va_arg(ap, GDTOA_Qtype);
+#else
+					u.ld = va_arg(ap, long double);
+#endif
+					fpbits(&u, &fpb);
+					RoundCheck
+					s = s0 = fpb.gdtoa(fpb.fpi, fpb.ex, fpb.bits,
+						&fpb.kind, 3, prec, &decpt, &se);
+					}
+#endif
+				if (decpt == 9999) {
+ fmt9999:
+					dot = prec = alt = 0;
+					if (*s == 'N')
+						goto have_s;
+					decpt = strlen(s);
+					}
+ f_fmt:
+				if (fpb.sign && (x||sign))
+					sign = '-';
+				if (prec > 0)
+					width -= prec;
+				if (width > 0) {
+					if (sign)
+						--width;
+					if (decpt <= 0) {
+						--width;
+						if (prec > 0)
+							--width;
+						}
+					else {
+						if (s == se)
+							decpt = 1;
+						width -= decpt;
+						if (prec > 0 || alt)
+							--width;
+						}
+					}
+				if (width > 0 && !left) {
+					if (lead0) {
+						if (sign)
+							put(sign)
+						sign = 0;
+						do put('0')
+							while(--width > 0);
+						}
+					else do put(' ')
+						while(--width > 0);
+					}
+				if (sign)
+					put(sign)
+				if (decpt <= 0) {
+					put('0')
+					if (prec > 0 || alt)
+						put('.')
+					while(decpt < 0) {
+						put('0')
+						prec--;
+						decpt++;
+						}
+					}
+				else {
+					do {
+						if ((c = *s))
+							s++;
+						else
+							c = '0';
+						put(c)
+						}
+						while(--decpt > 0);
+					if (prec > 0 || alt)
+						put('.')
+					}
+				while(--prec >= 0) {
+					if ((c = *s))
+						s++;
+					else
+						c = '0';
+					put(c)
+					}
+				while(--width >= 0)
+					put(' ')
+				if (s0)
+					freedtoa(s0);
+				continue;
+			case 'G':
+			case 'g':
+				if (!dot)
+					prec = 6;
+				if (prec < 0)
+					prec = 0;
+#ifdef GDTOA_H_INCLUDED
+				if (fpbits == dfpbits) {
+#endif
+					x = va_arg(ap, double);
+					s = s0 = dtoa(x, prec ? 2 : 0, prec, &decpt,
+						&fpb.sign, &se);
+#ifdef GDTOA_H_INCLUDED
+					}
+				else {
+#ifdef GDTOA_both
+					if (fpbits == GDTOA_LD_fpbits)
+						u.ld = va_arg(ap, long double);
+					else
+						u.Qd = va_arg(ap, GDTOA_Qtype);
+#else
+					u.ld = va_arg(ap, long double);
+#endif
+					fpbits(&u, &fpb);
+					RoundCheck
+					s = s0 = fpb.gdtoa(fpb.fpi, fpb.ex, fpb.bits,
+						&fpb.kind, prec ? 2 : 0, prec, &decpt, &se);
+					}
+#endif
+				if (decpt == 9999)
+					goto fmt9999;
+				c = se - s;
+				prec1 = prec;
+				if (!prec) {
+					prec = c;
+					prec1 = c + (s[1] || alt ? 5 : 4);
+					/* %.0g gives 10 rather than 1e1 */
+					}
+				if (decpt > -4 && decpt <= prec1) {
+					if (alt)
+						prec -= decpt;
+					else
+						prec = c - decpt;
+					if (prec < 0)
+						prec = 0;
+					goto f_fmt;
+					}
+				conv -= 2;
+				if (!alt && prec > c)
+					prec = c;
+				--prec;
+				goto e_fmt;
+			case 'e':
+			case 'E':
+				if (!dot)
+					prec = 6;
+				if (prec < 0)
+					prec = 0;
+#ifdef GDTOA_H_INCLUDED
+				if (fpbits == dfpbits) {
+#endif
+					x = va_arg(ap, double);
+					s = s0 = dtoa(x, prec ? 2 : 0, prec+1, &decpt,
+						&fpb.sign, &se);
+#ifdef GDTOA_H_INCLUDED
+					}
+				else {
+#ifdef GDTOA_both
+					if (fpbits == GDTOA_LD_fpbits)
+						u.ld = va_arg(ap, long double);
+					else
+						u.Qd = va_arg(ap, GDTOA_Qtype);
+#else
+					u.ld = va_arg(ap, long double);
+#endif
+					fpbits(&u, &fpb);
+					RoundCheck
+					s = s0 = fpb.gdtoa(fpb.fpi, fpb.ex, fpb.bits,
+						&fpb.kind, prec ? 2 : 0, prec, &decpt, &se);
+					}
+#endif
+				if (decpt == 9999)
+					goto fmt9999;
+ e_fmt:
+				if (fpb.sign && (x||sign))
+					sign = '-';
+				if ((width -= prec + 5) > 0) {
+					if (sign)
+						--width;
+					if (prec || alt)
+						--width;
+					}
+				if ((c = --decpt) < 0)
+					c = -c;
+				while(c >= 100) {
+					--width;
+					c /= 10;
+					}
+				if (width > 0 && !left) {
+					if (lead0) {
+						if (sign)
+							put(sign)
+						sign = 0;
+						do put('0')
+							while(--width > 0);
+						}
+					else do put(' ')
+						while(--width > 0);
+					}
+				if (sign)
+					put(sign)
+				put(*s++)
+				if (prec || alt)
+					put('.')
+				while(--prec >= 0) {
+					if ((c = *s))
+						s++;
+					else
+						c = '0';
+					put(c)
+					}
+				put(conv)
+				if (decpt < 0) {
+					put('-')
+					decpt = -decpt;
+					}
+				else
+					put('+')
+				for(c = 2, k = 10; 10*k <= decpt; c++, k *= 10);
+				for(;;) {
+					i1 = decpt / k;
+					put(i1 + '0')
+					if (--c <= 0)
+						break;
+					decpt -= i1*k;
+					decpt *= 10;
+					}
+				while(--width >= 0)
+					put(' ')
+				freedtoa(s0);
+				continue;
+#ifndef NO_PRINTF_A_FMT
+			case 'a':
+				digits = hex;
+				goto more_a;
+			case 'A':
+				digits = Hex;
+ more_a:
+#ifdef GDTOA_H_INCLUDED /*{{*/
+				if (fpbits == dfpbits)
+					u.d = va_arg(ap, double);
+#ifdef GDTOA_both /*{*/
+				else if (fpbits == GDTOA_LD_fpbits)
+					u.ld = va_arg(ap, long double);
+				else
+					u.Qd = va_arg(ap, GDTOA_Qtype);
+#else
+				else
+					u.ld = va_arg(ap, long double);
+#endif /*}*/
+#else /*}{*/
+				u.d = va_arg(ap, double);
+#endif /*}}*/
+				fpbits(&u, &fpb);
+				if (fpb.kind == STRTOG_Infinite) {
+					s = "Infinity";
+					s0 = 0;
+					goto fmt9999;
+					}
+				if (fpb.kind == STRTOG_NaN) {
+					s = "NaN";
+					s0 = 0;
+					goto fmt9999;
+					}
+				prec1 = fpiprec(&fpb);
+				if (dot && prec < prec1)
+					prec1 = bround(&fpb, prec, prec1);
+				bw = 1;
+				bex = fpb.ex + 4*prec1;
+				if (bex) {
+					if ((i1 = bex) < 0)
+						i1 = -i1;
+					while(i1 >= 10) {
+						++bw;
+						i1 /= 10;
+						}
+					}
+				if (fpb.sign && (sign || fpb.kind != STRTOG_Zero))
+					sign = '-';
+				if ((width -= bw + 5) > 0) {
+					if (sign)
+						--width;
+					if (prec1 || alt)
+						--width;
+					}
+				if ((width -= prec1) > 0 && !left && !lead0) {
+					do put(' ')
+						while(--width > 0);
+					}
+				if (sign)
+					put(sign)
+				put('0')
+				put(digits[17])
+				if (lead0 && width > 0 && !left) {
+					do put('0')
+						while(--width > 0);
+					}
+				i1 = prec1 & 7;
+				k = prec1 >> 3;
+				put(digits[(fpb.bits[k] >> 4*i1) & 0xf])
+				if (prec1 > 0 || alt)
+					put('.')
+				if (prec1 > 0) {
+					prec -= prec1;
+					while(prec1 > 0) {
+						if (--i1 < 0) {
+							if (--k < 0)
+								break;
+							i1 = 7;
+							}
+						put(digits[(fpb.bits[k] >> 4*i1) & 0xf])
+						--prec1;
+						}
+					if (alt && prec > 0)
+						do put(0)
+							while(--prec > 0);
+					}
+				put(digits[16])
+				if (bex < 0) {
+					put('-')
+					bex = -bex;
+					}
+				else
+					put('+')
+				for(c = 1; 10*c <= bex; c *= 10);
+				for(;;) {
+					i1 = bex / c;
+					put('0' + i1)
+					if (!--bw)
+						break;
+					bex -= i1 * c;
+					bex *= 10;
+					}
+				while(--width >= 0)
+					put(' ')
+				continue;
+#endif /* NO_PRINTF_A_FMT */
+			default:
+				put('%')
+				while(fmt0 < fmt)
+					put(*fmt0++)
+				continue;
+			}
+		}
+ done:
+	*outbuf = 0;
+	return (f->lastlen = outbuf - ob0) + rv;
+	}
+
+#define Bsize 256
+
+ int
+Printf
+#ifdef KR_headers
+	(va_alist)
+ va_dcl
+{
+	char *fmt;
+
+	va_list ap;
+	int rv;
+	Finfo f;
+	char buf[Bsize];
+
+	va_start(ap);
+	fmt = va_arg(ap, char*);
+	/*}*/
+#else
+	(const char *fmt, ...)
+{
+	va_list ap;
+	int rv;
+	Finfo f;
+	char buf[Bsize];
+
+	va_start(ap, fmt);
+#endif
+	f.u.cf = stdout;
+	f.ob0 = buf;
+	f.obe1 = buf + Bsize - 1;
+#ifdef _windows_
+	if (fileno(stdout) == stdout_fileno_ASL) {
+		rv = x_sprintf(f.obe1, Wput, &f, fmt, ap);
+		mwrite(buf, f.lastlen);
+		}
+	else
+#endif
+#ifdef PF_BUF
+	if (stdout == stderr_ASL) {
+		rv = x_sprintf(f.obe1, pfput, &f, fmt, ap);
+		pf_put(buf, f.lastlen);
+		}
+	else
+#endif
+		{
+		rv = x_sprintf(f.obe1, Fput, &f, fmt, ap);
+		fputs(buf, stdout);
+		}
+	va_end(ap);
+	return rv;
+	}
+
+ static char *
+Sput
+#ifdef KR_headers
+	(f, rvp) Finfo *f; int *rvp;
+#else
+	(Finfo *f, int *rvp)
+#endif
+{
+	if (Printf("\nBUG! Sput called!\n", f, rvp))
+		/* pass vp, rvp and return 0 to shut diagnostics off */
+		exit(250);
+	return 0;
+	}
+
+ int
+Sprintf
+#ifdef KR_headers
+	(va_alist)
+ va_dcl
+{
+	char *s, *fmt;
+	va_list ap;
+	int rv;
+	Finfo f;
+
+	va_start(ap);
+	s = va_arg(ap, char*);
+	fmt = va_arg(ap, char*);
+	/*}*/
+#else
+	(char *s, const char *fmt, ...)
+{
+	va_list ap;
+	int rv;
+	Finfo f;
+
+	va_start(ap, fmt);
+#endif
+	f.ob0 = s;
+	rv = x_sprintf(s, Sput, &f, fmt, ap);
+	va_end(ap);
+	return rv;
+	}
+
+ int
+Fprintf
+#ifdef KR_headers
+	(va_alist)
+ va_dcl
+{
+	FILE *F;
+	char *s, *fmt;
+	va_list ap;
+	int rv;
+	Finfo f;
+	char buf[Bsize];
+
+	va_start(ap);
+	F = va_arg(ap, FILE*);
+	fmt = va_arg(ap, char*);
+	/*}*/
+#else
+	(FILE *F, const char *fmt, ...)
+{
+	va_list ap;
+	int rv;
+	Finfo f;
+	char buf[Bsize];
+
+	va_start(ap, fmt);
+#endif
+	f.u.cf = F;
+	f.ob0 = buf;
+	f.obe1 = buf + Bsize - 1;
+#ifdef MESS
+	if (stdout_or_err(F)) {
+#ifdef _windows_
+		if (fileno(stdout) == stdout_fileno_ASL) {
+			rv = x_sprintf(f.obe1, Wput, &f, fmt, ap);
+			mwrite(buf, f.lastlen);
+			}
+		else
+#endif
+#ifdef PF_BUF
+		if (F == stderr_ASL) {
+			rv = x_sprintf(f.obe1, pfput, &f, fmt, ap);
+			pf_put(buf, f.lastlen);
+			}
+		else
+#endif
+			{
+			rv = x_sprintf(f.obe1, Fput, &f, fmt, ap);
+			fputs(buf, F);
+			}
+		}
+	else
+#endif /*MESS*/
+		{
+#ifdef PF_BUF
+		if (F == stderr_ASL) {
+			rv = x_sprintf(f.obe1, pfput, &f, fmt, ap);
+			pf_put(buf, f.lastlen);
+			}
+	else
+#endif
+			{
+			rv = x_sprintf(f.obe1, Fput, &f, fmt, ap);
+			fputs(buf, F);
+			}
+		}
+	va_end(ap);
+	return rv;
+	}
+
+ int
+Vsprintf
+#ifdef KR_headers
+	(s, fmt, ap) char *s, *fmt; va_list ap;
+#else
+	(char *s, const char *fmt, va_list ap)
+#endif
+{
+	Finfo f;
+	return x_sprintf(f.ob0 = s, Sput, &f, fmt, ap);
+	}
+
+ int
+Vfprintf
+#ifdef KR_headers
+	(F, fmt, ap) FILE *F; char *fmt; va_list ap;
+#else
+	(FILE *F, const char *fmt, va_list ap)
+#endif
+{
+	char buf[Bsize];
+	int rv;
+	Finfo f;
+
+	f.u.cf = F;
+	f.ob0 = buf;
+	f.obe1 = buf + Bsize - 1;
+#ifdef MESS
+	if (stdout_or_err(F)) {
+#ifdef _windows_
+		if (fileno(stdout) == stdout_fileno_ASL) {
+			rv = x_sprintf(f.obe1, Wput, &f, fmt, ap);
+			mwrite(buf, f.lastlen);
+			}
+		else
+#endif
+#ifdef PF_BUF
+		if (F == stderr_ASL) {
+			rv = x_sprintf(f.obe1, pfput, &f, fmt, ap);
+			pf_put(buf, f.lastlen);
+			}
+		else
+#endif
+			{
+			rv = x_sprintf(f.obe1, Fput, &f, fmt, ap);
+			fputs(buf, F);
+			}
+		}
+	else
+#endif /*MESS*/
+		{
+#ifdef PF_BUF
+		if (F == stderr_ASL) {
+			rv = x_sprintf(f.obe1, pfput, &f, fmt, ap);
+			pf_put(buf, f.lastlen);
+			}
+		else
+#endif
+			{
+			rv = x_sprintf(f.obe1, Fput, &f, fmt, ap);
+			fputs(buf, F);
+			}
+		}
+	va_end(ap);
+	return rv;
+	}
+
+ void
+Perror
+#ifdef KR_headers
+	(s) char *s;
+#else
+	(const char *s)
+#endif
+{
+	if (s && *s)
+		Fprintf(Stderr, "%s: ", s);
+	Fprintf(Stderr, "%s\n", strerror(errno));
+	}
+
+ static char *
+Snput
+#ifdef KR_headers
+	(f, rvp) Finfo *f; int *rvp;
+#else
+	(Finfo *f, int *rvp)
+#endif
+{
+	char *s, *s0;
+	size_t L;
+
+	*rvp += Bsize;
+	s0 = f->ob0;
+	s = f->u.sf;
+	if ((L = f->obe1 - s) > Bsize) {
+		L = Bsize;
+		goto copy;
+		}
+	if (L > 0) {
+ copy:
+		memcpy(s, s0, L);
+		f->u.sf = s + L;
+		}
+	return s0;
+	}
+
+ int
+Vsnprintf
+#ifdef KR_headers
+	(s, n, fmt, ap) char *s; size_t n; char *fmt; va_list ap;
+#else
+	(char *s, size_t n, const char *fmt, va_list ap)
+#endif
+{
+	Finfo f;
+	char buf[Bsize];
+	int L, rv;
+
+	if (n <= 0 || !s) {
+		n = 1;
+		s = buf;
+		}
+	f.u.sf = s;
+	f.ob0 = buf;
+	f.obe1 = s + n - 1;
+	rv = x_sprintf(buf + Bsize, Snput, &f, fmt, ap);
+	if (f.lastlen > (L = f.obe1 - f.u.sf))
+		f.lastlen = L;
+	if (f.lastlen > 0) {
+		memcpy(f.u.sf, buf, f.lastlen);
+		f.u.sf += f.lastlen;
+		}
+	*f.u.sf = 0;
+	return rv;
+	}
+ int
+Snprintf
+#ifdef KR_headers
+	(va_alist)
+ va_dcl
+{
+	char *s, *fmt;
+	int rv;
+	size_t n;
+	va_list ap;
+
+	va_start(ap);
+	s = va_arg(ap, char*);
+	n = va_arg(ap, size_t);
+	fmt = va_arg(ap, char*);
+	/*}*/
+#else
+	(char *s, size_t n, const char *fmt, ...)
+{
+	int rv;
+	va_list ap;
+
+	va_start(ap, fmt);
+#endif
+	rv = Vsnprintf(s, n, fmt, ap);
+	va_end(ap);
+	return rv;
+	}
+
+
+#ifdef __cplusplus
+}
+#endif
diff --git a/contrib/gdtoa/smisc.c b/contrib/gdtoa/smisc.c
index 50431fc..f4dbafb 100644
--- a/contrib/gdtoa/smisc.c
+++ b/contrib/gdtoa/smisc.c
@@ -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
diff --git a/contrib/gdtoa/stdio1.h b/contrib/gdtoa/stdio1.h
new file mode 100644
index 0000000..ef2bb63
--- /dev/null
+++ b/contrib/gdtoa/stdio1.h
@@ -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 */
diff --git a/contrib/gdtoa/strtoIg.c b/contrib/gdtoa/strtoIg.c
index 90c88db..6a17760 100644
--- a/contrib/gdtoa/strtoIg.c
+++ b/contrib/gdtoa/strtoIg.c
@@ -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);
diff --git a/contrib/gdtoa/strtod.c b/contrib/gdtoa/strtod.c
index 5550853..fe8cde8 100644
--- a/contrib/gdtoa/strtod.c
+++ b/contrib/gdtoa/strtod.c
@@ -59,6 +59,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
@@ -73,10 +95,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
@@ -90,7 +116,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;
 			}
@@ -117,7 +143,7 @@ strtod
 #endif /*}*/
 
 	sign = nz0 = nz = decpt = 0;
-	dval(rv) = 0.;
+	dval(&rv) = 0.;
 	for(s = s00;;s++) switch(*s) {
 		case '-':
 			sign = 1;
@@ -149,20 +175,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;
@@ -287,8 +305,8 @@ strtod
 					--s;
 					if (!match(&s,"inity"))
 						++s;
-					word0(rv) = 0x7ff00000;
-					word1(rv) = 0;
+					word0(&rv) = 0x7ff00000;
+					word1(&rv) = 0;
 					goto ret;
 					}
 				break;
@@ -299,13 +317,13 @@ strtod
 					if (*s == '(' /*)*/
 					 && hexnan(&s, &fpinan, bits)
 							== STRTOG_NaNbits) {
-						word0(rv) = 0x7ff80000 | bits[1];
-						word1(rv) = bits[0];
+						word0(&rv) = 0x7ff80000 | 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
@@ -329,13 +347,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
@@ -347,6 +365,7 @@ strtod
 			) {
 		if (!e)
 			goto ret;
+#ifndef ROUND_BIASED_without_Round_Up
 		if (e > 0) {
 			if (e <= Ten_pmax) {
 #ifdef VAX
@@ -355,11 +374,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
 				}
@@ -371,25 +390,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;
 				}
@@ -399,14 +418,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;
 
@@ -434,67 +454,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;
@@ -503,44 +529,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.
 				 */
@@ -558,7 +579,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) {
@@ -579,12 +600,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
@@ -594,7 +622,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
@@ -645,15 +673,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
@@ -662,66 +690,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*/
@@ -730,12 +758,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
 				) {
@@ -760,32 +788,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
@@ -800,7 +830,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 ==> */
@@ -811,10 +841,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
@@ -822,16 +852,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
 				}
@@ -843,14 +890,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 */
@@ -860,45 +907,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
@@ -907,58 +954,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)
@@ -968,7 +1015,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;
 				}
@@ -982,12 +1029,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)
@@ -995,36 +1047,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);
 	}
 
diff --git a/contrib/gdtoa/strtodI.c b/contrib/gdtoa/strtodI.c
index 98f8891..0b7b8a4 100644
--- a/contrib/gdtoa/strtodI.c
+++ b/contrib/gdtoa/strtodI.c
@@ -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;
 	}
diff --git a/contrib/gdtoa/strtodg.c b/contrib/gdtoa/strtodg.c
index a69ce30..5059869 100644
--- a/contrib/gdtoa/strtodg.c
+++ b/contrib/gdtoa/strtodg.c
@@ -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;
 				}
diff --git a/contrib/gdtoa/strtof.c b/contrib/gdtoa/strtof.c
index df9a75a..84bfe24 100644
--- a/contrib/gdtoa/strtof.c
+++ b/contrib/gdtoa/strtof.c
@@ -59,7 +59,7 @@ strtof(CONST char *s, char **sp)
 		break;
 
 	  case STRTOG_Normal:
-		u.L[0] = bits[0] & 0x7fffff | exp + 0x7f + 23 << 23;
+		u.L[0] = (bits[0] & 0x7fffff) | ((exp + 0x7f + 23) << 23);
 		break;
 
 	  case STRTOG_NaNbits:
diff --git a/contrib/gdtoa/strtopdd.c b/contrib/gdtoa/strtopdd.c
index c665976..738372d 100644
--- a/contrib/gdtoa/strtopdd.c
+++ b/contrib/gdtoa/strtopdd.c
@@ -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:
diff --git a/contrib/gdtoa/strtopf.c b/contrib/gdtoa/strtopf.c
index 74694a0..23ca5cb 100644
--- a/contrib/gdtoa/strtopf.c
+++ b/contrib/gdtoa/strtopf.c
@@ -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:
diff --git a/contrib/gdtoa/strtopx.c b/contrib/gdtoa/strtopx.c
index 07d0458..f7a25ff 100644
--- a/contrib/gdtoa/strtopx.c
+++ b/contrib/gdtoa/strtopx.c
@@ -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:
diff --git a/contrib/gdtoa/strtopxL.c b/contrib/gdtoa/strtopxL.c
index b0f5cdd..6519a41 100644
--- a/contrib/gdtoa/strtopxL.c
+++ b/contrib/gdtoa/strtopxL.c
@@ -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:
diff --git a/contrib/gdtoa/strtordd.c b/contrib/gdtoa/strtordd.c
index 9c2b46e..e0b8e6a 100644
--- a/contrib/gdtoa/strtordd.c
+++ b/contrib/gdtoa/strtordd.c
@@ -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:
diff --git a/contrib/gdtoa/strtorf.c b/contrib/gdtoa/strtorf.c
index 46b0ba2..21c6d28 100644
--- a/contrib/gdtoa/strtorf.c
+++ b/contrib/gdtoa/strtorf.c
@@ -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:
diff --git a/contrib/gdtoa/strtorx.c b/contrib/gdtoa/strtorx.c
index e9fd45f..cd938f1 100644
--- a/contrib/gdtoa/strtorx.c
+++ b/contrib/gdtoa/strtorx.c
@@ -89,7 +89,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:
diff --git a/contrib/gdtoa/strtorxL.c b/contrib/gdtoa/strtorxL.c
index ff62a61..58acc80 100644
--- a/contrib/gdtoa/strtorxL.c
+++ b/contrib/gdtoa/strtorxL.c
@@ -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:
diff --git a/contrib/gdtoa/ulp.c b/contrib/gdtoa/ulp.c
index 7810a5c..17e9f86 100644
--- a/contrib/gdtoa/ulp.c
+++ b/contrib/gdtoa/ulp.c
@@ -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);
 	}