summaryrefslogtreecommitdiffstats
path: root/lib/libc
diff options
context:
space:
mode:
Diffstat (limited to 'lib/libc')
-rw-r--r--lib/libc/gdtoa/Makefile.inc2
-rw-r--r--lib/libc/gdtoa/_hdtoa.c235
-rw-r--r--lib/libc/gdtoa/_hldtoa.c176
3 files changed, 206 insertions, 207 deletions
diff --git a/lib/libc/gdtoa/Makefile.inc b/lib/libc/gdtoa/Makefile.inc
index dce7dc3..0d6832f 100644
--- a/lib/libc/gdtoa/Makefile.inc
+++ b/lib/libc/gdtoa/Makefile.inc
@@ -3,7 +3,7 @@
# netlib gdtoa sources
.PATH: ${.CURDIR}/gdtoa
-MISRCS+=_hdtoa.c _ldtoa.c glue.c
+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
diff --git a/lib/libc/gdtoa/_hdtoa.c b/lib/libc/gdtoa/_hdtoa.c
index 5e8a1835..4ac6cb1 100644
--- a/lib/libc/gdtoa/_hdtoa.c
+++ b/lib/libc/gdtoa/_hdtoa.c
@@ -1,5 +1,5 @@
/*-
- * Copyright (c) 2004, 2005 David Schultz <das@FreeBSD.ORG>
+ * Copyright (c) 2004-2008 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@@ -37,63 +37,10 @@ __FBSDID("$FreeBSD$");
#define INFSTR "Infinity"
#define NANSTR "NaN"
-#define DBL_ADJ (DBL_MAX_EXP - 2 + ((DBL_MANT_DIG - 1) % 4))
-#define LDBL_ADJ (LDBL_MAX_EXP - 2 + ((LDBL_MANT_DIG - 1) % 4))
+#define DBL_ADJ (DBL_MAX_EXP - 2)
+#define SIGFIGS ((DBL_MANT_DIG + 3) / 4 + 1)
-/*
- * Round up the given digit string. If the digit string is fff...f,
- * this procedure sets it to 100...0 and returns 1 to indicate that
- * the exponent needs to be bumped. Otherwise, 0 is returned.
- */
-static int
-roundup(char *s0, int ndigits)
-{
- char *s;
-
- for (s = s0 + ndigits - 1; *s == 0xf; s--) {
- if (s == s0) {
- *s = 1;
- return (1);
- }
- *s = 0;
- }
- ++*s;
- return (0);
-}
-
-/*
- * Round the given digit string to ndigits digits according to the
- * current rounding mode. Note that this could produce a string whose
- * value is not representable in the corresponding floating-point
- * type. The exponent pointed to by decpt is adjusted if necessary.
- */
-static void
-dorounding(char *s0, int ndigits, int sign, int *decpt)
-{
- int adjust = 0; /* do we need to adjust the exponent? */
-
- switch (FLT_ROUNDS) {
- case 0: /* toward zero */
- default: /* implementation-defined */
- break;
- case 1: /* to nearest, halfway rounds to even */
- if ((s0[ndigits] > 8) ||
- (s0[ndigits] == 8 && s0[ndigits + 1] & 1))
- adjust = roundup(s0, ndigits);
- break;
- case 2: /* toward +inf */
- if (sign == 0)
- adjust = roundup(s0, ndigits);
- break;
- case 3: /* toward -inf */
- if (sign != 0)
- adjust = roundup(s0, ndigits);
- break;
- }
-
- if (adjust)
- *decpt += 4;
-}
+static const float one[] = { 1.0f, -1.0f };
/*
* This procedure converts a double-precision number in IEEE format
@@ -112,9 +59,9 @@ dorounding(char *s0, int ndigits, int sign, int *decpt)
*
* 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 chooses the
- * first digit so that subsequent digits are aligned on nibble
- * boundaries (before rounding).
+ * 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
@@ -123,10 +70,10 @@ char *
__hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign,
char **rve)
{
- static const int sigfigs = (DBL_MANT_DIG + 3) / 4;
union IEEEd2bits u;
char *s, *s0;
int bufsize;
+ uint32_t manh, manl;
u.d = d;
*sign = u.bits.sign;
@@ -145,11 +92,9 @@ __hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign,
case FP_INFINITE:
*decpt = INT_MAX;
return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
- case FP_NAN:
+ default: /* FP_NAN or unrecognized */
*decpt = INT_MAX;
return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
- default:
- abort();
}
/* FP_NORMAL or FP_SUBNORMAL */
@@ -158,162 +103,40 @@ __hdtoa(double d, const char *xdigs, int ndigits, int *decpt, int *sign,
ndigits = 1;
/*
- * For simplicity, we generate all the digits even if the
- * caller has requested fewer.
+ * If ndigits < 0, we are expected to auto-size, so we allocate
+ * enough space for all the digits.
*/
- bufsize = (sigfigs > ndigits) ? sigfigs : ndigits;
+ bufsize = (ndigits > 0) ? ndigits : SIGFIGS;
s0 = rv_alloc(bufsize);
- /*
- * We work from right to left, first adding any requested zero
- * padding, then the least significant portion of the
- * mantissa, followed by the most significant. The buffer is
- * filled with the byte values 0x0 through 0xf, which are
- * converted to xdigs[0x0] through xdigs[0xf] after the
- * rounding phase.
- */
- for (s = s0 + bufsize - 1; s > s0 + sigfigs - 1; s--)
- *s = 0;
- for (; s > s0 + sigfigs - (DBL_MANL_SIZE / 4) - 1 && s > s0; s--) {
- *s = u.bits.manl & 0xf;
- u.bits.manl >>= 4;
+ /* 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;
}
- for (; s > s0; s--) {
- *s = u.bits.manh & 0xf;
- u.bits.manh >>= 4;
- }
-
- /*
- * At this point, we have snarfed all the bits in the
- * mantissa, with the possible exception of the highest-order
- * (partial) nibble, which is dealt with by the next
- * statement. We also tack on the implicit normalization bit.
- */
- *s = u.bits.manh | (1U << ((DBL_MANT_DIG - 1) % 4));
- /* If ndigits < 0, we are expected to auto-size the precision. */
- if (ndigits < 0) {
- for (ndigits = sigfigs; s0[ndigits - 1] == 0; ndigits--)
- ;
+ 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 (sigfigs > ndigits && s0[ndigits] != 0)
- dorounding(s0, ndigits, u.bits.sign, decpt);
-
- s = s0 + ndigits;
- if (rve != NULL)
- *rve = s;
- *s-- = '\0';
- for (; s >= s0; s--)
- *s = xdigs[(unsigned int)*s];
-
- return (s0);
-}
-
-#if (LDBL_MANT_DIG > DBL_MANT_DIG)
-
-/*
- * This is the long double version of __hdtoa().
- */
-char *
-__hldtoa(long double e, const char *xdigs, int ndigits, int *decpt, int *sign,
- char **rve)
-{
- static const int sigfigs = (LDBL_MANT_DIG + 3) / 4;
- union IEEEl2bits u;
- char *s, *s0;
- int bufsize;
-
- 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:
- u.e *= 0x1p514L;
- *decpt = u.bits.exp - (514 + LDBL_ADJ);
- break;
- case FP_INFINITE:
- *decpt = INT_MAX;
- return (nrv_alloc(INFSTR, rve, sizeof(INFSTR) - 1));
- case FP_NAN:
- *decpt = INT_MAX;
- return (nrv_alloc(NANSTR, rve, sizeof(NANSTR) - 1));
- default:
- abort();
- }
-
- /* FP_NORMAL or FP_SUBNORMAL */
-
- if (ndigits == 0) /* dtoa() compatibility */
- ndigits = 1;
-
- /*
- * For simplicity, we generate all the digits even if the
- * caller has requested fewer.
- */
- bufsize = (sigfigs > ndigits) ? sigfigs : ndigits;
- s0 = rv_alloc(bufsize);
-
- /*
- * We work from right to left, first adding any requested zero
- * padding, then the least significant portion of the
- * mantissa, followed by the most significant. The buffer is
- * filled with the byte values 0x0 through 0xf, which are
- * converted to xdigs[0x0] through xdigs[0xf] after the
- * rounding phase.
- */
- for (s = s0 + bufsize - 1; s > s0 + sigfigs - 1; s--)
- *s = 0;
- for (; s > s0 + sigfigs - (LDBL_MANL_SIZE / 4) - 1 && s > s0; s--) {
- *s = u.bits.manl & 0xf;
- u.bits.manl >>= 4;
- }
- for (; s > s0; s--) {
- *s = u.bits.manh & 0xf;
- u.bits.manh >>= 4;
- }
-
- /*
- * At this point, we have snarfed all the bits in the
- * mantissa, with the possible exception of the highest-order
- * (partial) nibble, which is dealt with by the next
- * statement. We also tack on the implicit normalization bit.
- */
- *s = u.bits.manh | (1U << ((LDBL_MANT_DIG - 1) % 4));
-
/* If ndigits < 0, we are expected to auto-size the precision. */
if (ndigits < 0) {
- for (ndigits = sigfigs; s0[ndigits - 1] == 0; ndigits--)
+ for (ndigits = SIGFIGS; s0[ndigits - 1] == '0'; ndigits--)
;
}
- if (sigfigs > ndigits && s0[ndigits] != 0)
- dorounding(s0, ndigits, u.bits.sign, decpt);
-
s = s0 + ndigits;
+ *s = '\0';
if (rve != NULL)
*rve = s;
- *s-- = '\0';
- for (; s >= s0; s--)
- *s = xdigs[(unsigned int)*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/_hldtoa.c b/lib/libc/gdtoa/_hldtoa.c
new file mode 100644
index 0000000..24e42a2
--- /dev/null
+++ b/lib/libc/gdtoa/_hldtoa.c
@@ -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.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <float.h>
+#include <limits.h>
+#include <math.h>
+#include <stdint.h>
+
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#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) */
OpenPOWER on IntegriCloud