main/musl: include git commits up to 71d23b310383

multiple bug fixes, and sys/timeb.h with ftime() implemented increasing
compatibility with old applications.

7 files changed, 2352 insertions, 206 deletions

diff --git a/main/musl/0001-add-missing-i386-syscall-numbers.patch b/main/musl/0001-add-missing-i386-syscall-numbers.patch
deleted file mode 100644
index 8c0b0a892d..0000000000
--- a/main/musl/0001-add-missing-i386-syscall-numbers.patch
+++ /dev/null
@@ -1,41 +0,0 @@
-From 78178542e73e143bf44b3ba32cf0b58ced53f2d5 Mon Sep 17 00:00:00 2001
-From: Rich Felker <dalias@aerifal.cx>
-Date: Thu, 26 Sep 2013 14:17:36 -0400
-Subject: [PATCH] add missing i386 syscall numbers
-
-somehow the range 335-339 was missed when updating the file.
----
- arch/i386/bits/syscall.h | 10 ++++++++++
- 1 file changed, 10 insertions(+)
-
-diff --git a/arch/i386/bits/syscall.h b/arch/i386/bits/syscall.h
-index 800409a..2af242b 100644
---- a/arch/i386/bits/syscall.h
-+++ b/arch/i386/bits/syscall.h
-@@ -333,6 +333,11 @@
- #define __NR_inotify_init1	332
- #define __NR_preadv		333
- #define __NR_pwritev		334
-+#define __NR_rt_tgsigqueueinfo	335
-+#define __NR_perf_event_open	336
-+#define __NR_recvmmsg		337
-+#define __NR_fanotify_init	338
-+#define __NR_fanotify_mark	339
- #define __NR_prlimit64		340
- #define __NR_name_to_handle_at	341
- #define __NR_open_by_handle_at	342
-@@ -683,6 +688,11 @@
- #define SYS_inotify_init1	332
- #define SYS_preadv		333
- #define SYS_pwritev		334
-+#define SYS_rt_tgsigqueueinfo	335
-+#define SYS_perf_event_open	336
-+#define SYS_recvmmsg		337
-+#define SYS_fanotify_init	338
-+#define SYS_fanotify_mark	339
- #define SYS_prlimit64		340
- #define SYS_name_to_handle_at	341
- #define SYS_open_by_handle_at	342
--- 
-1.8.4
-
diff --git a/main/musl/0001-updates-from-git.patch b/main/musl/0001-updates-from-git.patch
new file mode 100644
index 0000000000..63d2468331
--- /dev/null
+++ b/main/musl/0001-updates-from-git.patch
@@ -0,0 +1,2347 @@
+git diff v0.9.14..71d23b310383
+
+diff --git a/arch/i386/bits/syscall.h b/arch/i386/bits/syscall.h
+index 800409a..2af242b 100644
+--- a/arch/i386/bits/syscall.h
++++ b/arch/i386/bits/syscall.h
+@@ -333,6 +333,11 @@
+ #define __NR_inotify_init1	332
+ #define __NR_preadv		333
+ #define __NR_pwritev		334
++#define __NR_rt_tgsigqueueinfo	335
++#define __NR_perf_event_open	336
++#define __NR_recvmmsg		337
++#define __NR_fanotify_init	338
++#define __NR_fanotify_mark	339
+ #define __NR_prlimit64		340
+ #define __NR_name_to_handle_at	341
+ #define __NR_open_by_handle_at	342
+@@ -683,6 +688,11 @@
+ #define SYS_inotify_init1	332
+ #define SYS_preadv		333
+ #define SYS_pwritev		334
++#define SYS_rt_tgsigqueueinfo	335
++#define SYS_perf_event_open	336
++#define SYS_recvmmsg		337
++#define SYS_fanotify_init	338
++#define SYS_fanotify_mark	339
+ #define SYS_prlimit64		340
+ #define SYS_name_to_handle_at	341
+ #define SYS_open_by_handle_at	342
+diff --git a/arch/mips/bits/syscall.h b/arch/mips/bits/syscall.h
+index 2726019..574d603 100644
+--- a/arch/mips/bits/syscall.h
++++ b/arch/mips/bits/syscall.h
+@@ -318,30 +318,6 @@
+ #define __NR_timerfd                 4318
+ #define __NR_eventfd                 4319
+ #define __NR_fallocate               4320
+-#define __NR_fallocate               4320
+-#define __NR_timerfd_create          4321
+-#define __NR_timerfd_gettime         4322
+-#define __NR_timerfd_settime         4323
+-#define __NR_signalfd4               4324
+-#define __NR_eventfd2                4325
+-#define __NR_epoll_create1           4326
+-#define __NR_dup3                    4327
+-#define __NR_pipe2                   4328
+-#define __NR_inotify_init1           4329
+-#define __NR_preadv                  4330
+-#define __NR_pwritev                 4331
+-#define __NR_rt_tgsigqueueinfo       4332
+-#define __NR_perf_event_open         4333
+-#define __NR_accept4                 4334
+-#define __NR_recvmmsg                4335
+-#define __NR_fanotify_init           4336
+-#define __NR_fanotify_mark           4337
+-#define __NR_prlimit64               4338
+-#define __NR_name_to_handle_at       4339
+-#define __NR_open_by_handle_at       4340
+-#define __NR_clock_adjtime           4341
+-#define __NR_syncfs                  4342
+-#define __NR_fallocate               4320
+ #define __NR_timerfd_create          4321
+ #define __NR_timerfd_gettime         4322
+ #define __NR_timerfd_settime         4323
+@@ -693,30 +669,6 @@
+ #define SYS_timerfd                 4318
+ #define SYS_eventfd                 4319
+ #define SYS_fallocate               4320
+-#define SYS_fallocate               4320
+-#define SYS_timerfd_create          4321
+-#define SYS_timerfd_gettime         4322
+-#define SYS_timerfd_settime         4323
+-#define SYS_signalfd4               4324
+-#define SYS_eventfd2                4325
+-#define SYS_epoll_create1           4326
+-#define SYS_dup3                    4327
+-#define SYS_pipe2                   4328
+-#define SYS_inotify_init1           4329
+-#define SYS_preadv                  4330
+-#define SYS_pwritev                 4331
+-#define SYS_rt_tgsigqueueinfo       4332
+-#define SYS_perf_event_open         4333
+-#define SYS_accept4                 4334
+-#define SYS_recvmmsg                4335
+-#define SYS_fanotify_init           4336
+-#define SYS_fanotify_mark           4337
+-#define SYS_prlimit64               4338
+-#define SYS_name_to_handle_at       4339
+-#define SYS_open_by_handle_at       4340
+-#define SYS_clock_adjtime           4341
+-#define SYS_syncfs                  4342
+-#define SYS_fallocate               4320
+ #define SYS_timerfd_create          4321
+ #define SYS_timerfd_gettime         4322
+ #define SYS_timerfd_settime         4323
+diff --git a/include/arpa/ftp.h b/include/arpa/ftp.h
+index 4041aeb..fb0a46f 100644
+--- a/include/arpa/ftp.h
++++ b/include/arpa/ftp.h
+@@ -1,5 +1,5 @@
+-#ifndef _ARPA_FTP_H_
+-#define _ARPA_FTP_H_
++#ifndef _ARPA_FTP_H
++#define _ARPA_FTP_H
+ #define PRELIM 1
+ #define COMPLETE 2
+ #define CONTINUE 3
+diff --git a/include/fcntl.h b/include/fcntl.h
+index b9bc269..55a89f9 100644
+--- a/include/fcntl.h
++++ b/include/fcntl.h
+@@ -59,8 +59,6 @@ int posix_fallocate(int, off_t, off_t);
+ #define AT_REMOVEDIR 0x200
+ #define AT_SYMLINK_FOLLOW 0x400
+ #define AT_EACCESS 0x200
+-#define AT_NO_AUTOMOUNT 0x800
+-#define AT_EMPTY_PATH 0x1000
+ 
+ #define POSIX_FADV_NORMAL     0
+ #define POSIX_FADV_RANDOM     1
+@@ -95,6 +93,9 @@ int posix_fallocate(int, off_t, off_t);
+ #endif
+ 
+ #if defined(_GNU_SOURCE) || defined(_BSD_SOURCE)
++#define AT_NO_AUTOMOUNT 0x800
++#define AT_EMPTY_PATH 0x1000
++
+ #define FAPPEND O_APPEND
+ #define FFSYNC O_FSYNC
+ #define FASYNC O_ASYNC
+diff --git a/include/langinfo.h b/include/langinfo.h
+index c6349ad..2153c42 100644
+--- a/include/langinfo.h
++++ b/include/langinfo.h
+@@ -5,6 +5,7 @@
+ extern "C" {
+ #endif
+ 
++#include <features.h>
+ #include <nl_types.h>
+ 
+ #define __NEED_locale_t
+@@ -75,8 +76,11 @@ extern "C" {
+ #define THOUSEP 0x10001
+ #define YESEXPR 0x50000
+ #define NOEXPR 0x50001
++
++#if defined(_GNU_SOURCE) || defined(_BSD_SOURCE)
+ #define YESSTR 0x50002
+ #define NOSTR 0x50003
++#endif
+ 
+ char *nl_langinfo(nl_item);
+ char *nl_langinfo_l(nl_item, locale_t);
+diff --git a/include/limits.h b/include/limits.h
+index a8460cc..574b406 100644
+--- a/include/limits.h
++++ b/include/limits.h
+@@ -84,12 +84,18 @@
+ #define NL_ARGMAX 9
+ #define NL_LANGMAX 32
+ #define NL_MSGMAX 32767
+-#define NL_NMAX (MB_LEN_MAX*4)
+ #define NL_SETMAX 255
+ #define NL_TEXTMAX 2048
+ 
+ #endif
+ 
++#if defined(_GNU_SOURCE) || defined(_BSD_SOURCE) \
++ || (defined(_XOPEN_SOURCE) && _XOPEN_SOURCE+0 < 700)
++
++#define NL_NMAX 16
++
++#endif
++
+ /* POSIX/SUS requirements follow. These numbers come directly
+  * from SUS and have nothing to do with the host system. */
+ 
+diff --git a/include/math.h b/include/math.h
+index c029156..dc17601 100644
+--- a/include/math.h
++++ b/include/math.h
+@@ -91,20 +91,20 @@ int __signbitl(long double);
+ static __inline int __is##rel(type __x, type __y) \
+ { return !isunordered(__x,__y) && __x op __y; }
+ 
+-__ISREL_DEF(lessf, <, float)
+-__ISREL_DEF(less, <, double)
++__ISREL_DEF(lessf, <, float_t)
++__ISREL_DEF(less, <, double_t)
+ __ISREL_DEF(lessl, <, long double)
+-__ISREL_DEF(lessequalf, <=, float)
+-__ISREL_DEF(lessequal, <=, double)
++__ISREL_DEF(lessequalf, <=, float_t)
++__ISREL_DEF(lessequal, <=, double_t)
+ __ISREL_DEF(lessequall, <=, long double)
+-__ISREL_DEF(lessgreaterf, !=, float)
+-__ISREL_DEF(lessgreater, !=, double)
++__ISREL_DEF(lessgreaterf, !=, float_t)
++__ISREL_DEF(lessgreater, !=, double_t)
+ __ISREL_DEF(lessgreaterl, !=, long double)
+-__ISREL_DEF(greaterf, >, float)
+-__ISREL_DEF(greater, >, double)
++__ISREL_DEF(greaterf, >, float_t)
++__ISREL_DEF(greater, >, double_t)
+ __ISREL_DEF(greaterl, >, long double)
+-__ISREL_DEF(greaterequalf, >=, float)
+-__ISREL_DEF(greaterequal, >=, double)
++__ISREL_DEF(greaterequalf, >=, float_t)
++__ISREL_DEF(greaterequal, >=, double_t)
+ __ISREL_DEF(greaterequall, >=, long double)
+ 
+ #define __tg_pred_2(x, y, p) ( \
+diff --git a/include/signal.h b/include/signal.h
+index e65a806..6f10a118 100644
+--- a/include/signal.h
++++ b/include/signal.h
+@@ -218,11 +218,8 @@ void (*sigset(int, void (*)(int)))(int);
+ #define SIGSTKSZ 8192
+ #endif
+ 
+-#if defined(_XOPEN_SOURCE) || defined(_GNU_SOURCE)
+-#define NSIG _NSIG
+-#endif
+-
+ #if defined(_BSD_SOURCE) || defined(_GNU_SOURCE)
++#define NSIG _NSIG
+ typedef void (*sig_t)(int);
+ #endif
+ 
+diff --git a/include/sys/socket.h b/include/sys/socket.h
+index 9fd2025..20eeee3 100644
+--- a/include/sys/socket.h
++++ b/include/sys/socket.h
+@@ -225,7 +225,7 @@ struct linger
+ #define MSG_EOR       0x0080
+ #define MSG_WAITALL   0x0100
+ #define MSG_FIN       0x0200
+-#define MSD_SYN       0x0400
++#define MSG_SYN       0x0400
+ #define MSG_CONFIRM   0x0800
+ #define MSG_RST       0x1000
+ #define MSG_ERRQUEUE  0x2000
+diff --git a/include/sys/timeb.h b/include/sys/timeb.h
+new file mode 100644
+index 0000000..108c1f5
+--- /dev/null
++++ b/include/sys/timeb.h
+@@ -0,0 +1,22 @@
++#ifndef _SYS_TIMEB_H
++#define _SYS_TIMEB_H
++#ifdef __cplusplus
++extern "C" {
++#endif
++
++#define __NEED_time_t
++
++#include <bits/alltypes.h>
++
++struct timeb {
++	time_t time;
++	unsigned short millitm;
++	short timezone, dstflag;
++};
++
++int ftime(struct timeb *);
++
++#ifdef __cplusplus
++}
++#endif
++#endif
+diff --git a/include/time.h b/include/time.h
+index 6b2a069..6e499ff 100644
+--- a/include/time.h
++++ b/include/time.h
+@@ -82,8 +82,8 @@ struct itimerspec
+ #define CLOCK_PROCESS_CPUTIME_ID 2
+ #define CLOCK_THREAD_CPUTIME_ID  3
+ #define CLOCK_MONOTONIC_RAW      4
+-#define CLOCK_REALTIME_COURSE    5
+-#define CLOCK_MONOTONIC_COURSE   6
++#define CLOCK_REALTIME_COARSE    5
++#define CLOCK_MONOTONIC_COARSE   6
+ #define CLOCK_BOOTTIME           7
+ #define CLOCK_REALTIME_ALARM     8
+ #define CLOCK_BOOTTIME_ALARM     9
+diff --git a/src/env/setenv.c b/src/env/setenv.c
+index c2c2544..76e8ee1 100644
+--- a/src/env/setenv.c
++++ b/src/env/setenv.c
+@@ -18,14 +18,13 @@ int setenv(const char *var, const char *value, int overwrite)
+ 	l1 = strlen(var);
+ 	l2 = strlen(value);
+ 	s = malloc(l1+l2+2);
+-	memcpy(s, var, l1);
+-	s[l1] = '=';
+-	memcpy(s+l1+1, value, l2);
+-	s[l1+l2+1] = 0;
+-	if (__putenv(s, 1)) {
+-		free(s);
+-		errno = ENOMEM;
+-		return -1;
++	if (s) {
++		memcpy(s, var, l1);
++		s[l1] = '=';
++		memcpy(s+l1+1, value, l2);
++		s[l1+l2+1] = 0;
++		if (!__putenv(s, 1)) return 0;
+ 	}
+-	return 0;
++	free(s);
++	return -1;
+ }
+diff --git a/src/ldso/dynlink.c b/src/ldso/dynlink.c
+index a525b3d..fc58832 100644
+--- a/src/ldso/dynlink.c
++++ b/src/ldso/dynlink.c
+@@ -1354,7 +1354,7 @@ int __dladdr(void *addr, Dl_info *info)
+ 		uint32_t *hashval;
+ 		buckets = p->ghashtab + 4 + (p->ghashtab[2]*sizeof(size_t)/4);
+ 		sym += p->ghashtab[1];
+-		for (i = 0; i < p->ghashtab[0]; i++) {
++		for (i = nsym = 0; i < p->ghashtab[0]; i++) {
+ 			if (buckets[i] > nsym)
+ 				nsym = buckets[i];
+ 		}
+diff --git a/src/malloc/__brk.c b/src/malloc/__brk.c
+index 0b561ea..4c9119b 100644
+--- a/src/malloc/__brk.c
++++ b/src/malloc/__brk.c
+@@ -3,5 +3,5 @@
+ 
+ uintptr_t __brk(uintptr_t newbrk)
+ {
+-	return syscall(SYS_brk, newbrk);
++	return __syscall(SYS_brk, newbrk);
+ }
+diff --git a/src/malloc/malloc.c b/src/malloc/malloc.c
+index fb65ab5..d6ad904 100644
+--- a/src/malloc/malloc.c
++++ b/src/malloc/malloc.c
+@@ -177,6 +177,7 @@ static struct chunk *expand_heap(size_t n)
+ 	return w;
+ fail:
+ 	unlock(mal.brk_lock);
++	errno = ENOMEM;
+ 	return 0;
+ }
+ 
+diff --git a/src/math/__log1p.h b/src/math/__log1p.h
+deleted file mode 100644
+index 5718711..0000000
+--- a/src/math/__log1p.h
++++ /dev/null
+@@ -1,94 +0,0 @@
+-/* origin: FreeBSD /usr/src/lib/msun/src/k_log.h */
+-/*
+- * ====================================================
+- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+- *
+- * Developed at SunSoft, a Sun Microsystems, Inc. business.
+- * Permission to use, copy, modify, and distribute this
+- * software is freely granted, provided that this notice
+- * is preserved.
+- * ====================================================
+- */
+-/*
+- * __log1p(f):
+- * Return log(1+f) - f for 1+f in ~[sqrt(2)/2, sqrt(2)].
+- *
+- * The following describes the overall strategy for computing
+- * logarithms in base e.  The argument reduction and adding the final
+- * term of the polynomial are done by the caller for increased accuracy
+- * when different bases are used.
+- *
+- * Method :
+- *   1. Argument Reduction: find k and f such that
+- *                      x = 2^k * (1+f),
+- *         where  sqrt(2)/2 < 1+f < sqrt(2) .
+- *
+- *   2. Approximation of log(1+f).
+- *      Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+- *               = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+- *               = 2s + s*R
+- *      We use a special Reme algorithm on [0,0.1716] to generate
+- *      a polynomial of degree 14 to approximate R The maximum error
+- *      of this polynomial approximation is bounded by 2**-58.45. In
+- *      other words,
+- *                      2      4      6      8      10      12      14
+- *          R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s
+- *      (the values of Lg1 to Lg7 are listed in the program)
+- *      and
+- *          |      2          14          |     -58.45
+- *          | Lg1*s +...+Lg7*s    -  R(z) | <= 2
+- *          |                             |
+- *      Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+- *      In order to guarantee error in log below 1ulp, we compute log
+- *      by
+- *              log(1+f) = f - s*(f - R)        (if f is not too large)
+- *              log(1+f) = f - (hfsq - s*(hfsq+R)).     (better accuracy)
+- *
+- *      3. Finally,  log(x) = k*ln2 + log(1+f).
+- *                          = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
+- *         Here ln2 is split into two floating point number:
+- *                      ln2_hi + ln2_lo,
+- *         where n*ln2_hi is always exact for |n| < 2000.
+- *
+- * Special cases:
+- *      log(x) is NaN with signal if x < 0 (including -INF) ;
+- *      log(+INF) is +INF; log(0) is -INF with signal;
+- *      log(NaN) is that NaN with no signal.
+- *
+- * Accuracy:
+- *      according to an error analysis, the error is always less than
+- *      1 ulp (unit in the last place).
+- *
+- * Constants:
+- * The hexadecimal values are the intended ones for the following
+- * constants. The decimal values may be used, provided that the
+- * compiler will convert from decimal to binary accurately enough
+- * to produce the hexadecimal values shown.
+- */
+-
+-static const double
+-Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
+-Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
+-Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
+-Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
+-Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
+-Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
+-Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
+-
+-/*
+- * We always inline __log1p(), since doing so produces a
+- * substantial performance improvement (~40% on amd64).
+- */
+-static inline double __log1p(double f)
+-{
+-	double_t hfsq,s,z,R,w,t1,t2;
+-
+-	s = f/(2.0+f);
+-	z = s*s;
+-	w = z*z;
+-	t1= w*(Lg2+w*(Lg4+w*Lg6));
+-	t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
+-	R = t2+t1;
+-	hfsq = 0.5*f*f;
+-	return s*(hfsq+R);
+-}
+diff --git a/src/math/__log1pf.h b/src/math/__log1pf.h
+deleted file mode 100644
+index f2fbef2..0000000
+--- a/src/math/__log1pf.h
++++ /dev/null
+@@ -1,35 +0,0 @@
+-/* origin: FreeBSD /usr/src/lib/msun/src/k_logf.h */
+-/*
+- * ====================================================
+- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+- *
+- * Developed at SunPro, a Sun Microsystems, Inc. business.
+- * Permission to use, copy, modify, and distribute this
+- * software is freely granted, provided that this notice
+- * is preserved.
+- * ====================================================
+- */
+-/*
+- * See comments in __log1p.h.
+- */
+-
+-/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
+-static const float
+-Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */
+-Lg2 = 0xccce13.0p-25, /* 0.40000972152 */
+-Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
+-Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
+-
+-static inline float __log1pf(float f)
+-{
+-	float_t hfsq,s,z,R,w,t1,t2;
+-
+-	s = f/(2.0f + f);
+-	z = s*s;
+-	w = z*z;
+-	t1 = w*(Lg2+w*Lg4);
+-	t2 = z*(Lg1+w*Lg3);
+-	R = t2+t1;
+-	hfsq = 0.5f * f * f;
+-	return s*(hfsq+R);
+-}
+diff --git a/src/math/acosh.c b/src/math/acosh.c
+index 4ce9b3d..badbf90 100644
+--- a/src/math/acosh.c
++++ b/src/math/acosh.c
+@@ -1,5 +1,10 @@
+ #include "libm.h"
+ 
++#if FLT_EVAL_METHOD==2
++#undef sqrt
++#define sqrt sqrtl
++#endif
++
+ /* acosh(x) = log(x + sqrt(x*x-1)) */
+ double acosh(double x)
+ {
+diff --git a/src/math/acoshf.c b/src/math/acoshf.c
+index 16550f1..8a4ec4d 100644
+--- a/src/math/acoshf.c
++++ b/src/math/acoshf.c
+@@ -1,5 +1,13 @@
+ #include "libm.h"
+ 
++#if FLT_EVAL_METHOD==2
++#undef sqrtf
++#define sqrtf sqrtl
++#elif FLT_EVAL_METHOD==1
++#undef sqrtf
++#define sqrtf sqrt
++#endif
++
+ /* acosh(x) = log(x + sqrt(x*x-1)) */
+ float acoshf(float x)
+ {
+diff --git a/src/math/erfl.c b/src/math/erfl.c
+index 42bb1a1..96b74de 100644
+--- a/src/math/erfl.c
++++ b/src/math/erfl.c
+@@ -266,23 +266,12 @@ static long double erfc2(uint32_t ix, long double x)
+ 		     s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
+ 		S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
+ 		     s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
+-	} else {  /* 2.857 <= |x| */
++	} else if (ix < 0x4001d555) {  /* 2.857 <= |x| < 6.6666259765625 */
+ 		R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
+ 		     s * (rb[5] + s * (rb[6] + s * rb[7]))))));
+ 		S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
+ 		     s * (sb[5] + s * (sb[6] + s))))));
+-	}
+-	if (ix < 0x4000b6db) {  /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */
+-		R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
+-		     s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
+-		S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
+-		     s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
+-	} else if (ix < 0x4001d555) {  /* 6.6666259765625 > |x| >= 1/.35 ~ 2.857143 */
+-		R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
+-		     s * (rb[5] + s * (rb[6] + s * rb[7]))))));
+-		S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
+-		     s * (sb[5] + s * (sb[6] + s))))));
+-	} else { /* 107 > |x| >= 6.666 */
++	} else { /* 6.666 <= |x| < 107 (erfc only) */
+ 		R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
+ 		     s * (rc[4] + s * rc[5]))));
+ 		S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
+diff --git a/src/math/fma.c b/src/math/fma.c
+index 84868be..02f5c86 100644
+--- a/src/math/fma.c
++++ b/src/math/fma.c
+@@ -431,12 +431,24 @@ double fma(double x, double y, double z)
+ 		/*
+ 		 * There is no need to worry about double rounding in directed
+ 		 * rounding modes.
+-		 * TODO: underflow is not raised properly, example in downward rounding:
++		 * But underflow may not be raised properly, example in downward rounding:
+ 		 * fma(0x1.000000001p-1000, 0x1.000000001p-30, -0x1p-1066)
+ 		 */
++		double ret;
++#if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
++		int e = fetestexcept(FE_INEXACT);
++		feclearexcept(FE_INEXACT);
++#endif
+ 		fesetround(oround);
+ 		adj = r.lo + xy.lo;
+-		return scalbn(r.hi + adj, spread);
++		ret = scalbn(r.hi + adj, spread);
++#if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
++		if (ilogb(ret) < -1022 && fetestexcept(FE_INEXACT))
++			feraiseexcept(FE_UNDERFLOW);
++		else if (e)
++			feraiseexcept(FE_INEXACT);
++#endif
++		return ret;
+ 	}
+ 
+ 	adj = add_adjusted(r.lo, xy.lo);
+diff --git a/src/math/fmaf.c b/src/math/fmaf.c
+index 745ee39..aa57feb 100644
+--- a/src/math/fmaf.c
++++ b/src/math/fmaf.c
+@@ -26,7 +26,8 @@
+  */
+ 
+ #include <fenv.h>
+-#include "libm.h"
++#include <math.h>
++#include <stdint.h>
+ 
+ /*
+  * Fused multiply-add: Compute x * y + z with a single rounding error.
+@@ -39,21 +40,35 @@ float fmaf(float x, float y, float z)
+ {
+ 	#pragma STDC FENV_ACCESS ON
+ 	double xy, result;
+-	uint32_t hr, lr;
++	union {double f; uint64_t i;} u;
++	int e;
+ 
+ 	xy = (double)x * y;
+ 	result = xy + z;
+-	EXTRACT_WORDS(hr, lr, result);
++	u.f = result;
++	e = u.i>>52 & 0x7ff;
+ 	/* Common case: The double precision result is fine. */
+-	if ((lr & 0x1fffffff) != 0x10000000 ||  /* not a halfway case */
+-		(hr & 0x7ff00000) == 0x7ff00000 ||  /* NaN */
++	if ((u.i & 0x1fffffff) != 0x10000000 || /* not a halfway case */
++		e == 0x7ff ||                   /* NaN */
+ 		result - xy == z ||                 /* exact */
+ 		fegetround() != FE_TONEAREST)       /* not round-to-nearest */
+ 	{
+ 		/*
+-		TODO: underflow is not raised correctly, example in
+-		downward rouding: fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f)
++		underflow may not be raised correctly, example:
++		fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f)
+ 		*/
++#if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
++		if (e < 0x3ff-126 && e >= 0x3ff-149 && fetestexcept(FE_INEXACT)) {
++			feclearexcept(FE_INEXACT);
++			/* TODO: gcc and clang bug workaround */
++			volatile float vz = z;
++			result = xy + vz;
++			if (fetestexcept(FE_INEXACT))
++				feraiseexcept(FE_UNDERFLOW);
++			else
++				feraiseexcept(FE_INEXACT);
++		}
++#endif
+ 		z = result;
+ 		return z;
+ 	}
+@@ -68,8 +83,11 @@ float fmaf(float x, float y, float z)
+ 	volatile double vxy = xy;  /* XXX work around gcc CSE bug */
+ 	double adjusted_result = vxy + z;
+ 	fesetround(FE_TONEAREST);
+-	if (result == adjusted_result)
+-		SET_LOW_WORD(adjusted_result, lr + 1);
++	if (result == adjusted_result) {
++		u.f = adjusted_result;
++		u.i++;
++		adjusted_result = u.f;
++	}
+ 	z = adjusted_result;
+ 	return z;
+ }
+diff --git a/src/math/fmal.c b/src/math/fmal.c
+index c68db25..4506aac 100644
+--- a/src/math/fmal.c
++++ b/src/math/fmal.c
+@@ -264,12 +264,24 @@ long double fmal(long double x, long double y, long double z)
+ 		/*
+ 		 * There is no need to worry about double rounding in directed
+ 		 * rounding modes.
+-		 * TODO: underflow is not raised correctly, example in downward rounding:
++		 * But underflow may not be raised correctly, example in downward rounding:
+ 		 * fmal(0x1.0000000001p-16000L, 0x1.0000000001p-400L, -0x1p-16440L)
+ 		 */
++		long double ret;
++#if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
++		int e = fetestexcept(FE_INEXACT);
++		feclearexcept(FE_INEXACT);
++#endif
+ 		fesetround(oround);
+ 		adj = r.lo + xy.lo;
+-		return scalbnl(r.hi + adj, spread);
++		ret = scalbnl(r.hi + adj, spread);
++#if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
++		if (ilogbl(ret) < -16382 && fetestexcept(FE_INEXACT))
++			feraiseexcept(FE_UNDERFLOW);
++		else if (e)
++			feraiseexcept(FE_INEXACT);
++#endif
++		return ret;
+ 	}
+ 
+ 	adj = add_adjusted(r.lo, xy.lo);
+diff --git a/src/math/lgammal.c b/src/math/lgammal.c
+index cc4895e..58054e5 100644
+--- a/src/math/lgammal.c
++++ b/src/math/lgammal.c
+@@ -341,12 +341,12 @@ long double __lgammal_r(long double x, int *sg) {
+ 	} else if (ix < 0x40028000) {  /* 8.0 */
+ 		/* x < 8.0 */
+ 		i = (int)x;
+-		t = 0.0;
+ 		y = x - (double)i;
+ 		p = y * (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))))));
+ 		q = r0 + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * (r6 + y))))));
+ 		r = 0.5 * y + p / q;
+-		z = 1.0;/* lgamma(1+s) = log(s) + lgamma(s) */
++		z = 1.0;
++		/* lgamma(1+s) = log(s) + lgamma(s) */
+ 		switch (i) {
+ 		case 7:
+ 			z *= (y + 6.0); /* FALLTHRU */
+diff --git a/src/math/log.c b/src/math/log.c
+index 9805120..e61e113 100644
+--- a/src/math/log.c
++++ b/src/math/log.c
+@@ -10,7 +10,7 @@
+  * ====================================================
+  */
+ /* log(x)
+- * Return the logrithm of x
++ * Return the logarithm of x
+  *
+  * Method :
+  *   1. Argument Reduction: find k and f such that
+@@ -60,12 +60,12 @@
+  * to produce the hexadecimal values shown.
+  */
+ 
+-#include "libm.h"
++#include <math.h>
++#include <stdint.h>
+ 
+ static const double
+ ln2_hi = 6.93147180369123816490e-01,  /* 3fe62e42 fee00000 */
+ ln2_lo = 1.90821492927058770002e-10,  /* 3dea39ef 35793c76 */
+-two54  = 1.80143985094819840000e+16,  /* 43500000 00000000 */
+ Lg1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
+ Lg2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
+ Lg3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */
+@@ -76,63 +76,43 @@ Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
+ 
+ double log(double x)
+ {
+-	double hfsq,f,s,z,R,w,t1,t2,dk;
+-	int32_t k,hx,i,j;
+-	uint32_t lx;
+-
+-	EXTRACT_WORDS(hx, lx, x);
++	union {double f; uint64_t i;} u = {x};
++	double_t hfsq,f,s,z,R,w,t1,t2,dk;
++	uint32_t hx;
++	int k;
+ 
++	hx = u.i>>32;
+ 	k = 0;
+-	if (hx < 0x00100000) {  /* x < 2**-1022  */
+-		if (((hx&0x7fffffff)|lx) == 0)
+-			return -two54/0.0;  /* log(+-0)=-inf */
+-		if (hx < 0)
+-			return (x-x)/0.0;   /* log(-#) = NaN */
+-		/* subnormal number, scale up x */
++	if (hx < 0x00100000 || hx>>31) {
++		if (u.i<<1 == 0)
++			return -1/(x*x);  /* log(+-0)=-inf */
++		if (hx>>31)
++			return (x-x)/0.0; /* log(-#) = NaN */
++		/* subnormal number, scale x up */
+ 		k -= 54;
+-		x *= two54;
+-		GET_HIGH_WORD(hx,x);
+-	}
+-	if (hx >= 0x7ff00000)
+-		return x+x;
+-	k += (hx>>20) - 1023;
+-	hx &= 0x000fffff;
+-	i = (hx+0x95f64)&0x100000;
+-	SET_HIGH_WORD(x, hx|(i^0x3ff00000));  /* normalize x or x/2 */
+-	k += i>>20;
++		x *= 0x1p54;
++		u.f = x;
++		hx = u.i>>32;
++	} else if (hx >= 0x7ff00000) {
++		return x;
++	} else if (hx == 0x3ff00000 && u.i<<32 == 0)
++		return 0;
++
++	/* reduce x into [sqrt(2)/2, sqrt(2)] */
++	hx += 0x3ff00000 - 0x3fe6a09e;
++	k += (int)(hx>>20) - 0x3ff;
++	hx = (hx&0x000fffff) + 0x3fe6a09e;
++	u.i = (uint64_t)hx<<32 | (u.i&0xffffffff);
++	x = u.f;
++
+ 	f = x - 1.0;
+-	if ((0x000fffff&(2+hx)) < 3) {  /* -2**-20 <= f < 2**-20 */
+-		if (f == 0.0) {
+-			if (k == 0) {
+-				return 0.0;
+-			}
+-			dk = (double)k;
+-			return dk*ln2_hi + dk*ln2_lo;
+-		}
+-		R = f*f*(0.5-0.33333333333333333*f);
+-		if (k == 0)
+-			return f - R;
+-		dk = (double)k;
+-		return dk*ln2_hi - ((R-dk*ln2_lo)-f);
+-	}
++	hfsq = 0.5*f*f;
+ 	s = f/(2.0+f);
+-	dk = (double)k;
+ 	z = s*s;
+-	i = hx - 0x6147a;
+ 	w = z*z;
+-	j = 0x6b851 - hx;
+ 	t1 = w*(Lg2+w*(Lg4+w*Lg6));
+ 	t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
+-	i |= j;
+ 	R = t2 + t1;
+-	if (i > 0) {
+-		hfsq = 0.5*f*f;
+-		if (k == 0)
+-			return f - (hfsq-s*(hfsq+R));
+-		return dk*ln2_hi - ((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
+-	} else {
+-		if (k == 0)
+-			return f - s*(f-R);
+-		return dk*ln2_hi - ((s*(f-R)-dk*ln2_lo)-f);
+-	}
++	dk = k;
++	return s*(hfsq+R) + dk*ln2_lo - hfsq + f + dk*ln2_hi;
+ }
+diff --git a/src/math/log10.c b/src/math/log10.c
+index ed65d9b..8102687 100644
+--- a/src/math/log10.c
++++ b/src/math/log10.c
+@@ -10,72 +10,91 @@
+  * ====================================================
+  */
+ /*
+- * Return the base 10 logarithm of x.  See e_log.c and k_log.h for most
+- * comments.
++ * Return the base 10 logarithm of x.  See log.c for most comments.
+  *
+- *    log10(x) = (f - 0.5*f*f + k_log1p(f)) / ln10 + k * log10(2)
+- * in not-quite-routine extra precision.
++ * Reduce x to 2^k (1+f) and calculate r = log(1+f) - f + f*f/2
++ * as in log.c, then combine and scale in extra precision:
++ *    log10(x) = (f - f*f/2 + r)/log(10) + k*log10(2)
+  */
+ 
+-#include "libm.h"
+-#include "__log1p.h"
++#include <math.h>
++#include <stdint.h>
+ 
+ static const double
+-two54     = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+ ivln10hi  = 4.34294481878168880939e-01, /* 0x3fdbcb7b, 0x15200000 */
+ ivln10lo  = 2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */
+ log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
+-log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
++log10_2lo = 3.69423907715893078616e-13, /* 0x3D59FEF3, 0x11F12B36 */
++Lg1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
++Lg2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
++Lg3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */
++Lg4 = 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */
++Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
++Lg6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
++Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
+ 
+ double log10(double x)
+ {
+-	double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2;
+-	int32_t i,k,hx;
+-	uint32_t lx;
+-
+-	EXTRACT_WORDS(hx, lx, x);
++	union {double f; uint64_t i;} u = {x};
++	double_t hfsq,f,s,z,R,w,t1,t2,dk,y,hi,lo,val_hi,val_lo;
++	uint32_t hx;
++	int k;
+ 
++	hx = u.i>>32;
+ 	k = 0;
+-	if (hx < 0x00100000) {  /* x < 2**-1022  */
+-		if (((hx&0x7fffffff)|lx) == 0)
+-			return -two54/0.0;  /* log(+-0)=-inf */
+-		if (hx<0)
+-			return (x-x)/0.0;   /* log(-#) = NaN */
+-		/* subnormal number, scale up x */
++	if (hx < 0x00100000 || hx>>31) {
++		if (u.i<<1 == 0)
++			return -1/(x*x);  /* log(+-0)=-inf */
++		if (hx>>31)
++			return (x-x)/0.0; /* log(-#) = NaN */
++		/* subnormal number, scale x up */
+ 		k -= 54;
+-		x *= two54;
+-		GET_HIGH_WORD(hx, x);
+-	}
+-	if (hx >= 0x7ff00000)
+-		return x+x;
+-	if (hx == 0x3ff00000 && lx == 0)
+-		return 0.0;  /* log(1) = +0 */
+-	k += (hx>>20) - 1023;
+-	hx &= 0x000fffff;
+-	i = (hx+0x95f64)&0x100000;
+-	SET_HIGH_WORD(x, hx|(i^0x3ff00000));  /* normalize x or x/2 */
+-	k += i>>20;
+-	y = (double)k;
++		x *= 0x1p54;
++		u.f = x;
++		hx = u.i>>32;
++	} else if (hx >= 0x7ff00000) {
++		return x;
++	} else if (hx == 0x3ff00000 && u.i<<32 == 0)
++		return 0;
++
++	/* reduce x into [sqrt(2)/2, sqrt(2)] */
++	hx += 0x3ff00000 - 0x3fe6a09e;
++	k += (int)(hx>>20) - 0x3ff;
++	hx = (hx&0x000fffff) + 0x3fe6a09e;
++	u.i = (uint64_t)hx<<32 | (u.i&0xffffffff);
++	x = u.f;
++
+ 	f = x - 1.0;
+ 	hfsq = 0.5*f*f;
+-	r = __log1p(f);
++	s = f/(2.0+f);
++	z = s*s;
++	w = z*z;
++	t1 = w*(Lg2+w*(Lg4+w*Lg6));
++	t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
++	R = t2 + t1;
+ 
+ 	/* See log2.c for details. */
++	/* hi+lo = f - hfsq + s*(hfsq+R) ~ log(1+f) */
+ 	hi = f - hfsq;
+-	SET_LOW_WORD(hi, 0);
+-	lo = (f - hi) - hfsq + r;
++	u.f = hi;
++	u.i &= (uint64_t)-1<<32;
++	hi = u.f;
++	lo = f - hi - hfsq + s*(hfsq+R);
++
++	/* val_hi+val_lo ~ log10(1+f) + k*log10(2) */
+ 	val_hi = hi*ivln10hi;
+-	y2 = y*log10_2hi;
+-	val_lo = y*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi;
++	dk = k;
++	y = dk*log10_2hi;
++	val_lo = dk*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi;
+ 
+ 	/*
+-	 * Extra precision in for adding y*log10_2hi is not strictly needed
++	 * Extra precision in for adding y is not strictly needed
+ 	 * since there is no very large cancellation near x = sqrt(2) or
+ 	 * x = 1/sqrt(2), but we do it anyway since it costs little on CPUs
+ 	 * with some parallelism and it reduces the error for many args.
+ 	 */
+-	w = y2 + val_hi;
+-	val_lo += (y2 - w) + val_hi;
++	w = y + val_hi;
++	val_lo += (y - w) + val_hi;
+ 	val_hi = w;
+ 
+ 	return val_lo + val_hi;
+diff --git a/src/math/log10f.c b/src/math/log10f.c
+index e10749b..9ca2f01 100644
+--- a/src/math/log10f.c
++++ b/src/math/log10f.c
+@@ -13,57 +13,65 @@
+  * See comments in log10.c.
+  */
+ 
+-#include "libm.h"
+-#include "__log1pf.h"
++#include <math.h>
++#include <stdint.h>
+ 
+ static const float
+-two25     =  3.3554432000e+07, /* 0x4c000000 */
+ ivln10hi  =  4.3432617188e-01, /* 0x3ede6000 */
+ ivln10lo  = -3.1689971365e-05, /* 0xb804ead9 */
+ log10_2hi =  3.0102920532e-01, /* 0x3e9a2080 */
+-log10_2lo =  7.9034151668e-07; /* 0x355427db */
++log10_2lo =  7.9034151668e-07, /* 0x355427db */
++/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
++Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */
++Lg2 = 0xccce13.0p-25, /* 0.40000972152 */
++Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
++Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
+ 
+ float log10f(float x)
+ {
+-	float f,hfsq,hi,lo,r,y;
+-	int32_t i,k,hx;
+-
+-	GET_FLOAT_WORD(hx, x);
++	union {float f; uint32_t i;} u = {x};
++	float_t hfsq,f,s,z,R,w,t1,t2,dk,hi,lo;
++	uint32_t ix;
++	int k;
+ 
++	ix = u.i;
+ 	k = 0;
+-	if (hx < 0x00800000) {  /* x < 2**-126  */
+-		if ((hx&0x7fffffff) == 0)
+-			return -two25/0.0f;  /* log(+-0)=-inf */
+-		if (hx < 0)
+-			return (x-x)/0.0f;   /* log(-#) = NaN */
++	if (ix < 0x00800000 || ix>>31) {  /* x < 2**-126  */
++		if (ix<<1 == 0)
++			return -1/(x*x);  /* log(+-0)=-inf */
++		if (ix>>31)
++			return (x-x)/0.0f; /* log(-#) = NaN */
+ 		/* subnormal number, scale up x */
+ 		k -= 25;
+-		x *= two25;
+-		GET_FLOAT_WORD(hx, x);
+-	}
+-	if (hx >= 0x7f800000)
+-		return x+x;
+-	if (hx == 0x3f800000)
+-		return 0.0f;  /* log(1) = +0 */
+-	k += (hx>>23) - 127;
+-	hx &= 0x007fffff;
+-	i = (hx+(0x4afb0d))&0x800000;
+-	SET_FLOAT_WORD(x, hx|(i^0x3f800000));  /* normalize x or x/2 */
+-	k += i>>23;
+-	y = (float)k;
++		x *= 0x1p25f;
++		u.f = x;
++		ix = u.i;
++	} else if (ix >= 0x7f800000) {
++		return x;
++	} else if (ix == 0x3f800000)
++		return 0;
++
++	/* reduce x into [sqrt(2)/2, sqrt(2)] */
++	ix += 0x3f800000 - 0x3f3504f3;
++	k += (int)(ix>>23) - 0x7f;
++	ix = (ix&0x007fffff) + 0x3f3504f3;
++	u.i = ix;
++	x = u.f;
++
+ 	f = x - 1.0f;
+-	hfsq = 0.5f * f * f;
+-	r = __log1pf(f);
++	s = f/(2.0f + f);
++	z = s*s;
++	w = z*z;
++	t1= w*(Lg2+w*Lg4);
++	t2= z*(Lg1+w*Lg3);
++	R = t2 + t1;
++	hfsq = 0.5f*f*f;
+ 
+-// FIXME
+-//      /* See log2f.c and log2.c for details. */
+-//      if (sizeof(float_t) > sizeof(float))
+-//              return (r - hfsq + f) * ((float_t)ivln10lo + ivln10hi) +
+-//                  y * ((float_t)log10_2lo + log10_2hi);
+ 	hi = f - hfsq;
+-	GET_FLOAT_WORD(hx, hi);
+-	SET_FLOAT_WORD(hi, hx&0xfffff000);
+-	lo = (f - hi) - hfsq + r;
+-	return y*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi +
+-	        hi*ivln10hi + y*log10_2hi;
++	u.f = hi;
++	u.i &= 0xfffff000;
++	hi = u.f;
++	lo = f - hi - hfsq + s*(hfsq+R);
++	dk = k;
++	return dk*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi + hi*ivln10hi + dk*log10_2hi;
+ }
+diff --git a/src/math/log10l.c b/src/math/log10l.c
+index f0eeeaf..c7aacf9 100644
+--- a/src/math/log10l.c
++++ b/src/math/log10l.c
+@@ -117,16 +117,15 @@ static const long double S[4] = {
+ 
+ long double log10l(long double x)
+ {
+-	long double y;
+-	volatile long double z;
++	long double y, z;
+ 	int e;
+ 
+ 	if (isnan(x))
+ 		return x;
+ 	if(x <= 0.0) {
+ 		if(x == 0.0)
+-			return -1.0 / (x - x);
+-		return (x - x) / (x - x);
++			return -1.0 / (x*x);
++		return (x - x) / 0.0;
+ 	}
+ 	if (x == INFINITY)
+ 		return INFINITY;
+diff --git a/src/math/log1p.c b/src/math/log1p.c
+index a71ac42..0097134 100644
+--- a/src/math/log1p.c
++++ b/src/math/log1p.c
+@@ -10,6 +10,7 @@
+  * ====================================================
+  */
+ /* double log1p(double x)
++ * Return the natural logarithm of 1+x.
+  *
+  * Method :
+  *   1. Argument Reduction: find k and f such that
+@@ -23,31 +24,9 @@
+  *      and add back the correction term c/u.
+  *      (Note: when x > 2**53, one can simply return log(x))
+  *
+- *   2. Approximation of log1p(f).
+- *      Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+- *               = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+- *               = 2s + s*R
+- *      We use a special Reme algorithm on [0,0.1716] to generate
+- *      a polynomial of degree 14 to approximate R The maximum error
+- *      of this polynomial approximation is bounded by 2**-58.45. In
+- *      other words,
+- *                      2      4      6      8      10      12      14
+- *          R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s  +Lp6*s  +Lp7*s
+- *      (the values of Lp1 to Lp7 are listed in the program)
+- *      and
+- *          |      2          14          |     -58.45
+- *          | Lp1*s +...+Lp7*s    -  R(z) | <= 2
+- *          |                             |
+- *      Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+- *      In order to guarantee error in log below 1ulp, we compute log
+- *      by
+- *              log1p(f) = f - (hfsq - s*(hfsq+R)).
++ *   2. Approximation of log(1+f): See log.c
+  *
+- *      3. Finally, log1p(x) = k*ln2 + log1p(f).
+- *                           = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
+- *         Here ln2 is split into two floating point number:
+- *                      ln2_hi + ln2_lo,
+- *         where n*ln2_hi is always exact for |n| < 2000.
++ *   3. Finally, log1p(x) = k*ln2 + log(1+f) + c/u. See log.c
+  *
+  * Special cases:
+  *      log1p(x) is NaN with signal if x < -1 (including -INF) ;
+@@ -79,94 +58,65 @@
+ static const double
+ ln2_hi = 6.93147180369123816490e-01,  /* 3fe62e42 fee00000 */
+ ln2_lo = 1.90821492927058770002e-10,  /* 3dea39ef 35793c76 */
+-two54  = 1.80143985094819840000e+16,  /* 43500000 00000000 */
+-Lp1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
+-Lp2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
+-Lp3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */
+-Lp4 = 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */
+-Lp5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
+-Lp6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
+-Lp7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
++Lg1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
++Lg2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
++Lg3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */
++Lg4 = 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */
++Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
++Lg6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
++Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
+ 
+ double log1p(double x)
+ {
+-	double hfsq,f,c,s,z,R,u;
+-	int32_t k,hx,hu,ax;
+-
+-	GET_HIGH_WORD(hx, x);
+-	ax = hx & 0x7fffffff;
++	union {double f; uint64_t i;} u = {x};
++	double_t hfsq,f,c,s,z,R,w,t1,t2,dk;
++	uint32_t hx,hu;
++	int k;
+ 
++	hx = u.i>>32;
+ 	k = 1;
+-	if (hx < 0x3FDA827A) {  /* 1+x < sqrt(2)+ */
+-		if (ax >= 0x3ff00000) {  /* x <= -1.0 */
+-			if (x == -1.0)
+-				return -two54/0.0; /* log1p(-1)=+inf */
+-			return (x-x)/(x-x);         /* log1p(x<-1)=NaN */
++	if (hx < 0x3fda827a || hx>>31) {  /* 1+x < sqrt(2)+ */
++		if (hx >= 0xbff00000) {  /* x <= -1.0 */
++			if (x == -1)
++				return x/0.0; /* log1p(-1) = -inf */
++			return (x-x)/0.0;     /* log1p(x<-1) = NaN */
+ 		}
+-		if (ax < 0x3e200000) {   /* |x| < 2**-29 */
+-			/* if 0x1p-1022 <= |x| < 0x1p-54, avoid raising underflow */
+-			if (ax < 0x3c900000 && ax >= 0x00100000)
+-				return x;
+-#if FLT_EVAL_METHOD != 0
+-			FORCE_EVAL((float)x);
+-#endif
+-			return x - x*x*0.5;
++		if (hx<<1 < 0x3ca00000<<1) {  /* |x| < 2**-53 */
++			/* underflow if subnormal */
++			if ((hx&0x7ff00000) == 0)
++				FORCE_EVAL((float)x);
++			return x;
+ 		}
+-		if (hx > 0 || hx <= (int32_t)0xbfd2bec4) {  /* sqrt(2)/2- <= 1+x < sqrt(2)+ */
++		if (hx <= 0xbfd2bec4) {  /* sqrt(2)/2- <= 1+x < sqrt(2)+ */
+ 			k = 0;
++			c = 0;
+ 			f = x;
+-			hu = 1;
+ 		}
+-	}
+-	if (hx >= 0x7ff00000)
+-		return x+x;
+-	if (k != 0) {
+-		if (hx < 0x43400000) {
+-			u = 1 + x;
+-			GET_HIGH_WORD(hu, u);
+-			k = (hu>>20) - 1023;
+-			c = k > 0 ? 1.0-(u-x) : x-(u-1.0); /* correction term */
+-			c /= u;
+-		} else {
+-			u = x;
+-			GET_HIGH_WORD(hu,u);
+-			k = (hu>>20) - 1023;
++	} else if (hx >= 0x7ff00000)
++		return x;
++	if (k) {
++		u.f = 1 + x;
++		hu = u.i>>32;
++		hu += 0x3ff00000 - 0x3fe6a09e;
++		k = (int)(hu>>20) - 0x3ff;
++		/* correction term ~ log(1+x)-log(u), avoid underflow in c/u */
++		if (k < 54) {
++			c = k >= 2 ? 1-(u.f-x) : x-(u.f-1);
++			c /= u.f;
++		} else
+ 			c = 0;
+-		}
+-		hu &= 0x000fffff;
+-		/*
+-		 * The approximation to sqrt(2) used in thresholds is not
+-		 * critical.  However, the ones used above must give less
+-		 * strict bounds than the one here so that the k==0 case is
+-		 * never reached from here, since here we have committed to
+-		 * using the correction term but don't use it if k==0.
+-		 */
+-		if (hu < 0x6a09e) {  /* u ~< sqrt(2) */
+-			SET_HIGH_WORD(u, hu|0x3ff00000); /* normalize u */
+-		} else {
+-			k += 1;
+-			SET_HIGH_WORD(u, hu|0x3fe00000); /* normalize u/2 */
+-			hu = (0x00100000-hu)>>2;
+-		}
+-		f = u - 1.0;
++		/* reduce u into [sqrt(2)/2, sqrt(2)] */
++		hu = (hu&0x000fffff) + 0x3fe6a09e;
++		u.i = (uint64_t)hu<<32 | (u.i&0xffffffff);
++		f = u.f - 1;
+ 	}
+ 	hfsq = 0.5*f*f;
+-	if (hu == 0) {   /* |f| < 2**-20 */
+-		if (f == 0.0) {
+-			if(k == 0)
+-				return 0.0;
+-			c += k*ln2_lo;
+-			return k*ln2_hi + c;
+-		}
+-		R = hfsq*(1.0 - 0.66666666666666666*f);
+-		if (k == 0)
+-			return f - R;
+-		return k*ln2_hi - ((R-(k*ln2_lo+c))-f);
+-	}
+ 	s = f/(2.0+f);
+ 	z = s*s;
+-	R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))));
+-	if (k == 0)
+-		return f - (hfsq-s*(hfsq+R));
+-	return k*ln2_hi - ((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
++	w = z*z;
++	t1 = w*(Lg2+w*(Lg4+w*Lg6));
++	t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
++	R = t2 + t1;
++	dk = k;
++	return s*(hfsq+R) + (dk*ln2_lo+c) - hfsq + f + dk*ln2_hi;
+ }
+diff --git a/src/math/log1pf.c b/src/math/log1pf.c
+index e6940d2..23985c3 100644
+--- a/src/math/log1pf.c
++++ b/src/math/log1pf.c
+@@ -1,8 +1,5 @@
+ /* origin: FreeBSD /usr/src/lib/msun/src/s_log1pf.c */
+ /*
+- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+- */
+-/*
+  * ====================================================
+  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+  *
+@@ -18,95 +15,63 @@
+ static const float
+ ln2_hi = 6.9313812256e-01, /* 0x3f317180 */
+ ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */
+-two25  = 3.355443200e+07,  /* 0x4c000000 */
+-Lp1 = 6.6666668653e-01, /* 3F2AAAAB */
+-Lp2 = 4.0000000596e-01, /* 3ECCCCCD */
+-Lp3 = 2.8571429849e-01, /* 3E924925 */
+-Lp4 = 2.2222198546e-01, /* 3E638E29 */
+-Lp5 = 1.8183572590e-01, /* 3E3A3325 */
+-Lp6 = 1.5313838422e-01, /* 3E1CD04F */
+-Lp7 = 1.4798198640e-01; /* 3E178897 */
++/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
++Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */
++Lg2 = 0xccce13.0p-25, /* 0.40000972152 */
++Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
++Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
+ 
+ float log1pf(float x)
+ {
+-	float hfsq,f,c,s,z,R,u;
+-	int32_t k,hx,hu,ax;
+-
+-	GET_FLOAT_WORD(hx, x);
+-	ax = hx & 0x7fffffff;
++	union {float f; uint32_t i;} u = {x};
++	float_t hfsq,f,c,s,z,R,w,t1,t2,dk;
++	uint32_t ix,iu;
++	int k;
+ 
++	ix = u.i;
+ 	k = 1;
+-	if (hx < 0x3ed413d0) {  /* 1+x < sqrt(2)+  */
+-		if (ax >= 0x3f800000) {  /* x <= -1.0 */
+-			if (x == -1.0f)
+-				return -two25/0.0f; /* log1p(-1)=+inf */
+-			return (x-x)/(x-x);         /* log1p(x<-1)=NaN */
++	if (ix < 0x3ed413d0 || ix>>31) {  /* 1+x < sqrt(2)+  */
++		if (ix >= 0xbf800000) {  /* x <= -1.0 */
++			if (x == -1)
++				return x/0.0f; /* log1p(-1)=+inf */
++			return (x-x)/0.0f;     /* log1p(x<-1)=NaN */
+ 		}
+-		if (ax < 0x38000000) {   /* |x| < 2**-15 */
+-			/* if 0x1p-126 <= |x| < 0x1p-24, avoid raising underflow */
+-			if (ax < 0x33800000 && ax >= 0x00800000)
+-				return x;
+-#if FLT_EVAL_METHOD != 0
+-			FORCE_EVAL(x*x);
+-#endif
+-			return x - x*x*0.5f;
++		if (ix<<1 < 0x33800000<<1) {   /* |x| < 2**-24 */
++			/* underflow if subnormal */
++			if ((ix&0x7f800000) == 0)
++				FORCE_EVAL(x*x);
++			return x;
+ 		}
+-		if (hx > 0 || hx <= (int32_t)0xbe95f619) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */
++		if (ix <= 0xbe95f619) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */
+ 			k = 0;
++			c = 0;
+ 			f = x;
+-			hu = 1;
+ 		}
+-	}
+-	if (hx >= 0x7f800000)
+-		return x+x;
+-	if (k != 0) {
+-		if (hx < 0x5a000000) {
+-			u = 1 + x;
+-			GET_FLOAT_WORD(hu, u);
+-			k = (hu>>23) - 127;
+-			/* correction term */
+-			c = k > 0 ? 1.0f-(u-x) : x-(u-1.0f);
+-			c /= u;
+-		} else {
+-			u = x;
+-			GET_FLOAT_WORD(hu,u);
+-			k = (hu>>23) - 127;
++	} else if (ix >= 0x7f800000)
++		return x;
++	if (k) {
++		u.f = 1 + x;
++		iu = u.i;
++		iu += 0x3f800000 - 0x3f3504f3;
++		k = (int)(iu>>23) - 0x7f;
++		/* correction term ~ log(1+x)-log(u), avoid underflow in c/u */
++		if (k < 25) {
++			c = k >= 2 ? 1-(u.f-x) : x-(u.f-1);
++			c /= u.f;
++		} else
+ 			c = 0;
+-		}
+-		hu &= 0x007fffff;
+-		/*
+-		 * The approximation to sqrt(2) used in thresholds is not
+-		 * critical.  However, the ones used above must give less
+-		 * strict bounds than the one here so that the k==0 case is
+-		 * never reached from here, since here we have committed to
+-		 * using the correction term but don't use it if k==0.
+-		 */
+-		if (hu < 0x3504f4) {  /* u < sqrt(2) */
+-			SET_FLOAT_WORD(u, hu|0x3f800000);  /* normalize u */
+-		} else {
+-			k += 1;
+-			SET_FLOAT_WORD(u, hu|0x3f000000);  /* normalize u/2 */
+-			hu = (0x00800000-hu)>>2;
+-		}
+-		f = u - 1.0f;
+-	}
+-	hfsq = 0.5f * f * f;
+-	if (hu == 0) {  /* |f| < 2**-20 */
+-		if (f == 0.0f) {
+-			if (k == 0)
+-				return 0.0f;
+-			c += k*ln2_lo;
+-			return k*ln2_hi+c;
+-		}
+-		R = hfsq*(1.0f - 0.66666666666666666f * f);
+-		if (k == 0)
+-			return f - R;
+-		return k*ln2_hi - ((R-(k*ln2_lo+c))-f);
++		/* reduce u into [sqrt(2)/2, sqrt(2)] */
++		iu = (iu&0x007fffff) + 0x3f3504f3;
++		u.i = iu;
++		f = u.f - 1;
+ 	}
+ 	s = f/(2.0f + f);
+ 	z = s*s;
+-	R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))));
+-	if (k == 0)
+-		return f - (hfsq-s*(hfsq+R));
+-	return k*ln2_hi - ((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
++	w = z*z;
++	t1= w*(Lg2+w*Lg4);
++	t2= z*(Lg1+w*Lg3);
++	R = t2 + t1;
++	hfsq = 0.5f*f*f;
++	dk = k;
++	return s*(hfsq+R) + (dk*ln2_lo+c) - hfsq + f + dk*ln2_hi;
+ }
+diff --git a/src/math/log1pl.c b/src/math/log1pl.c
+index edb48df..37da46d 100644
+--- a/src/math/log1pl.c
++++ b/src/math/log1pl.c
+@@ -118,7 +118,7 @@ long double log1pl(long double xm1)
+ 	/* Test for domain errors.  */
+ 	if (x <= 0.0) {
+ 		if (x == 0.0)
+-			return -1/x; /* -inf with divbyzero */
++			return -1/(x*x); /* -inf with divbyzero */
+ 		return 0/0.0f; /* nan with invalid */
+ 	}
+ 
+diff --git a/src/math/log2.c b/src/math/log2.c
+index 1974215..0aafad4 100644
+--- a/src/math/log2.c
++++ b/src/math/log2.c
+@@ -10,55 +10,66 @@
+  * ====================================================
+  */
+ /*
+- * Return the base 2 logarithm of x.  See log.c and __log1p.h for most
+- * comments.
++ * Return the base 2 logarithm of x.  See log.c for most comments.
+  *
+- * This reduces x to {k, 1+f} exactly as in e_log.c, then calls the kernel,
+- * then does the combining and scaling steps
+- *    log2(x) = (f - 0.5*f*f + k_log1p(f)) / ln2 + k
+- * in not-quite-routine extra precision.
++ * Reduce x to 2^k (1+f) and calculate r = log(1+f) - f + f*f/2
++ * as in log.c, then combine and scale in extra precision:
++ *    log2(x) = (f - f*f/2 + r)/log(2) + k
+  */
+ 
+-#include "libm.h"
+-#include "__log1p.h"
++#include <math.h>
++#include <stdint.h>
+ 
+ static const double
+-two54   = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+ ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
+-ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
++ivln2lo = 1.67517131648865118353e-10, /* 0x3de705fc, 0x2eefa200 */
++Lg1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
++Lg2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
++Lg3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */
++Lg4 = 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */
++Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
++Lg6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
++Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
+ 
+ double log2(double x)
+ {
+-	double f,hfsq,hi,lo,r,val_hi,val_lo,w,y;
+-	int32_t i,k,hx;
+-	uint32_t lx;
+-
+-	EXTRACT_WORDS(hx, lx, x);
++	union {double f; uint64_t i;} u = {x};
++	double_t hfsq,f,s,z,R,w,t1,t2,y,hi,lo,val_hi,val_lo;
++	uint32_t hx;
++	int k;
+ 
++	hx = u.i>>32;
+ 	k = 0;
+-	if (hx < 0x00100000) {  /* x < 2**-1022  */
+-		if (((hx&0x7fffffff)|lx) == 0)
+-			return -two54/0.0;  /* log(+-0)=-inf */
+-		if (hx < 0)
+-			return (x-x)/0.0;   /* log(-#) = NaN */
+-		/* subnormal number, scale up x */
++	if (hx < 0x00100000 || hx>>31) {
++		if (u.i<<1 == 0)
++			return -1/(x*x);  /* log(+-0)=-inf */
++		if (hx>>31)
++			return (x-x)/0.0; /* log(-#) = NaN */
++		/* subnormal number, scale x up */
+ 		k -= 54;
+-		x *= two54;
+-		GET_HIGH_WORD(hx, x);
+-	}
+-	if (hx >= 0x7ff00000)
+-		return x+x;
+-	if (hx == 0x3ff00000 && lx == 0)
+-		return 0.0;  /* log(1) = +0 */
+-	k += (hx>>20) - 1023;
+-	hx &= 0x000fffff;
+-	i = (hx+0x95f64) & 0x100000;
+-	SET_HIGH_WORD(x, hx|(i^0x3ff00000));  /* normalize x or x/2 */
+-	k += i>>20;
+-	y = (double)k;
++		x *= 0x1p54;
++		u.f = x;
++		hx = u.i>>32;
++	} else if (hx >= 0x7ff00000) {
++		return x;
++	} else if (hx == 0x3ff00000 && u.i<<32 == 0)
++		return 0;
++
++	/* reduce x into [sqrt(2)/2, sqrt(2)] */
++	hx += 0x3ff00000 - 0x3fe6a09e;
++	k += (int)(hx>>20) - 0x3ff;
++	hx = (hx&0x000fffff) + 0x3fe6a09e;
++	u.i = (uint64_t)hx<<32 | (u.i&0xffffffff);
++	x = u.f;
++
+ 	f = x - 1.0;
+ 	hfsq = 0.5*f*f;
+-	r = __log1p(f);
++	s = f/(2.0+f);
++	z = s*s;
++	w = z*z;
++	t1 = w*(Lg2+w*(Lg4+w*Lg6));
++	t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
++	R = t2 + t1;
+ 
+ 	/*
+ 	 * f-hfsq must (for args near 1) be evaluated in extra precision
+@@ -90,13 +101,19 @@ double log2(double x)
+ 	 * The multi-precision calculations for the multiplications are
+ 	 * routine.
+ 	 */
++
++	/* hi+lo = f - hfsq + s*(hfsq+R) ~ log(1+f) */
+ 	hi = f - hfsq;
+-	SET_LOW_WORD(hi, 0);
+-	lo = (f - hi) - hfsq + r;
++	u.f = hi;
++	u.i &= (uint64_t)-1<<32;
++	hi = u.f;
++	lo = f - hi - hfsq + s*(hfsq+R);
++
+ 	val_hi = hi*ivln2hi;
+ 	val_lo = (lo+hi)*ivln2lo + lo*ivln2hi;
+ 
+ 	/* spadd(val_hi, val_lo, y), except for not using double_t: */
++	y = k;
+ 	w = y + val_hi;
+ 	val_lo += (y - w) + val_hi;
+ 	val_hi = w;
+diff --git a/src/math/log2f.c b/src/math/log2f.c
+index e0d6a9e..b3e305f 100644
+--- a/src/math/log2f.c
++++ b/src/math/log2f.c
+@@ -13,67 +13,62 @@
+  * See comments in log2.c.
+  */
+ 
+-#include "libm.h"
+-#include "__log1pf.h"
++#include <math.h>
++#include <stdint.h>
+ 
+ static const float
+-two25   =  3.3554432000e+07, /* 0x4c000000 */
+ ivln2hi =  1.4428710938e+00, /* 0x3fb8b000 */
+-ivln2lo = -1.7605285393e-04; /* 0xb9389ad4 */
++ivln2lo = -1.7605285393e-04, /* 0xb9389ad4 */
++/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
++Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */
++Lg2 = 0xccce13.0p-25, /* 0.40000972152 */
++Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
++Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
+ 
+ float log2f(float x)
+ {
+-	float f,hfsq,hi,lo,r,y;
+-	int32_t i,k,hx;
+-
+-	GET_FLOAT_WORD(hx, x);
++	union {float f; uint32_t i;} u = {x};
++	float_t hfsq,f,s,z,R,w,t1,t2,hi,lo;
++	uint32_t ix;
++	int k;
+ 
++	ix = u.i;
+ 	k = 0;
+-	if (hx < 0x00800000) {  /* x < 2**-126  */
+-		if ((hx&0x7fffffff) == 0)
+-			return -two25/0.0f;  /* log(+-0)=-inf */
+-		if (hx < 0)
+-			return (x-x)/0.0f;   /* log(-#) = NaN */
++	if (ix < 0x00800000 || ix>>31) {  /* x < 2**-126  */
++		if (ix<<1 == 0)
++			return -1/(x*x);  /* log(+-0)=-inf */
++		if (ix>>31)
++			return (x-x)/0.0f; /* log(-#) = NaN */
+ 		/* subnormal number, scale up x */
+ 		k -= 25;
+-		x *= two25;
+-		GET_FLOAT_WORD(hx, x);
+-	}
+-	if (hx >= 0x7f800000)
+-		return x+x;
+-	if (hx == 0x3f800000)
+-		return 0.0f;  /* log(1) = +0 */
+-	k += (hx>>23) - 127;
+-	hx &= 0x007fffff;
+-	i = (hx+(0x4afb0d))&0x800000;
+-	SET_FLOAT_WORD(x, hx|(i^0x3f800000));  /* normalize x or x/2 */
+-	k += i>>23;
+-	y = (float)k;
+-	f = x - 1.0f;
+-	hfsq = 0.5f * f * f;
+-	r = __log1pf(f);
++		x *= 0x1p25f;
++		u.f = x;
++		ix = u.i;
++	} else if (ix >= 0x7f800000) {
++		return x;
++	} else if (ix == 0x3f800000)
++		return 0;
+ 
+-	/*
+-	 * We no longer need to avoid falling into the multi-precision
+-	 * calculations due to compiler bugs breaking Dekker's theorem.
+-	 * Keep avoiding this as an optimization.  See log2.c for more
+-	 * details (some details are here only because the optimization
+-	 * is not yet available in double precision).
+-	 *
+-	 * Another compiler bug turned up.  With gcc on i386,
+-	 * (ivln2lo + ivln2hi) would be evaluated in float precision
+-	 * despite runtime evaluations using double precision.  So we
+-	 * must cast one of its terms to float_t.  This makes the whole
+-	 * expression have type float_t, so return is forced to waste
+-	 * time clobbering its extra precision.
+-	 */
+-// FIXME
+-//      if (sizeof(float_t) > sizeof(float))
+-//              return (r - hfsq + f) * ((float_t)ivln2lo + ivln2hi) + y;
++	/* reduce x into [sqrt(2)/2, sqrt(2)] */
++	ix += 0x3f800000 - 0x3f3504f3;
++	k += (int)(ix>>23) - 0x7f;
++	ix = (ix&0x007fffff) + 0x3f3504f3;
++	u.i = ix;
++	x = u.f;
++
++	f = x - 1.0f;
++	s = f/(2.0f + f);
++	z = s*s;
++	w = z*z;
++	t1= w*(Lg2+w*Lg4);
++	t2= z*(Lg1+w*Lg3);
++	R = t2 + t1;
++	hfsq = 0.5f*f*f;
+ 
+ 	hi = f - hfsq;
+-	GET_FLOAT_WORD(hx,hi);
+-	SET_FLOAT_WORD(hi,hx&0xfffff000);
+-	lo = (f - hi) - hfsq + r;
+-	return (lo+hi)*ivln2lo + lo*ivln2hi + hi*ivln2hi + y;
++	u.f = hi;
++	u.i &= 0xfffff000;
++	hi = u.f;
++	lo = f - hi - hfsq + s*(hfsq+R);
++	return (lo+hi)*ivln2lo + lo*ivln2hi + hi*ivln2hi + k;
+ }
+diff --git a/src/math/log2l.c b/src/math/log2l.c
+index 345b395..d00531d 100644
+--- a/src/math/log2l.c
++++ b/src/math/log2l.c
+@@ -117,7 +117,7 @@ long double log2l(long double x)
+ 		return x;
+ 	if (x <= 0.0) {
+ 		if (x == 0.0)
+-			return -1/(x+0); /* -inf with divbyzero */
++			return -1/(x*x); /* -inf with divbyzero */
+ 		return 0/0.0f; /* nan with invalid */
+ 	}
+ 
+diff --git a/src/math/logf.c b/src/math/logf.c
+index c7f7dbe..52230a1 100644
+--- a/src/math/logf.c
++++ b/src/math/logf.c
+@@ -13,12 +13,12 @@
+  * ====================================================
+  */
+ 
+-#include "libm.h"
++#include <math.h>
++#include <stdint.h>
+ 
+ static const float
+ ln2_hi = 6.9313812256e-01, /* 0x3f317180 */
+ ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */
+-two25  = 3.355443200e+07,  /* 0x4c000000 */
+ /* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
+ Lg1 = 0xaaaaaa.0p-24, /* 0.66666662693 */
+ Lg2 = 0xccce13.0p-25, /* 0.40000972152 */
+@@ -27,61 +27,43 @@ Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
+ 
+ float logf(float x)
+ {
+-	float hfsq,f,s,z,R,w,t1,t2,dk;
+-	int32_t k,ix,i,j;
+-
+-	GET_FLOAT_WORD(ix, x);
++	union {float f; uint32_t i;} u = {x};
++	float_t hfsq,f,s,z,R,w,t1,t2,dk;
++	uint32_t ix;
++	int k;
+ 
++	ix = u.i;
+ 	k = 0;
+-	if (ix < 0x00800000) {  /* x < 2**-126  */
+-		if ((ix & 0x7fffffff) == 0)
+-			return -two25/0.0f;  /* log(+-0)=-inf */
+-		if (ix < 0)
+-			return (x-x)/0.0f;   /* log(-#) = NaN */
++	if (ix < 0x00800000 || ix>>31) {  /* x < 2**-126  */
++		if (ix<<1 == 0)
++			return -1/(x*x);  /* log(+-0)=-inf */
++		if (ix>>31)
++			return (x-x)/0.0f; /* log(-#) = NaN */
+ 		/* subnormal number, scale up x */
+ 		k -= 25;
+-		x *= two25;
+-		GET_FLOAT_WORD(ix, x);
+-	}
+-	if (ix >= 0x7f800000)
+-		return x+x;
+-	k += (ix>>23) - 127;
+-	ix &= 0x007fffff;
+-	i = (ix + (0x95f64<<3)) & 0x800000;
+-	SET_FLOAT_WORD(x, ix|(i^0x3f800000));  /* normalize x or x/2 */
+-	k += i>>23;
++		x *= 0x1p25f;
++		u.f = x;
++		ix = u.i;
++	} else if (ix >= 0x7f800000) {
++		return x;
++	} else if (ix == 0x3f800000)
++		return 0;
++
++	/* reduce x into [sqrt(2)/2, sqrt(2)] */
++	ix += 0x3f800000 - 0x3f3504f3;
++	k += (int)(ix>>23) - 0x7f;
++	ix = (ix&0x007fffff) + 0x3f3504f3;
++	u.i = ix;
++	x = u.f;
++
+ 	f = x - 1.0f;
+-	if ((0x007fffff & (0x8000 + ix)) < 0xc000) {  /* -2**-9 <= f < 2**-9 */
+-		if (f == 0.0f) {
+-			if (k == 0)
+-				return 0.0f;
+-			dk = (float)k;
+-			return dk*ln2_hi + dk*ln2_lo;
+-		}
+-		R = f*f*(0.5f - 0.33333333333333333f*f);
+-		if (k == 0)
+-			return f-R;
+-		dk = (float)k;
+-		return dk*ln2_hi - ((R-dk*ln2_lo)-f);
+-	}
+ 	s = f/(2.0f + f);
+-	dk = (float)k;
+ 	z = s*s;
+-	i = ix-(0x6147a<<3);
+ 	w = z*z;
+-	j = (0x6b851<<3)-ix;
+ 	t1= w*(Lg2+w*Lg4);
+ 	t2= z*(Lg1+w*Lg3);
+-	i |= j;
+ 	R = t2 + t1;
+-	if (i > 0) {
+-		hfsq = 0.5f * f * f;
+-		if (k == 0)
+-			return f - (hfsq-s*(hfsq+R));
+-		return dk*ln2_hi - ((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
+-	} else {
+-		if (k == 0)
+-			return f - s*(f-R);
+-		return dk*ln2_hi - ((s*(f-R)-dk*ln2_lo)-f);
+-	}
++	hfsq = 0.5f*f*f;
++	dk = k;
++	return s*(hfsq+R) + dk*ln2_lo - hfsq + f + dk*ln2_hi;
+ }
+diff --git a/src/math/logl.c b/src/math/logl.c
+index ef2b551..03c5188 100644
+--- a/src/math/logl.c
++++ b/src/math/logl.c
+@@ -35,9 +35,9 @@
+  *
+  *     log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
+  *
+- * Otherwise, setting  z = 2(x-1)/x+1),
++ * Otherwise, setting  z = 2(x-1)/(x+1),
+  *
+- *     log(x) = z + z**3 P(z)/Q(z).
++ *     log(x) = log(1+z/2) - log(1-z/2) = z + z**3 P(z)/Q(z).
+  *
+  *
+  * ACCURACY:
+@@ -116,7 +116,7 @@ long double logl(long double x)
+ 		return x;
+ 	if (x <= 0.0) {
+ 		if (x == 0.0)
+-			return -1/(x+0); /* -inf with divbyzero */
++			return -1/(x*x); /* -inf with divbyzero */
+ 		return 0/0.0f; /* nan with invalid */
+ 	}
+ 
+@@ -127,7 +127,7 @@ long double logl(long double x)
+ 	x = frexpl(x, &e);
+ 
+ 	/* logarithm using log(x) = z + z**3 P(z)/Q(z),
+-	 * where z = 2(x-1)/x+1)
++	 * where z = 2(x-1)/(x+1)
+ 	 */
+ 	if (e > 2 || e < -2) {
+ 		if (x < SQRTH) {  /* 2(2x-1)/(2x+1) */
+diff --git a/src/math/modfl.c b/src/math/modfl.c
+index fc85bb5..f736bba 100644
+--- a/src/math/modfl.c
++++ b/src/math/modfl.c
+@@ -14,7 +14,6 @@ long double modfl(long double x, long double *iptr)
+ long double modfl(long double x, long double *iptr)
+ {
+ 	union ldshape u = {x};
+-	uint64_t mask;
+ 	int e = (u.i.se & 0x7fff) - 0x3fff;
+ 	int s = u.i.se >> 15;
+ 	long double absx;
+diff --git a/src/multibyte/mbsrtowcs.c b/src/multibyte/mbsrtowcs.c
+index b9bbc33..066cce6 100644
+--- a/src/multibyte/mbsrtowcs.c
++++ b/src/multibyte/mbsrtowcs.c
+@@ -59,7 +59,7 @@ resume0:
+ 			return wn0;
+ 		}
+ 		if (*s-1u < 0x7f && (uintptr_t)s%4 == 0) {
+-			while (wn>=4 && !(( *(uint32_t*)s | *(uint32_t*)s-0x01010101) & 0x80808080)) {
++			while (wn>=5 && !(( *(uint32_t*)s | *(uint32_t*)s-0x01010101) & 0x80808080)) {
+ 				*ws++ = *s++;
+ 				*ws++ = *s++;
+ 				*ws++ = *s++;
+diff --git a/src/network/__ipparse.c b/src/network/__ipparse.c
+index b0647aa..2480265 100644
+--- a/src/network/__ipparse.c
++++ b/src/network/__ipparse.c
+@@ -1,27 +1,31 @@
+ #include <string.h>
+ #include <stdlib.h>
++#include <ctype.h>
+ #include <sys/socket.h>
+ #include <netinet/in.h>
+ #include <arpa/inet.h>
+ #include "__dns.h"
+-#include <stdio.h>
+ 
+ int __ipparse(void *dest, int family, const char *s0)
+ {
+ 	const char *s = s0;
+ 	unsigned char *d = dest;
+ 	unsigned long a[16] = { 0 };
+-	const char *z;
++	char *z;
+ 	int i;
+ 
+ 	if (family == AF_INET6) goto not_v4;
+ 
+ 	for (i=0; i<4; i++) {
+-		a[i] = strtoul(s, (char **)&z, 0);
+-		if (z==s || (*z && *z != '.')) goto not_v4;
++		a[i] = strtoul(s, &z, 0);
++		if (z==s || (*z && *z != '.') || !isdigit(*s)) {
++			if (family == AF_INET) return -1;
++			goto not_v4;
++		}
+ 		if (!*z) break;
+ 		s=z+1;
+ 	}
++	if (i==4) return -1;
+ 	switch (i) {
+ 	case 0:
+ 		a[1] = a[0] & 0xffffff;
+@@ -35,7 +39,10 @@ int __ipparse(void *dest, int family, const char *s0)
+ 	}
+ 	((struct sockaddr_in *)d)->sin_family = AF_INET;
+ 	d = (void *)&((struct sockaddr_in *)d)->sin_addr;
+-	for (i=0; i<4; i++) d[i] = a[i];
++	for (i=0; i<4; i++) {
++		if (a[i] > 255) return -1;
++		d[i] = a[i];
++	}
+ 	return 0;
+ 
+ not_v4:
+diff --git a/src/network/inet_addr.c b/src/network/inet_addr.c
+new file mode 100644
+index 0000000..8413728
+--- /dev/null
++++ b/src/network/inet_addr.c
+@@ -0,0 +1,11 @@
++#include <sys/socket.h>
++#include <netinet/in.h>
++#include <arpa/inet.h>
++#include "__dns.h"
++
++in_addr_t inet_addr(const char *p)
++{
++	struct sockaddr_in sin;
++	if (__ipparse(&sin, AF_INET, p)) return -1;
++	return sin.sin_addr.s_addr;
++}
+diff --git a/src/network/inet_legacy.c b/src/network/inet_legacy.c
+index e802557..dd75420 100644
+--- a/src/network/inet_legacy.c
++++ b/src/network/inet_legacy.c
+@@ -1,16 +1,8 @@
+ #include <sys/socket.h>
+ #include <netinet/in.h>
+ #include <arpa/inet.h>
+-#include <stdio.h>
+ #include "__dns.h"
+ 
+-in_addr_t inet_addr(const char *p)
+-{
+-	struct sockaddr_in sin;
+-	if (__ipparse(&sin, AF_INET, p)) return -1;
+-	return sin.sin_addr.s_addr;
+-}
+-
+ in_addr_t inet_network(const char *p)
+ {
+ 	return ntohl(inet_addr(p));
+@@ -18,15 +10,10 @@ in_addr_t inet_network(const char *p)
+ 
+ int inet_aton(const char *cp, struct in_addr *inp)
+ {
+-	return inet_pton(AF_INET, cp, (void *)inp) > 0;
+-}
+-
+-char *inet_ntoa(struct in_addr in)
+-{
+-	static char buf[16];
+-	unsigned char *a = (void *)&in;
+-	snprintf(buf, sizeof buf, "%d.%d.%d.%d", a[0], a[1], a[2], a[3]);
+-	return buf;
++	struct sockaddr_in sin;
++	int r = __ipparse(&sin, AF_INET, cp);
++	*inp = sin.sin_addr;
++	return r;
+ }
+ 
+ struct in_addr inet_makeaddr(int net, int host)
+diff --git a/src/network/inet_ntoa.c b/src/network/inet_ntoa.c
+new file mode 100644
+index 0000000..71411e0
+--- /dev/null
++++ b/src/network/inet_ntoa.c
+@@ -0,0 +1,10 @@
++#include <arpa/inet.h>
++#include <stdio.h>
++
++char *inet_ntoa(struct in_addr in)
++{
++	static char buf[16];
++	unsigned char *a = (void *)&in;
++	snprintf(buf, sizeof buf, "%d.%d.%d.%d", a[0], a[1], a[2], a[3]);
++	return buf;
++}
+diff --git a/src/network/inet_pton.c b/src/network/inet_pton.c
+index 5c4850a..f840dd4 100644
+--- a/src/network/inet_pton.c
++++ b/src/network/inet_pton.c
+@@ -1,7 +1,6 @@
+ #include <sys/socket.h>
+ #include <netdb.h>
+ #include <arpa/inet.h>
+-#include <stdlib.h>
+ #include <ctype.h>
+ #include <errno.h>
+ #include <string.h>
+@@ -18,52 +17,46 @@ int inet_pton(int af, const char *restrict s, void *restrict a0)
+ {
+ 	uint16_t ip[8];
+ 	unsigned char *a = a0;
+-	const char *z;
+-	unsigned long x;
+ 	int i, j, v, d, brk=-1, need_v4=0;
+ 
+-	/* Reimplement this because inet_pton cannot accept special v4 forms */
+ 	if (af==AF_INET) {
+-		for (i=0; i<4 && *s; i++) {
+-			a[i] = x = strtoul(s, (char **)&z, 10);
+-			if (!isdigit(*s) || z==s || (*z && *z != '.') || x>255)
+-				return 0;
+-			s=z+1;
++		for (i=0; i<4; i++) {
++			for (v=j=0; j<3 && isdigit(s[j]); j++)
++				v = 10*v + s[j]-'0';
++			if (j==0 || (j>1 && s[0]=='0') || v>255) return 0;
++			a[i] = v;
++			if (s[j]==0 && i==3) return 1;
++			if (s[j]!='.') return 0;
++			s += j+1;
+ 		}
+-		return 1;
++		return 0;
+ 	} else if (af!=AF_INET6) {
+ 		errno = EAFNOSUPPORT;
+ 		return -1;
+ 	}
+ 
+-	if (s[0]==':' && s[1]==':') s++;
++	if (*s==':' && *++s!=':') return 0;
+ 
+-	for (i=0; ; i++, s+=j+1) {
++	for (i=0; ; i++) {
+ 		if (s[0]==':' && brk<0) {
+ 			brk=i;
+-			j=0;
+ 			ip[i]=0;
+-			if (!s[1]) break;
++			if (!*++s) break;
+ 			continue;
+ 		}
+-		if (hexval(s[0])<0) return -1;
+-		while (s[0]=='0' && s[1]=='0') s++;
+-		for (v=j=0; j<5 && (d=hexval(s[j]))>=0; j++)
++		for (v=j=0; j<4 && (d=hexval(s[j]))>=0; j++)
+ 			v=16*v+d;
+-		if (v > 65535) return -1;
++		if (j==0 || v > 65535) return 0;
+ 		ip[i] = v;
+-		if (!s[j]) {
+-			if (brk<0 && i!=7) return -1;
+-			break;
+-		}
+-		if (i<7) {
+-			if (s[j]==':') continue;
+-			if (s[j]!='.') return -1;
++		if (!s[j] && (brk>=0 || i==7)) break;
++		if (i==7) return 0;
++		if (s[j]!=':') {
++			if (s[j]!='.' || (i<6 && brk<0)) return 0;
+ 			need_v4=1;
+ 			i++;
+ 			break;
+ 		}
+-		return -1;
++		s += j+1;
+ 	}
+ 	if (brk>=0) {
+ 		memmove(ip+brk+7-i, ip+brk, 2*(i+1-brk));
+@@ -73,6 +66,6 @@ int inet_pton(int af, const char *restrict s, void *restrict a0)
+ 		*a++ = ip[j]>>8;
+ 		*a++ = ip[j];
+ 	}
+-	if (need_v4 &&inet_pton(AF_INET, (void *)s, a-4) <= 0) return -1;
++	if (need_v4 && inet_pton(AF_INET, (void *)s, a-4) <= 0) return 0;
+ 	return 1;
+ }
+diff --git a/src/passwd/getgr_r.c b/src/passwd/getgr_r.c
+index 234c901..3fe2e2b 100644
+--- a/src/passwd/getgr_r.c
++++ b/src/passwd/getgr_r.c
+@@ -26,14 +26,14 @@ static int getgr_r(const char *name, gid_t gid, struct group *gr, char *buf, siz
+ 	while (__getgrent_a(f, gr, &line, &len, &mem, &nmem)) {
+ 		if (name && !strcmp(name, gr->gr_name)
+ 		|| !name && gr->gr_gid == gid) {
+-			if (size < len + nmem*sizeof(char *) + 32) {
++			if (size < len + (nmem+1)*sizeof(char *) + 32) {
+ 				rv = ERANGE;
+ 				break;
+ 			}
+ 			*res = gr;
+ 			buf += (16-(uintptr_t)buf)%16;
+ 			gr->gr_mem = (void *)buf;
+-			buf += nmem*sizeof(char *);
++			buf += (nmem+1)*sizeof(char *);
+ 			memcpy(buf, line, len);
+ 			FIX(name);
+ 			FIX(passwd);
+diff --git a/src/process/execl.c b/src/process/execl.c
+index 327d78b..5ee5c81 100644
+--- a/src/process/execl.c
++++ b/src/process/execl.c
+@@ -16,6 +16,7 @@ int execl(const char *path, const char *argv0, ...)
+ 		for (i=1; i<argc; i++)
+ 			argv[i] = va_arg(ap, char *);
+ 		argv[i] = NULL;
++		va_end(ap);
+ 		return execv(path, argv);
+ 	}
+ }
+diff --git a/src/process/execle.c b/src/process/execle.c
+index 6490836..37ca503 100644
+--- a/src/process/execle.c
++++ b/src/process/execle.c
+@@ -14,10 +14,10 @@ int execle(const char *path, const char *argv0, ...)
+ 		char **envp;
+ 		va_start(ap, argv0);
+ 		argv[0] = (char *)argv0;
+-		for (i=1; i<argc; i++)
++		for (i=1; i<=argc; i++)
+ 			argv[i] = va_arg(ap, char *);
+-		argv[i] = NULL;
+ 		envp = va_arg(ap, char **);
++		va_end(ap);
+ 		return execve(path, argv, envp);
+ 	}
+ }
+diff --git a/src/process/execlp.c b/src/process/execlp.c
+index c7b3f51..5eed886 100644
+--- a/src/process/execlp.c
++++ b/src/process/execlp.c
+@@ -16,6 +16,7 @@ int execlp(const char *file, const char *argv0, ...)
+ 		for (i=1; i<argc; i++)
+ 			argv[i] = va_arg(ap, char *);
+ 		argv[i] = NULL;
++		va_end(ap);
+ 		return execvp(file, argv);
+ 	}
+ }
+diff --git a/src/regex/regcomp.c b/src/regex/regcomp.c
+index 5cedfd5..1633754 100644
+--- a/src/regex/regcomp.c
++++ b/src/regex/regcomp.c
+@@ -516,7 +516,7 @@ tre_new_item(tre_mem_t mem, int min, int max, int *i, int *max_i,
+       if (*max_i > 1024)
+ 	return REG_ESPACE;
+       *max_i *= 2;
+-      new_items = xrealloc(array, sizeof(*items) * *max_i);
++      new_items = xrealloc(array, sizeof(*array) * *max_i);
+       if (new_items == NULL)
+ 	return REG_ESPACE;
+       *items = array = new_items;
+@@ -765,7 +765,7 @@ tre_parse_bracket(tre_parse_ctx_t *ctx, tre_ast_node_t **result)
+ 	  if (num_neg_classes > 0)
+ 	    {
+ 	      l->neg_classes = tre_mem_alloc(ctx->mem,
+-					     (sizeof(l->neg_classes)
++					     (sizeof(*l->neg_classes)
+ 					      * (num_neg_classes + 1)));
+ 	      if (l->neg_classes == NULL)
+ 		{
+@@ -805,7 +805,7 @@ tre_parse_bracket(tre_parse_ctx_t *ctx, tre_ast_node_t **result)
+ 	  if (num_neg_classes > 0)
+ 	    {
+ 	      l->neg_classes = tre_mem_alloc(ctx->mem,
+-					     (sizeof(l->neg_classes)
++					     (sizeof(*l->neg_classes)
+ 					      * (num_neg_classes + 1)));
+ 	      if (l->neg_classes == NULL)
+ 		{
+@@ -3167,7 +3167,7 @@ regcomp(regex_t *restrict preg, const char *restrict regex, int cflags)
+ 		 sizeof(*tag_directions) * (tnfa->num_tags + 1));
+ 	}
+       tnfa->minimal_tags = xcalloc((unsigned)tnfa->num_tags * 2 + 1,
+-				   sizeof(tnfa->minimal_tags));
++				   sizeof(*tnfa->minimal_tags));
+       if (tnfa->minimal_tags == NULL)
+ 	ERROR_EXIT(REG_ESPACE);
+ 
+diff --git a/src/stdio/vfprintf.c b/src/stdio/vfprintf.c
+index a2b287b..b5948bd 100644
+--- a/src/stdio/vfprintf.c
++++ b/src/stdio/vfprintf.c
+@@ -530,7 +530,6 @@ static int printf_core(FILE *f, const char *fmt, va_list *ap, union arg *nl_arg,
+ 		/* Check validity of argument type (nl/normal) */
+ 		if (st==NOARG) {
+ 			if (argpos>=0) return -1;
+-			else if (!f) continue;
+ 		} else {
+ 			if (argpos>=0) nl_type[argpos]=st, arg=nl_arg[argpos];
+ 			else if (f) pop_arg(&arg, st, ap);
+@@ -660,8 +659,12 @@ int vfprintf(FILE *restrict f, const char *restrict fmt, va_list ap)
+ 	unsigned char internal_buf[80], *saved_buf = 0;
+ 	int ret;
+ 
++	/* the copy allows passing va_list* even if va_list is an array */
+ 	va_copy(ap2, ap);
+-	if (printf_core(0, fmt, &ap2, nl_arg, nl_type) < 0) return -1;
++	if (printf_core(0, fmt, &ap2, nl_arg, nl_type) < 0) {
++		va_end(ap2);
++		return -1;
++	}
+ 
+ 	FLOCK(f);
+ 	if (!f->buf_size) {
+diff --git a/src/stdio/vfwprintf.c b/src/stdio/vfwprintf.c
+index eb07931..984ff7b 100644
+--- a/src/stdio/vfwprintf.c
++++ b/src/stdio/vfwprintf.c
+@@ -167,7 +167,7 @@ static const char sizeprefix['y'-'a'] = {
+ 
+ static int wprintf_core(FILE *f, const wchar_t *fmt, va_list *ap, union arg *nl_arg, int *nl_type)
+ {
+-	wchar_t *a, *z, *s=(wchar_t *)fmt, *s0;
++	wchar_t *a, *z, *s=(wchar_t *)fmt;
+ 	unsigned l10n=0, litpct, fl;
+ 	int w, p;
+ 	union arg arg;
+@@ -242,7 +242,6 @@ static int wprintf_core(FILE *f, const wchar_t *fmt, va_list *ap, union arg *nl_
+ 		} else p = -1;
+ 
+ 		/* Format specifier state machine */
+-		s0=s;
+ 		st=0;
+ 		do {
+ 			if (OOB(*s)) return -1;
+@@ -254,7 +253,6 @@ static int wprintf_core(FILE *f, const wchar_t *fmt, va_list *ap, union arg *nl_
+ 		/* Check validity of argument type (nl/normal) */
+ 		if (st==NOARG) {
+ 			if (argpos>=0) return -1;
+-			else if (!f) continue;
+ 		} else {
+ 			if (argpos>=0) nl_type[argpos]=st, arg=nl_arg[argpos];
+ 			else if (f) pop_arg(&arg, st, ap);
+@@ -288,8 +286,7 @@ static int wprintf_core(FILE *f, const wchar_t *fmt, va_list *ap, union arg *nl_
+ 		case 'S':
+ 			a = arg.p;
+ 			z = wmemchr(a, 0, p);
+-			if (!z) z=a+p;
+-			else p=z-a;
++			if (z) p=z-a;
+ 			if (w<p) w=p;
+ 			if (!(fl&LEFT_ADJ)) fprintf(f, "%.*s", w-p, "");
+ 			out(f, a, p);
+@@ -350,8 +347,12 @@ int vfwprintf(FILE *restrict f, const wchar_t *restrict fmt, va_list ap)
+ 	union arg nl_arg[NL_ARGMAX];
+ 	int ret;
+ 
++	/* the copy allows passing va_list* even if va_list is an array */
+ 	va_copy(ap2, ap);
+-	if (wprintf_core(0, fmt, &ap2, nl_arg, nl_type) < 0) return -1;
++	if (wprintf_core(0, fmt, &ap2, nl_arg, nl_type) < 0) {
++		va_end(ap2);
++		return -1;
++	}
+ 
+ 	FLOCK(f);
+ 	ret = wprintf_core(f, fmt, &ap2, nl_arg, nl_type);
+diff --git a/src/thread/pthread_key_create.c b/src/thread/pthread_key_create.c
+index e51cb02..c29935c 100644
+--- a/src/thread/pthread_key_create.c
++++ b/src/thread/pthread_key_create.c
+@@ -17,7 +17,7 @@ int pthread_key_create(pthread_key_t *k, void (*dtor)(void *))
+ 	__pthread_self_init();
+ 	if (!dtor) dtor = nodtor;
+ 	do {
+-		if (!a_cas_p(keys+j, 0, dtor)) {
++		if (!a_cas_p(keys+j, 0, (void *)dtor)) {
+ 			*k = j;
+ 			return 0;
+ 		}
+diff --git a/src/time/ftime.c b/src/time/ftime.c
+new file mode 100644
+index 0000000..a1734d0
+--- /dev/null
++++ b/src/time/ftime.c
+@@ -0,0 +1,12 @@
++#include <sys/timeb.h>
++#include <time.h>
++
++int ftime(struct timeb *tp)
++{
++	struct timespec ts;
++	clock_gettime(CLOCK_REALTIME, &ts);
++	tp->time = ts.tv_sec;
++	tp->millitm = ts.tv_nsec / 1000000;
++	tp->timezone = tp->dstflag = 0;
++	return 0;
++}
+diff --git a/src/unistd/faccessat.c b/src/unistd/faccessat.c
+index 5b2c5e3..821e13f 100644
+--- a/src/unistd/faccessat.c
++++ b/src/unistd/faccessat.c
+@@ -14,8 +14,8 @@ static int checker(void *p)
+ {
+ 	struct ctx *c = p;
+ 	int ret;
+-	if (__syscall(SYS_setgid, __syscall(SYS_getegid))
+-	    || __syscall(SYS_setuid, __syscall(SYS_geteuid)))
++	if (__syscall(SYS_setregid, __syscall(SYS_getegid), -1)
++	    || __syscall(SYS_setreuid, __syscall(SYS_geteuid), -1))
+ 		__syscall(SYS_exit, 1);
+ 	ret = __syscall(SYS_faccessat, c->fd, c->filename, c->amode, 0);
+ 	__syscall(SYS_write, c->p, &ret, sizeof ret);
+@@ -34,7 +34,7 @@ int faccessat(int fd, const char *filename, int amode, int flag)
+ 	sigset_t set;
+ 	int ret, p[2];
+ 
+-	if (pipe(p)) return __syscall_ret(-EBUSY);
++	if (pipe2(p, O_CLOEXEC)) return __syscall_ret(-EBUSY);
+ 	struct ctx c = { .fd = fd, .filename = filename, .amode = amode, .p = p[1] };
+ 
+ 	__block_all_sigs(&set);
+diff --git a/src/unistd/getcwd.c b/src/unistd/getcwd.c
+index 2e540cd..a7b925d 100644
+--- a/src/unistd/getcwd.c
++++ b/src/unistd/getcwd.c
+@@ -7,7 +7,13 @@
+ char *getcwd(char *buf, size_t size)
+ {
+ 	char tmp[PATH_MAX];
+-	if (!buf) buf = tmp, size = PATH_MAX;
++	if (!buf) {
++		buf = tmp;
++		size = PATH_MAX;
++	} else if (!size) {
++		errno = EINVAL;
++		return 0;
++	}
+ 	if (syscall(SYS_getcwd, buf, size) < 0) return 0;
+ 	return buf == tmp ? strdup(buf) : buf;
+ }
diff --git a/main/musl/0002-fix-buffer-overflow-in-mbsrtowcs.patch b/main/musl/0002-fix-buffer-overflow-in-mbsrtowcs.patch
deleted file mode 100644
index 8b68cc8777..0000000000
--- a/main/musl/0002-fix-buffer-overflow-in-mbsrtowcs.patch
+++ /dev/null
@@ -1,41 +0,0 @@
-From 211264e46a2f1bc382a84435e904d1548de672b0 Mon Sep 17 00:00:00 2001
-From: Rich Felker <dalias@aerifal.cx>
-Date: Fri, 27 Sep 2013 13:09:46 -0400
-Subject: [PATCH] fix buffer overflow in mbsrtowcs
-
-issue reported by Michael Forney:
-
-"If wn becomes 0 after processing a chunk of 4, mbsrtowcs currently
-continues on, wrapping wn around to -1, causing the rest of the string
-to be processed.
-
-This resulted in buffer overruns if there was only space in ws for wn
-wide characters."
-
-the original patch submitted added an additional check for !wn after
-the loop; to avoid extra branching, I instead just changed the wn>=4
-check to wn>=5 to ensure that at least one slot remains after the
-word-at-a-time loop runs. this should not slow down the tail
-processing on real-world usage, since an extra slot that can't be
-processed in the word-at-a-time loop is needed for the null
-termination anyway.
----
- src/multibyte/mbsrtowcs.c | 2 +-
- 1 file changed, 1 insertion(+), 1 deletion(-)
-
-diff --git a/src/multibyte/mbsrtowcs.c b/src/multibyte/mbsrtowcs.c
-index b9bbc33..066cce6 100644
---- a/src/multibyte/mbsrtowcs.c
-+++ b/src/multibyte/mbsrtowcs.c
-@@ -59,7 +59,7 @@ resume0:
- 			return wn0;
- 		}
- 		if (*s-1u < 0x7f && (uintptr_t)s%4 == 0) {
--			while (wn>=4 && !(( *(uint32_t*)s | *(uint32_t*)s-0x01010101) & 0x80808080)) {
-+			while (wn>=5 && !(( *(uint32_t*)s | *(uint32_t*)s-0x01010101) & 0x80808080)) {
- 				*ws++ = *s++;
- 				*ws++ = *s++;
- 				*ws++ = *s++;
--- 
-1.8.4
-
diff --git a/main/musl/0003-fix-off-by-one-error-in-getgrnam_r-and-getgrgid_r-cl.patch b/main/musl/0003-fix-off-by-one-error-in-getgrnam_r-and-getgrgid_r-cl.patch
deleted file mode 100644
index ebc1910c7f..0000000000
--- a/main/musl/0003-fix-off-by-one-error-in-getgrnam_r-and-getgrgid_r-cl.patch
+++ /dev/null
@@ -1,38 +0,0 @@
-From 23b8e3bc95620b0bd90a78ce0d926942c12b45da Mon Sep 17 00:00:00 2001
-From: Rich Felker <dalias@aerifal.cx>
-Date: Sun, 29 Sep 2013 02:52:33 -0400
-Subject: [PATCH] fix off-by-one error in getgrnam_r and getgrgid_r, clobbering
- gr_name
-
-bug report and patch by Michael Forney. the terminating null pointer
-at the end of the gr_mem array was overwriting the beginning of the
-string data, causing the gr_name member to always be a zero-length
-string.
----
- src/passwd/getgr_r.c | 4 ++--
- 1 file changed, 2 insertions(+), 2 deletions(-)
-
-diff --git a/src/passwd/getgr_r.c b/src/passwd/getgr_r.c
-index 234c901..3fe2e2b 100644
---- a/src/passwd/getgr_r.c
-+++ b/src/passwd/getgr_r.c
-@@ -26,14 +26,14 @@ static int getgr_r(const char *name, gid_t gid, struct group *gr, char *buf, siz
- 	while (__getgrent_a(f, gr, &line, &len, &mem, &nmem)) {
- 		if (name && !strcmp(name, gr->gr_name)
- 		|| !name && gr->gr_gid == gid) {
--			if (size < len + nmem*sizeof(char *) + 32) {
-+			if (size < len + (nmem+1)*sizeof(char *) + 32) {
- 				rv = ERANGE;
- 				break;
- 			}
- 			*res = gr;
- 			buf += (16-(uintptr_t)buf)%16;
- 			gr->gr_mem = (void *)buf;
--			buf += nmem*sizeof(char *);
-+			buf += (nmem+1)*sizeof(char *);
- 			memcpy(buf, line, len);
- 			FIX(name);
- 			FIX(passwd);
--- 
-1.8.4
-
diff --git a/main/musl/0004-fix-new-environment-always-being-null-with-execle.patch b/main/musl/0004-fix-new-environment-always-being-null-with-execle.patch
deleted file mode 100644
index 31a1672e22..0000000000
--- a/main/musl/0004-fix-new-environment-always-being-null-with-execle.patch
+++ /dev/null
@@ -1,36 +0,0 @@
-From 2b2aff37aced66e4a50a38a14607a9b1dc0ee001 Mon Sep 17 00:00:00 2001
-From: Rich Felker <dalias@aerifal.cx>
-Date: Thu, 3 Oct 2013 10:16:01 -0400
-Subject: [PATCH] fix new environment always being null with execle
-MIME-Version: 1.0
-Content-Type: text/plain; charset=UTF-8
-Content-Transfer-Encoding: 8bit
-
-the va_arg call for the argv[]-terminating null pointer was missing,
-so this pointer was being wrongly used as the environment pointer.
-
-issue reported by Timo Teräs. proposed patch slightly modified to
-simplify the resulting code.
----
- src/process/execle.c | 3 +--
- 1 file changed, 1 insertion(+), 2 deletions(-)
-
-diff --git a/src/process/execle.c b/src/process/execle.c
-index 6490836..6f2ec05 100644
---- a/src/process/execle.c
-+++ b/src/process/execle.c
-@@ -14,9 +14,8 @@ int execle(const char *path, const char *argv0, ...)
- 		char **envp;
- 		va_start(ap, argv0);
- 		argv[0] = (char *)argv0;
--		for (i=1; i<argc; i++)
-+		for (i=1; i<=argc; i++)
- 			argv[i] = va_arg(ap, char *);
--		argv[i] = NULL;
- 		envp = va_arg(ap, char **);
- 		return execve(path, argv, envp);
- 	}
--- 
-1.8.4
-
-
diff --git a/main/musl/0005-fix-clockid-macro-names-in-time.h-reported-by-Paul-S.patch b/main/musl/0005-fix-clockid-macro-names-in-time.h-reported-by-Paul-S.patch
deleted file mode 100644
index 3314319346..0000000000
--- a/main/musl/0005-fix-clockid-macro-names-in-time.h-reported-by-Paul-S.patch
+++ /dev/null
@@ -1,28 +0,0 @@
-From 551da2c1916e3166a0562ac6f2e76fd46781f3d7 Mon Sep 17 00:00:00 2001
-From: Szabolcs Nagy <nsz@port70.net>
-Date: Mon, 7 Oct 2013 13:28:03 +0000
-Subject: [PATCH] fix clockid macro names in time.h, reported by Paul Schutte
-
----
- include/time.h | 4 ++--
- 1 file changed, 2 insertions(+), 2 deletions(-)
-
-diff --git a/include/time.h b/include/time.h
-index 6b2a069..6e499ff 100644
---- a/include/time.h
-+++ b/include/time.h
-@@ -82,8 +82,8 @@ struct itimerspec
- #define CLOCK_PROCESS_CPUTIME_ID 2
- #define CLOCK_THREAD_CPUTIME_ID  3
- #define CLOCK_MONOTONIC_RAW      4
--#define CLOCK_REALTIME_COURSE    5
--#define CLOCK_MONOTONIC_COURSE   6
-+#define CLOCK_REALTIME_COARSE    5
-+#define CLOCK_MONOTONIC_COARSE   6
- #define CLOCK_BOOTTIME           7
- #define CLOCK_REALTIME_ALARM     8
- #define CLOCK_BOOTTIME_ALARM     9
--- 
-1.8.4
-
-
diff --git a/main/musl/APKBUILD b/main/musl/APKBUILD
index 0d3e44f8f3..3b9e818395 100644
--- a/main/musl/APKBUILD
+++ b/main/musl/APKBUILD
@@ -2,7 +2,7 @@
 # Maintainer: Timo Teräs <timo.teras@iki.fi>
 pkgname=musl
 pkgver=0.9.14
-pkgrel=5
+pkgrel=6
 pkgdesc="the musl c library (libc) implementation"
 url="http://www.musl-libc.org/"
 arch="all"
@@ -14,12 +14,7 @@ install=""
 subpackages="$pkgname-dev $pkgname-utils"
 [ "${CTARGET#*musl}" = "$CTARGET" ] && subpackages="$subpackages musl-gcc:crosstool"
 source="http://www.musl-libc.org/releases/musl-$pkgver.tar.gz
-	0001-add-missing-i386-syscall-numbers.patch
-	0002-fix-buffer-overflow-in-mbsrtowcs.patch
-	0003-fix-off-by-one-error-in-getgrnam_r-and-getgrgid_r-cl.patch
-	0004-fix-new-environment-always-being-null-with-execle.patch
-	0005-fix-clockid-macro-names-in-time.h-reported-by-Paul-S.patch
-
+	0001-updates-from-git.patch
 	1001-add-rfc3678-mcast-structs.patch
 	1002-add-linux-tcp-state-enums.patch
 	1003-add-basic-dns-record-parsing-functions.patch
@@ -115,11 +110,7 @@ crosstool() {
 }
 
 md5sums="bfb685695aa942e64c63170589e575b2  musl-0.9.14.tar.gz
-f116cf69bcbcb7080ef3aa521acce8b8  0001-add-missing-i386-syscall-numbers.patch
-bfefbd099f555fe8fd22e7ffc3accbef  0002-fix-buffer-overflow-in-mbsrtowcs.patch
-5d722e38a7ca2032c9f202db8ff7e369  0003-fix-off-by-one-error-in-getgrnam_r-and-getgrgid_r-cl.patch
-216d6915ea8397102c22a7aeaafdd98f  0004-fix-new-environment-always-being-null-with-execle.patch
-aa4c3ccf3e074fcac609420573ce466a  0005-fix-clockid-macro-names-in-time.h-reported-by-Paul-S.patch
+2ef3f1267630bd7ec0dbf4171f4ccea8  0001-updates-from-git.patch
 b28080b5c8c1d44521277aa3255d280a  1001-add-rfc3678-mcast-structs.patch
 34044ab59029e9510a2760d669e1a377  1002-add-linux-tcp-state-enums.patch
 6cdf1c56450d59f3a3acf452b2db4c2e  1003-add-basic-dns-record-parsing-functions.patch
@@ -129,11 +120,7 @@ b28080b5c8c1d44521277aa3255d280a  1001-add-rfc3678-mcast-structs.patch
 ef81489a6258501cf45db58dfc6d5211  getent
 33e4fd94e2560e008e2c3b431d0e3419  ldconfig"
 sha256sums="982e9de1287cf95f9aa526adba008660d8885bfccc41faf5c613ea47f1922872  musl-0.9.14.tar.gz
-84886493008bdca79ad223708e5568baeb948a520499b9f4eea2f1526aefb304  0001-add-missing-i386-syscall-numbers.patch
-a6cb8b279e5b737d43c2de1bd6229f5e6599e9514bcb41ddaa411cd44dc65ba7  0002-fix-buffer-overflow-in-mbsrtowcs.patch
-b6b161383b287505eecb53595039b8fe26c622508e783c86ee82d38c1ea582f4  0003-fix-off-by-one-error-in-getgrnam_r-and-getgrgid_r-cl.patch
-2e9f262bf9c117f1ca59afd9615daa6f8d71aaddb719f5259218b088244622fb  0004-fix-new-environment-always-being-null-with-execle.patch
-f2fd3b921bd57190d0bf19c4a4beb5f808a5d19c71092e2962a3e7b1ca08c67e  0005-fix-clockid-macro-names-in-time.h-reported-by-Paul-S.patch
+79d286c3cd31bd9710784388bc7fa7e3cedf0d025e2b6a83d30346be86febac3  0001-updates-from-git.patch
 720cb88dd1ef57fc806a22b46b3c47d7a0a38a34d31edb583e97bfa7a47eb44c  1001-add-rfc3678-mcast-structs.patch
 53637d1dfdbff7131277252d63cba7c3fc1f7b61c7b735e503fbccaa6dcdd887  1002-add-linux-tcp-state-enums.patch
 54686df1392c52f4e9c62648dcb544f4bd48111be8d9734b7f65d8452b7ead12  1003-add-basic-dns-record-parsing-functions.patch
@@ -143,11 +130,7 @@ d9b644ec20bc33e81a7c52b9fcf7973d835923a69faf50f03db45534b811bd96  getopt_long.c
 d6996273f5aaaed429058257e4646b243d9e3a4d8609522f802762453f5be4cb  getent
 306c6ca7407560340797866e077e053627ad409277d1b9da58106fce4cf717cb  ldconfig"
 sha512sums="e5c3f7b1549dc2f9cbd3359cc413f761d5967607c23705f651c33d0ae93f00582193a41fe1f87158467d58d8eba2d7c09e0fe2f2b2c02c1dda78eee1a4cecff6  musl-0.9.14.tar.gz
-e6cdb6b4c87e2488d31ac46898010dc9d41f31f2ed9c6f2f8a763b76e63587a37a54a7557cd7db1c524b1fcbd17e5418ec1058b30dc17cc91c3fb3ac7cd47fc8  0001-add-missing-i386-syscall-numbers.patch
-bb1a0025675fc0241a0ae15e04e12cfbb3c12604ccdb8ac9a2a192d76d346c2804d051f6698a636671d9b4154132b1fe4270bc464bdd370a9b2070c4573361e0  0002-fix-buffer-overflow-in-mbsrtowcs.patch
-60fd9640ea6a8c46e8724ec80d228bcabb9ae3f0f366d2a21c7e15da4a31395aa6dfadeb480f22c720cbd773b73b245174adc3031ee04e69efecf4c45af8538f  0003-fix-off-by-one-error-in-getgrnam_r-and-getgrgid_r-cl.patch
-f37b4bb6e15c0fa02cb7643a3f291158901a13f5c39651d3feda371f9de9a05ce82ce7e81a45440422a4bf9198cd2baa1cb8fbebe267b6726e8d293aa5b836d3  0004-fix-new-environment-always-being-null-with-execle.patch
-f3f2f0cad04b58891e42d0846d929c05d793c1c7ab61aa0ac1377f38c1291bce45e39a7acdc62749804859a5de1b20f9ec0839c58a7015cb12c0b7c85c7a6194  0005-fix-clockid-macro-names-in-time.h-reported-by-Paul-S.patch
+ec170e0680627426bf4fd53f7b09fdd5aad2546905875a5db58414ef8006ca81ceca0f69719500005961accb1263b2627b1246e02ccc35dc8bf3b557d64f71fd  0001-updates-from-git.patch
 d63d3ea0c59dbff0cb24cc30ba0b013a33e81bc91de5f6c8b0082fa3e261b0c0dac3f2dc30c3b0370f7d1c4536e2042cd85ecb8e0357a805bd93855e42396d92  1001-add-rfc3678-mcast-structs.patch
 53249ab17705190b17905623cc22b42616e1648f91f8f129bd1a5299b0231ce216154e009a0c8f62444f12fa6e51655cefeade7d383ed0427d879d508f338209  1002-add-linux-tcp-state-enums.patch
 335ec63cfbc7f348f33cfba1238c069fed4c8c51a51d1ea39eff4b7dfcceeae7ff4afb1038fa7c5545a42854ed553a936c0b1ff6c4795fa25b982398a2cc02bd  1003-add-basic-dns-record-parsing-functions.patch