Replace our strtod with the gdtoa package.
authorPeter Avalos <pavalos@theshell.com>
Mon, 19 Jan 2009 19:48:45 +0000 (14:48 -0500)
committerPeter Avalos <pavalos@theshell.com>
Tue, 7 Apr 2009 07:09:01 +0000 (21:09 -1000)
Obtained-from: FreeBSD

22 files changed:
include/stdlib.h
lib/libc/Makefile
lib/libc/Makefile.inc
lib/libc/amd64/Makefile.inc
lib/libc/amd64/gd_qnan.h [new file with mode: 0644]
lib/libc/amd64/stdlib/gdtoa.mk [deleted file]
lib/libc/gdtoa/Makefile.inc [new file with mode: 0644]
lib/libc/gdtoa/_hdtoa.c [new file with mode: 0644]
lib/libc/gdtoa/_hldtoa.c [new file with mode: 0644]
lib/libc/gdtoa/_ldtoa.c [new file with mode: 0644]
lib/libc/gdtoa/glue.c [new file with mode: 0644]
lib/libc/gdtoa/machdep_ldisx.c [new file with mode: 0644]
lib/libc/i386/Makefile.inc [new file with mode: 0644]
lib/libc/i386/_fpmath.h [new file with mode: 0644]
lib/libc/i386/arith.h [new file with mode: 0644]
lib/libc/i386/gd_qnan.h [new file with mode: 0644]
lib/libc/include/fpmath.h [new file with mode: 0644]
lib/libc/stdio/vfprintf.c
lib/libc/stdlib/Makefile.inc
lib/libc/stdlib/strtod.3
lib/libc/stdlib/strtod.c [deleted file]
lib/libc_rtld/Makefile

index 7f717c0..dee562f 100644 (file)
@@ -98,12 +98,12 @@ int  rand(void);
 void   *realloc(void *, size_t);
 void    srand(unsigned);
 double  strtod(const char * __restrict, char ** __restrict);
-/* float        strtof(const char * __restrict, char ** __restrict); */
+float   strtof(const char * __restrict, char ** __restrict);
 #if !defined(_KERNEL_VIRTUAL)
 long    strtol(const char * __restrict, char ** __restrict, int);
 #endif
-/* long double
-        strtold(const char * __restrict, char ** __restrict); */
+long double
+        strtold(const char * __restrict, char ** __restrict);
 #if !defined(_KERNEL_VIRTUAL)
 unsigned long
         strtoul(const char * __restrict, char ** __restrict, int);
index d8751a4..d4b6cab 100644 (file)
@@ -9,7 +9,7 @@
 # stubs, remove just -DSYSLIBC_RCS from CFLAGS.
 LIB=c
 SHLIB_MAJOR= 6
-CFLAGS+=-DLIBC_RCS -DSYSLIBC_RCS -I${.CURDIR}/include
+CFLAGS+=-DLIBC_RCS -DSYSLIBC_RCS -I${.CURDIR}/include ${AINC}
 AINC=  -I${.OBJDIR} -I${.CURDIR}/${MACHINE_ARCH}
 CLEANFILES+=tags
 PRECIOUSLIB=   yes
index ae9e9b5..1a83c29 100644 (file)
@@ -13,23 +13,19 @@ NOASM=
 
 #WARNS=6
 
-#
-# If there is a machine dependent makefile, use it:
-#
-.if exists(${.CURDIR}/../libc/${MACHINE_ARCH}/Makefile.inc)
 .include "${.CURDIR}/../libc/${MACHINE_ARCH}/Makefile.inc"
-.endif
 
 .include "${.CURDIR}/../libc/citrus/Makefile.inc"
 .include "${.CURDIR}/../libc/db/Makefile.inc"
 .include "${.CURDIR}/../libc/compat-43/Makefile.inc"
+.include "${.CURDIR}/../libc/gdtoa/Makefile.inc"
 .include "${.CURDIR}/../libc/gen/Makefile.inc"
 .include "${.CURDIR}/../libc/gmon/Makefile.inc"
 .include "${.CURDIR}/../libc/iconv/Makefile.inc"
 .include "${.CURDIR}/../libc/locale/Makefile.inc"
 .include "${.CURDIR}/../libc/net/Makefile.inc"
 .include "${.CURDIR}/../libc/nls/Makefile.inc"
-.if !defined(NO_QUAD)
+.if ${MACHINE_ARCH} != "amd64"
 .include "${.CURDIR}/../libc/quad/Makefile.inc"
 .endif
 .include "${.CURDIR}/../libc/regex/Makefile.inc"
@@ -67,3 +63,6 @@ SRCS+=        ${_src}
 .endif
 .endfor
 .endif
+
+# Disable warnings in contributed sources.
+CWARNFLAGS:=   ${.IMPSRC:Ngdtoa_*.c:C/^.+$/${CWARNFLAGS}/}
index 6ecb550..1fb1310 100644 (file)
@@ -1,10 +1,9 @@
-# $FreeBSD: src/lib/libc/amd64/Makefile.inc,v 1.1 2003/07/22 06:34:57 peter Exp $
+# $FreeBSD: src/lib/libc/amd64/Makefile.inc,v 1.6 2007/12/03 07:17:32 das Exp $
 # $DragonFly: src/lib/libc/amd64/Makefile.inc,v 1.1 2004/02/02 05:43:14 dillon Exp $
 #
 # Machine dependent definitions for the amd64 architecture.
 #
 
