gawk-diffs
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[gawk-diffs] [SCM] gawk branch, long-double, updated. a37db774455e596fa9


From: John Haque
Subject: [gawk-diffs] [SCM] gawk branch, long-double, updated. a37db774455e596fa9c964a2e0f0fe73af47e221
Date: Sun, 27 Jan 2013 10:57:21 +0000

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "gawk".

The branch, long-double has been updated
       via  a37db774455e596fa9c964a2e0f0fe73af47e221 (commit)
      from  26bfdfdc9bb065dfc1ea7f2494070f826e81b874 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
http://git.sv.gnu.org/cgit/gawk.git/commit/?id=a37db774455e596fa9c964a2e0f0fe73af47e221

commit a37db774455e596fa9c964a2e0f0fe73af47e221
Author: John Haque <address@hidden>
Date:   Sun Jan 27 04:53:25 2013 -0600

    Formatting and strtold emulation for __float128, organize misc dir.

diff --git a/ChangeLog b/ChangeLog
index d0daef3..4caad77 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,10 @@
+2013-01-27         John Haque      <address@hidden>
+
+       * long_double.c (GAWK_FORMAT_INT): New define.
+       (format_float_1, format_uint_1, double_to_int): Move definitions here
+       from long_double.h.
+       * long_double.h (format_float_1, format_uint_1, double_to_int): Adjust.
+
 2013-01-23         John Haque      <address@hidden>
 
        * long_double.c (format_uint_1): New inline routines to format integers.
diff --git a/Makefile.am b/Makefile.am
index 00d1b2c..5aec821 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -68,7 +68,10 @@ distcleancheck_listfiles = \
 # The order to do things in.
 # Build explicitly in "." in order to build gawk first, so
 # that `make check' without a prior `make' works.
-SUBDIRS = \
+
+EXTRADIRS = misc
+
+BASEDIRS = \
        . \
        awklib \
        doc \
@@ -76,6 +79,8 @@ SUBDIRS = \
        extension \
        test
 
+SUBDIRS = $(EXTRADIRS) $(BASEDIRS)
+
 # what to make and install
 bin_PROGRAMS = gawk
 include_HEADERS = gawkapi.h
@@ -132,10 +137,12 @@ base_sources = \
        version.c \
        xalloc.h
 
-gawk_SOURCES = $(base_sources) float128.c
+gawk_SOURCES = $(base_sources)
+
+EXTRA_LIBS = misc/libmisc.a
 
 # Get extra libs as needed, Automake will supply LIBINTL and SOCKET_LIBS.
-LDADD = $(LIBSIGSEGV) $(LIBINTL) $(SOCKET_LIBS) @LIBREADLINE@ @LIBMPFR@
+LDADD = $(EXTRA_LIBS) $(LIBSIGSEGV) $(LIBINTL) $(SOCKET_LIBS) @LIBREADLINE@ 
@LIBMPFR@
 
 # Directory for gawk's data files. Automake supplies datadir.
 pkgdatadir = $(datadir)/awk
diff --git a/Makefile.in b/Makefile.in
index 243af73..007f0fa 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -115,12 +115,12 @@ am__objects_1 = array.$(OBJEXT) awkgram.$(OBJEXT) 
builtin.$(OBJEXT) \
        random.$(OBJEXT) re.$(OBJEXT) regex.$(OBJEXT) \
        replace.$(OBJEXT) str_array.$(OBJEXT) symbol.$(OBJEXT) \
        version.$(OBJEXT)
-am_gawk_OBJECTS = $(am__objects_1) float128.$(OBJEXT)
+am_gawk_OBJECTS = $(am__objects_1)
 gawk_OBJECTS = $(am_gawk_OBJECTS)
 gawk_LDADD = $(LDADD)
 am__DEPENDENCIES_1 =
-gawk_DEPENDENCIES = $(am__DEPENDENCIES_1) $(am__DEPENDENCIES_1) \
-       $(am__DEPENDENCIES_1)
+gawk_DEPENDENCIES = $(EXTRA_LIBS) $(am__DEPENDENCIES_1) \
+       $(am__DEPENDENCIES_1) $(am__DEPENDENCIES_1)
 DEFAULT_INCLUDES = address@hidden@
 depcomp = $(SHELL) $(top_srcdir)/depcomp
 am__depfiles_maybe = depfiles
@@ -403,7 +403,8 @@ distcleancheck_listfiles = \
 # The order to do things in.
 # Build explicitly in "." in order to build gawk first, so
 # that `make check' without a prior `make' works.
-SUBDIRS = \
+EXTRADIRS = misc
+BASEDIRS = \
        . \
        awklib \
        doc \
@@ -411,6 +412,7 @@ SUBDIRS = \
        extension \
        test
 
+SUBDIRS = $(EXTRADIRS) $(BASEDIRS)
 include_HEADERS = gawkapi.h
 
 # sources for both gawk and dgawk
@@ -465,10 +467,11 @@ base_sources = \
        version.c \
        xalloc.h
 
-gawk_SOURCES = $(base_sources) float128.c
+gawk_SOURCES = $(base_sources)
+EXTRA_LIBS = misc/libmisc.a
 
 # Get extra libs as needed, Automake will supply LIBINTL and SOCKET_LIBS.
-LDADD = $(LIBSIGSEGV) $(LIBINTL) $(SOCKET_LIBS) @LIBREADLINE@ @LIBMPFR@
+LDADD = $(EXTRA_LIBS) $(LIBSIGSEGV) $(LIBINTL) $(SOCKET_LIBS) @LIBREADLINE@ 
@LIBMPFR@
 
 # stuff for compiling gawk/pgawk
 DEFPATH = '".$(PATH_SEPARATOR)$(pkgdatadir)"'
@@ -598,7 +601,6 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @address@hidden/$(DEPDIR)/address@hidden@
 @AMDEP_TRUE@@am__include@ @address@hidden/$(DEPDIR)/address@hidden@
 @AMDEP_TRUE@@am__include@ @address@hidden/$(DEPDIR)/address@hidden@
address@hidden@@am__include@ @address@hidden/$(DEPDIR)/address@hidden@
 @AMDEP_TRUE@@am__include@ @address@hidden/$(DEPDIR)/address@hidden@
 @AMDEP_TRUE@@am__include@ @address@hidden/$(DEPDIR)/address@hidden@
 @AMDEP_TRUE@@am__include@ @address@hidden/$(DEPDIR)/address@hidden@
diff --git a/configure b/configure
index 3d80b0d..e3b43ce 100755
--- a/configure
+++ b/configure
@@ -5558,8 +5558,6 @@ $as_echo "no" >&6; }
        if test -f $srcdir/.dev
        then
                CFLAGS="$CFLAGS -DLDBLTEST=1"   # DO NOT TURN OFF ASSERTIONS
-#              cp $srcdir/misc/float128.c $srcdir/
-#              sed -i 's/^gawk_SOURCES = $(base_sources)$/& float128.c/' 
$srcdir/Makefile.am
        else
                CFLAGS="$CFLAGS -DNDEBUG"       # turn off assertions
        fi
diff --git a/configure.ac b/configure.ac
index 2a13506..98e3c94 100644
--- a/configure.ac
+++ b/configure.ac
@@ -98,8 +98,6 @@ else
        if test -f $srcdir/.dev
        then
                CFLAGS="$CFLAGS -DLDBLTEST=1"   # DO NOT TURN OFF ASSERTIONS
-#              cp $srcdir/misc/float128.c $srcdir/
-#              sed -i 's/^gawk_SOURCES = $(base_sources)$/& float128.c/' 
$srcdir/Makefile.am
        else
                CFLAGS="$CFLAGS -DNDEBUG"       # turn off assertions
        fi
