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 #include #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 + +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 -#include "libm.h" +#include +#include /* * 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 +#include 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 +#include 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 +#include 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 +#include 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 +#include 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 +#include 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 #include +#include #include #include #include #include "__dns.h" -#include 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 +#include +#include +#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 #include #include -#include #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 *)∈ - 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 +#include + +char *inet_ntoa(struct in_addr in) +{ + static char buf[16]; + unsigned char *a = (void *)∈ + 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 #include #include -#include #include #include #include @@ -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 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 +#include + +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; }