-#
-# AMD64 is 64-bit, so it doesn't need quad functions:
-#
-NO_QUAD=1
+# Long double is 80 bits
+GDTOASRCS+=strtorx.c
+MDSRCS+=machdep_ldisx.c
diff --git a/lib/libc/amd64/gd_qnan.h b/lib/libc/amd64/gd_qnan.h
new file mode 100644 (file)
index 0000000..065ab0a
--- /dev/null
@@ -0,0 +1,21 @@
+/*
+ * MD header for contrib/gdtoa
+ *
+ * This file can be generated by compiling and running contrib/gdtoa/qnan.c
+ * on the target architecture after arith.h has been generated.
+ *
+ * $FreeBSD: src/lib/libc/amd64/gd_qnan.h,v 1.2 2007/12/16 21:15:08 das Exp $
+ */
+
+#define f_QNAN 0x7fc00000
+#define d_QNAN0 0x0
+#define d_QNAN1 0x7ff80000
+#define ld_QNAN0 0x0
+#define ld_QNAN1 0xc0000000
+#define ld_QNAN2 0x7fff
+#define ld_QNAN3 0x0
+#define ldus_QNAN0 0x0
+#define ldus_QNAN1 0x0
+#define ldus_QNAN2 0x0
+#define ldus_QNAN3 0xc000
+#define ldus_QNAN4 0x7fff
diff --git a/lib/libc/amd64/stdlib/gdtoa.mk b/lib/libc/amd64/stdlib/gdtoa.mk
deleted file mode 100644 (file)
index 4147ac3..0000000
+++ /dev/null
@@ -1,6 +0,0 @@
-# $FreeBSD: src/lib/libc/amd64/stdlib/gdtoa.mk,v 1.1 2003/03/12 20:29:59 das Exp $
-# $DragonFly: src/lib/libc/amd64/stdlib/gdtoa.mk,v 1.1 2004/02/02 05:43:14 dillon Exp $
-
-# Long double is 80 bits
-GDTOASRCS+=strtopx.c
-MDSRCS+=machdep_ldisx.c
diff --git a/lib/libc/gdtoa/Makefile.inc b/lib/libc/gdtoa/Makefile.inc
new file mode 100644 (file)
index 0000000..02df92d
--- /dev/null
@@ -0,0 +1,18 @@
+# $FreeBSD: src/lib/libc/gdtoa/Makefile.inc,v 1.10 2008/04/12 03:11:36 das Exp $
+
+# netlib gdtoa sources
+.PATH: ${.CURDIR}/../libc/gdtoa
+
+MISRCS+=_hdtoa.c _hldtoa.c _ldtoa.c glue.c
+GDTOASRCS+=dmisc.c dtoa.c gdtoa.c gethex.c gmisc.c \
+       hd_init.c hexnan.c misc.c smisc.c \
+       strtod.c strtodg.c strtof.c strtord.c sum.c ulp.c
+
+CFLAGS+=-I${.CURDIR}/../../contrib/gdtoa
+
+.for src in ${GDTOASRCS}
+MISRCS+=gdtoa_${src}
+CLEANFILES+=gdtoa_${src}
+gdtoa_${src}:
+       ln -sf ${.CURDIR}/../../contrib/gdtoa/${src} ${.TARGET}
+.endfor
diff --git a/lib/libc/gdtoa/_hdtoa.c b/lib/libc/gdtoa/_hdtoa.c
new file mode 100644 (file)
index 0000000..0c20616
--- /dev/null
@@ -0,0 +1,143 @@
+/*-
+ * Copyright (c) 2004-2008 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * $FreeBSD: src/lib/libc/gdtoa/_hdtoa.c,v 1.7 2008/04/12 14:53:52 das Exp $
+ */
+
+#include <float.h>
+#include <limits.h>
+#include <math.h>
+
+#include "../stdio/floatio.h"
+#include "fpmath.h"
+#include "gdtoaimp.h"
+
+/* Strings values used by dtoa() */
+#define        INFSTR  "Infinity"
+#define        NANSTR  "NaN"
+
+#define        DBL_ADJ (DBL_MAX_EXP - 2)
+#define        SIGFIGS ((DBL_MANT_DIG + 3) / 4 + 1)
+
+static const float one[] = { 1.0f, -1.0f };
+
+/*
+ * This procedure converts a double-precision number in IEEE format
+ * into a string of hexadecimal digits and an exponent of 2.  Its
+ * behavior is bug-for-bug compatible with dtoa() in mode 2, with the
+ * following exceptions:
+ *
+ * - An ndigits < 0 causes it to use as many digits as necessary to
+ *   represent the number exactly.
+ * - The additional xdigs argument should point to either the string
+ *   "0123456789ABCDEF" or the string "0123456789abcdef", depending on
+ *   which case is desired.
+ * - This routine does not repeat dtoa's mistake of setting decpt
+ *   to 9999 in the case of an infinity or NaN.  INT_MAX is used
+ *   for this purpose instead.
+ *
+ * Note that the C99 standard does not specify what the leading digit
+ * should be for non-zero numbers.  For instance, 0x1.3p3 is the same
+ * as 0x2.6p2 is the same as 0x4.cp3.  This implementation always makes
+ * the leading digit a 1. This ensures that the exponent printed is the
+ * actual base-2 exponent, i.e., ilogb(d).
+ *
+ * Inputs:     d, xdigs, ndigits
+ * Outputs:    decpt, sign, rve
+ */
+char *
+__hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign,
+    char **rve)
+{
+       union IEEEd2bits u;
+       char *s, *s0;
+       int bufsize;
+       uint32_t manh, manl;
+
+       u.d = d;
+       *sign = u.bits.sign;
+
+       switch (fpclassify(d)) {
+       case FP_NORMAL:
+               *decpt = u.bits.exp - DBL_ADJ;
+               break;
+       case FP_ZERO:
+               *decpt = 1;
+               return (nrv_alloc("0", rve, 1));
+       case FP_SUBNORMAL:
+               u.d *= 0x1p514;
+               *decpt = u.bits.exp - (514 + DBL_ADJ);
+               break;
+       case FP_INFINITE:
+               *decpt = INT_MAX;
+               return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
+       default:        /* FP_NAN or unrecognized */
+               *decpt = INT_MAX;
+               return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
+       }
+
+       /* FP_NORMAL or FP_SUBNORMAL */
+
+       if (ndigits == 0)               /* dtoa() compatibility */
+               ndigits = 1;
+
+       /*
+        * If ndigits < 0, we are expected to auto-size, so we allocate
+        * enough space for all the digits.
+        */
+       bufsize = (ndigits > 0) ? ndigits : SIGFIGS;
+       s0 = rv_alloc(bufsize);
+
+       /* Round to the desired number of digits. */
+       if (SIGFIGS > ndigits && ndigits > 0) {
+               float redux = one[u.bits.sign];
+               int offset = 4 * ndigits + DBL_MAX_EXP - 4 - DBL_MANT_DIG;
+               u.bits.exp = offset;
+               u.d += redux;
+               u.d -= redux;
+               *decpt += u.bits.exp - offset;
+       }
+
+       manh = u.bits.manh;
+       manl = u.bits.manl;
+       *s0 = '1';
+       for (s = s0 + 1; s < s0 + bufsize; s++) {
+               *s = xdigs[(manh >> (DBL_MANH_SIZE - 4)) & 0xf];
+               manh = (manh << 4) | (manl >> (DBL_MANL_SIZE - 4));
+               manl <<= 4;
+       }
+
+       /* If ndigits < 0, we are expected to auto-size the precision. */
+       if (ndigits < 0) {
+               for (ndigits = SIGFIGS; s0[ndigits - 1] == '0'; ndigits--)
+                       ;
+       }
+
+       s = s0 + ndigits;
+       *s = '\0';
+       if (rve != NULL)
+               *rve = s;
+       return (s0);
+}
diff --git a/lib/libc/gdtoa/_hldtoa.c b/lib/libc/gdtoa/_hldtoa.c
new file mode 100644 (file)
index 0000000..5efd007
--- /dev/null
@@ -0,0 +1,176 @@
+/*-
+ * Copyright (c) 2004-2008 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * $FreeBSD: src/lib/libc/gdtoa/_hldtoa.c,v 1.2 2008/04/12 14:53:52 das Exp $
+ */
+
+#include <float.h>
+#include <limits.h>
+#include <math.h>
+#include <stdint.h>
+
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "../stdio/floatio.h"
+#include "fpmath.h"
+#include "gdtoaimp.h"
+
+#if (LDBL_MANT_DIG > DBL_MANT_DIG)
+
+/* Strings values used by dtoa() */
+#define        INFSTR  "Infinity"
+#define        NANSTR  "NaN"
+
+#ifdef LDBL_IMPLICIT_NBIT
+#define        MANH_SIZE       LDBL_MANH_SIZE
+#else
+#define        MANH_SIZE       (LDBL_MANH_SIZE - 1)
+#endif
+
+#if MANH_SIZE > 32
+typedef uint64_t manh_t;
+#else
+typedef uint32_t manh_t;
+#endif
+
+#if LDBL_MANL_SIZE > 32
+typedef uint64_t manl_t;
+#else
+typedef uint32_t manl_t;
+#endif
+
+#define        LDBL_ADJ        (LDBL_MAX_EXP - 2)
+#define        SIGFIGS         ((LDBL_MANT_DIG + 3) / 4 + 1)
+
+static const float one[] = { 1.0f, -1.0f };
+
+/*
+ * This is the long double version of __hdtoa().
+ */
+char *
+__hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign,
+    char **rve)
+{
+       union IEEEl2bits u;
+       char *s, *s0;
+       manh_t manh;
+       manl_t manl;
+       int bufsize;
+#ifdef __i386__
+       fp_prec_t oldprec;
+#endif
+
+       u.e = e;
+       *sign = u.bits.sign;
+
+       switch (fpclassify(e)) {
+       case FP_NORMAL:
+               *decpt = u.bits.exp - LDBL_ADJ;
+               break;
+       case FP_ZERO:
+               *decpt = 1;
+               return (nrv_alloc("0", rve, 1));
+       case FP_SUBNORMAL:
+#ifdef __i386__
+               oldprec = fpsetprec(FP_PE);
+#endif
+               u.e *= 0x1p514L;
+               *decpt = u.bits.exp - (514 + LDBL_ADJ);
+#ifdef __i386__
+               fpsetprec(oldprec);
+#endif
+               break;
+       case FP_INFINITE:
+               *decpt = INT_MAX;
+               return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
+       default:        /* FP_NAN or unrecognized */
+               *decpt = INT_MAX;
+               return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
+       }
+
+       /* FP_NORMAL or FP_SUBNORMAL */
+
+       if (ndigits == 0)               /* dtoa() compatibility */
+               ndigits = 1;
+
+       /*
+        * If ndigits < 0, we are expected to auto-size, so we allocate
+        * enough space for all the digits.
+        */
+       bufsize = (ndigits > 0) ? ndigits : SIGFIGS;
+       s0 = rv_alloc(bufsize);
+
+       /* Round to the desired number of digits. */
+       if (SIGFIGS > ndigits && ndigits > 0) {
+               float redux = one[u.bits.sign];
+               int offset = 4 * ndigits + LDBL_MAX_EXP - 4 - LDBL_MANT_DIG;
+#ifdef __i386__
+               oldprec = fpsetprec(FP_PE);
+#endif
+               u.bits.exp = offset;
+               u.e += redux;
+               u.e -= redux;
+               *decpt += u.bits.exp - offset;
+#ifdef __i386__
+               fpsetprec(oldprec);
+#endif
+       }
+
+       mask_nbit_l(u);
+       manh = u.bits.manh;
+       manl = u.bits.manl;
+       *s0 = '1';
+       for (s = s0 + 1; s < s0 + bufsize; s++) {
+               *s = xdigs[(manh >> (MANH_SIZE - 4)) & 0xf];
+               manh = (manh << 4) | (manl >> (LDBL_MANL_SIZE - 4));
+               manl <<= 4;
+       }
+
+       /* If ndigits < 0, we are expected to auto-size the precision. */
+       if (ndigits < 0) {
+               for (ndigits = SIGFIGS; s0[ndigits - 1] == '0'; ndigits--)
+                       ;
+       }
+
+       s = s0 + ndigits;
+       *s = '\0';
+       if (rve != NULL)
+               *rve = s;
+       return (s0);
+}
+
+#else  /* (LDBL_MANT_DIG == DBL_MANT_DIG) */
+
+char *
+__hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign,
+    char **rve)
+{
+
+       return (__hdtoa((double)e, xdigs, ndigits, decpt, sign, rve));
+}
+
+#endif /* (LDBL_MANT_DIG == DBL_MANT_DIG) */
diff --git a/lib/libc/gdtoa/_ldtoa.c b/lib/libc/gdtoa/_ldtoa.c
new file mode 100644 (file)
index 0000000..05aa83f
--- /dev/null
@@ -0,0 +1,106 @@
+/*-
+ * Copyright (c) 2003 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * $FreeBSD: src/lib/libc/gdtoa/_ldtoa.c,v 1.5 2007/12/09 19:48:57 das Exp $
+ */
+
+#include <float.h>
+#include <inttypes.h>
+#include <limits.h>
+#include <math.h>
+#include <stdlib.h>
+#include "fpmath.h"
+#include "gdtoaimp.h"
+
+/*
+ * ldtoa() is a wrapper for gdtoa() that makes it smell like dtoa(),
+ * except that the floating point argument is passed by reference.
+ * When dtoa() is passed a NaN or infinity, it sets expt to 9999.
+ * However, a long double could have a valid exponent of 9999, so we
+ * use INT_MAX in ldtoa() instead.
+ */
+char *
+__ldtoa(long double *ld, int mode, int ndigits, int *decpt, int *sign,
+    char **rve)
+{
+       FPI fpi = {
+               LDBL_MANT_DIG,                  /* nbits */
+               LDBL_MIN_EXP - LDBL_MANT_DIG,   /* emin */
+               LDBL_MAX_EXP - LDBL_MANT_DIG,   /* emax */
+               FLT_ROUNDS,                     /* rounding */
+#ifdef Sudden_Underflow        /* unused, but correct anyway */
+               1
+#else
+               0
+#endif
+       };
+       int be, kind;
+       char *ret;
+       union IEEEl2bits u;
+       uint32_t bits[(LDBL_MANT_DIG + 31) / 32];
+       void *vbits = bits;
+
+       u.e = *ld;
+
+       /*
+        * gdtoa doesn't know anything about the sign of the number, so
+        * if the number is negative, we need to swap rounding modes of
+        * 2 (upwards) and 3 (downwards).
+        */
+       *sign = u.bits.sign;
+       fpi.rounding ^= (fpi.rounding >> 1) & u.bits.sign;
+
+       be = u.bits.exp - (LDBL_MAX_EXP - 1) - (LDBL_MANT_DIG - 1);
+       LDBL_TO_ARRAY32(u, bits);
+
+       switch (fpclassify(u.e)) {
+       case FP_NORMAL:
+               kind = STRTOG_Normal;
+#ifdef LDBL_IMPLICIT_NBIT
+               bits[LDBL_MANT_DIG / 32] |= 1 << ((LDBL_MANT_DIG - 1) % 32);
+#endif /* LDBL_IMPLICIT_NBIT */
+               break;
+       case FP_ZERO:
+               kind = STRTOG_Zero;
+               break;
+       case FP_SUBNORMAL:
+               kind = STRTOG_Denormal;
+               be++;
+               break;
+       case FP_INFINITE:
+               kind = STRTOG_Infinite;
+               break;
+       case FP_NAN:
+               kind = STRTOG_NaN;
+               break;
+       default:
+               abort();
+       }
+
+       ret = gdtoa(&fpi, be, vbits, &kind, mode, ndigits, decpt, rve);
+       if (*decpt == -32768)
+               *decpt = INT_MAX;
+       return ret;
+}
diff --git a/lib/libc/gdtoa/glue.c b/lib/libc/gdtoa/glue.c
new file mode 100644 (file)
index 0000000..5178f59
--- /dev/null
@@ -0,0 +1,13 @@
+/*
+ * Machine-independent glue to integrate David Gay's gdtoa
+ * package into libc.
+ *
+ * $FreeBSD: src/lib/libc/gdtoa/glue.c,v 1.2 2003/06/21 08:20:14 das Exp $
+ */
+
+#include <pthread.h>
+
+pthread_mutex_t __gdtoa_locks[] = {
+       PTHREAD_MUTEX_INITIALIZER,
+       PTHREAD_MUTEX_INITIALIZER
+};
diff --git a/lib/libc/gdtoa/machdep_ldisx.c b/lib/libc/gdtoa/machdep_ldisx.c
new file mode 100644 (file)
index 0000000..aa49110
--- /dev/null
@@ -0,0 +1,46 @@
+/*-
+ * Copyright (c) 2003 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * $FreeBSD: src/lib/libc/gdtoa/machdep_ldisx.c,v 1.3 2007/12/03 07:17:32 das Exp $
+ */
+
+/*
+ * Machine-dependent glue to integrate David Gay's gdtoa
+ * package into libc for architectures where a long double
+ * is an IEEE extended precision number.
+ */
+
+#include <float.h>
+
+#include "gdtoaimp.h"
+
+long double
+strtold(const char * __restrict s, char ** __restrict sp)
+{
+       long double result;
+
+       strtorx(s, sp, FLT_ROUNDS, &result);
+       return result;
+}
diff --git a/lib/libc/i386/Makefile.inc b/lib/libc/i386/Makefile.inc
new file mode 100644 (file)
index 0000000..d96e1a5
--- /dev/null
@@ -0,0 +1,5 @@
+# $FreeBSD: src/lib/libc/i386/Makefile.inc,v 1.3 2007/12/03 07:17:32 das Exp $
+
+# Long double is 80 bits
+GDTOASRCS+=strtorx.c
+MDSRCS+=machdep_ldisx.c
diff --git a/lib/libc/i386/_fpmath.h b/lib/libc/i386/_fpmath.h
new file mode 100644 (file)
index 0000000..455631c
--- /dev/null
@@ -0,0 +1,54 @@
+/*-
+ * Copyright (c) 2002, 2003 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * $FreeBSD: src/lib/libc/i386/_fpmath.h,v 1.6 2008/01/17 16:39:06 bde Exp $
+ */
+
+union IEEEl2bits {
+       long double     e;
+       struct {
+               unsigned int    manl    :32;
+               unsigned int    manh    :32;
+               unsigned int    exp     :15;
+               unsigned int    sign    :1;
+               unsigned int    junk    :16;
+       } bits;
+       struct {
+               unsigned long long man  :64;
+               unsigned int    expsign :16;
+               unsigned int    junk    :16;
+       } xbits;
+};
+
+#define        LDBL_NBIT       0x80000000
+#define        mask_nbit_l(u)  ((u).bits.manh &= ~LDBL_NBIT)
+
+#define        LDBL_MANH_SIZE  32
+#define        LDBL_MANL_SIZE  32
+
+#define        LDBL_TO_ARRAY32(u, a) do {                      \
+       (a)[0] = (uint32_t)(u).bits.manl;               \
+       (a)[1] = (uint32_t)(u).bits.manh;               \
+} while (0)
diff --git a/lib/libc/i386/arith.h b/lib/libc/i386/arith.h
new file mode 100644 (file)
index 0000000..c55a4df
--- /dev/null
@@ -0,0 +1,15 @@
+/*
+ * MD header for contrib/gdtoa
+ *
+ * $FreeBSD: src/lib/libc/i386/arith.h,v 1.2 2003/05/08 13:50:43 das Exp $
+ */
+
+/*
+ * NOTE: The definitions in this file must be correct or strtod(3) and
+ * floating point formats in printf(3) will break!  The file can be
+ * generated by running contrib/gdtoa/arithchk.c on the target
+ * architecture.  See contrib/gdtoa/gdtoaimp.h for details.
+ */
+
+#define IEEE_8087
+#define Arith_Kind_ASL 1
diff --git a/lib/libc/i386/gd_qnan.h b/lib/libc/i386/gd_qnan.h
new file mode 100644 (file)
index 0000000..8273014
--- /dev/null
@@ -0,0 +1,21 @@
+/*
+ * MD header for contrib/gdtoa
+ *
+ * This file can be generated by compiling and running contrib/gdtoa/qnan.c
+ * on the target architecture after arith.h has been generated.
+ *
+ * $FreeBSD: src/lib/libc/i386/gd_qnan.h,v 1.2 2007/12/16 21:15:08 das Exp $
+ */
+
+#define f_QNAN 0x7fc00000
+#define d_QNAN0 0x0
+#define d_QNAN1 0x7ff80000
+#define ld_QNAN0 0x0
+#define ld_QNAN1 0xc0000000
+#define ld_QNAN2 0x7fff
+#define ld_QNAN3 0x0
+#define ldus_QNAN0 0x0
+#define ldus_QNAN1 0x0
+#define ldus_QNAN2 0x0
+#define ldus_QNAN3 0xc000
+#define ldus_QNAN4 0x7fff
diff --git a/lib/libc/include/fpmath.h b/lib/libc/include/fpmath.h
new file mode 100644 (file)
index 0000000..b9c13a9
--- /dev/null
@@ -0,0 +1,75 @@
+/*-
+ * Copyright (c) 2003 Mike Barcroft <mike@FreeBSD.org>
+ * Copyright (c) 2002 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * $FreeBSD: src/lib/libc/include/fpmath.h,v 1.4 2008/12/23 22:20:59 marcel Exp $
+ */
+
+#include <sys/endian.h>
+#include "_fpmath.h"
+
+#ifndef _IEEE_WORD_ORDER
+#define        _IEEE_WORD_ORDER        _BYTE_ORDER
+#endif
+
+union IEEEf2bits {
+       float   f;
+       struct {
+#if _BYTE_ORDER == _LITTLE_ENDIAN
+               unsigned int    man     :23;
+               unsigned int    exp     :8;
+               unsigned int    sign    :1;
+#else /* _BIG_ENDIAN */
+               unsigned int    sign    :1;
+               unsigned int    exp     :8;
+               unsigned int    man     :23;
+#endif
+       } bits;
+};
+
+#define        DBL_MANH_SIZE   20
+#define        DBL_MANL_SIZE   32
+
+union IEEEd2bits {
+       double  d;
+       struct {
+#if _BYTE_ORDER == _LITTLE_ENDIAN
+#if _IEEE_WORD_ORDER == _LITTLE_ENDIAN
+               unsigned int    manl    :32;
+#endif
+               unsigned int    manh    :20;
+               unsigned int    exp     :11;
+               unsigned int    sign    :1;
+#if _IEEE_WORD_ORDER == _BIG_ENDIAN
+               unsigned int    manl    :32;
+#endif
+#else /* _BIG_ENDIAN */
+               unsigned int    sign    :1;
+               unsigned int    exp     :11;
+               unsigned int    manh    :20;
+               unsigned int    manl    :32;
+#endif
+       } bits;
+};
index c66bb9f..f8ba337 100644 (file)
@@ -419,7 +419,10 @@ vfprintf(FILE *fp, const char *fmt0, va_list ap)
 #define        BUF             ((MAXEXPDIG * 2) + MAXFRACT + 1) /* + decimal point */
 #define        DEFPREC         6
 