diff --git a/float128.c b/float128.c
deleted file mode 100644
index 15582b8..0000000
--- a/float128.c
+++ /dev/null
@@ -1,128 +0,0 @@
-#include "awk.h"
-
-#if defined(LDBLTEST) && LDBLTEST == 1
-#include <math.h>
-#include "random.h"
-#if 0
-#include "floatmagic.h"        /* definition of isnan -- XXX: only for float 
or double or long double */
-#endif
-
-#include "format.h"
-
-/* XXX: No libquadmath, and don't want or need it anyway. */
-
-/*
- * Macros copied from quadmath.h
- *     https://github.com/mirrors/gcc/blob/master/libquadmath/quadmath.h
- */
-
-#define FLT128_MAX 1.18973149535723176508575932662800702e4932Q
-#define FLT128_MIN 3.36210314311209350626267781732175260e-4932Q
-#define FLT128_EPSILON 1.92592994438723585305597794258492732e-34Q
-#define FLT128_DENORM_MIN 6.475175119438025110924438958227646552e-4966Q
-#define FLT128_MANT_DIG 113
-#define FLT128_MIN_EXP (-16381)
-#define FLT128_MAX_EXP 16384
-#define FLT128_DIG 33
-#define FLT128_MIN_10_EXP (-4931)
-#define FLT128_MAX_10_EXP 4932
-
-/* #define HUGE_VALQ __builtin_huge_valq() */
-/* The following alternative is valid, but brings the warning:
-   (floating constant exceeds range of ‘__float128’)  */
-#define HUGE_VALQ (__extension__ 0x1.0p32767Q)
-
-/* end of Macros from quadmath.h */
-
-
-#define LDBL_FRAC_BITS FLT128_MANT_DIG
-#define        LDBL_INT_BITS   128
-
-#define AWKLDBL        __float128
-#define LDBL_VAL(n)    (*((AWKLDBL *) (n)->qnumbr))
-#define LDC(x)         x##Q
-
-#define get_long_double(d)     emalloc(d, void *, sizeof(AWKLDBL), "float128")
-#define free_long_double(d)    efree(d)
-
-static int format_float_1(char *str, size_t size, const char *format, int fw, 
int prec, AWKLDBL x);
-static int format_uint_finite_p(char *str, size_t size, AWKLDBL x);
-
-/*
- * format_uint_1 --- format a AWKLDBL as an unsigned integer. The double value
- *     must be finite and >= 0.
- */
-
-static inline int
-format_uint_1(char *str, size_t size, AWKLDBL x)
-{
-       int ret;
-       if ((ret = format_uint_finite_p(str, size, x)) < 0)
-               return snprintf(str, size, "%.0f", (double) x);
-       return ret;
-}
-
-static AWKLDBL double_to_int(AWKLDBL);
-
-/* if have __float128, also have strtold. No ? */
-
-static inline AWKLDBL
-gawk_strtold(const char *str, char **endptr)
-{
-       return strtold(str, endptr);
-}
-
-/* Define isnan() and isinf() before including the two files */
-
-static inline int isnan_awkldbl(AWKLDBL x) { return x != x; }
-#ifdef isnan
-#undef isnan
-#endif
-#define isnan isnan_awkldbl
-
-static inline int isinf_awkldbl(AWKLDBL x) { return isnan(x - x); }
-#ifdef isinf
-#undef isinf
-#endif
-#define isinf isinf_awkldbl
-
-#define GAWK_INFINITY  HUGE_VALQ
-#define GAWK_NAN       (LDC(0.0) / LDC(0.0))
-
-
-/* XXX: Don't want floorl() or ceill() */ 
-
-#ifdef HAVE_FLOORL
-#undef HAVE_FLOORL
-#endif
-#ifdef HAVE_CEILL
-#undef HAVE_CEILL
-#endif
-#ifdef PRINTF_HAS_LF_FORMAT
-#undef PRINTF_HAS_LF_FORMAT
-#endif
-
-numbr_handler_t float128_hndlr;
-
-#define awkldbl_hndlr float128_hndlr
-
-#define gawk_int_t long long
-#define gawk_uint_t unsigned long long
-#ifdef SIZEOF_GAWK_INT
-#undef SIZEOF_GAWK_INT
-#undef GAWK_INT_MAX
-#undef GAWK_INT_MIN
-#undef GAWK_UINT_MAX
-#endif
-
-#define SIZEOF_GAWK_INT        8
-#define        GAWK_INT_MAX    LLONG_MAX
-#define        GAWK_INT_MIN    LLONG_MIN
-#define        GAWK_UINT_MAX   ULLONG_MAX
-
-#define TEST_NUMBR 1
-
-#include "misc/gawk_math.h"
-#include "long_double.h"
-
-#endif /* LDBLTEST == 1 */
diff --git a/long_double.c b/long_double.c
index 593cd97..c2a5c58 100644
--- a/long_double.c
+++ b/long_double.c
@@ -156,36 +156,53 @@ gawk_sqrtl(AWKLDBL x)
 #endif
 
 
-#if ! defined(PRINTF_HAS_LF_FORMAT)
-#ifdef HAVE_FLOORL
-#undef HAVE_FLOORL
-#endif
-#ifdef HAVE_CEILL
-#undef HAVE_CEILL
+/* N.B: if floorl() or ceill() or "%Lf" is missing format integers ourself */
+
+#if ! (defined(HAVE_FLOORL) && defined(HAVE_CEILL) && 
defined(PRINTF_HAS_LF_FORMAT)) 
+#define GAWK_FMT_INT 1
 #endif
-#else  /* PRINTF_HAS_LF_FORMAT */
+
+
+#if ! defined(PRINTF_HAS_LF_FORMAT)
 
 /*
- * XXX: have "%Lf" but no floorl() and/or ceill(). ceill() (or ceil() for
- * double) isn't really needed.
+ * format_float_1 --- format a single AWKLDBL value according to FORMAT.
+ *     The value must be finite.       
  */
 
-#if ! (defined(HAVE_FLOORL) && defined(HAVE_CEILL))
-#ifdef HAVE_FLOORL
-#undef HAVE_FLOORL
-#endif
-#ifdef HAVE_CEILL
-#undef HAVE_CEILL
+static int
+format_float_1(char *str, size_t size, const char *format, int fw, int prec, 
AWKLDBL x)
+{
+       char alt_format[16];
+       size_t len;
+
+       len = strlen(format);
+
+       /* expect %Lf, %LF, %Le, %LE, %Lg or %LG */
+       assert(len >= 2 && format[len - 2] == 'L');
+
+       memcpy(alt_format, format, len - 2);
+       alt_format[len - 2] = format[len - 1];  /* skip the `L' */ 
+       alt_format[len - 1] = '\0';
+       return snprintf(str, size, alt_format, fw, prec, (double) x);
+}
+#else
+
+static inline int
+format_float_1(char *str, size_t size, const char *format, int fw, int prec, 
AWKLDBL x)
+{
+       return snprintf(str, size, format, fw, prec, x);
+}
 #endif
-#endif /* ! (HAVE_FLOORL && HAVE_CEILL) */ 
-#endif /* PRINTF_HAS_LF_FORMAT */
 
-#if ! defined(PRINTF_HAS_LF_FORMAT)
-static int format_float_1(char *str, size_t size, const char *format, int fw, 
int prec, AWKLDBL x);
+
+#if    defined(GAWK_FMT_INT)
+
 static int format_uint_finite_p(char *str, size_t size, AWKLDBL x);
+static AWKLDBL gawk_floorl_finite_p(AWKLDBL x, gawk_uint_t *chunk);
 
 /*
- * format_uint_finite_p --- format a long double as an unsigned integer. The 
double value
+ * format_uint_1 --- format a long double as an unsigned integer. The double 
value
  *     must be finite and >= 0.
  */
 
@@ -193,33 +210,41 @@ static inline int
 format_uint_1(char *str, size_t size, AWKLDBL x)
 {
        int ret;
-       if ((ret = format_uint_finite_p(str, size, x)) < 0)
+       if ((ret = format_uint_finite_p(str, size, x)) < 0) {
+               /* outside useful range */
                return snprintf(str, size, "%.0f", (double) x);
+       }
        return ret;
 }
-#else
 
-/*
- * format_float_1 --- format a single AWKLDBL value according to FORMAT.
- *     The value must be finite.       
- */
+/* double_to_int --- convert double to int, used in several places */
 
-static inline int
-format_float_1(char *str, size_t size, const char *format, int fw, int prec, 
AWKLDBL x)
+static AWKLDBL
+double_to_int(AWKLDBL x)
 {
-       return snprintf(str, size, format, fw, prec, x);
+       AWKLDBL intval;
+       int sign = 1;
+
+       if (isnan(x) || isinf(x) || x == LDC(0.0))
+               return x;
+       if (x < LDC(0.0)) {
+               sign = -1;
+               x = -x;
+       }
+       if ((intval = gawk_floorl_finite_p(x, NULL)) < LDC(0.0)) {
+               /* outside range, use floor() for C double */
+               intval = (AWKLDBL) Floor((double) x);
+       }
+       return intval * ((AWKLDBL) sign);
 }
 
+#else
+
 static inline int
 format_uint_1(char *str, size_t size, AWKLDBL x)
 {
        return snprintf(str, size, "%.0Lf", x);
 }
-#endif
-
-#if ! (defined(HAVE_FLOORL) && defined(HAVE_CEILL))
-static AWKLDBL double_to_int(AWKLDBL);
-#else
 
 /* double_to_int --- convert double to int, used in several places */
 
@@ -230,6 +255,7 @@ double_to_int(AWKLDBL x)
 }
 #endif
 
+
 #include "long_double.h"
 
 #else  /* ! USE_LONG_DOUBLE */
diff --git a/long_double.h b/long_double.h
index 023ece1..203e04b 100644
--- a/long_double.h
+++ b/long_double.h
@@ -30,11 +30,15 @@
  */
 
 
-static AWKLDBL *pow2d_table;   /* 2^n table for 0 <= n <= LDBL_INT_BITS */
-static AWKLDBL init_pow2d_table(unsigned int);
-
-#define POW2LD(n)      (pow2d_table != NULL ? pow2d_table[n] : 
init_pow2d_table(n))
+static AWKLDBL *pow2d_table;   /* 2^n table for 0 <= n <= 128 */
+static void init_pow2d_table(void);
+static AWKLDBL pow2ld__p(int n);
 
+static inline AWKLDBL
+pow2ld(unsigned int n)
+{
+       return n <= 128 ? pow2d_table[n] : pow2ld__p(n);
+}
 
 /* Can declare these, since we always use the random shipped with gawk */
 extern char *initstate(unsigned long seed, char *state, long n);
@@ -177,6 +181,8 @@ awkldbl_init(bltin_t **numbr_bltins)
        false_node->flags |= NUMINT;
        true_node->flags |= NUMINT;
 
+       init_pow2d_table();     /* FIXME -- initialize only if needed ? */
+
        *numbr_bltins = awkldbl_bltins;
        return true;
 }
@@ -1609,19 +1615,28 @@ fmt1:
 
 /* init_pow2d_table --- populate powers of 2 table */
 
-static AWKLDBL
-init_pow2d_table(unsigned int n)
+static void
+init_pow2d_table()
 {
        unsigned i;
        emalloc(pow2d_table, AWKLDBL *, 129 * sizeof (AWKLDBL), 
"init_pow2d_table");
        pow2d_table[0] = LDC(1.0);
        for (i = 1; i <= 128; i++)
                pow2d_table[i] = pow2d_table[i - 1] * LDC(2.0);
-       return pow2d_table[n];
 }
 
+/* pow2ld__p --- return 2^n for some integer n >= 0 */
+
+static inline AWKLDBL
+pow2ld__p(int n)
+{
+       AWKLDBL dval = LDC(1.0);
+        for (; n > 128; n -= 128)
+               dval *= pow2d_table[128];
+        return dval * pow2d_table[n];
+}
 
-#if ! (defined(HAVE_FLOORL) && defined(HAVE_CEILL) && 
defined(PRINTF_HAS_LF_FORMAT))
+#if defined(GAWK_FMT_INT)
 
 /*
  * gawk_floorl_finite_p --- provide floor() for long double. The double value
@@ -1636,9 +1651,9 @@ gawk_floorl_finite_p(AWKLDBL x, gawk_uint_t *chunk)
        AWKLDBL intval = LDC(0.0);
 
 #if    LDBL_INT_BITS == 128
-       if (x >= POW2LD(113))
+       if (x >= pow2ld(113))
 #else
-       if (x >= POW2LD(64))
+       if (x >= pow2ld(64))
 #endif
                return -LDC(1.0);
 
@@ -1651,15 +1666,15 @@ gawk_floorl_finite_p(AWKLDBL x, gawk_uint_t *chunk)
                low = 0;
                while (low <= high) {
                        mid = (low + high) / 2;
-                       if (x < POW2LD(mid))
+                       if (x < pow2ld(mid))
                                high = mid - 1;
                        else
                                low = mid + 1;
                }
 
-               if (x < POW2LD(low)) {
-                       x -= POW2LD(low - 1);
-                       intval += POW2LD(low - 1);
+               if (x < pow2ld(low)) {
+                       x -= pow2ld(low - 1);
+                       intval += pow2ld(low - 1);
 
                        if (chunk != NULL) {
                                /*
@@ -1669,15 +1684,15 @@ gawk_floorl_finite_p(AWKLDBL x, gawk_uint_t *chunk)
                                 */
 
 #if    LDBL_INT_BITS == 128
-                               if (low <= 32) chunk[0] += (gawk_uint_t) 
POW2LD(low - 1);
-                               else if (low <= 64) chunk[1] += (gawk_uint_t) 
POW2LD(low - 33);
-                               else if (low <= 96) chunk[2] += (gawk_uint_t) 
POW2LD(low - 65);
-                               else chunk[3] += (gawk_uint_t) POW2LD(low - 97);
+                               if (low <= 32) chunk[0] += (gawk_uint_t) 
pow2ld(low - 1);
+                               else if (low <= 64) chunk[1] += (gawk_uint_t) 
pow2ld(low - 33);
+                               else if (low <= 96) chunk[2] += (gawk_uint_t) 
pow2ld(low - 65);
+                               else chunk[3] += (gawk_uint_t) pow2ld(low - 97);
 #else
-                               if (low <= 16) chunk[0] += (gawk_uint_t) 
POW2LD(low - 1);
-                               else if (low <= 32) chunk[1] += (gawk_uint_t) 
POW2LD(low - 17);
-                               else if (low <= 48) chunk[2] += (gawk_uint_t) 
POW2LD(low - 33);
-                               else chunk[3] += (gawk_uint_t) POW2LD(low - 49);
+                               if (low <= 16) chunk[0] += (gawk_uint_t) 
pow2ld(low - 1);
+                               else if (low <= 32) chunk[1] += (gawk_uint_t) 
pow2ld(low - 17);
+                               else if (low <= 48) chunk[2] += (gawk_uint_t) 
pow2ld(low - 33);
+                               else chunk[3] += (gawk_uint_t) pow2ld(low - 49);
 #endif
                        }
                }