-static char *cvt (double, int, int, char *, int *, int, int *, char **);
+extern char *__dtoa(double, int, int, int *, int *, char **);
+extern void __freedtoa(char *s);
+
+static char *cvt(double, int, int, char *, int *, int, int *);
 static int exponent (char *, int, int);
 
 #else /* no NO_FLOATING_POINT */
@@ -854,11 +857,11 @@ fp_begin:         if (prec == -1)
                        }
                        flags |= FPT;
                        if (dtoaresult != NULL) {
-                               free(dtoaresult);
+                               __freedtoa(dtoaresult);
                                dtoaresult = NULL;
                        }
-                       cp = cvt(_double, prec, flags, &softsign,
-                               &expt, ch, &ndig, &dtoaresult);
+                       dtoaresult = cp = cvt(_double, prec, flags, &softsign,
+                               &expt, ch, &ndig);
                        if (ch == 'g' || ch == 'G') {
                                if (expt <= -4 || expt > prec)
                                        ch = (ch == 'g') ? 'e' : 'E';
@@ -1152,7 +1155,7 @@ done:
 error:
 #ifndef NO_FLOATING_POINT
        if (dtoaresult != NULL)
-               free(dtoaresult);
+               __freedtoa(dtoaresult);
 #endif
        if (convbuf != NULL)
                free(convbuf);
@@ -1518,11 +1521,9 @@ __grow_type_table(int nextarg, enum typeid **typetable, int *tablesize)
 
 #ifndef NO_FLOATING_POINT
 
-extern char *__dtoa (double, int, int, int *, int *, char **, char **);
-
 static char *
 cvt(double value, int ndigits, int flags, char *sign, int *decpt,
-    int ch, int *length, char **dtoaresultp)
+    int ch, int *length)
 {
        int mode, dsgn;
        char *digits, *bp, *rve;
@@ -1539,7 +1540,7 @@ cvt(double value, int ndigits, int flags, char *sign, int *decpt,
                        ndigits++;
                mode = 2;               /* ndigits significant digits */
        }
-       digits = __dtoa(value, mode, ndigits, decpt, &dsgn, &rve, dtoaresultp);
+       digits = __dtoa(value, mode, ndigits, decpt, &dsgn, &rve);
        *sign = dsgn != 0;
        if ((ch != 'g' && ch != 'G') || flags & ALT) {
                /* print trailing zeros */
index eb3f737..a6da14a 100644 (file)
@@ -10,7 +10,7 @@ MISRCS+=a64l.c abort.c abs.c atexit.c atof.c atoi.c atol.c atoll.c \
        getsubopt.c hcreate.c heapsort.c imaxabs.c imaxdiv.c \
        insque.c l64a.c labs.c ldiv.c llabs.c lldiv.c lsearch.c malloc.c \
        merge.c qsort.c qsort_r.c radixsort.c rand.c random.c \
-       reallocf.c realpath.c remque.c strfmon.c strtod.c strtoimax.c \
+       reallocf.c realpath.c remque.c strfmon.c strtoimax.c \
        strtol.c strtoll.c strtonum.c strtoq.c strtoul.c strtoull.c \
        strtoumax.c strtouq.c system.c tdelete.c tfind.c tsearch.c twalk.c
 
index b7375c9..02f8556 100644 (file)
 .\" 2. Redistributions in binary form must reproduce the above copyright
 .\"    notice, this list of conditions and the following disclaimer in the
 .\"    documentation and/or other materials provided with the distribution.
-.\" 3. All advertising materials mentioning features or use of this software
-.\"    must display the following acknowledgement:
-.\"    This product includes software developed by the University of
-.\"    California, Berkeley and its contributors.
 .\" 4. Neither the name of the University nor the names of its contributors
 .\"    may be used to endorse or promote products derived from this software
 .\"    without specific prior written permission.
 .\" SUCH DAMAGE.
 .\"
 .\"     @(#)strtod.3   8.1 (Berkeley) 6/4/93
-.\" $FreeBSD: src/lib/libc/stdlib/strtod.3,v 1.4.2.7 2001/12/14 18:33:58 ru Exp $
+.\" $FreeBSD: src/lib/libc/stdlib/strtod.3,v 1.21 2007/01/09 00:28:10 imp Exp $
 .\" $DragonFly: src/lib/libc/stdlib/strtod.3,v 1.2 2003/06/17 04:26:46 dillon Exp $
 .\"
-.Dd June 4, 1993
+.Dd March 2, 2003
 .Dt STRTOD 3
 .Os
 .Sh NAME
-.Nm strtod
+.Nm strtod ,
+.Nm strtof ,
+.Nm strtold
 .Nd convert
 .Tn ASCII
-string to double
+string to floating point
 .Sh LIBRARY
 .Lb libc
 .Sh SYNOPSIS
 .In stdlib.h
 .Ft double
-.Fn strtod "const char *nptr" "char **endptr"
+.Fn strtod "const char * restrict nptr" "char ** restrict endptr"
+.Ft float
+.Fn strtof "const char * restrict nptr" "char ** restrict endptr"
+.Ft "long double"
+.Fn strtold "const char * restrict nptr" "char ** restrict endptr"
 .Sh DESCRIPTION
-The
-.Fn strtod
-function converts the initial portion of the string
+These conversion
+functions convert the initial portion of the string
 pointed to by
 .Fa nptr
 to
-.Em double
-representation.
+.Vt double ,
+.Vt float ,
+and
+.Vt "long double"
+representation, respectively.
 .Pp
 The expected form of the string is an optional plus (``+'') or minus
-sign (``\-'') followed by a sequence of digits optionally containing
-a decimal-point character, optionally followed by an exponent.
-An exponent consists of an ``E'' or ``e'', followed by an optional plus
-or minus sign, followed by a sequence of digits.
+sign (``\-'') followed by either:
+.Bl -bullet
+.It
+a decimal significand consisting of a sequence of decimal digits
+optionally containing a decimal-point character, or
+.It
+a hexadecimal significand consisting of a ``0X'' or ``0x'' followed
+by a sequence of hexadecimal digits optionally containing a
+decimal-point character.
+.El
+.Pp
+In both cases, the significand may be optionally followed by an
+exponent.
+An exponent consists of an ``E'' or ``e'' (for decimal
+constants) or a ``P'' or ``p'' (for hexadecimal constants),
+followed by an optional plus or minus sign, followed by a
+sequence of decimal digits.
+For decimal constants, the exponent indicates the power of 10 by
+which the significand should be scaled.
+For hexadecimal constants, the scaling is instead done by powers
+of 2.
+.Pp
+Alternatively, if the portion of the string following the optional
+plus or minus sign begins with ``INFINITY'' or ``NAN'', ignoring
+case, it is interpreted as an infinity or a quiet NaN, respectively.
 .Pp
-Leading white-space characters in the string (as defined by the
+In any of the above cases, leading white-space characters in the
+string (as defined by the
 .Xr isspace 3
 function) are skipped.
+The decimal point
+character is defined in the program's locale (category
+.Dv LC_NUMERIC ) .
 .Sh RETURN VALUES
 The
-.Fn strtod
-function returns the converted value, if any.
+.Fn strtod ,
+.Fn strtof ,
+and
+.Fn strtold
+functions return the converted value, if any.
 .Pp
 If
 .Fa endptr
@@ -89,8 +121,11 @@ is stored in the location referenced by
 .Fa endptr .
 .Pp
 If the correct value would cause overflow, plus or minus
-.Dv HUGE_VAL
-is returned (according to the sign of the value), and
+.Dv HUGE_VAL ,
+.Dv HUGE_VALF ,
+or
+.Dv HUGE_VALL
+is returned (according to the sign and type of the return value), and
 .Er ERANGE
 is stored in
 .Va errno .
@@ -109,30 +144,41 @@ Overflow or underflow occurred.
 .Xr atoi 3 ,
 .Xr atol 3 ,
 .Xr strtol 3 ,
-.Xr strtoul 3
+.Xr strtoul 3 ,
+.Xr wcstod 3
 .Sh STANDARDS
 The
 .Fn strtod
 function
 conforms to
-.St -isoC .
+.St -isoC-99 ,
+with the exception of the bug noted below.
 .Sh AUTHORS
 The author of this software is
 .An David M. Gay .
 .Pp
-Copyright (c) 1991 by AT&T.
-.Pp
-Permission to use, copy, modify, and distribute this software for any
-purpose without fee is hereby granted, provided that this entire notice
-is included in all copies of any software which is or includes a copy
-or modification of this software and in all copies of the supporting
-documentation for such software.
-.Pp
-THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
-WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
-REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
-OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
-.Pp
-Contact your vendor for a free copy of the source code to
-.Fn strtod
-and accompanying functions.
+.Bd -literal
+Copyright (c) 1998 by 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.
+.Ed
+.Sh BUGS
+These routines do not recognize the C99 ``NaN(...)'' syntax.
diff --git a/lib/libc/stdlib/strtod.c b/lib/libc/stdlib/strtod.c
deleted file mode 100644 (file)
index 6334a5b..0000000
+++ /dev/null
@@ -1,2098 +0,0 @@
-/*-
- * Copyright (c) 1993
- *     The Regents of the University of California.  All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice, this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- * 3. All advertising materials mentioning features or use of this software
- *    must display the following acknowledgement:
- *     This product includes software developed by the University of
- *     California, Berkeley and its contributors.
- * 4. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
- * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
- * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
- * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
- * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
- * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
- * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
- * SUCH DAMAGE.
- *
- * $FreeBSD: src/lib/libc/stdlib/strtod.c,v 1.3.8.4 2002/08/31 22:26:35 dwmalone Exp $
- * $DragonFly: src/lib/libc/stdlib/strtod.c,v 1.6 2006/09/10 21:22:32 swildner Exp $
- *
- * @(#)strtod.c        8.1 (Berkeley) 6/4/93
- */
-
-/****************************************************************
- *
- * The author of this software is David M. Gay.
- *
- * Copyright (c) 1991 by AT&T.
- *
- * Permission to use, copy, modify, and distribute this software for any
- * purpose without fee is hereby granted, provided that this entire notice
- * is included in all copies of any software which is or includes a copy
- * or modification of this software and in all copies of the supporting
- * documentation for such software.
- *
- * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
- * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
- * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
- * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
- *
- ***************************************************************/
-
-/* Please send bug reports to
-       David M. Gay
-       AT&T Bell Laboratories, Room 2C-463
-       600 Mountain Avenue
-       Murray Hill, NJ 07974-2070
-       U.S.A.
-       dmg@research.att.com or research!dmg
- */
-
-/* strtod for IEEE--arithmetic machines.
- *
- * This strtod returns a nearest machine number to the input decimal
- * string (or sets errno to ERANGE).  With IEEE arithmetic, ties are
- * broken by the IEEE round-even rule.  Otherwise ties are broken by
- * biased rounding (add half and chop).
- *
- * Inspired loosely by William D. Clinger's paper "How to Read Floating
- * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
- *
- * Modifications:
- *
- *     1. We only require IEEE double-precision arithmetic
- *             (not IEEE double-extended).
- *     2. We get by with floating-point arithmetic in a case that
- *             Clinger missed -- when we're computing d * 10^n
- *             for a small integer d and the integer n is not too
- *             much larger than 22 (the maximum integer k for which
- *             we can represent 10^k exactly), we may be able to
- *             compute (d*10^k) * 10^(e-k) with just one roundoff.
- *     3. Rather than a bit-at-a-time adjustment of the binary
- *             result in the hard case, we use floating-point
- *             arithmetic to determine the adjustment to within
- *             one bit; only in really hard cases do we need to
- *             compute a second residual.
- *     4. Because of 3., we don't need a large table of powers of 10
- *             for ten-to-e (just some small tables, e.g. of 10^k
- *             for 0 <= k <= 22).
- */
-
-/*
- * #define IEEE_8087 for IEEE-arithmetic machines where the least
- *     significant byte has the lowest address.
- * #define IEEE_MC68k for IEEE-arithmetic machines where the most
- *     significant byte has the lowest address.
- * #define Sudden_Underflow for IEEE-format machines without gradual
- *     underflow (i.e., that flush to zero on underflow).
- * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
- * #define No_leftright to omit left-right logic in fast floating-point
- *     computation of dtoa.
- * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
- * #define ROUND_BIASED for IEEE-format with biased rounding.
- * #define Inaccurate_Divide for IEEE-format with correctly rounded
- *     products but inaccurate quotients, e.g., for Intel i860.
- * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
- *     integer arithmetic.  Whether this speeds things up or slows things
- *     down depends on the machine and the number being converted.
- */
-
-#if defined(i386) || defined(mips) && defined(MIPSEL)
-#define IEEE_8087
-#else
-#define IEEE_MC68k
-#endif
-
-#ifdef DEBUG
-#include "stdio.h"
-#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
-#endif
-
-#include <locale.h>
-#ifdef __cplusplus
-#include "malloc.h"
-#include "memory.h"
-#else
-#include "stdlib.h"
-#include "string.h"
-#endif
-
-#include "errno.h"
-#include <ctype.h>
-#include "float.h"
-
-#ifndef __MATH_H__
-#include "math.h"
-#endif
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-#ifndef CONST
-#define CONST const
-#endif
-
-#ifdef Unsigned_Shifts
-#define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
-#else
-#define Sign_Extend(a,b) /*no-op*/
-#endif
-
-#if defined(IEEE_8087) + defined(IEEE_MC68k) != 1
-Exactly one of IEEE_8087 or IEEE_MC68k should be defined.
-#endif
-
-union doubleasulongs {
-       double x;
-       unsigned long w[2];
-};
-#ifdef IEEE_8087
-#define word0(x) (((union doubleasulongs *)&x)->w)[1]
-#define word1(x) (((union doubleasulongs *)&x)->w)[0]
-#else
-#define word0(x) (((union doubleasulongs *)&x)->w)[0]
-#define word1(x) (((union doubleasulongs *)&x)->w)[1]
-#endif
-
-/* The following definition of Storeinc is appropriate for MIPS processors.
- * An alternative that might be better on some machines is
- * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
- */
-#if defined(IEEE_8087)
-#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
-((unsigned short *)a)[0] = (unsigned short)c, a++)
-#else
-#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
-((unsigned short *)a)[1] = (unsigned short)c, a++)
-#endif
-
-/* #define P DBL_MANT_DIG */
-/* Ten_pmax = floor(P*log(2)/log(5)) */
-/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
-/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
-/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
-
-#if defined(IEEE_8087) + defined(IEEE_MC68k)
-#define Exp_shift  20
-#define Exp_shift1 20
-#define Exp_msk1    0x100000
-#define Exp_msk11   0x100000
-#define Exp_mask  0x7ff00000
-#define P 53
-#define Bias 1023
-#define IEEE_Arith
-#define Emin (-1022)
-#define Exp_1  0x3ff00000
-#define Exp_11 0x3ff00000
-#define Ebits 11
-#define Frac_mask  0xfffff
-#define Frac_mask1 0xfffff
-#define Ten_pmax 22
-#define Bletch 0x10
-#define Bndry_mask  0xfffff
-#define Bndry_mask1 0xfffff
-#define LSB 1
-#define Sign_bit 0x80000000
-#define Log2P 1
-#define Tiny0 0
-#define Tiny1 1
-#define Quick_max 14
-#define Int_max 14
-#define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
-#endif
-
-#ifndef IEEE_Arith
-#define ROUND_BIASED
-#endif
-
-#define rounded_product(a,b) a *= b
-#define rounded_quotient(a,b) a /= b
-
-#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
-#define Big1 0xffffffff
-
-#ifndef Just_16
-/* When Pack_32 is not defined, we store 16 bits per 32-bit long.
- * This makes some inner loops simpler and sometimes saves work
- * during multiplications, but it often seems to make things slightly
- * slower.  Hence the default is now to store 32 bits per long.
- */
-#ifndef Pack_32
-#define Pack_32
-#endif
-#endif
-
-#define Kmax 15
-
-#ifdef __cplusplus
-extern "C" double strtod(const char *s00, char **se);
-extern "C" char *__dtoa(double d, int mode, int ndigits,
-                       int *decpt, int *sign, char **rve, char **resultp);
-#endif
-
-struct
-Bigint {
-       struct Bigint *next;
-       int k, maxwds, sign, wds;
-       unsigned long x[1];
-};
-
-typedef struct Bigint Bigint;
-
-static Bigint *
-Balloc(int k)
-{
-       int x;
-       Bigint *rv;
-
-       x = 1 << k;
-       rv = (Bigint *)malloc(sizeof(Bigint) + (x-1)*sizeof(long));
-       rv->k = k;
-       rv->maxwds = x;
-       rv->sign = rv->wds = 0;
-       return rv;
-}
-
-static void
-Bfree(Bigint *v)
-{
-       free(v);
-}
-
-#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
-y->wds*sizeof(long) + 2*sizeof(int))
-
-static Bigint *
-multadd(Bigint *b, int m, int a)       /* multiply by m and add a */
-{
-       int i, wds;
-       unsigned long *x, y;
-#ifdef Pack_32
-       unsigned long xi, z;
-#endif
-       Bigint *b1;
-
-       wds = b->wds;
-       x = b->x;
-       i = 0;
-       do {
-#ifdef Pack_32
-               xi = *x;
-               y = (xi & 0xffff) * m + a;
-               z = (xi >> 16) * m + (y >> 16);
-               a = (int)(z >> 16);
-               *x++ = (z << 16) + (y & 0xffff);
-#else
-               y = *x * m + a;
-               a = (int)(y >> 16);
-               *x++ = y & 0xffff;
-#endif
-       } while (++i < wds);
-       if (a) {
-               if (wds >= b->maxwds) {
-                       b1 = Balloc(b->k+1);
-                       Bcopy(b1, b);
-                       Bfree(b);
-                       b = b1;
-                       }
-               b->x[wds++] = a;
-               b->wds = wds;
-       }
-       return b;
-}
-
-static Bigint *
-s2b(CONST char *s, int nd0, int nd, unsigned long y9)
-{
-       Bigint *b;
-       int i, k;
-       long x, y;
-
-       x = (nd + 8) / 9;
-       for (k = 0, y = 1; x > y; y <<= 1, k++) ;
-#ifdef Pack_32
-       b = Balloc(k);
-       b->x[0] = y9;
-       b->wds = 1;
-#else
-       b = Balloc(k+1);
-       b->x[0] = y9 & 0xffff;
-       b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
-#endif
-
-       i = 9;
-       if (9 < nd0) {
-               s += 9;
-               do
-                       b = multadd(b, 10, *s++ - '0');
-               while (++i < nd0);
-               s++;
-       } else
-               s += 10;
-       for (; i < nd; i++)
-               b = multadd(b, 10, *s++ - '0');
-       return b;
-}
-
-static int
-hi0bits(unsigned long x)
-{
-       int k = 0;
-
-       if (!(x & 0xffff0000)) {
-               k = 16;
-               x <<= 16;
-       }
-       if (!(x & 0xff000000)) {
-               k += 8;
-               x <<= 8;
-       }
-       if (!(x & 0xf0000000)) {
-               k += 4;
-               x <<= 4;
-       }
-       if (!(x & 0xc0000000)) {
-               k += 2;
-               x <<= 2;
-       }
-       if (!(x & 0x80000000)) {
-               k++;
-               if (!(x & 0x40000000))
-                       return 32;
-       }
-       return k;
-}
-
-static int
-lo0bits(unsigned long *y)
-{
-       int k;
-       unsigned long x = *y;
-
-       if (x & 7) {
-               if (x & 1)
-                       return 0;
-               if (x & 2) {
-                       *y = x >> 1;
-                       return 1;
-               }
-               *y = x >> 2;
-               return 2;
-       }
-       k = 0;
-       if (!(x & 0xffff)) {
-               k = 16;
-               x >>= 16;
-       }
-       if (!(x & 0xff)) {
-               k += 8;
-               x >>= 8;
-       }
-       if (!(x & 0xf)) {
-               k += 4;
-               x >>= 4;
-       }
-       if (!(x & 0x3)) {
-               k += 2;
-               x >>= 2;
-       }
-       if (!(x & 1)) {
-               k++;
-               x >>= 1;
-               if (!x & 1)
-                       return 32;
-       }
-       *y = x;
-       return k;
-}
-
-static Bigint *
-i2b(int i)
-{
-       Bigint *b;
-
-       b = Balloc(1);
-       b->x[0] = i;
-       b->wds = 1;
-       return b;
-       }
-
-static Bigint *
-mult(Bigint *a, Bigint *b)
-{
-       Bigint *c;
-       int k, wa, wb, wc;
-       unsigned long carry, y, z;
-       unsigned long *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
-#ifdef Pack_32
-       unsigned long z2;
-#endif
-
-       if (a->wds < b->wds) {
-               c = a;
-               a = b;
-               b = c;
-       }
-       k = a->k;
-       wa = a->wds;
-       wb = b->wds;
-       wc = wa + wb;
-       if (wc > a->maxwds)
-               k++;
-       c = Balloc(k);
-       for (x = c->x, xa = x + wc; x < xa; x++)
-               *x = 0;
-       xa = a->x;
-       xae = xa + wa;
-       xb = b->x;
-       xbe = xb + wb;
-       xc0 = c->x;
-#ifdef Pack_32
-       for (; xb < xbe; xb++, xc0++) {
-               if ( (y = *xb & 0xffff) ) {
-                       x = xa;
-                       xc = xc0;
-                       carry = 0;
-                       do {
-                               z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
-                               carry = z >> 16;
-                               z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
-                               carry = z2 >> 16;
-                               Storeinc(xc, z2, z);
-                       } while (x < xae);
-                       *xc = carry;
-               }
-               if ( (y = *xb >> 16) ) {
-                       x = xa;
-                       xc = xc0;
-                       carry = 0;
-                       z2 = *xc;
-                       do {
-                               z = (*x & 0xffff) * y + (*xc >> 16) + carry;
-                               carry = z >> 16;
-                               Storeinc(xc, z, z2);
-                               z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
-                               carry = z2 >> 16;
-                       } while (x < xae);
-                       *xc = z2;
-               }
-       }
-#else
-       for (; xb < xbe; xc0++) {
-               if (y = *xb++) {
-                       x = xa;
-                       xc = xc0;
-                       carry = 0;
-                       do {
-                               z = *x++ * y + *xc + carry;
-                               carry = z >> 16;
-                               *xc++ = z & 0xffff;
-                       } while (x < xae);
-                       *xc = carry;
-               }
-       }
-#endif
-       for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
-       c->wds = wc;
-       return c;
-}
-
-static Bigint *p5s;
-
-static Bigint *
-pow5mult(Bigint *b, int k)
-{
-       Bigint *b1, *p5, *p51;
-       int i;
-       static int p05[3] = { 5, 25, 125 };
-
-       if ( (i = k & 3) )
-               b = multadd(b, p05[i-1], 0);
-
-       if (!(k >>= 2))
-               return b;
-       if (!(p5 = p5s)) {
-               /* first time */
-               p5 = p5s = i2b(625);
-               p5->next = 0;
-       }
-       for (;;) {
-               if (k & 1) {
-                       b1 = mult(b, p5);
-                       Bfree(b);
-                       b = b1;
-               }
-               if (!(k >>= 1))
-                       break;
-               if (!(p51 = p5->next)) {
-                       p51 = p5->next = mult(p5,p5);
-                       p51->next = 0;
-               }
-               p5 = p51;
-       }
-       return b;
-}
-
-static Bigint *
-lshift(Bigint *b, int k)
-{
-       int i, k1, n, n1;
-       Bigint *b1;
-       unsigned long *x, *x1, *xe, z;
-
-#ifdef Pack_32
-       n = k >> 5;
-#else
-       n = k >> 4;
-#endif
-       k1 = b->k;
-       n1 = n + b->wds + 1;
-       for (i = b->maxwds; n1 > i; i <<= 1)
-               k1++;
-       b1 = Balloc(k1);
-       x1 = b1->x;
-       for (i = 0; i < n; i++)
-               *x1++ = 0;
-       x = b->x;
-       xe = x + b->wds;
-#ifdef Pack_32
-       if (k &= 0x1f) {
-               k1 = 32 - k;
-               z = 0;
-               do {
-                       *x1++ = *x << k | z;
-                       z = *x++ >> k1;
-               } while (x < xe);
-               if ( (*x1 = z) )
-                       ++n1;
-       }
-#else
-       if (k &= 0xf) {
-               k1 = 16 - k;
-               z = 0;
-               do {
-                       *x1++ = *x << k  & 0xffff | z;
-                       z = *x++ >> k1;
-               } while (x < xe);
-               if (*x1 = z)
-                       ++n1;
-       }
-#endif
-       else
-               do
-                       *x1++ = *x++;
-               while (x < xe);
-       b1->wds = n1 - 1;
-       Bfree(b);
-       return b1;
-}
-
-static int
-cmp(Bigint *a, Bigint *b)
-{
-       unsigned long *xa, *xa0, *xb, *xb0;
-       int i, j;
-
-       i = a->wds;
-       j = b->wds;
-#ifdef DEBUG
-       if (i > 1 && !a->x[i-1])
-               Bug("cmp called with a->x[a->wds-1] == 0");
-       if (j > 1 && !b->x[j-1])
-               Bug("cmp called with b->x[b->wds-1] == 0");
-#endif
-       if (i -= j)
-               return i;
-       xa0 = a->x;
-       xa = xa0 + j;
-       xb0 = b->x;
-       xb = xb0 + j;
-       for (;;) {
-               if (*--xa != *--xb)
-                       return *xa < *xb ? -1 : 1;
-               if (xa <= xa0)
-                       break;
-       }
-       return 0;
-}
-
-static Bigint *
-diff(Bigint *a, Bigint *b)
-{
-       Bigint *c;
-       int i, wa, wb;
-       long borrow, y; /* We need signed shifts here. */
-       unsigned long *xa, *xae, *xb, *xbe, *xc;
-#ifdef Pack_32
-       long z;
-#endif
-
-       i = cmp(a,b);
-       if (!i) {
-               c = Balloc(0);
-               c->wds = 1;
-               c->x[0] = 0;
-               return c;
-       }
-       if (i < 0) {
-               c = a;
-               a = b;
-               b = c;
-               i = 1;
-       } else
-               i = 0;
-       c = Balloc(a->k);
-       c->sign = i;
-       wa = a->wds;
-       xa = a->x;
-       xae = xa + wa;
-       wb = b->wds;
-       xb = b->x;
-       xbe = xb + wb;
-       xc = c->x;
-       borrow = 0;
-#ifdef Pack_32
-       do {
-               y = (*xa & 0xffff) - (*xb & 0xffff) + borrow;
-               borrow = y >> 16;
-               Sign_Extend(borrow, y);
-               z = (*xa++ >> 16) - (*xb++ >> 16) + borrow;
-               borrow = z >> 16;
-               Sign_Extend(borrow, z);
-               Storeinc(xc, z, y);
-       } while (xb < xbe);
-       while (xa < xae) {
-               y = (*xa & 0xffff) + borrow;
-               borrow = y >> 16;
-               Sign_Extend(borrow, y);
-               z = (*xa++ >> 16) + borrow;
-               borrow = z >> 16;
-               Sign_Extend(borrow, z);
-               Storeinc(xc, z, y);
-       }
-#else
-       do {
-               y = *xa++ - *xb++ + borrow;
-               borrow = y >> 16;
-               Sign_Extend(borrow, y);
-               *xc++ = y & 0xffff;
-       } while (xb < xbe);
-       while (xa < xae) {
-               y = *xa++ + borrow;
-               borrow = y >> 16;
-               Sign_Extend(borrow, y);
-               *xc++ = y & 0xffff;
-       }
-#endif
-       while (!*--xc)
-               wa--;
-       c->wds = wa;
-       return c;
-}
-
-static double
-ulp(double x)
-{
-       long L;
-       double a;
-
-       L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
-#ifndef Sudden_Underflow
-       if (L > 0) {
-#endif
-               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;
-               } else {
-                       word0(a) = 0;
-                       L -= Exp_shift;
-                       word1(a) = L >= 31 ? 1 : 1 << (31 - L);
-               }
-       }
-#endif
-       return a;
-}
-
-static double
-b2d(Bigint *a, int *e)
-{
-       unsigned long *xa, *xa0, w, y, z;
-       int k;
-       double d;
-#define d0 word0(d)
-#define d1 word1(d)
-
-       xa0 = a->x;
-       xa = xa0 + a->wds;
-       y = *--xa;
-#ifdef DEBUG
-       if (!y) Bug("zero y in b2d");
-#endif
-       k = hi0bits(y);
-       *e = 32 - k;
-#ifdef Pack_32
-       if (k < Ebits) {
-               d0 = Exp_1 | (y >> (Ebits - k));
-               w = xa > xa0 ? *--xa : 0;
-               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));
-               y = xa > xa0 ? *--xa : 0;
-               d1 = (z << k) | (y >> (32 - k));
-       } else {
-               d0 = Exp_1 | y;
-               d1 = z;
-       }
-#else
-       if (k < Ebits + 16) {
-               z = xa > xa0 ? *--xa : 0;
-               d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
-               w = xa > xa0 ? *--xa : 0;
-               y = xa > xa0 ? *--xa : 0;
-               d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
-               goto ret_d;
-       }
-       z = xa > xa0 ? *--xa : 0;
-       w = xa > xa0 ? *--xa : 0;
-       k -= Ebits + 16;
-       d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
-       y = xa > xa0 ? *--xa : 0;
-       d1 = w << k + 16 | y << k;
-#endif
-ret_d:
-#undef d0
-#undef d1
-       return d;
-}
-
-static Bigint *
-d2b(double d, int *e, int *bits)
-{
-       Bigint *b;
-       int de, i, k;
-       unsigned long *x, y, z;
-#define d0 word0(d)
-#define d1 word1(d)
-
-#ifdef Pack_32
-       b = Balloc(1);
-#else
-       b = Balloc(2);
-#endif
-       x = b->x;
-
-       z = d0 & Frac_mask;
-       d0 &= 0x7fffffff;       /* clear sign bit, which we ignore */
-#ifdef Sudden_Underflow
-       de = (int)(d0 >> Exp_shift);
-       z |= Exp_msk11;
-#else
-       if ( (de = (int)(d0 >> Exp_shift)) )
-               z |= Exp_msk1;
-#endif
-#ifdef Pack_32
-       if ( (y = d1) ) {
-               if ( (k = lo0bits(&y)) ) {
-                       x[0] = y | (z << (32 - k));
-                       z >>= k;
-                       }
-               else
-                       x[0] = y;
-               i = b->wds = (x[1] = z) ? 2 : 1;
-       } else {
-#ifdef DEBUG
-               if (!z)
-                       Bug("Zero passed to d2b");
-#endif
-               k = lo0bits(&z);
-               x[0] = z;
-               i = b->wds = 1;
-               k += 32;
-       }
-#else
-       if (y = d1) {
-               if (k = lo0bits(&y))
-                       if (k >= 16) {
-                               x[0] = y | z << 32 - k & 0xffff;
-                               x[1] = z >> k - 16 & 0xffff;
-                               x[2] = z >> k;
-                               i = 2;
-                       } else {
-                               x[0] = y & 0xffff;
-                               x[1] = y >> 16 | z << 16 - k & 0xffff;
-                               x[2] = z >> k & 0xffff;
-                               x[3] = z >> k+16;
-                               i = 3;
-                       }
-               else {
-                       x[0] = y & 0xffff;
-                       x[1] = y >> 16;
-                       x[2] = z & 0xffff;
-                       x[3] = z >> 16;
-                       i = 3;
-               }
-       } else {
-#ifdef DEBUG
-               if (!z)
-                       Bug("Zero passed to d2b");
-#endif
-               k = lo0bits(&z);
-               if (k >= 16) {
-                       x[0] = z;
-                       i = 0;
-               } else {
-                       x[0] = z & 0xffff;
-                       x[1] = z >> 16;
-                       i = 1;
-               }
-               k += 32;
-       }
-       while (!x[i])
-               --i;
-       b->wds = i + 1;
-#endif
-#ifndef Sudden_Underflow
-       if (de) {
-#endif
-               *e = de - Bias - (P-1) + k;
-               *bits = P - k;
-#ifndef Sudden_Underflow
-       } else {
-               *e = de - Bias - (P-1) + 1 + k;
-#ifdef Pack_32
-               *bits = 32*i - hi0bits(x[i-1]);
-#else
-               *bits = (i+2)*16 - hi0bits(x[i]);
-#endif
-       }
-#endif
-       return b;
-}
-#undef d0
-#undef d1
-
-static double
-ratio(Bigint *a, Bigint *b)
-{
-       double da, db;
-       int k, ka, kb;
-
-       da = b2d(a, &ka);
-       db = b2d(b, &kb);
-#ifdef Pack_32
-       k = ka - kb + 32*(a->wds - b->wds);
-#else
-       k = ka - kb + 16*(a->wds - b->wds);
-#endif
-       if (k > 0)
-               word0(da) += k*Exp_msk1;
-       else {
-               k = -k;
-               word0(db) += k*Exp_msk1;
-       }
-       return da / db;
-}
-
-static double
-tens[] = {
-               1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
-               1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
-               1e20, 1e21, 1e22
-               };
-
-static double
-#ifdef IEEE_Arith
-bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
-static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
-#define n_bigtens 5
-#else
-bigtens[] = { 1e16, 1e32 };
-static double tinytens[] = { 1e-16, 1e-32 };
-#define n_bigtens 2
-#endif
-
-double
-strtod(CONST char *s00, char **se)
-{
-       int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
-                e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
-       CONST char *s, *s0, *s1;
-       double aadj, aadj1, adj, rv, rv0;
-       long L;
-       unsigned long y, z;
-       Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
-       char decimal_point = localeconv()->decimal_point[0];
-
-       sign = nz0 = nz = 0;
-       rv = 0.;
-       for (s = s00;;s++) switch(*s) {
-               case '-':
-                       sign = 1;
-                       /* no break */
-               case '+':
-                       if (*++s)
-                               goto break2;
-                       /* no break */
-               case 0:
-                       s = s00;
-                       goto ret;
-               default:
-                       if (isspace((unsigned char)*s))
-                               continue;
-                       goto break2;
-       }
-break2:
-       if (*s == '0') {
-               nz0 = 1;
-               while (*++s == '0') ;
-               if (!*s)
-                       goto ret;
-       }
-       s0 = s;
-       y = z = 0;
-       for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
-               if (nd < 9)
-                       y = 10*y + c - '0';
-               else if (nd < 16)
-                       z = 10*z + c - '0';
-       nd0 = nd;
-       if ((char)c == decimal_point) {
-               c = *++s;
-               if (!nd) {
-                       for (; c == '0'; c = *++s)
-                               nz++;
-                       if (c > '0' && c <= '9') {
-                               s0 = s;
-                               nf += nz;
-                               nz = 0;
-                               goto have_dig;
-                       }
-                       goto dig_done;
-               }
-               for (; c >= '0' && c <= '9'; c = *++s) {
-have_dig:
-                       nz++;
-                       if (c -= '0') {
-                               nf += nz;
-                               for (i = 1; i < nz; i++)
-                                       if (nd++ < 9)
-                                               y *= 10;
-                                       else if (nd <= DBL_DIG + 1)
-                                               z *= 10;
-                               if (nd++ < 9)
-                                       y = 10*y + c;
-                               else if (nd <= DBL_DIG + 1)
-                                       z = 10*z + c;
-                               nz = 0;
-                       }
-               }
-       }
-dig_done:
-       e = 0;
-       if (c == 'e' || c == 'E') {
-               if (!nd && !nz && !nz0) {
-                       s = s00;
-                       goto ret;
-               }
-               s00 = s;
-               esign = 0;
-               switch(c = *++s) {
-                       case '-':
-                               esign = 1;
-                       case '+':
-                               c = *++s;
-               }
-               if (c >= '0' && c <= '9') {
-                       while (c == '0')
-                               c = *++s;
-                       if (c > '0' && c <= '9') {
-                               L = c - '0';
-                               s1 = s;
-                               while ((c = *++s) >= '0' && c <= '9')
-                                       L = 10*L + c - '0';
-                               if (s - s1 > 8 || L > 19999)
-                                       /* Avoid confusion from exponents
-                                        * so large that e might overflow.
-                                        */
-                                       e = 19999; /* safe for 16 bit ints */
-                               else
-                                       e = (int)L;
-                               if (esign)
-                                       e = -e;
-                       } else
-                               e = 0;
-               } else
-                       s = s00;
-       }
-       if (!nd) {
-               if (!nz && !nz0)
-                       s = s00;
-               goto ret;
-       }
-       e1 = e -= nf;
-
-       /* Now we have nd0 digits, starting at s0, followed by a
-        * decimal point, followed by nd-nd0 digits.  The number we're
-        * after is the integer represented by those digits times
-        * 10**e */
-
-       if (!nd0)
-               nd0 = nd;
-       k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
-       rv = y;
-       if (k > 9)
-               rv = tens[k - 9] * rv + z;
-       if (nd <= DBL_DIG) {
-               if (!e)
-                       goto ret;
-               if (e > 0) {
-                       if (e <= Ten_pmax) {
-                               /* rv = */ rounded_product(rv, tens[e]);
-                               goto ret;
-                               }
-                       i = DBL_DIG - nd;
-                       if (e <= Ten_pmax + i) {
-                               /* A fancier test would sometimes let us do
-                                * this for larger i values.
-                                */
-                               e -= i;
-                               rv *= tens[i];
-                               /* rv = */ rounded_product(rv, tens[e]);
-                               goto ret;
-                       }
-               }
-#ifndef Inaccurate_Divide
-               else if (e >= -Ten_pmax) {
-                       /* rv = */ rounded_quotient(rv, tens[-e]);
-                       goto ret;
-               }
-#endif
-       }
-       e1 += nd - k;
-
-       /* Get starting approximation = rv * 10**e1 */
-
-       if (e1 > 0) {
-               if ( (i = e1 & 15) )
-                       rv *= tens[i];
-               if ( (e1 &= ~15) ) {
-                       if (e1 > DBL_MAX_10_EXP) {
-ovfl:
-                               errno = ERANGE;
-                               rv = HUGE_VAL;
-                               goto ret;
-                       }
-                       if (e1 >>= 4) {
-                               for (j = 0; e1 > 1; j++, e1 >>= 1)
-                                       if (e1 & 1)
-                                               rv *= bigtens[j];
-                       /* The last multiplication could overflow. */
-                               word0(rv) -= P*Exp_msk1;
-                               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;
-                                       }
-                               else
-                                       word0(rv) += P*Exp_msk1;
-                       }
-               }
-       } else if (e1 < 0) {
-               e1 = -e1;
-               if ( (i = e1 & 15) )
-                       rv /= tens[i];
-               if ( (e1 &= ~15) ) {
-                       e1 >>= 4;
-                       for (j = 0; e1 > 1; j++, e1 >>= 1)
-                               if (e1 & 1)
-                                       rv *= tinytens[j];
-                       /* The last multiplication could underflow. */
-                       rv0 = rv;
-                       rv *= tinytens[j];
-                       if (!rv) {
-                               rv = 2.*rv0;
-                               rv *= tinytens[j];
-                               if (!rv) {
-undfl:
-                                       rv = 0.;
-                                       errno = ERANGE;
-                                       goto ret;
-                                       }
-                               word0(rv) = Tiny0;
-                               word1(rv) = Tiny1;
-                               /* The refinement below will clean
-                                * this approximation up.
-                                */
-                       }
-               }
-       }
-
-       /* Now the hard part -- adjusting rv to the correct value.*/
-
-       /* Put digits into bd: true value = bd * 10^e */
-
-       bd0 = s2b(s0, nd0, nd, y);
-
-       for (;;) {
-               bd = Balloc(bd0->k);
-               Bcopy(bd, bd0);
-               bb = d2b(rv, &bbe, &bbbits);    /* rv = bb * 2^bbe */
-               bs = i2b(1);
-
-               if (e >= 0) {
-                       bb2 = bb5 = 0;
-                       bd2 = bd5 = e;
-               } else {
-                       bb2 = bb5 = -e;
-                       bd2 = bd5 = 0;
-               }
-               if (bbe >= 0)
-                       bb2 += bbe;
-               else
-                       bd2 -= bbe;
-               bs2 = bb2;
-#ifdef Sudden_Underflow
-               j = P + 1 - bbbits;
-#else
-               i = bbe + bbbits - 1;   /* logb(rv) */
-               if (i < Emin)   /* denormal */
-                       j = bbe + (P-Emin);
-               else
-                       j = P + 1 - bbbits;
-#endif
-               bb2 += j;
-               bd2 += j;
-               i = bb2 < bd2 ? bb2 : bd2;
-               if (i > bs2)
-                       i = bs2;
-               if (i > 0) {
-                       bb2 -= i;
-                       bd2 -= i;
-                       bs2 -= i;
-                       }
-               if (bb5 > 0) {
-                       bs = pow5mult(bs, bb5);
-                       bb1 = mult(bs, bb);
-                       Bfree(bb);
-                       bb = bb1;
-                       }
-               if (bb2 > 0)
-                       bb = lshift(bb, bb2);
-               if (bd5 > 0)
-                       bd = pow5mult(bd, bd5);
-               if (bd2 > 0)
-                       bd = lshift(bd, bd2);
-               if (bs2 > 0)
-                       bs = lshift(bs, bs2);
-               delta = diff(bb, bd);
-               dsign = delta->sign;
-               delta->sign = 0;
-               i = cmp(delta, bs);
-               if (i < 0) {
-                       /* 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)
-                               break;
-                       delta = lshift(delta,Log2P);
-                       if (cmp(delta, bs) > 0)
-                               goto drop_down;
-                       break;
-               }
-               if (i == 0) {
-                       /* exactly half-way between */
-                       if (dsign) {
-                               if ((word0(rv) & Bndry_mask1) == Bndry_mask1
-                                &&  word1(rv) == 0xffffffff) {
-                                       /*boundary case -- increment exponent*/
-                                       word0(rv) = (word0(rv) & Exp_mask)
-                                               + Exp_msk1;
-                                       word1(rv) = 0;
-                                       break;
-                               }
-                       } else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
-drop_down:
-                               /* boundary case -- decrement exponent */
-#ifdef Sudden_Underflow
-                               L = word0(rv) & Exp_mask;
-                               if (L <= Exp_msk1)
-                                       goto undfl;
-                               L -= Exp_msk1;
-#else
-                               L = (word0(rv) & Exp_mask) - Exp_msk1;
-#endif
-                               word0(rv) = L | Bndry_mask1;
-                               word1(rv) = 0xffffffff;
-                               break;
-                       }
-#ifndef ROUND_BIASED
-                       if (!(word1(rv) & LSB))
-                               break;
-#endif
-                       if (dsign)
-                               rv += ulp(rv);
-#ifndef ROUND_BIASED
-                       else {
-                               rv -= ulp(rv);
-#ifndef Sudden_Underflow
-                               if (!rv)
-                                       goto undfl;
-#endif
-                       }
-#endif
-                       break;
-               }
-               if ((aadj = ratio(delta, bs)) <= 2.) {
-                       if (dsign)
-                               aadj = aadj1 = 1.;
-                       else if (word1(rv) || word0(rv) & Bndry_mask) {
-#ifndef Sudden_Underflow
-                               if (word1(rv) == Tiny1 && !word0(rv))
-                                       goto undfl;
-#endif
-                               aadj = 1.;
-                               aadj1 = -1.;
-                       } else {
-                               /* special case -- power of FLT_RADIX to be */
-                               /* rounded down... */
-
-                               if (aadj < 2./FLT_RADIX)
-                                       aadj = 1./FLT_RADIX;
-                               else
-                                       aadj *= 0.5;
-                               aadj1 = -aadj;
-                       }
-               } else {
-                       aadj *= 0.5;
-                       aadj1 = dsign ? aadj : -aadj;
-#ifdef Check_FLT_ROUNDS
-                       switch(FLT_ROUNDS) {
-                               case 2: /* towards +infinity */
-                                       aadj1 -= 0.5;
-                                       break;
-                               case 0: /* towards 0 */
-                               case 3: /* towards -infinity */
-                                       aadj1 += 0.5;
-                       }
-#else
-                       if (FLT_ROUNDS == 0)
-                               aadj1 += 0.5;
-#endif
-               }
-               y = word0(rv) & Exp_mask;
-
-               /* Check for overflow */
-
-               if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
-                       rv0 = rv;
-                       word0(rv) -= P*Exp_msk1;
-                       adj = aadj1 * ulp(rv);
-                       rv += adj;
-                       if ((word0(rv) & Exp_mask) >=
-                                       Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
-                               if (word0(rv0) == Big0 && word1(rv0) == Big1)
-                                       goto ovfl;
-                               word0(rv) = Big0;
-                               word1(rv) = Big1;
-                               goto cont;
-                       } else
-                               word0(rv) += P*Exp_msk1;
-               } else {
-#ifdef Sudden_Underflow
-                       if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
-                               rv0 = rv;
-                               word0(rv) += P*Exp_msk1;
-                               adj = aadj1 * ulp(rv);
-                               rv += adj;
-                               if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
-                               {
-                                       if (word0(rv0) == Tiny0
-                                        && word1(rv0) == Tiny1)
-                                               goto undfl;
-                                       word0(rv) = Tiny0;
-                                       word1(rv) = Tiny1;
-                                       goto cont;
-                               } else
-                                       word0(rv) -= P*Exp_msk1;
-                       } else {
-                               adj = aadj1 * ulp(rv);
-                               rv += adj;
-                       }
-#else
-                       /* 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.,
-                        * 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);
-                               if (!dsign)
-                                       aadj1 = -aadj1;
-                       }
-                       adj = aadj1 * ulp(rv);
-                       rv += adj;
-#endif
-               }
-               z = word0(rv) & Exp_mask;
-               if (y == z) {
-                       /* Can we stop now? */
-                       L = aadj;
-                       aadj -= L;
-                       /* The tolerances below are conservative. */
-                       if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
-                               if (aadj < .4999999 || aadj > .5000001)
-                                       break;
-                       } else if (aadj < .4999999/FLT_RADIX)
-                               break;
-               }
-cont:
-               Bfree(bb);
-               Bfree(bd);
-               Bfree(bs);
-               Bfree(delta);
-       }
-       Bfree(bb);
-       Bfree(bd);
-       Bfree(bs);
-       Bfree(bd0);
-       Bfree(delta);
-ret:
-       if (se)
-               *se = (char *)s;
-       return sign ? -rv : rv;
-}
-
-static int
-quorem(Bigint *b, Bigint *S)
-{
-       int n;
-       long borrow, y;
-       unsigned long carry, q, ys;
-       unsigned long *bx, *bxe, *sx, *sxe;
-#ifdef Pack_32
-       long z;
-       unsigned long si, zs;
-#endif
-
-       n = S->wds;
-#ifdef DEBUG
-       /*debug*/ if (b->wds > n)
-       /*debug*/       Bug("oversize b in quorem");
-#endif
-       if (b->wds < n)
-               return 0;
-       sx = S->x;
-       sxe = sx + --n;
-       bx = b->x;
-       bxe = bx + n;
-       q = *bxe / (*sxe + 1);  /* ensure q <= true quotient */
-#ifdef DEBUG
-       /*debug*/ if (q > 9)
-       /*debug*/       Bug("oversized quotient in quorem");
-#endif
-       if (q) {
-               borrow = 0;
-               carry = 0;
-               do {
-#ifdef Pack_32
-                       si = *sx++;
-                       ys = (si & 0xffff) * q + carry;
-                       zs = (si >> 16) * q + (ys >> 16);
-                       carry = zs >> 16;
-                       y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
-                       borrow = y >> 16;
-                       Sign_Extend(borrow, y);
-                       z = (*bx >> 16) - (zs & 0xffff) + borrow;
-                       borrow = z >> 16;
-                       Sign_Extend(borrow, z);
-                       Storeinc(bx, z, y);
-#else
-                       ys = *sx++ * q + carry;
-                       carry = ys >> 16;
-                       y = *bx - (ys & 0xffff) + borrow;
-                       borrow = y >> 16;
-                       Sign_Extend(borrow, y);
-                       *bx++ = y & 0xffff;
-#endif
-               } while (sx <= sxe);
-               if (!*bxe) {
-                       bx = b->x;
-                       while (--bxe > bx && !*bxe)
-                               --n;
-                       b->wds = n;
-               }
-       }
-       if (cmp(b, S) >= 0) {
-               q++;
-               borrow = 0;
-               carry = 0;
-               bx = b->x;
-               sx = S->x;
-               do {
-#ifdef Pack_32
-                       si = *sx++;
-                       ys = (si & 0xffff) + carry;
-                       zs = (si >> 16) + (ys >> 16);
-                       carry = zs >> 16;
-                       y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
-                       borrow = y >> 16;
-                       Sign_Extend(borrow, y);
-                       z = (*bx >> 16) - (zs & 0xffff) + borrow;
-                       borrow = z >> 16;
-                       Sign_Extend(borrow, z);
-                       Storeinc(bx, z, y);
-#else
-                       ys = *sx++ + carry;
-                       carry = ys >> 16;
-                       y = *bx - (ys & 0xffff) + borrow;
-                       borrow = y >> 16;
-                       Sign_Extend(borrow, y);
-                       *bx++ = y & 0xffff;
-#endif
-               } while (sx <= sxe);
-               bx = b->x;
-               bxe = bx + n;
-               if (!*bxe) {
-                       while (--bxe > bx && !*bxe)
-                               --n;
-                       b->wds = n;
-               }
-       }
-       return q;
-}
-
-/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
- *
- * Inspired by "How to Print Floating-Point Numbers Accurately" by
- * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
- *
- * Modifications:
- *     1. Rather than iterating, we use a simple numeric overestimate
- *        to determine k = floor(log10(d)).  We scale relevant
- *        quantities using O(log2(k)) rather than O(k) multiplications.
- *     2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
- *        try to generate digits strictly left to right.  Instead, we
- *        compute with fewer bits and propagate the carry if necessary
- *        when rounding the final digit up.  This is often faster.
- *     3. Under the assumption that input will be rounded nearest,
- *        mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
- *        That is, we allow equality in stopping tests when the
- *        round-nearest rule will give the same floating-point value
- *        as would satisfaction of the stopping test with strict
- *        inequality.
- *     4. We remove common factors of powers of 2 from relevant
- *        quantities.
- *     5. When converting floating-point integers less than 1e16,
- *        we use floating-point arithmetic rather than resorting
- *        to multiple-precision integers.
- *     6. When asked to produce fewer than 15 digits, we first try
- *        to get by with floating-point arithmetic; we resort to
- *        multiple-precision integer arithmetic only if we cannot
- *        guarantee that the floating-point calculation has given
- *        the correctly rounded result.  For k requested digits and
- *        "uniformly" distributed input, the probability is
- *        something like 10^(k-15) that we must resort to the long
- *        calculation.
- */
-
-char *
-__dtoa(double d, int mode, int ndigits, int *decpt, int *sign, char **rve,
-       char **resultp)
-{
- /*    Arguments ndigits, decpt, sign are similar to those
-       of ecvt and fcvt; trailing zeros are suppressed from
-       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.
-
-       mode:
-               0 ==> shortest string that yields d when read in
-                       and rounded to nearest.
-               1 ==> like 0, but with Steele & White stopping rule;
-                       e.g. with IEEE P754 arithmetic , mode 0 gives
-                       1e23 whereas mode 1 gives 9.999999999999999e22.
-               2 ==> max(1,ndigits) significant digits.  This gives a
-                       return value similar to that of ecvt, except
-                       that trailing zeros are suppressed.
-               3 ==> through ndigits past the decimal point.  This
-                       gives a return value similar to that from fcvt,
-                       except that trailing zeros are suppressed, and
-                       ndigits can be negative.
-               4-9 should give the same return values as 2-3, i.e.,
-                       4 <= mode <= 9 ==> same return as mode
-                       2 + (mode & 1).  These modes are mainly for
-                       debugging; often they run slower but sometimes
-                       faster than modes 2-3.
-               4,5,8,9 ==> left-to-right digit generation.
-               6-9 ==> don't try fast floating-point estimate
-                       (if applicable).
-
-               Values of mode other than 0-9 are treated as mode 0.
-
-               Sufficient space is allocated to the return value
-               to hold the suppressed trailing zeros.
-       */
-
-       int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
-               j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
-               spec_case, try_quick;
-       long L;
-#ifndef Sudden_Underflow
-       int denorm;
-       unsigned long x;
-#endif
-       Bigint *b, *b1, *delta, *mlo, *mhi, *S;
-       double d2, ds, eps;
-       char *s, *s0;
-
-       /*
-        * XXX initialize to silence GCC warnings
-        */
-       ilim = 0;
-       ilim1 = 0;
-       mlo = NULL;
-       spec_case = 0;
-
-       if (word0(d) & Sign_bit) {
-               /* set sign for everything, including 0's and NaNs */
-               *sign = 1;
-               word0(d) &= ~Sign_bit;  /* clear sign bit */
-       }
-       else
-               *sign = 0;
-
-#ifdef IEEE_Arith
-       if ((word0(d) & Exp_mask) == Exp_mask)
-       {
-               /* Infinity or NaN */
-               *decpt = 9999;
-               s = !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" : "NaN";
-               if (rve)
-                       *rve = s[3] ? s + 8 : s + 3;
-               return s;
-       }
-#endif
-       if (!d) {
-               *decpt = 1;
-               s = "0";
-               if (rve)
-                       *rve = s + 1;
-               return s;
-       }
-
-       b = d2b(d, &be, &bbits);
-#ifdef Sudden_Underflow
-       i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
-#else
-       if ( (i = (int)((word0(d) >> Exp_shift1) & (Exp_mask>>Exp_shift1))) ) {
-#endif
-               d2 = d;
-               word0(d2) &= Frac_mask1;
-               word0(d2) |= Exp_11;
-
-               /* 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)
-                *
-                * This suggests computing an approximation k to log10(d) by
-                *
-                * k = (i - Bias)*0.301029995663981
-                *      + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
-                *
-                * We want k to be too large rather than too small.
-                * The error in the first-order Taylor series approximation
-                * is in our favor, so we just round up the constant enough
-                * to compensate for any error in the multiplication of
-                * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
-                * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
-                * adding 1e-13 to the constant term more than suffices.
-                * Hence we adjust the constant term to 0.1760912590558.
-                * (We could get a more accurate k by invoking log10,
-                *  but this is probably not worthwhile.)
-                */
-
-               i -= Bias;
-#ifndef Sudden_Underflow
-               denorm = 0;
-       } else {
-               /* 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));
-               d2 = x;
-               word0(d2) -= 31*Exp_msk1; /* adjust exponent */
-               i -= (Bias + (P-1) - 1) + 1;
-               denorm = 1;
-       }
-#endif
-       ds = (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 (d < tens[k])
-                       k--;
-               k_check = 0;
-       }
-       j = bbits - i - 1;
-       if (j >= 0) {
-               b2 = 0;
-               s2 = j;
-       } else {
-               b2 = -j;
-               s2 = 0;
-       }
-       if (k >= 0) {
-               b5 = 0;
-               s5 = k;
-               s2 += k;
-       } else {
-               b2 -= k;
-               b5 = -k;
-               s5 = 0;
-       }
-       if (mode < 0 || mode > 9)
-               mode = 0;
-       try_quick = 1;
-       if (mode > 5) {
-               mode -= 4;
-               try_quick = 0;
-       }
-       leftright = 1;
-       switch(mode) {
-               case 0:
-               case 1:
-                       ilim = ilim1 = -1;
-                       i = 18;
-                       ndigits = 0;
-                       break;
-               case 2:
-                       leftright = 0;
-                       /* no break */
-               case 4:
-                       if (ndigits <= 0)
-                               ndigits = 1;
-                       ilim = ilim1 = i = ndigits;
-                       break;
-               case 3:
-                       leftright = 0;
-                       /* no break */
-               case 5:
-                       i = ndigits + k + 1;
-                       ilim = i;
-                       ilim1 = i - 1;
-                       if (i <= 0)
-                               i = 1;
-       }
-       *resultp = (char *) malloc(i + 1);
-       s = s0 = *resultp;
-
-       if (ilim >= 0 && ilim <= Quick_max && try_quick) {
-
-               /* Try to get by with floating-point arithmetic. */
-
-               i = 0;
-               d2 = d;
-               k0 = k;
-               ilim0 = ilim;
-               ieps = 2; /* conservative */
-               if (k > 0) {
-                       ds = tens[k&0xf];
-                       j = k >> 4;
-                       if (j & Bletch) {
-                               /* prevent overflows */
-                               j &= Bletch - 1;
-                               d /= bigtens[n_bigtens-1];
-                               ieps++;
-                       }
-                       for (; j; j >>= 1, i++)
-                               if (j & 1) {
-                                       ieps++;
-                                       ds *= bigtens[i];
-                               }
-                       d /= ds;
-               } else if ( (j1 = -k) ) {
-                       d *= tens[j1 & 0xf];
-                       for (j = j1 >> 4; j; j >>= 1, i++)
-                               if (j & 1) {
-                                       ieps++;
-                                       d *= bigtens[i];
-                               }
-               }
-               if (k_check && d < 1. && ilim > 0) {
-                       if (ilim1 <= 0)
-                               goto fast_failed;
-                       ilim = ilim1;
-                       k--;
-                       d *= 10.;
-                       ieps++;
-               }
-               eps = ieps*d + 7.;
-               word0(eps) -= (P-1)*Exp_msk1;
-               if (ilim == 0) {
-                       S = mhi = 0;
-                       d -= 5.;
-                       if (d > eps)
-                               goto one_digit;
-                       if (d < -eps)
-                               goto no_digits;
-                       goto fast_failed;
-               }
-#ifndef No_leftright
-               if (leftright) {
-                       /* Use Steele & White method of only
-                        * generating digits needed.
-                        */
-                       eps = 0.5/tens[ilim-1] - eps;
-                       for (i = 0;;) {
-                               L = d;
-                               d -= L;
-                               *s++ = '0' + (int)L;
-                               if (d < eps)
-                                       goto ret1;
-                               if (1. - d < eps)
-                                       goto bump_up;
-                               if (++i >= ilim)
-                                       break;
-                               eps *= 10.;
-                               d *= 10.;
-                       }
-               } else {
-#endif
-                       /* Generate ilim digits, then fix them up. */
-                       eps *= tens[ilim-1];
-                       for (i = 1;; i++, d *= 10.) {
-                               L = d;
-                               d -= L;
-                               *s++ = '0' + (int)L;
-                               if (i == ilim) {
-                                       if (d > 0.5 + eps)
-                                               goto bump_up;
-                                       else if (d < 0.5 - eps) {
-                                               while (*--s == '0');
-                                               s++;
-                                               goto ret1;
-                                       }
-                                       break;
-                               }
-                       }
-#ifndef No_leftright
-               }
-#endif
-fast_failed:
-               s = s0;
-               d = d2;
-               k = k0;
-               ilim = ilim0;
-       }
-
-       /* Do we have a "small" integer? */
-
-       if (be >= 0 && k <= Int_max) {
-               /* Yes. */
-               ds = tens[k];
-               if (ndigits < 0 && ilim <= 0) {
-                       S = mhi = 0;
-                       if (ilim < 0 || d <= 5*ds)
-                               goto no_digits;
-                       goto one_digit;
-               }
-               for (i = 1;; i++) {
-                       L = d / ds;
-                       d -= L*ds;
-#ifdef Check_FLT_ROUNDS
-                       /* If FLT_ROUNDS == 2, L will usually be high by 1 */
-                       if (d < 0) {
-                               L--;
-                               d += ds;
-                       }
-#endif
-                       *s++ = '0' + (int)L;
-                       if (i == ilim) {
-                               d += d;
-                               if (d > ds || (d == ds && L & 1)) {
-bump_up:
-                                       while (*--s == '9')
-                                               if (s == s0) {
-                                                       k++;
-                                                       *s = '0';
-                                                       break;
-                                               }
-                                       ++*s++;
-                               }
-                               break;
-                       }
-                       if (!(d *= 10.))
-                               break;
-               }
-               goto ret1;
-       }
-
-       m2 = b2;
-       m5 = b5;
-       mhi = mlo = 0;
-       if (leftright) {
-               if (mode < 2) {
-                       i =
-#ifndef Sudden_Underflow
-                               denorm ? be + (Bias + (P-1) - 1 + 1) :
-#endif
-                               1 + P - bbits;
-               } else {
-                       j = ilim - 1;
-                       if (m5 >= j)
-                               m5 -= j;
-                       else {
-                               s5 += j -= m5;
-                               b5 += j;
-                               m5 = 0;
-                       }
-                       if ((i = ilim) < 0) {
-                               m2 -= i;
-                               i = 0;
-                       }
-               }
-               b2 += i;
-               s2 += i;
-               mhi = i2b(1);
-       }
-       if (m2 > 0 && s2 > 0) {
-               i = m2 < s2 ? m2 : s2;
-               b2 -= i;
-               m2 -= i;
-               s2 -= i;
-       }
-       if (b5 > 0) {
-               if (leftright) {
-                       if (m5 > 0) {
-                               mhi = pow5mult(mhi, m5);
-                               b1 = mult(mhi, b);
-                               Bfree(b);
-                               b = b1;
-                               }
-                       if ( (j = b5 - m5) )
-                               b = pow5mult(b, j);
-               } else
-                       b = pow5mult(b, b5);
-       }
-       S = i2b(1);
-       if (s5 > 0)
-               S = pow5mult(S, s5);
-
-       /* Check for special case that d is a normalized power of 2. */
-
-       if (mode < 2) {
-               if (!word1(d) && !(word0(d) & Bndry_mask)
-#ifndef Sudden_Underflow
-                && word0(d) & Exp_mask
-#endif
-                               ) {
-                       /* The special case */
-                       b2 += Log2P;
-                       s2 += Log2P;
-                       spec_case = 1;
-               } else
-                       spec_case = 0;
-       }
-
-       /* Arrange for convenient computation of quotients:
-        * shift left if necessary so divisor has 4 leading 0 bits.
-        *
-        * Perhaps we should just compute leading 28 bits of S once
-        * 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) )
-               i = 32 - i;
-#else
-       if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) )
-               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)
-               b = lshift(b, b2);
-       if (s2 > 0)
-               S = lshift(S, s2);
-       if (k_check) {
-               if (cmp(b,S) < 0) {
-                       k--;
-                       b = multadd(b, 10, 0);  /* we botched the k estimate */
-                       if (leftright)
-                               mhi = multadd(mhi, 10, 0);
-                       ilim = ilim1;
-               }
-       }
-       if (ilim <= 0 && mode > 2) {
-               if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
-                       /* no digits, fcvt style */
-no_digits:
-                       k = -1 - ndigits;
-                       goto ret;
-               }
-one_digit:
-               *s++ = '1';
-               k++;
-               goto ret;
-       }
-       if (leftright) {
-               if (m2 > 0)
-                       mhi = lshift(mhi, m2);
-
-               /* Compute mlo -- check for special case
-                * that d is a normalized power of 2.
-                */
-
-               mlo = mhi;
-               if (spec_case) {
-                       mhi = Balloc(mhi->k);
-                       Bcopy(mhi, mlo);
-                       mhi = lshift(mhi, Log2P);
-               }
-
-               for (i = 1;;i++) {
-                       dig = quorem(b,S) + '0';
-                       /* Do we yet have the shortest decimal string
-                        * that will round to d?
-                        */
-                       j = cmp(b, mlo);
-                       delta = diff(S, mhi);
-                       j1 = delta->sign ? 1 : cmp(b, delta);
-                       Bfree(delta);
-#ifndef ROUND_BIASED
-                       if (j1 == 0 && !mode && !(word1(d) & 1)) {
-                               if (dig == '9')
-                                       goto round_9_up;
-                               if (j > 0)
-                                       dig++;
-                               *s++ = dig;
-                               goto ret;
-                       }
-#endif
-                       if (j < 0 || (j == 0 && !mode
-#ifndef ROUND_BIASED
-                                                       && !(word1(d) & 1)
-#endif
-                                       )) {
-                               if (j1 > 0) {
-                                       b = lshift(b, 1);
-                                       j1 = cmp(b, S);
-                                       if ((j1 > 0 || (j1 == 0 && dig & 1))
-                                       && dig++ == '9')
-                                               goto round_9_up;
-                               }
-                               *s++ = dig;
-                               goto ret;
-                       }
-                       if (j1 > 0) {
-                               if (dig == '9') { /* possible if i == 1 */
-round_9_up:
-                                       *s++ = '9';
-                                       goto roundoff;
-                               }
-                               *s++ = dig + 1;
-                               goto ret;
-                       }
-                       *s++ = dig;
-                       if (i == ilim)
-                               break;
-                       b = multadd(b, 10, 0);
-                       if (mlo == mhi)
-                               mlo = mhi = multadd(mhi, 10, 0);
-                       else {
-                               mlo = multadd(mlo, 10, 0);
-                               mhi = multadd(mhi, 10, 0);
-                       }
-               }
-       } else
-               for (i = 1;; i++) {
-                       *s++ = dig = quorem(b,S) + '0';
-                       if (i >= ilim)
-                               break;
-                       b = multadd(b, 10, 0);
-               }
-
-       /* Round off last digit */
-
-       b = lshift(b, 1);
-       j = cmp(b, S);
-       if (j > 0 || (j == 0 && dig & 1)) {
-roundoff:
-               while (*--s == '9')
-                       if (s == s0) {
-                               k++;
-                               *s++ = '1';
-                               goto ret;
-                       }
-               ++*s++;
-       } else {
-               while (*--s == '0');
-               s++;
-       }
-ret:
-       Bfree(S);
-       if (mhi) {
-               if (mlo && mlo != mhi)
-                       Bfree(mlo);
-               Bfree(mhi);
-       }
-ret1:
-       Bfree(b);
-       if (s == s0) {  /* don't return empty string */
-               *s++ = '0';
-               k = 0;
-       }
-       *s = 0;
-       *decpt = k + 1;
-       if (rve)
-               *rve = s;
-       return s0;
-       }
-#ifdef __cplusplus
-}
-#endif
index 6584c5c..9922735 100644 (file)
@@ -7,6 +7,7 @@ NOPROFILE=
 CFLAGS+=-I${.CURDIR}/../libc/include
 CFLAGS+=-I${.CURDIR}/../libc/rpc
 CFLAGS+=-I${.CURDIR}/../../include -D__thread=
+CFLAGS+=${AINC}
 
 AINC=  -I${.OBJDIR} -I${.CURDIR}/../libc/${MACHINE_ARCH}
 PRECIOUSLIB=   yes