@@ -1693,36 +1708,6 @@ gawk_floorl_finite_p(AWKLDBL x, gawk_uint_t *chunk)
        }
        return intval;  /* >= 0.0 */
 }
-#endif
-
-
-#if ! (defined(HAVE_FLOORL) && defined(HAVE_CEILL))
-
-/* double_to_int --- convert double to int, used in several places */
-
-static AWKLDBL
-double_to_int(AWKLDBL x)
-{
-       AWKLDBL intval;
-       int sign = 1;
-
-       if (isnan(x) || isinf(x) || x == LDC(0.0))
-               return x;
-
-       if (x < LDC(0.0)) {
-               sign = -1;
-               x = -x;
-       }
-
-       if ((intval = gawk_floorl_finite_p(x, NULL)) < LDC(0.0))
-               intval = (AWKLDBL) Floor((double) x);   /* outside range, use 
floor() for C double */
-       return intval * ((AWKLDBL) sign);
-}
-
-#endif
-
-
-#if ! defined(PRINTF_HAS_LF_FORMAT)
 
 /*
  * format_uint_finite_p --- format a AWKLDBL as an unsigned integer. The 
AWKLDBL value
@@ -1810,25 +1795,4 @@ format_uint_finite_p(char *str, size_t size, AWKLDBL x)
        return len;
 }
 
-/*
- * format_float_1 --- format a single AWKLDBL value according to FORMAT.
- *     The value must be finite.       
- */
-
-static int
-format_float_1(char *str, size_t size, const char *format, int fw, int prec, 
AWKLDBL x)
-{
-       char alt_format[16];
-       size_t len;
-
-       len = strlen(format);
-
-       /* expect %Lf, %LF, %Le, %LE, %Lg or %LG */
-       assert(len >= 2 && format[len - 2] == 'L');
-
-       memcpy(alt_format, format, len - 2);
-       alt_format[len - 2] = format[len - 1];  /* skip the `L' */ 
-       alt_format[len - 1] = '\0';
-       return snprintf(str, size, alt_format, fw, prec, (double) x);
-}
 #endif
diff --git a/misc/.gitignore b/misc/.gitignore
new file mode 100644
index 0000000..2460008
--- /dev/null
+++ b/misc/.gitignore
@@ -0,0 +1 @@
+!Makefile
diff --git a/misc/Makefile b/misc/Makefile
new file mode 100644
index 0000000..8aef648
--- /dev/null
+++ b/misc/Makefile
@@ -0,0 +1,56 @@
+LIBSTATIC = libmisc.a
+
+CC = gcc
+LD = gcc
+AR = ar
+AWKPROG = ../gawk
+CMP = cmp
+
+AWK = LC_ALL=$${GAWKLOCALE:-C} LANG=$${GAWKLOCALE:-C} $(AWKPROG)
+
+ALL_CFLAGS = $(CFLAGS) -fPIC -DGAWK -DHAVE_CONFIG_H -c -I.. -I.
+
+all: $(LIBSTATIC)
+
+OBJS := float128.o
+
+$(LIBSTATIC): $(OBJS)
+       $(AR) rc $(LIBSTATIC) $(OBJS)
+       ranlib $(LIBSTATIC)
+
+float128.o: float128.c ../awk.h ../long_double.h ./gawk_math.h ./gawk_math.c
+
+.c.o:
+       $(CC) $(ALL_CFLAGS) $< -o $@
+
+clean:
+       rm -f $(OBJS) $(LIBSTATIC)
+
+distclean: clean
+
+check:
+
+LDBL_TESTS = \
+       ldblint64 \
+       ldblint128
+
+# An attempt to print something that can be grepped for in build logs
+pass-fail:
+       @COUNT=`ls "$(TESTDIR)"/_* 2>/dev/null | wc -l` ; \
+       if test $$COUNT = 0 ; \
+       then    echo ALL TESTS PASSED ; \
+       else    echo $$COUNT TESTS FAILED ; \
+       fi
+
+ldbl-tests: $(LDBL_TESTS)
+       @$(MAKE) pass-fail TESTDIR=ldbl-tests
+
+ldblint64:
+       @echo $@
+       @$(AWK) -B -f ldbl-tests/address@hidden > ldbl-tests/_$@ 2>&1
+       @-$(CMP) ldbl-tests/address@hidden ldbl-tests/_$@ && rm -f 
ldbl-tests/_$@
+
+ldblint128:
+       @echo $@
+       @$(AWK) -B1 -f ldbl-tests/address@hidden > ldbl-tests/_$@ 2>&1
+       @-$(CMP) ldbl-tests/address@hidden ldbl-tests/_$@ && rm -f 
ldbl-tests/_$@
diff --git a/misc/float128.c b/misc/float128.c
new file mode 100644
index 0000000..3f6ae12
--- /dev/null
+++ b/misc/float128.c
@@ -0,0 +1,447 @@
+/*
+ * float128.c - routines for a quad-precision (binary 128) floating-point.
+ */
+
+/* 
+ * Copyright (C) 2013 the Free Software Foundation, Inc.
+ * 
+ * This file is part of GAWK, the GNU implementation of the
+ * AWK Programming Language.
+ * 
+ * GAWK is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * GAWK is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, 
USA
+ */
+
+#include "awk.h"
+
+#if defined(LDBLTEST) && LDBLTEST == 1
+#include <math.h>
+#include "random.h"
+#if 0
+#include "floatmagic.h"        /* definition of isnan -- XXX: only for float 
or double or long double */
+#endif
+
+#include "format.h"
+
+#include <gmp.h>
+#include <mpfr.h>
+#ifndef MPFR_RNDN
+/* for compatibility with MPFR 2.X */
+#define MPFR_RNDN GMP_RNDN
+#define MPFR_RNDZ GMP_RNDZ
+#define MPFR_RNDU GMP_RNDU
+#define MPFR_RNDD GMP_RNDD
+#endif
+
+/* XXX: No libquadmath, and don't want or need it anyway. */
+
+/*
+ * Macros copied from quadmath.h
+ *     https://github.com/mirrors/gcc/blob/master/libquadmath/quadmath.h
+ */
+
+#define FLT128_MAX 1.18973149535723176508575932662800702e4932Q
+#define FLT128_MIN 3.36210314311209350626267781732175260e-4932Q
+#define FLT128_EPSILON 1.92592994438723585305597794258492732e-34Q
+#define FLT128_DENORM_MIN 6.475175119438025110924438958227646552e-4966Q
+#define FLT128_MANT_DIG 113
+#define FLT128_MIN_EXP (-16381)
+#define FLT128_MAX_EXP 16384
+#define FLT128_DIG 33
+#define FLT128_MIN_10_EXP (-4931)
+#define FLT128_MAX_10_EXP 4932
+
+/* #define HUGE_VALQ __builtin_huge_valq() */
+/* The following alternative is valid, but brings the warning:
+   (floating constant exceeds range of ‘__float128’)  */
+#define HUGE_VALQ (__extension__ 0x1.0p32767Q)
+
+#define IEEE_FLOAT128_BIAS 0x3fff
+
+/* end of Macros from quadmath.h */
+
+
+#define LDBL_FRAC_BITS FLT128_MANT_DIG
+#define        LDBL_INT_BITS   128
+
+#define AWKLDBL        __float128
+#define LDBL_VAL(n)    (*((AWKLDBL *) (n)->qnumbr))
+#define LDC(x)         x##Q
+
+#define get_long_double(d)     emalloc(d, void *, sizeof(AWKLDBL), "float128")
+#define free_long_double(d)    efree(d)
+
+/* we want to format integers ourself */
+#define GAWK_FMT_INT 1
+
+#define gawk_int_t long long
+#define gawk_uint_t unsigned long long
+#ifdef SIZEOF_GAWK_INT
+#undef SIZEOF_GAWK_INT
+#undef GAWK_INT_MAX
+#undef GAWK_INT_MIN
+#undef GAWK_UINT_MAX
+#endif
+
+#define SIZEOF_GAWK_INT        8
+#define        GAWK_INT_MAX    LLONG_MAX
+#define        GAWK_INT_MIN    LLONG_MIN
+#define        GAWK_UINT_MAX   ULLONG_MAX
+
+static int format_uint_finite_p(char *str, size_t size, AWKLDBL x);
+static AWKLDBL gawk_floorl_finite_p(AWKLDBL x, gawk_uint_t *chunk);
+static int format_float_1(char *str, size_t size, const char *format, int fw, 
int prec, AWKLDBL x);
+static inline int format_uint_1(char *str, size_t size, AWKLDBL x);
+static AWKLDBL double_to_int(AWKLDBL x);
+
+static __float128 strtofloat128(const char *str, char **endptr);
+#define gawk_strtold strtofloat128
+
+#if 0
+/* if have __float128, also have strtold. No ? */
+
+static inline AWKLDBL
+gawk_strtold(const char *str, char **endptr)
+{
+       return strtold(str, endptr);
+}
+#endif
+
+#define GAWK_INFINITY  HUGE_VALQ
+#define GAWK_NAN       (LDC(0.0) / LDC(0.0))
+
+/* Define isnan() and isinf() before including the two .h files */
+
+static inline int isnan_awkldbl(AWKLDBL x) { return x != x; }
+#ifdef isnan
+#undef isnan
+#endif
+#define isnan isnan_awkldbl
+
+static inline int isinf_awkldbl(AWKLDBL x) { return isnan(x - x); }
+#ifdef isinf
+#undef isinf
+#endif
+#define isinf isinf_awkldbl
+
+numbr_handler_t float128_hndlr;
+
+#define awkldbl_hndlr float128_hndlr
+
+#define TEST_NUMBR 1
+
+#include "misc/gawk_math.h"
+#include "long_double.h"
+#include "misc/gawk_math.c"
+
+
+/* XXX: This is for little-endian architecture only */
+typedef union {
+       __float128  val;
+
+       /* 64 mantissa LSBs */
+       /* 48 mantissa MSBs */
+       /* 15 exponent bits */
+       /* 1 sign bit */
+
+       struct {
+               uint8_t si[14];
+               uint16_t se;
+       } w1;
+       struct {
+               uint64_t low;
+               uint64_t high;
+       } w2;
+} gfloat128_t;
+
+
+/* hex_to_float128 --- hexadecimal string to __float128 conversion */
+
+static __float128
+hex_to_float128(const char *hexstr)
+{
+       int sign, leading_bit, exp_sign;
+       char *p;
+       uint8_t ch, leading_hex;
+       int exponent;
+       int i, k;
+       gfloat128_t fv;
+
+       memset(& fv, '\0', sizeof (gfloat128_t));
+       
+       p = (char *) hexstr;
+       /* optional sign */
+       sign = 0;
+       if (*p == '-') {
+               sign = 1;
+               p++;
+       } else if (*p == '+')
+               p++;
+
+#define XX (('a' - '0') - 10)
+
+       /* optional "0x" prefix */
+       if (*p == '0' && *(p + 1) == 'x')
+               p += 2;
+
+       leading_bit = 0;
+       k = 0;
+       if (*p != '0') {
+               /* normalized number */
+               leading_hex = ch = (uint8_t)((*p - '0') % XX);
+               while (ch >= 2) {
+                       ch /= 2;
+                       k++;
+               }               
+               leading_bit = 1;
+       }
+       p++;
+
+       if (*p == '.')  /* may not have a decimal point: 0xfp+0 */
+               *p++;
+
+       /* maximum of 28 hex digits */
+       for (i = 13; i >= 0; i--) {
+               if (*p == 'p')
+                       break;
+               ch = (uint8_t) ((*p++ - '0') % XX) ;
+               fv.w1.si[i] |= (ch << 4);
+               if (*p == 'p')
+                       break;
+               ch = (uint8_t) ((*p++ - '0') % XX);
+               fv.w1.si[i] |= ch;
+       }
+
+       if (k > 0) {
+               /* normalized number leading digit not 1; 1 <= k <= 3 */
+               /* shift mantissa to right by k; adjust 4 MSBs */
+               fv.w2.low = fv.w2.low >> k;     /* drop k LSBs from low part */
+               fv.w2.low |= (fv.w2.high << (64 - k));  /* move k LSBs from low 
to MSBs of high */
+               fv.w2.high = fv.w2.high >> k;   /* drop k LSBs from high */
+               fv.w1.si[13] |= leading_hex << (8 - k); /* move k LSBs from 
leading hex digit ... */
+               fv.w1.se = 0;   /* clear sign + exponent part */
+       }
+
+       *p++;   /* skip `p' */
+
+       /* sign of exponent */
+       exp_sign = 0;
+       if (*p == '-') {
+               exp_sign = 1;
+               p++;
+       } else if (*p == '+')
+               p++;
+       exponent = (int) strtol(p, NULL, 10);
+
+       /* adjust exponent for 1 <= significand < 2 */
+       if (exp_sign == 0)
+               exponent += k;
+       else if (exponent > k)
+               exponent -= k;
+       else {
+               exponent = k - exponent;
+               exp_sign = 0;
+       }
+       assert(exponent >= 0);
+
+       /* calculate biased exponent */
+       if (exponent == IEEE_FLOAT128_BIAS - 1 && exp_sign && ! leading_bit)    
/* de-normal number */
+               exponent = 0;
+       else if (leading_bit) {
+               if (exp_sign)
+                       exponent = IEEE_FLOAT128_BIAS - exponent;
+               else
+                       exponent += IEEE_FLOAT128_BIAS;
+       }
+#if 0
+       fprintf(stderr, "======= hex2float=======\n%s\n", hexstr);
+       fprintf(stderr, "exponent (biased) = %d\n", exponent);
+       fprintf(stderr, "sign = %d\n", sign);
+       fprintf(stderr, "significand = %012llx %016llx\n", fv.w2.high, 
fv.w2.low);
+       fprintf(stderr, "========================\n");
+#endif
+       /* set exponent */
+       fv.w1.se |= (uint16_t) exponent;
+       /* set sign */
+       fv.w1.se |= (uint16_t) (sign << 15);
+       return fv.val;
+}
+
+
+/* float128_to_hex --- __float128 to hexadecimal string conversion */
+
+static const char *
+float128_to_hex(__float128 x)
+{
+       static char buf[64];
+       char *p; 
+       int i, exponent, sign, exp_sign;
+       int extra_bit;
+       gfloat128_t fv;
+
+       /* XXX: we don't need to handle infinity and nan here. */
+
+       fv.val = x;
+       p = buf;
+
+       sign = (fv.w1.se & (1 << 15)) != 0;
+       exponent = (fv.w1.se & ((1 << 15) - 1));
+       extra_bit = 1;
+       exp_sign = 0;
+
+       /* calculate unbiased exponent */
+       if (exponent == 0) {
+               int nonzero_significand = ((fv.w2.high & ((1ULL << 48) - 1)) | 
fv.w2.low) != 0;
+
+               extra_bit = 0;
+               if (nonzero_significand) {
+                       exp_sign = 1;
+                       exponent = IEEE_FLOAT128_BIAS - 1;
+               }
+       } else if (exponent >= IEEE_FLOAT128_BIAS)
+               exponent -= IEEE_FLOAT128_BIAS;
+       else {
+               exp_sign = 1;
+               exponent = -(exponent - IEEE_FLOAT128_BIAS);    
+       }
+
+       if (sign)
+               *p++ = '-';
+       *p++ = extra_bit ? '1' : '0';
+       *p++ = '.';
+
+       /* all 28 hex digits -- XXX: quit early if last 64-bits all 0s ? */
+       for (i = 13; i >= 0; i--) {
+               *p++ = lchbuf[(fv.w1.si[i] & 0xF0) >> 4];
+               *p++ = lchbuf[(fv.w1.si[i] & 0x0F)];
+       }
+
+       *p++ = 'p';
+       *p++ = exp_sign ? '-' : '+';
+       sprintf(p, "%d", exponent);     /* decimal exponent */
+
+       /* assert(hex_to_float128(buf) == x); */
+       return buf;
+}
+
+/*
+ * format_float_1 --- format a single AWKLDBL value according to FORMAT.
+ *     The value must be finite.
+ */
+
+static int
+format_float_1(char *str, size_t size, const char *format, int fw, int prec, 
AWKLDBL x)
+{
+       char alt_format[32];
+       size_t len;
+       int ret = -1;
+       const char *hexstr;
+       mpfr_t ap_float;
+
+       if (isnan(x) || isinf(x)) {
+               /* FIXME: ignoring sign */
+               if (size < 4)
+                       return 3;
+               strcpy(str, isnan(x) ? "nan" : "inf");
+               return 3;
+       }
+
+       len = strlen(format);
+
+       /* expect %Lf, %LF, %Le, %LE, %Lg or %LG */
+       assert(len >= 2 && format[len - 2] == 'L');
+       
+       memcpy(alt_format, format, len + 1);
+       alt_format[len - 2] = 'R';      /* replace `L' with `R' */ 
+       alt_format[len] = alt_format[len - 1];
+       alt_format[len - 1] = 'N';      /* rounding mode = MPFR_RNDN */
+       alt_format[len + 1] = '\0';
+
+       /* loss-less conversion to hexadecimal string */
+       hexstr = float128_to_hex(x);
+
+       /* fprintf(stderr, "hexfloat = %s\n", hexstr); */
+
+       mpfr_init2(ap_float, 113);
+       if (mpfr_set_str(ap_float, hexstr, 16, MPFR_RNDN) < 0) {
+               /* XXX: Invalid string? can it be? */
+               return -1;
+       }
+
+       ret = mpfr_snprintf(str, size, alt_format, fw, prec, ap_float);
+       mpfr_clear(ap_float);
+       return ret;
+}
+
+/*
+ * format_uint_1 --- format a long double as an unsigned integer. The double 
value
+ *     must be finite and >= 0.
+ */
+
+static inline int
+format_uint_1(char *str, size_t size, AWKLDBL x)
+{
+       int ret;
+       if ((ret = format_uint_finite_p(str, size, x)) < 0) {
+               /* outside useful range */
+               /* FIXME -- use format_float_1 (MPFR) */
+               return snprintf(str, size, "%.0Lf", (long double) x);
+       }
+       return ret;
+}
+
+/* double_to_int --- convert double to int, used in several places */
+
+static AWKLDBL
+double_to_int(AWKLDBL x)
+{
+       AWKLDBL intval;
+       int sign = 1;
+
+       if (isnan(x) || isinf(x) || x == LDC(0.0))
+               return x;
+       if (x < LDC(0.0)) {
+               sign = -1;
+               x = -x;
+       }
+       if ((intval = gawk_floorl_finite_p(x, NULL)) < LDC(0.0)) {
+               /* outside range */
+               /* FIXME -- use format_float_1 with "%.0Lf" (MPFR) */
+               intval = (AWKLDBL) floorl((long double) x);
+       }
+       return intval * ((AWKLDBL) sign);
+}
+
+/* strtofloat128 --- convert string to __float128 */
+
+static __float128
+strtofloat128(const char *str, char **endptr)
+{
+       static char hexstr[64];
+       mpfr_t ap_float;
+
+       /* store as MPFR number */
+       mpfr_init2(ap_float, 113);
+       (void) mpfr_strtofr(ap_float, str, endptr, 10, MPFR_RNDN);
+
+       /* get back a hexadecimal string representation */
+       (void) mpfr_snprintf(hexstr, 64, "%RNa", ap_float);
+       mpfr_clear(ap_float);
+
+       /* fprintf(stderr, "hexfloat (mpfr): %s\n", hexstr); */
+
+       return hex_to_float128(hexstr); 
+}
+
+#endif /* LDBLTEST == 1 */
diff --git a/misc/fp_math.awk b/misc/fp_math.awk
index 80fe7d8..e2847e6 100644
--- a/misc/fp_math.awk
+++ b/misc/fp_math.awk
@@ -4,7 +4,6 @@
 #      TODO:   * finish fmod().
 #              * replace all usages of %, and int() or any other math builtin; 
and() is ok!.
 #                implement double_to_int() and remove floor/ceil.
-#              * fix sqrt(x) for x > DBL_MAX or x < -DBL_MIN.
 #
 
 
@@ -337,19 +336,37 @@ function fp_sqrt(x,       \
        if (x == 0)
                return x        # return +0 or -0
 
-       yn = sqrt(x)    # double-precision sqrt
+       if (x >= DBL_MIN && x <= DBL_MAX)
+               yn = sqrt(x)    # double-precision sqrt
+       else {
+               frac = frexpl(x, e);    # in [1, 2)
+               exponent = e[0]
+               if (and(exponent, 1) != 0) {
+                       frac /= 2
+                       exponent++
+               }
+               yn = sqrt(frac)
+               exponent /= 2
+               if (exponent >= 0)
+                       yn *= pow2(exponent)
+               else
+                       yn /= pow2(-exponent)
+       }
+
        do {
                last = yn
                yn = (yn + x / yn) \
                                / 2
                err = (yn - last) / yn
+               if (err < 0)
+                       err = -err
        } while (err > __REL_ERROR__)
        return yn
 }
 
-# __log2 -- compute log(2)
+# klog2 -- compute k * log2 without loss of precision 
 
-function __log2( \
+function klog2(k, \
                i, three_pow, sum, err, term)
 {
        #
@@ -360,18 +377,18 @@ function __log2( \
        # y = (x - 1) / (x + 1) => y = 1 / 3
        i = 1
        three_pow = 3
-       sum = 1 / 3
-
+       sum = (2 * k) / 3
        do {
                three_pow *= 9
                i += 2
-               term = 1 / (three_pow * i)
+               term = (2 * k) / (three_pow * i)
                sum += term
                err = term / sum
        } while (err > __REL_ERROR__)
-       return 2 * sum
+       return sum
 }
 
+
 function pow2(n, k)
 {
        if (n <= 64)
@@ -381,28 +398,17 @@ function pow2(n, k)
        return k * __POW2__[n]
 }
 
-function fp_log(x,     \
-               k, i, y, ypow2, ypow_odd, sum, err, term, sign, high, low, mid)
+
+# the returned fraction is normalized in [1, 2) and NOT [0.5, 1)  
+# exponent is in e[0]
+
+function frexpl(x, e, \
+               y, low, high, mid)
 {
-       #
-       #       ln(x) = 2 * arctanh(y)
-       #             = 2 * (y + y^3 / 3 + y^5 /5 + ..) where y = (x - 1) / (x 
+ 1) 
-       #
-       x = +x
-       if (isnan(x) || x == __PLUS_INF__)
-               return x
-       if (x == __MINUS_INF__ || x < 0)        # errno EDOM ?
-               return __PLUS_NAN__
-       if (x == 0)     # errno ERANGE ?
-               return __MINUS_INF__
-       
-       if (x == 1)     # special case
-               return 0
-       if (x == 2)     # special case
-               return __LOG2__
-       k = 0
+       # XXX -- (isnormal(x) && x > 0) is assumed to be true
+
+       low = 0
        if (x > 2) {
-               low = 0
                high = LDBL_MAX_EXP - 1         # XXX: should be 4 * 
LDBL_MAX_EXP - 1 if FLT_RADIX = 16
                while (low <= high) {
                        mid = int ((low + high) / 2)
@@ -413,11 +419,8 @@ function fp_log(x, \
                                high = mid - 1
                }
                x /= pow2(low)
-               k = low
-#              printf("x = %0.18e, k = %d\n", x / pow2(low), low)
-
+               e[0] = low
        } else if (x < 1) {
-               low = 0
                high = LDBL_MAX_EXP - 1         # could be -LDBL_MIN_EXP, but 
no harm in using LDBL_MAX_EXP
                while (low <= high) {
                        mid = int ((low + high) / 2)
@@ -428,17 +431,64 @@ function fp_log(x,        \
                                high = mid - 1
                }
                x *= pow2(low)
-               k = -low
+               e[0] = -low
+       } else
+               e[0] = 0
+
+       if (x == 2) {
+               x = 1
+               e[0]++
        }
+       return x        # x in [1, 2)
+}
+
+
+#
+# $ ./gawk 'BEGIN { printf("%.18e\n", log(7.12111e900))}'
+# inf
+# $ ./gawk -B 'BEGIN { printf("%.18e\n", log(7.12111e900))}'           # with 
logl(), without it inf
+# 2.074289647306790437e+03
+# $ ./gawk -B -f misc/fp_math.awk -e 'BEGIN { printf("%.18e\n", 
fp_log(7.12111e900))}'
+# 2.074289647306790437e+03
+# $ ./gawk -M -vPREC=quad 'BEGIN { printf("%.18e\n", log(7.12111e900))}'
+# 2.074289647306790438e+03
+#
+
+
+function fp_log(x,     \
+               e, exponent, i, ypow2, ypow_odd, sum, err, term, sign)
+{
+       #
+       #       ln(x) = 2 * arctanh(y)
+       #             = 2 * (y + y^3 / 3 + y^5 /5 + ..) where y = (x - 1) / (x 
+ 1) 
+       #
+       x = +x
+       if (isnan(x) || x == __PLUS_INF__)
+               return x
+       if (x == __MINUS_INF__ || x < 0)        # errno EDOM ?
+               return __PLUS_NAN__
+       if (x == 0)     # errno ERANGE ?
+               return __MINUS_INF__
+       
+       if (x == 1)     # special case
+               return 0
+       if (x == 2)     # special case
+               return __LOG2__
+
+       x = frexpl(x, e)
+       exponent = e[0]
 
        # arctanh(x) series has faster convergence when x is close to 1
        # range reduction -- 1/sqrt(2) <= x <= sqrt(2)
        if (x > __SQRT_2__ || x < __1_OVER_SQRT_2__) {
                x *= __1_OVER_SQRT_2__
-               k += 0.5
+               exponent += 0.5
        }
 
        y = (x - 1) / (x + 1)
+       if (y == 0)     # tricky special case 
+               return exponent * __LOG2__;
+
        sign = 1
        if (y < 0) {
                sign = -1
@@ -458,7 +508,7 @@ function fp_log(x,  \
                err = term / sum
        } while (err > __REL_ERROR__)
 
-       return sign * (sum = 2 * sum + k * __LOG2__)
+       return sign * (sum = 2 * sum + exponent * __LOG2__)
 }
 
 
@@ -470,6 +520,8 @@ function taylor_exp(x,      \
        # sinh(x) = sqrt(cosh(x) * cosh(x) - 1)
        # exp(x) = cosh(x) + sinh(x)
 
+       if (x == 0)
+               return 1;
        xpow2 = x * x
        xpow_even = 1
        i = 0
@@ -485,12 +537,15 @@ function taylor_exp(x,    \
                err  = term / cosh
        } while (err > __REL_ERROR__)
 
-       sinh = fp_sqrt(cosh * cosh - 1) # sqrt is cheap
+#      sinh = fp_sqrt(cosh * cosh - 1) # sqrt is cheap
+       z = cosh - 1.0
+       sinh = z * fp_sqrt(1 + 2 / z)   
        return (sinh + cosh)
 }
 
+
 function fp_exp(x,     \
-               exp_val, k, sign)
+               y, exp_val, k, sign)
 {
        x = +x
        if (isnan(x) || x == __PLUS_INF__)
@@ -505,19 +560,18 @@ function fp_exp(x,        \
                sign = -1
                x = -x
        }
+
        k = fp_floor(x / __LOG2__)
 
        # range reduction -- 0 < x < log(2) (0.693...)
        if (k == 0)
                exp_val = taylor_exp(x)
        else {
-               exp_val = taylor_exp(x - k * __LOG2__)
-               while (k >= 64) {
-                       exp_val *= __POW2__[64]
-                       k -= 64
-               }
-               exp_val *= __POW2__[k]
+               y = x - klog2(k)
+               exp_val = taylor_exp(y)
+               exp_val *= pow2(k)
        }
+
        return sign < 0 ? (1 / exp_val) : exp_val
 }
 
@@ -791,6 +845,8 @@ BEGIN {
        __REL_ERROR__ = 5.0e-35
        LDBL_MAX_EXP = 16384
        LDBL_MANT_DIG = 113
+       DBL_MAX = 1.7976931348623157e+308
+       DBL_MIN = 2.2250738585072014E-308
 
        # We can read hex/octal numbers without using strtod/strtodl
        if (0x10000000000000000 == 0x10000000000000001) {
@@ -812,7 +868,7 @@ BEGIN {
 
        __SQRT_2__ = fp_sqrt(2)
        __1_OVER_SQRT_2__ = 1 / __SQRT_2__
-       __LOG2__ = __log2()     # cannot use fp_log()
+       __LOG2__ = klog2(1)     # cannot use fp_log()
        __PI_OVER_4__ = 4 * euler_atan_one_over(5) - euler_atan_one_over(239)
 
        # pre-calculate 2^0 - 2^64
diff --git a/misc/gawk_math.c b/misc/gawk_math.c
new file mode 100644
index 0000000..3e8e6ef
--- /dev/null
+++ b/misc/gawk_math.c
@@ -0,0 +1,72 @@
+/*
+ * gawk_math.c - routines for replacement AWKLDBL math functions.
+ */
+
+/* 
+ * Copyright (C) 2013 the Free Software Foundation, Inc.
+ * 
+ * This file is part of GAWK, the GNU implementation of the
+ * AWK Programming Language.
+ * 
+ * GAWK is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * GAWK is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, 
USA
+ */
+
+static AWKLDBL
+gawk_sinl(AWKLDBL x)
+{
+       return sin( (double) x);
+}
+
+static AWKLDBL
+gawk_cosl(AWKLDBL x)
+{
+       return cos( (double) x);
+}
+
+static AWKLDBL
+gawk_atan2l(AWKLDBL y, AWKLDBL x)
+{
+       return atan2( (double) y, (double) x);
+}
+
+static AWKLDBL
+gawk_logl(AWKLDBL x)
+{
+       return log( (double) x);
+}
+
+static AWKLDBL
+gawk_expl(AWKLDBL x)
+{
+       return exp( (double) x);
+}
+
+static AWKLDBL
+gawk_fmodl(AWKLDBL x, AWKLDBL y)
+{
+       return fmod( (double) x, (double) y);
+}
+
+static AWKLDBL
+gawk_powl(AWKLDBL x, AWKLDBL y)
+{
+       return pow( (double) x, (double) y);
+}
+
+static AWKLDBL
+gawk_sqrtl(AWKLDBL x)
+{
+       return sqrt( (double) x);
+}
diff --git a/misc/gawk_math.h b/misc/gawk_math.h
index 909bd69..5109d55 100644
--- a/misc/gawk_math.h
+++ b/misc/gawk_math.h
@@ -1,3 +1,28 @@
+/*
+ * gawk_math.h - replacement AWKLDBL math functions.
+ */
+
+/* 
+ * Copyright (C) 2013 the Free Software Foundation, Inc.
+ * 
+ * This file is part of GAWK, the GNU implementation of the
+ * AWK Programming Language.
+ * 
+ * GAWK is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * GAWK is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, 
USA
+ */
+
 static AWKLDBL gawk_sinl(AWKLDBL x);
 static AWKLDBL gawk_cosl(AWKLDBL x);
 static AWKLDBL gawk_atan2l(AWKLDBL y, AWKLDBL x);
@@ -7,82 +32,3 @@ static AWKLDBL gawk_fmodl(AWKLDBL x, AWKLDBL y);
 static AWKLDBL gawk_powl(AWKLDBL x, AWKLDBL y);
 static AWKLDBL gawk_sqrtl(AWKLDBL x);
 
-static AWKLDBL
-gawk_sinl(AWKLDBL x)
-{
-#ifdef HAVE_SINL
-       return sinl( (long double) x);
-#else
-       return sin( (double) x);
-#endif
-}
-
-static AWKLDBL
-gawk_cosl(AWKLDBL x)
-{
-#ifdef HAVE_COSL
-       return cosl( (long double) x);
-#else
-       return cos( (double) x);
-#endif
-}
-
-static AWKLDBL
-gawk_atan2l(AWKLDBL y, AWKLDBL x)
-{
-#ifdef HAVE_ATAN2L
-       return atan2l( (long double) y, (long double) x);
-#else
-       return atan2( (double) y, (double) x);
-#endif
-}
-
-static AWKLDBL
-gawk_logl(AWKLDBL x)
-{
-#ifdef HAVE_LOGL
-       return logl( (long double) x);
-#else
-       return log( (double) x);
-#endif
-}
-
-static AWKLDBL
-gawk_expl(AWKLDBL x)
-{
-#ifdef HAVE_EXPL
-       return expl( (long double) x);
-#else
-       return exp( (double) x);
-#endif
-}
-
-static AWKLDBL
-gawk_fmodl(AWKLDBL x, AWKLDBL y)
-{
-#ifdef HAVE_FMODL
-       return fmodl( (long double) x, (long double) y);
-#else
-       return fmod( (double) x, (double) y);
-#endif
-}
-
-static AWKLDBL
-gawk_powl(AWKLDBL x, AWKLDBL y)
-{
-#ifdef HAVE_POWL
-       return powl( (long double) x, (long double) y);
-#else
-       return pow( (double) x, (double) y);
-#endif
-}
-
-static AWKLDBL
-gawk_sqrtl(AWKLDBL x)
-{
-#ifdef HAVE_SQRTL
-       return sqrtl( (long double) x);
-#else
-       return sqrt( (double) x);
-#endif 
-}
diff --git a/misc/ldblint64.awk b/misc/ldbl-tests/ldblint128.awk
similarity index 62%
copy from misc/ldblint64.awk
copy to misc/ldbl-tests/ldblint128.awk
index 65f3939..ad2b083 100644
--- a/misc/ldblint64.awk
+++ b/misc/ldbl-tests/ldblint128.awk
@@ -1,4 +1,4 @@
-# Usage gawk -B -f ldblint54.awk
+# Usage gawk -B -f ldblint128.awk
 
 BEGIN {
        for (i = -9.0; i < 2.0; i++) {
@@ -9,6 +9,6 @@ BEGIN {
                        j = -j
                        j = " - " j
                }
-               printf("2^64%s = %d\n", j, 2^64 + i)
+               printf("2^113%s = %d\n", j, 2^113 + i)
        }
 }
diff --git a/misc/ldbl-tests/ldblint128.ok b/misc/ldbl-tests/ldblint128.ok
new file mode 100644
index 0000000..0526d3b
--- /dev/null
+++ b/misc/ldbl-tests/ldblint128.ok
@@ -0,0 +1,11 @@
+2^113 - 9 = 10384593717069655257060992658440183
+2^113 - 8 = 10384593717069655257060992658440184
+2^113 - 7 = 10384593717069655257060992658440185
+2^113 - 6 = 10384593717069655257060992658440186
+2^113 - 5 = 10384593717069655257060992658440187
+2^113 - 4 = 10384593717069655257060992658440188
+2^113 - 3 = 10384593717069655257060992658440189
+2^113 - 2 = 10384593717069655257060992658440190
+2^113 - 1 = 10384593717069655257060992658440191
+2^113 + 0 = 10384593717069655257060992658440192
+2^113 + 1 = 10384593717069655257060992658440192
diff --git a/misc/ldblint64.awk b/misc/ldbl-tests/ldblint64.awk
similarity index 100%
rename from misc/ldblint64.awk
rename to misc/ldbl-tests/ldblint64.awk
diff --git a/misc/ldblint64.ok b/misc/ldbl-tests/ldblint64.ok
similarity index 100%
rename from misc/ldblint64.ok
rename to misc/ldbl-tests/ldblint64.ok

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog                                         |    7 +
 Makefile.am                                       |   13 +-
 Makefile.in                                       |   16 +-
 configure                                         |    2 -
 configure.ac                                      |    2 -
 float128.c                                        |  128 ------
 long_double.c                                     |   94 +++--
 long_double.h                                     |  110 ++----
 misc/.gitignore                                   |    1 +
 misc/Makefile                                     |   56 +++
 misc/float128.c                                   |  447 +++++++++++++++++++++
 misc/fp_math.awk                                  |  144 +++++--
 random.h => misc/gawk_math.c                      |   76 +++--
 misc/gawk_math.h                                  |  104 ++----
 misc/{ldblint64.awk => ldbl-tests/ldblint128.awk} |    4 +-
 misc/ldbl-tests/ldblint128.ok                     |   11 +
 misc/{ => ldbl-tests}/ldblint64.awk               |    0
 misc/{ => ldbl-tests}/ldblint64.ok                |    0
 18 files changed, 814 insertions(+), 401 deletions(-)
 delete mode 100644 float128.c
 create mode 100644 misc/.gitignore
 create mode 100644 misc/Makefile
 create mode 100644 misc/float128.c
 copy random.h => misc/gawk_math.c (52%)
 copy misc/{ldblint64.awk => ldbl-tests/ldblint128.awk} (62%)
 create mode 100644 misc/ldbl-tests/ldblint128.ok
 rename misc/{ => ldbl-tests}/ldblint64.awk (100%)
 rename misc/{ => ldbl-tests}/ldblint64.ok (100%)


hooks/post-receive
-- 
gawk



reply via email to

[Prev in Thread] Current Thread [Next in Thread]