diff options
| author | Samuel Thibault <samuel.thibault@ens-lyon.org> | 2016-10-09 19:04:57 +0200 |
|---|---|---|
| committer | Samuel Thibault <samuel.thibault@ens-lyon.org> | 2016-10-09 19:04:57 +0200 |
| commit | 7bb5f8a836b916d6ebf7b6921b136e99cea2442d (patch) | |
| tree | 28a7ed786dae726ad14f100e8626eee872b1ba11 /sysdeps/ieee754/ldbl-96 | |
| parent | f76453c31593957fec1a99b986bfa5506618b79c (diff) | |
| parent | ab30899d880f9741a409cbc0d7a28399bdac21bf (diff) | |
Merge tag 'glibc-2.23' into baseline
The GNU C Library
=================
The GNU C Library version 2.23 is now available.
The GNU C Library is used as *the* C library in the GNU system and
in GNU/Linux systems, as well as many other systems that use Linux
as the kernel.
The GNU C Library is primarily designed to be a portable
and high performance C library. It follows all relevant
standards including ISO C11 and POSIX.1-2008. It is also
internationalized and has one of the most complete
internationalization interfaces known.
The GNU C Library webpage is at http://www.gnu.org/software/libc/
Packages for the 2.23 release may be downloaded from:
http://ftpmirror.gnu.org/libc/
http://ftp.gnu.org/gnu/libc/
The mirror list is at http://www.gnu.org/order/ftp.html
NEWS for version 2.23
=====================
* Unicode 8.0.0 Support: Character encoding, character type info, and
transliteration tables are all updated to Unicode 8.0.0, using new
and/or improved generator scripts contributed by Mike FABIAN (Red Hat).
These updates cause user visible changes, such as the fixes for bugs
89, 16061, and 18568.
* sched_setaffinity, pthread_setaffinity_np no longer attempt to guess the
kernel-internal CPU set size. This means that requests that change the
CPU affinity which failed before (for example, an all-ones CPU mask) will
now succeed. Applications that need to determine the effective CPU
affinities need to call sched_getaffinity or pthread_getaffinity_np after
setting it because the kernel can adjust it (and the previous size check
would not detect this in the majority of cases).
* The fts.h header can now be used with -D_FILE_OFFSET_BITS=64. With LFS
the following new symbols are used: fts64_children, fts64_close,
fts64_open, fts64_read and fts64_set.
* getaddrinfo now detects certain invalid responses on an internal netlink
socket. If such responses are received, an affected process will
terminate with an error message of "Unexpected error <number> on netlink
descriptor <number>" or "Unexpected netlink response of size <number> on
descriptor <number>". The most likely cause for these errors is a
multi-threaded application which erroneously closes and reuses the netlink
file descriptor while it is used by getaddrinfo.
* A defect in the malloc implementation, present since glibc 2.15 (2012) or
glibc 2.10 via --enable-experimental-malloc (2009), could result in the
unnecessary serialization of memory allocation requests across threads.
The defect is now corrected. Users should see a substantial increase in
the concurent throughput of allocation requests for applications which
trigger this bug. Affected applications typically create create and
destroy threads frequently. (Bug 19048 was reported and analyzed by
Ericsson.)
* There is now a --disable-timezone-tools configure option for disabling the
building and installing of the timezone related utilities (zic, zdump, and
tzselect). This is useful for people who build the timezone data and code
independent of the GNU C Library.
* The obsolete header <regexp.h> has been removed. Programs that require
this header must be updated to use <regex.h> instead.
* The obsolete functions bdflush, create_module, get_kernel_syms,
query_module and uselib are no longer available to newly linked binaries;
the header <sys/kdaemon.h> has been removed. These functions and header
were specific to systems using the Linux kernel and could not usefully be
used with the GNU C Library on systems with version 2.6 or later of the
Linux kernel.
* Optimized string, wcsmbs and memory functions for IBM z13.
Implemented by Stefan Liebler.
* Newly linked programs that define a variable called signgam will no longer
have it set by the lgamma, lgammaf and lgammal functions. Programs that
require signgam to be set by those functions must ensure that they use the
variable provided by the GNU C Library and declared in <math.h>, without
defining their own copy.
* The minimum GCC version that can be used to build this version of the GNU
C Library is GCC 4.7. Older GCC versions, and non-GNU compilers, can
still be used to compile programs using the GNU C Library.
Security related changes:
* An out-of-bounds value in a broken-out struct tm argument to strftime no
longer causes a crash. Reported by Adam Nielsen. (CVE-2015-8776)
* The LD_POINTER_GUARD environment variable can no longer be used to disable
the pointer guard feature. It is always enabled. Previously,
LD_POINTER_GUARD could be used to disable security hardening in binaries
running in privileged AT_SECURE mode. Reported by Hector Marco-Gisbert.
(CVE-2015-8777)
* An integer overflow in hcreate and hcreate_r could lead to an
out-of-bounds memory access. Reported by Szabolcs Nagy. (CVE-2015-8778)
* The catopen function no longer has unbounded stack usage. Reported by
Max. (CVE-2015-8779)
* The nan, nanf and nanl functions no longer have unbounded stack usage
depending on the length of the string passed as an argument to the
functions. Reported by Joseph Myers. (CVE-2014-9761)
* A stack-based buffer overflow was found in libresolv when invoked from
libnss_dns, allowing specially crafted DNS responses to seize control
of execution flow in the DNS client. The buffer overflow occurs in
the functions send_dg (send datagram) and send_vc (send TCP) for the
NSS module libnss_dns.so.2 when calling getaddrinfo with AF_UNSPEC
family. The use of AF_UNSPEC triggers the low-level resolver code to
send out two parallel queries for A and AAAA. A mismanagement of the
buffers used for those queries could result in the response of a query
writing beyond the alloca allocated buffer created by
_nss_dns_gethostbyname4_r. Buffer management is simplified to remove
the overflow. Thanks to the Google Security Team and Red Hat for
reporting the security impact of this issue, and Robert Holiday of
Ciena for reporting the related bug 18665. (CVE-2015-7547)
The following bugs are resolved with this release:
[89] localedata: Locales nb_NO and nn_NO should transliterate æøå
[887] math: Math library function "logb" and "nextafter" inconsistent
[2542] math: Incorrect return from float gamma (-0X1.FA471547C2FE5P+1)
[2543] math: Incorrect return from float gamma (-0X1.9260DCP+1)
[2558] math: Incorrect return from double gamma (-0X1.FA471547C2FE5P+1)
[2898] libc: [improve] warning: the use of `mktemp' is dangerous, better
use `mkstemp'
[4404] localedata: German translation of "Alarm clock" is misleading
[6799] math: nextafter() and nexttoward() doen't set errno on
overflow/underflow errors
[6803] math: scalb(), scalbln(), scalbn() do not set errno on
overflow/underflow
[10432] nis: _nss_nis_setnetgrent assertion failure
[11460] libc: fts has no LFS support
[12926] network: getaddrinfo()/make_request() may spin forever
[13065] nptl: Race condition in pthread barriers
[13690] nptl: pthread_mutex_unlock potentially cause invalid access
[14341] dynamic-link: Dynamic linker crash when DT_JMPREL and DT_REL{,A}
are not contiguous
[14551] math: [ldbl-128ibm] strtold overflow handling for IBM long double
[14912] libc: Rename non-installed bits/*.h headers
[15002] libc: Avoid undefined behavior in posix_fallocate overflow check
[15367] math: Let gcc use __builtin_isinf
[15384] math: One constant fewer in ieee754/dbl-64/wordsize-64/s_finite.c
[15421] math: lgamma wrongly sets signgam for ISO C
[15470] math: [arm] On ARM llrintl() and llroundl() do not raise
FE_INVALID with argument out of range
[15491] math: [i386/x86_64] x86 nearbyint implementations wrongly clear
all exceptions
[15786] dynamic-link: ifunc resolver functions can smash function
arguments
[15918] math: Unnecessary check for equality in hypotf()
[16061] localedata: Review / update transliteration data
[16068] math: [i386/x86_64] x86 and x86_64 fesetenv exclude state they
should include
[16141] time: strptime %z offset restriction
[16171] math: drem should be alias of remainder
[16296] math: fegetround is pure?
[16347] math: [ldbl-128ibm] ldbl-128/e_lgammal_r.c may not be suitable.
[16364] libc: sleep may leave SIGCHLD blocked on sync cancellation on
GNU/Linux
[16399] math: [mips] lrint / llrint / lround / llround missing exceptions
[16415] math: Clean up ldbl-128 / ldbl-128ibm expm1l for large positive
arguments
[16422] math: [powerpc] math-float, math-double failing llrint tests with
"Exception "Inexact" set" on ppc32
[16495] localedata: nl_NL: date_fmt: shuffle year/month around
[16517] math: Missing underflow exception from tanf/tan/tanl
[16519] math: Missing underflow exception from sinhf
[16520] math: Missing underflow exception from tanhf
[16521] math: Missing underflow exception from exp2
[16620] math: [ldbl-128ibm] exp10l spurious overflows / bad directed
rounding results
[16734] stdio: fopen calls mmap to allocate its buffer
[16961] math: nan function incorrect handling of bad sequences
[16962] math: nan function unbounded stack allocation (CVE-2014-9761)
[16973] localedata: Fix lang_lib/lang_term as per ISO 639-2
[16985] locale: localedef: confusing error message when opening output
fails
[17118] math: ctanh(INFINITY + 2 * I) returns incorrect value
[17197] locale: Redundant shift character in iconv conversion output at
block boundary
[17243] libc: trunk/posix/execl.c:53: va_args problem ?
[17244] libc: trunk/sysdeps/unix/sysv/linux/semctl.c:116: va_args muxup ?
[17250] dynamic-link: static linking breaks nss loading
(getaddrinfo/getpwnam/etc...)
[17404] libc: atomic_exchange_rel lacking a barrier on MIPS16, GCC before
4.7?
[17441] math: isnan() should use __builtin_isnan() in GCC
[17514] nptl: Assert failure unlocking ERRORCHECK mutex after timedlock
(related to lock elision)
[17787] manual: Exponent on page 324 of the PDF ends prematurely
[17886] time: strptime should be able to parse "Z" as a timezone with %z
[17887] time: strptime should be able to parse "+01:00" style timezones
[17905] libc: catopen() Multiple unbounded stack allocations
(CVE-2015-8779)
[18084] libc: backtrace (..., 0) dumps core on x86
[18086] libc: nice() sets errno to 0 on success
[18240] libc: hcreate, hcreate_r should fail with ENOMEM if element count
is too large (CVE-2015-8778)
[18251] dynamic-link: SONAME missing when audit modules provides path
[18265] libc: add attributes for wchar string and memory functions
[18370] math: csqrt missing underflows
[18421] libc: [hppa] read-only segment has dynamic relocations
[18472] libc: Obsolete syscall wrappers should be compat symbols
[18480] libc: hppa glibc miscompilation in sched_setaffinity()
[18491] localedata: Update tr_TR LC_CTYPE as part of Unicode updates
[18525] localedata: Remove locale timezone information
[18560] libc: [powerpc] spurious bits/ipc.h definitions
[18568] localedata: Update locale data to Unicode 8.0
[18589] locale: sort-test.sh fails at random
[18595] math: ctan, ctanh missing underflows
[18604] libc: assert macro-expands its argument
[18610] math: S390: fetestexcept() reports any exception if DXC-code
contains a vector instruction exception.
[18611] math: j1, jn missing errno setting on underflow
[18618] localedata: sync Chechen locale definitions with other *_RU
locales
[18647] math: powf(-0x1.000002p0, 0x1p30) returns 0 instead of +inf
[18661] libc: Some x86-64 assembly codes don't align stack to 16 bytes
[18665] network: In send_dg, the recvfrom function is NOT always using the
buffer size of a newly created buffer (CVE-2015-7547)
[18674] libc: [i386] trunk/sysdeps/i386/tst-auditmod3b.c:84: possible
missing break ?
[18675] libc: fpathconf(_PC_NAME_MAX) fails against large filesystems for
32bit processes
[18681] libc: regexp.h is obsolete and buggy, and should be desupported
[18699] math: tilegx cproj() for various complex infinities does not yield
infinity
[18724] libc: Harden put*ent functions against data injection
[18743] nptl: PowerPC: findutils testcase fails with --enable-lock-elision
[18755] build: build errors with -DNDEBUG
[18757] stdio: fmemopen fails to set errno on failure
[18778] dynamic-link: ld.so crashes if failed dlopen causes libpthread to
be forced unloaded
[18781] libc: openat64 lacks O_LARGEFILE
[18787] libc: [hppa] sysdeps/unix/sysv/linux/hppa/bits/atomic.h:71:6:
error: can’t find a register in class ‘R1_REGS’ while reloading ‘asm’
[18789] math: [ldbl-128ibm] sinhl inaccurate near 0
[18790] math: [ldbl-128ibm] tanhl inaccurate
[18795] libc: stpncpy fortification misses buffer lengths that are
statically too large
[18796] build: build fails for --disable-mathvec
[18803] math: hypot missing underflows
[18820] stdio: fmemopen may leak memory on failure
[18823] math: csqrt spurious underflows
[18824] math: fma spurious underflows
[18825] math: pow missing underflows
[18857] math: [ldbl-128ibm] nearbyintl wrongly uses signaling comparisons
[18868] nptl: pthread_barrier_init typo has in-theory-undefined behavior
[18870] build: sem_open.c fails to compile with missing symbol
FUTEX_SHARED
[18872] stdio: Fix memory leak in printf_positional
[18873] libc: posix_fallocate overflow check ineffective
[18875] math: Excess precision leads incorrect libm
[18877] libc: arm: mmap offset regression
[18887] libc: memory corruption when using getmntent on blank lines
[18918] localedata: hu_HU: change time to HH:MM:SS format
[18921] libc: Regression: extraneous stat() and fstat() performed by
opendir()
[18928] dynamic-link: LD_POINTER_GUARD is not ignored for privileged
binaries (CVE-2015-8777)
[18951] math: tgamma missing underflows
[18952] math: [ldbl-128/ldbl-128ibm] lgammal spurious "invalid", incorrect
signgam
[18953] localedata: lt_LT: change currency symbol to the euro
[18956] math: powf inaccuracy
[18961] math: [i386] exp missing underflows
[18966] math: [i386] exp10 missing underflows
[18967] math: math.h XSI POSIX namespace (gamma, isnan, scalb)
[18969] build: multiple string test failures due to missing locale
dependencies
[18970] libc: Reference of pthread_setcancelstate in libc.a
[18977] math: float / long double Bessel functions not in XSI POSIX
[18980] math: i386 libm functions return with excess range and precision
[18981] math: i386 scalb*, ldexp return with excess range and precision
[18982] stdio: va_list and vprintf
[18985] time: Passing out of range data to strftime() causes a segfault
(CVE-2015-8776)
[19003] math: [x86_64] fma4 version of pow inappropriate contraction
[19007] libc: FAIL: elf/check-localplt with -z now and binutils 2.26
[19012] locale: iconv_open leaks memory on error path
[19016] math: clog, clog10 inaccuracy
[19018] nptl: Mangle function pointers in tls_dtor_list
[19032] math: [i386] acosh (-qNaN) spurious "invalid" exception
[19046] math: ldbl-128 / ldbl-128ibm lgamma bad overflow handling
[19048] malloc: malloc: arena free list can become cyclic, increasing
contention
[19049] math: [powerpc] erfc incorrect zero sign
[19050] math: [powerpc] log* incorrect zero sign
[19058] math: [x86_64] Link fail with -fopenmp and -flto
[19059] math: nexttoward overflow incorrect in non-default rounding modes
[19071] math: ldbl-96 lroundl incorrect just below powers of 2
[19074] network: Data race in _res_hconf_reorder_addrs
[19076] math: [ldbl-128ibm] log1pl (-1) wrong sign of infinity
[19077] math: [ldbl-128ibm] logl (1) incorrect sign of zero result
[19078] math: [ldbl-128ibm] expl overflow incorrect in non-default
rounding modes
[19079] math: dbl-64/wordsize-64 lround based on llround incorrect for
ILP32
[19085] math: ldbl-128 lrintl, lroundl missing exceptions for 32-bit long
[19086] manual: posix_fallocate64 documented argument order is wrong.
[19088] math: lround, llround missing exceptions close to overflow
threshold
[19094] math: lrint, llrint missing exceptions close to overflow threshold
[19095] math: dbl-64 lrint incorrect for 64-bit long
[19122] dynamic-link: Unnecessary PLT relocations in librtld.os
[19124] dynamic-link: ld.so failed to build with older assmebler
[19125] math: [powerpc32] llroundf, llround incorrect exceptions
[19129] dynamic-link: [arm] Concurrent lazy TLSDESC resolution can crash
[19134] math: [powerpc32] lround, lroundf spurious exceptions
[19137] libc: i386/epoll_pwait.S doesn't support cancellation
[19143] nptl: Remove CPU set size checking from sched_setaffinity,
pthread_setaffinity_np
[19156] math: [ldbl-128] j0l spurious underflows
[19164] nptl: tst-getcpu fails with many possible CPUs
[19168] math: math/test-ildoubl and math/test-ldouble failure
[19174] nptl: PowerPC: TLE enabled pthread mutex performs poorly.
[19178] dynamic-link: ELF_RTYPE_CLASS_EXTERN_PROTECTED_DATA confuses
prelink
[19181] math: [i386/x86_64] fesetenv (FE_DFL_ENV), fesetenv
(FE_NOMASK_ENV) do not clear SSE exceptions
[19182] malloc: malloc deadlock between ptmalloc_lock_all and
_int_new_arena/reused_arena
[19189] math: [ldbl-128] log1pl (-qNaN) spurious "invalid" exception
[19201] math: dbl-64 remainder incorrect sign of zero result
[19205] math: bits/math-finite.h conditions do not match math.h and
bits/mathcalls.h
[19209] math: bits/math-finite.h wrongly maps ldexp to scalbn
[19211] math: lgamma functions do not set signgam for -ffinite-math-only
for C99-based standards
[19212] libc: features.h not -Wundef clean
[19213] math: [i386/x86_64] log* (1) incorrect zero sign for -ffinite-
math-only
[19214] libc: Family and model identification for AMD CPU's are incorrect.
[19219] libc: GLIBC build fails for ia64 with missing __nearbyintl
[19228] math: [powerpc] nearbyint wrongly clears "inexact", leaves traps
disabled
[19235] math: [powerpc64] lround, lroundf, llround, llroundf spurious
"inexact" exceptions
[19238] math: [powerpc] round, roundf spurious "inexact" for integer
arguments
[19242] libc: strtol incorrect in Turkish locales
[19243] malloc: reused_arena can pick an arena on the free list, leading
to an assertion failure and reference count corruption
[19253] time: tzset() ineffective when temporary TZ did not include DST
rules
[19266] math: strtod ("NAN(I)") incorrect in Turkish locales
[19270] math: [hppa] Shared libm missing __isnanl
[19285] libc: [hppa] sysdeps/unix/sysv/linux/hppa/bits/mman.h: missing
MAP_HUGETLB and MAP_STACK defines
[19313] nptl: Wrong __cpu_mask for x32
[19347] libc: grantpt: try to force a specific gid even without pt_chown
[19349] math: [ldbl-128ibm] tanhl inaccurate for small arguments
[19350] math: [ldbl-128ibm] sinhl spurious overflows
[19351] math: [ldbl-128ibm] logl inaccurate near 1
[19363] time: x32: times() return value wrongly truncates/sign extends
from 32bit
[19367] dynamic-link: Improve branch prediction on Silvermont
[19369] network: Default domain name not reset by res_ninit when "search"
/ "domain" entry is removed from resolv.conf
[19375] math: powerpc: incorrect results for POWER7 logb with negative
subnormals
[19385] localedata: bg_BG: time separator should be colon, not comma
[19408] libc: linux personality syscall wrapper may erroneously return an
error on 32-bit architectures
[19415] libc: dladdr returns wrong names on hppa
[19432] libc: iconv rejects redundant escape sequences in IBM900, IBM903,
IBM905, IBM907, and IBM909
[19439] math: Unix98 isinf and isnan functions conflict with C++11
[19443] build: build failures with -DDEBUG
[19451] build: Make check fails on test-double-vlen2
[19462] libc: Glibc failed to build with -Os
[19465] math: Wrong code with -Os
[19466] time: time/tst-mktime2.c is compiled into an infinite loop with
-Os
[19467] string: Fast_Unaligned_Load needs to be enabled for Excavator core
CPU's.
[19475] libc: Glibc 2.22 doesn't build on sparc [PATCH]
[19486] math: S390: Math tests fail with "Exception Inexact set".
[19529] libc: [ARM]: FAIL: stdlib/tst-makecontext
[19550] libc: [mips] mmap negative offset handling inconsistent with other
architectures
[19590] math: Fail to build shared objects that use libmvec.so functions.
Contributors
============
This release was made possible by the contributions of many people.
The maintainers are grateful to everyone who has contributed
changes or bug reports. These include:
Adhemerval Zanella
Alan Modra
Amit Pawar
Andreas Schwab
Andrew Bennett
Andrew Senkevich
Andrew Stubbs
Anton Blanchard
Arjun Shankar
Arslanbek Astemirov
Aurelien Jarno
Brett Neumeier
Carlos Eduardo Seo
Carlos O'Donell
Chris Metcalf
Chung-Lin Tang
Damyan Ivanov
Daniel Marjamäki
David Kastrup
David Lamparter
David S. Miller
Dmitry V. Levin
Egmont Koblinger
Evert
Flavio Cruz
Florian Weimer
Gabriel F. T. Gomes
Geoffrey Thomas
Gleb Fotengauer-Malinovskiy
Gunnar Hjalmarsson
H.J. Lu
Helge Deller
James Perkins
John David Anglin
Joseph Myers
Justus Winter
Khem Raj
Ludovic Courtès
Maciej W. Rozycki
Manolis Ragkousis
Marcin Kościelnicki
Mark Wielaard
Marko Myllynen
Martin Sebor
Maxim Ostapenko
Mike FABIAN
Mike Frysinger
Namhyung Kim
Ondrej Bilka
Ondřej Bílka
Paul E. Murphy
Paul Eggert
Paul Murphy
Paul Pluzhnikov
Petar Jovanovic
Phil Blundell
Rajalakshmi Srinivasaraghavan
Rasmus Villemoes
Richard Henderson
Rob Wu
Roland McGrath
Samuel Thibault
Siddhesh Poyarekar
Stan Shebs
Stefan Liebler
Steve Ellcey
Szabolcs Nagy
Thomas Schwinge
Torvald Riegel
Tulio Magno Quites Machado Filho
Vincent Bernat
Wilco Dijkstra
Zack Weinberg
Diffstat (limited to 'sysdeps/ieee754/ldbl-96')
44 files changed, 785 insertions, 152 deletions
diff --git a/sysdeps/ieee754/ldbl-96/e_asinl.c b/sysdeps/ieee754/ldbl-96/e_asinl.c index 2973bf071b..f52b931459 100644 --- a/sysdeps/ieee754/ldbl-96/e_asinl.c +++ b/sysdeps/ieee754/ldbl-96/e_asinl.c @@ -112,11 +112,7 @@ __ieee754_asinl (long double x) { /* |x|<0.5 */ if (ix < 0x3fde8000) { /* if |x| < 2**-33 */ - if (fabsl (x) < LDBL_MIN) - { - long double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); if (huge + x > one) return x; /* return x with inexact if x!=0 */ } diff --git a/sysdeps/ieee754/ldbl-96/e_atanhl.c b/sysdeps/ieee754/ldbl-96/e_atanhl.c index 9a957c9065..b99a83c6ee 100644 --- a/sysdeps/ieee754/ldbl-96/e_atanhl.c +++ b/sysdeps/ieee754/ldbl-96/e_atanhl.c @@ -55,11 +55,7 @@ __ieee754_atanhl(long double x) return x/zero; if(ix<0x3fdf) { math_force_eval(huge+x); - if (fabsl (x) < LDBL_MIN) - { - long double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); return x; /* x<2**-32 */ } SET_LDOUBLE_EXP(x,ix); diff --git a/sysdeps/ieee754/ldbl-96/e_gammal_r.c b/sysdeps/ieee754/ldbl-96/e_gammal_r.c index 9da5db33f0..8dd7a03918 100644 --- a/sysdeps/ieee754/ldbl-96/e_gammal_r.c +++ b/sysdeps/ieee754/ldbl-96/e_gammal_r.c @@ -1,5 +1,5 @@ /* Implementation of gamma function according to ISO C. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -186,6 +186,7 @@ __ieee754_gammal_r (long double x, int *signgamp) ret = M_PIl / (-x * sinpix * gammal_positive (-x, &exp2_adj)); ret = __scalbnl (ret, -exp2_adj); + math_check_force_underflow_nonneg (ret); } } } diff --git a/sysdeps/ieee754/ldbl-96/e_hypotl.c b/sysdeps/ieee754/ldbl-96/e_hypotl.c index d3152f91e5..ee3a07055b 100644 --- a/sysdeps/ieee754/ldbl-96/e_hypotl.c +++ b/sysdeps/ieee754/ldbl-96/e_hypotl.c @@ -132,7 +132,9 @@ long double __ieee754_hypotl(long double x, long double y) t1 = 1.0; GET_LDOUBLE_EXP(exp,t1); SET_LDOUBLE_EXP(t1,exp+k); - return t1*w; + w *= t1; + math_check_force_underflow_nonneg (w); + return w; } else return w; } strong_alias (__ieee754_hypotl, __hypotl_finite) diff --git a/sysdeps/ieee754/ldbl-96/e_j1l.c b/sysdeps/ieee754/ldbl-96/e_j1l.c index 46e28dfe11..e8a7349cf4 100644 --- a/sysdeps/ieee754/ldbl-96/e_j1l.c +++ b/sysdeps/ieee754/ldbl-96/e_j1l.c @@ -154,11 +154,9 @@ __ieee754_j1l (long double x) if (huge + x > one) /* inexact if x!=0 necessary */ { long double ret = 0.5 * x; - if (fabsl (ret) < LDBL_MIN) - { - long double force_underflow = ret * ret; - math_force_eval (force_underflow); - } + math_check_force_underflow (ret); + if (ret == 0 && x != 0) + __set_errno (ERANGE); return ret; } } diff --git a/sysdeps/ieee754/ldbl-96/e_jnl.c b/sysdeps/ieee754/ldbl-96/e_jnl.c index 2f3a452f55..92f96921a7 100644 --- a/sysdeps/ieee754/ldbl-96/e_jnl.c +++ b/sysdeps/ieee754/ldbl-96/e_jnl.c @@ -289,12 +289,12 @@ __ieee754_jnl (int n, long double x) ret = b; } if (ret == 0) - ret = __copysignl (LDBL_MIN, ret) * LDBL_MIN; - else if (fabsl (ret) < LDBL_MIN) { - long double force_underflow = ret * ret; - math_force_eval (force_underflow); + ret = __copysignl (LDBL_MIN, ret) * LDBL_MIN; + __set_errno (ERANGE); } + else + math_check_force_underflow (ret); return ret; } strong_alias (__ieee754_jnl, __jnl_finite) diff --git a/sysdeps/ieee754/ldbl-96/e_lgammal_r.c b/sysdeps/ieee754/ldbl-96/e_lgammal_r.c index 0cc35f9252..9862fe8d5c 100644 --- a/sysdeps/ieee754/ldbl-96/e_lgammal_r.c +++ b/sysdeps/ieee754/ldbl-96/e_lgammal_r.c @@ -306,6 +306,8 @@ __ieee754_lgammal_r (long double x, int *signgamp) } if (se & 0x8000) { + if (x < -2.0L && x > -33.0L) + return __lgamma_negl (x, signgamp); t = sin_pi (x); if (t == zero) return one / fabsl (t); /* -integer */ @@ -428,11 +430,7 @@ __ieee754_lgammal_r (long double x, int *signgamp) in warnings that it may be used uninitialized although in the cases where it is used it has always been set. */ DIAG_PUSH_NEEDS_COMMENT; -#if __GNUC_PREREQ (4, 7) DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wmaybe-uninitialized"); -#else - DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wuninitialized"); -#endif if (se & 0x8000) r = nadj - r; DIAG_POP_NEEDS_COMMENT; diff --git a/sysdeps/ieee754/ldbl-96/e_rem_pio2l.c b/sysdeps/ieee754/ldbl-96/e_rem_pio2l.c index d05afce623..0d8e64675c 100644 --- a/sysdeps/ieee754/ldbl-96/e_rem_pio2l.c +++ b/sysdeps/ieee754/ldbl-96/e_rem_pio2l.c @@ -1,5 +1,5 @@ /* Extended-precision floating point argument reduction. - Copyright (C) 1999-2015 Free Software Foundation, Inc. + Copyright (C) 1999-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Based on quad-precision code by Jakub Jelinek <jj@ultra.linux.cz> diff --git a/sysdeps/ieee754/ldbl-96/e_sinhl.c b/sysdeps/ieee754/ldbl-96/e_sinhl.c index 4978f348bb..095b142621 100644 --- a/sysdeps/ieee754/ldbl-96/e_sinhl.c +++ b/sysdeps/ieee754/ldbl-96/e_sinhl.c @@ -36,6 +36,7 @@ static char rcsid[] = "$NetBSD: $"; * only sinhl(0)=0 is exact for finite x. */ +#include <float.h> #include <math.h> #include <math_private.h> @@ -58,8 +59,10 @@ __ieee754_sinhl(long double x) if (jx & 0x8000) h = -h; /* |x| in [0,25], return sign(x)*0.5*(E+E/(E+1))) */ if (ix < 0x4003 || (ix == 0x4003 && i0 <= 0xc8000000)) { /* |x|<25 */ - if (ix<0x3fdf) /* |x|<2**-32 */ + if (ix<0x3fdf) { /* |x|<2**-32 */ + math_check_force_underflow (x); if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */ + } t = __expm1l(fabsl(x)); if(ix<0x3fff) return h*(2.0*t-t*t/(t+one)); return h*(t+t/(t+one)); diff --git a/sysdeps/ieee754/ldbl-96/gamma_product.c b/sysdeps/ieee754/ldbl-96/gamma_product.c index e113e23723..419d11598f 100644 --- a/sysdeps/ieee754/ldbl-96/gamma_product.c +++ b/sysdeps/ieee754/ldbl-96/gamma_product.c @@ -1,5 +1,5 @@ /* Compute a product of X, X+1, ..., with an error estimate. - Copyright (C) 2013-2015 Free Software Foundation, Inc. + Copyright (C) 2013-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -36,10 +36,7 @@ __gamma_product (double x, double x_eps, int n, double *eps) for (int i = 1; i < n; i++) ret *= x_full + i; -#if FLT_EVAL_METHOD != 0 - volatile -#endif - double fret = ret; + double fret = math_narrow_eval ((double) ret); *eps = (ret - fret) / fret; return fret; diff --git a/sysdeps/ieee754/ldbl-96/gamma_productl.c b/sysdeps/ieee754/ldbl-96/gamma_productl.c index 32990b8f1c..849b57d95d 100644 --- a/sysdeps/ieee754/ldbl-96/gamma_productl.c +++ b/sysdeps/ieee754/ldbl-96/gamma_productl.c @@ -1,5 +1,5 @@ /* Compute a product of X, X+1, ..., with an error estimate. - Copyright (C) 2013-2015 Free Software Foundation, Inc. + Copyright (C) 2013-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or diff --git a/sysdeps/ieee754/ldbl-96/k_cosl.c b/sysdeps/ieee754/ldbl-96/k_cosl.c index 95d64bd014..08b11b3733 100644 --- a/sysdeps/ieee754/ldbl-96/k_cosl.c +++ b/sysdeps/ieee754/ldbl-96/k_cosl.c @@ -1,5 +1,5 @@ /* Extended-precision floating point cosine on <-pi/4,pi/4>. - Copyright (C) 1999-2015 Free Software Foundation, Inc. + Copyright (C) 1999-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Based on quad-precision cosine by Jakub Jelinek <jj@ultra.linux.cz> diff --git a/sysdeps/ieee754/ldbl-96/k_sinl.c b/sysdeps/ieee754/ldbl-96/k_sinl.c index b7b5ae359d..6ba7ceddc0 100644 --- a/sysdeps/ieee754/ldbl-96/k_sinl.c +++ b/sysdeps/ieee754/ldbl-96/k_sinl.c @@ -1,5 +1,5 @@ /* Quad-precision floating point sine on <-pi/4,pi/4>. - Copyright (C) 1999-2015 Free Software Foundation, Inc. + Copyright (C) 1999-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Based on quad-precision sine by Jakub Jelinek <jj@ultra.linux.cz> @@ -96,11 +96,7 @@ __kernel_sinl(long double x, long double y, int iy) polynomial of degree 17. */ if (absx < 0x1p-33L) { - if (fabsl (x) < LDBL_MIN) - { - long double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); if (!((int)x)) return x; /* generate inexact */ } z = x * x; diff --git a/sysdeps/ieee754/ldbl-96/k_tanl.c b/sysdeps/ieee754/ldbl-96/k_tanl.c index 31cd236aa2..0c050c112f 100644 --- a/sysdeps/ieee754/ldbl-96/k_tanl.c +++ b/sysdeps/ieee754/ldbl-96/k_tanl.c @@ -56,6 +56,8 @@ * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y))) */ +#include <float.h> +#include <libc-internal.h> #include <math.h> #include <math_private.h> static const long double @@ -94,8 +96,13 @@ __kernel_tanl (long double x, long double y, int iy) { /* generate inexact */ if (x == 0 && iy == -1) return one / fabsl (x); + else if (iy == 1) + { + math_check_force_underflow_nonneg (absx); + return x; + } else - return (iy == 1) ? x : -one / x; + return -one / x; } } if (absx >= 0.6743316650390625L) @@ -126,8 +133,15 @@ __kernel_tanl (long double x, long double y, int iy) { v = (long double) iy; w = (v - 2.0 * (x - (w * w / (w + v) - r))); + /* SIGN is set for arguments that reach this code, but not + otherwise, resulting in warnings that it may be used + uninitialized although in the cases where it is used it has + always been set. */ + DIAG_PUSH_NEEDS_COMMENT; + DIAG_IGNORE_NEEDS_COMMENT (4.8, "-Wmaybe-uninitialized"); if (sign < 0) w = -w; + DIAG_POP_NEEDS_COMMENT; return w; } if (iy == 1) diff --git a/sysdeps/ieee754/ldbl-96/ldbl2mpn.c b/sysdeps/ieee754/ldbl-96/ldbl2mpn.c index 26efea18f6..fe7002f640 100644 --- a/sysdeps/ieee754/ldbl-96/ldbl2mpn.c +++ b/sysdeps/ieee754/ldbl-96/ldbl2mpn.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1995-2015 Free Software Foundation, Inc. +/* Copyright (C) 1995-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or diff --git a/sysdeps/ieee754/ldbl-96/lgamma_negl.c b/sysdeps/ieee754/ldbl-96/lgamma_negl.c new file mode 100644 index 0000000000..99ecf7e85f --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/lgamma_negl.c @@ -0,0 +1,418 @@ +/* lgammal expanding around zeros. + Copyright (C) 2015-2016 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <float.h> +#include <math.h> +#include <math_private.h> + +static const long double lgamma_zeros[][2] = + { + { -0x2.74ff92c01f0d82acp+0L, 0x1.360cea0e5f8ed3ccp-68L }, + { -0x2.bf6821437b201978p+0L, -0x1.95a4b4641eaebf4cp-64L }, + { -0x3.24c1b793cb35efb8p+0L, -0xb.e699ad3d9ba6545p-68L }, + { -0x3.f48e2a8f85fca17p+0L, -0xd.4561291236cc321p-68L }, + { -0x4.0a139e16656030cp+0L, -0x3.9f0b0de18112ac18p-64L }, + { -0x4.fdd5de9bbabf351p+0L, -0xd.0aa4076988501d8p-68L }, + { -0x5.021a95fc2db64328p+0L, -0x2.4c56e595394decc8p-64L }, + { -0x5.ffa4bd647d0357ep+0L, 0x2.b129d342ce12071cp-64L }, + { -0x6.005ac9625f233b6p+0L, -0x7.c2d96d16385cb868p-68L }, + { -0x6.fff2fddae1bbff4p+0L, 0x2.9d949a3dc02de0cp-64L }, + { -0x7.000cff7b7f87adf8p+0L, 0x3.b7d23246787d54d8p-64L }, + { -0x7.fffe5fe05673c3c8p+0L, -0x2.9e82b522b0ca9d3p-64L }, + { -0x8.0001a01459fc9f6p+0L, -0xc.b3cec1cec857667p-68L }, + { -0x8.ffffd1c425e81p+0L, 0x3.79b16a8b6da6181cp-64L }, + { -0x9.00002e3bb47d86dp+0L, -0x6.d843fedc351deb78p-64L }, + { -0x9.fffffb606bdfdcdp+0L, -0x6.2ae77a50547c69dp-68L }, + { -0xa.0000049f93bb992p+0L, -0x7.b45d95e15441e03p-64L }, + { -0xa.ffffff9466e9f1bp+0L, -0x3.6dacd2adbd18d05cp-64L }, + { -0xb.0000006b9915316p+0L, 0x2.69a590015bf1b414p-64L }, + { -0xb.fffffff70893874p+0L, 0x7.821be533c2c36878p-64L }, + { -0xc.00000008f76c773p+0L, -0x1.567c0f0250f38792p-64L }, + { -0xc.ffffffff4f6dcf6p+0L, -0x1.7f97a5ffc757d548p-64L }, + { -0xd.00000000b09230ap+0L, 0x3.f997c22e46fc1c9p-64L }, + { -0xd.fffffffff36345bp+0L, 0x4.61e7b5c1f62ee89p-64L }, + { -0xe.000000000c9cba5p+0L, -0x4.5e94e75ec5718f78p-64L }, + { -0xe.ffffffffff28c06p+0L, -0xc.6604ef30371f89dp-68L }, + { -0xf.0000000000d73fap+0L, 0xc.6642f1bdf07a161p-68L }, + { -0xf.fffffffffff28cp+0L, -0x6.0c6621f512e72e5p-64L }, + { -0x1.000000000000d74p+4L, 0x6.0c6625ebdb406c48p-64L }, + { -0x1.0ffffffffffff356p+4L, -0x9.c47e7a93e1c46a1p-64L }, + { -0x1.1000000000000caap+4L, 0x9.c47e7a97778935ap-64L }, + { -0x1.1fffffffffffff4cp+4L, 0x1.3c31dcbecd2f74d4p-64L }, + { -0x1.20000000000000b4p+4L, -0x1.3c31dcbeca4c3b3p-64L }, + { -0x1.2ffffffffffffff6p+4L, -0x8.5b25cbf5f545ceep-64L }, + { -0x1.300000000000000ap+4L, 0x8.5b25cbf5f547e48p-64L }, + { -0x1.4p+4L, 0x7.950ae90080894298p-64L }, + { -0x1.4p+4L, -0x7.950ae9008089414p-64L }, + { -0x1.5p+4L, 0x5.c6e3bdb73d5c63p-68L }, + { -0x1.5p+4L, -0x5.c6e3bdb73d5c62f8p-68L }, + { -0x1.6p+4L, 0x4.338e5b6dfe14a518p-72L }, + { -0x1.6p+4L, -0x4.338e5b6dfe14a51p-72L }, + { -0x1.7p+4L, 0x2.ec368262c7033b3p-76L }, + { -0x1.7p+4L, -0x2.ec368262c7033b3p-76L }, + { -0x1.8p+4L, 0x1.f2cf01972f577ccap-80L }, + { -0x1.8p+4L, -0x1.f2cf01972f577ccap-80L }, + { -0x1.9p+4L, 0x1.3f3ccdd165fa8d4ep-84L }, + { -0x1.9p+4L, -0x1.3f3ccdd165fa8d4ep-84L }, + { -0x1.ap+4L, 0xc.4742fe35272cd1cp-92L }, + { -0x1.ap+4L, -0xc.4742fe35272cd1cp-92L }, + { -0x1.bp+4L, 0x7.46ac70b733a8c828p-96L }, + { -0x1.bp+4L, -0x7.46ac70b733a8c828p-96L }, + { -0x1.cp+4L, 0x4.2862898d42174ddp-100L }, + { -0x1.cp+4L, -0x4.2862898d42174ddp-100L }, + { -0x1.dp+4L, 0x2.4b3f31686b15af58p-104L }, + { -0x1.dp+4L, -0x2.4b3f31686b15af58p-104L }, + { -0x1.ep+4L, 0x1.3932c5047d60e60cp-108L }, + { -0x1.ep+4L, -0x1.3932c5047d60e60cp-108L }, + { -0x1.fp+4L, 0xa.1a6973c1fade217p-116L }, + { -0x1.fp+4L, -0xa.1a6973c1fade217p-116L }, + { -0x2p+4L, 0x5.0d34b9e0fd6f10b8p-120L }, + { -0x2p+4L, -0x5.0d34b9e0fd6f10b8p-120L }, + { -0x2.1p+4L, 0x2.73024a9ba1aa36a8p-124L }, + }; + +static const long double e_hi = 0x2.b7e151628aed2a6cp+0L; +static const long double e_lo = -0x1.408ea77f630b0c38p-64L; + +/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's + approximation to lgamma function. */ + +static const long double lgamma_coeff[] = + { + 0x1.5555555555555556p-4L, + -0xb.60b60b60b60b60bp-12L, + 0x3.4034034034034034p-12L, + -0x2.7027027027027028p-12L, + 0x3.72a3c5631fe46aep-12L, + -0x7.daac36664f1f208p-12L, + 0x1.a41a41a41a41a41ap-8L, + -0x7.90a1b2c3d4e5f708p-8L, + 0x2.dfd2c703c0cfff44p-4L, + -0x1.6476701181f39edcp+0L, + 0xd.672219167002d3ap+0L, + -0x9.cd9292e6660d55bp+4L, + 0x8.911a740da740da7p+8L, + -0x8.d0cc570e255bf5ap+12L, + 0xa.8d1044d3708d1c2p+16L, + -0xe.8844d8a169abbc4p+20L, + }; + +#define NCOEFF (sizeof (lgamma_coeff) / sizeof (lgamma_coeff[0])) + +/* Polynomial approximations to (|gamma(x)|-1)(x-n)/(x-x0), where n is + the integer end-point of the half-integer interval containing x and + x0 is the zero of lgamma in that half-integer interval. Each + polynomial is expressed in terms of x-xm, where xm is the midpoint + of the interval for which the polynomial applies. */ + +static const long double poly_coeff[] = + { + /* Interval [-2.125, -2] (polynomial degree 13). */ + -0x1.0b71c5c54d42eb6cp+0L, + -0xc.73a1dc05f349517p-4L, + -0x1.ec841408528b6baep-4L, + -0xe.37c9da26fc3b492p-4L, + -0x1.03cd87c5178991ap-4L, + -0xe.ae9ada65ece2f39p-4L, + 0x9.b1185505edac18dp-8L, + -0xe.f28c130b54d3cb2p-4L, + 0x2.6ec1666cf44a63bp-4L, + -0xf.57cb2774193bbd5p-4L, + 0x4.5ae64671a41b1c4p-4L, + -0xf.f48ea8b5bd3a7cep-4L, + 0x6.7d73788a8d30ef58p-4L, + -0x1.11e0e4b506bd272ep+0L, + /* Interval [-2.25, -2.125] (polynomial degree 13). */ + -0xf.2930890d7d675a8p-4L, + -0xc.a5cfde054eab5cdp-4L, + 0x3.9c9e0fdebb0676e4p-4L, + -0x1.02a5ad35605f0d8cp+0L, + 0x9.6e9b1185d0b92edp-4L, + -0x1.4d8332f3d6a3959p+0L, + 0x1.1c0c8cacd0ced3eap+0L, + -0x1.c9a6f592a67b1628p+0L, + 0x1.d7e9476f96aa4bd6p+0L, + -0x2.921cedb488bb3318p+0L, + 0x2.e8b3fd6ca193e4c8p+0L, + -0x3.cb69d9d6628e4a2p+0L, + 0x4.95f12c73b558638p+0L, + -0x5.d392d0b97c02ab6p+0L, + /* Interval [-2.375, -2.25] (polynomial degree 14). */ + -0xd.7d28d505d618122p-4L, + -0xe.69649a304098532p-4L, + 0xb.0d74a2827d055c5p-4L, + -0x1.924b09228531c00ep+0L, + 0x1.d49b12bccee4f888p+0L, + -0x3.0898bb7dbb21e458p+0L, + 0x4.207a6cad6fa10a2p+0L, + -0x6.39ee630b46093ad8p+0L, + 0x8.e2e25211a3fb5ccp+0L, + -0xd.0e85ccd8e79c08p+0L, + 0x1.2e45882bc17f9e16p+4L, + -0x1.b8b6e841815ff314p+4L, + 0x2.7ff8bf7504fa04dcp+4L, + -0x3.c192e9c903352974p+4L, + 0x5.8040b75f4ef07f98p+4L, + /* Interval [-2.5, -2.375] (polynomial degree 15). */ + -0xb.74ea1bcfff94b2cp-4L, + -0x1.2a82bd590c375384p+0L, + 0x1.88020f828b968634p+0L, + -0x3.32279f040eb80fa4p+0L, + 0x5.57ac825175943188p+0L, + -0x9.c2aedcfe10f129ep+0L, + 0x1.12c132f2df02881ep+4L, + -0x1.ea94e26c0b6ffa6p+4L, + 0x3.66b4a8bb0290013p+4L, + -0x6.0cf735e01f5990bp+4L, + 0xa.c10a8db7ae99343p+4L, + -0x1.31edb212b315feeap+8L, + 0x2.1f478592298b3ebp+8L, + -0x3.c546da5957ace6ccp+8L, + 0x7.0e3d2a02579ba4bp+8L, + -0xc.b1ea961c39302f8p+8L, + /* Interval [-2.625, -2.5] (polynomial degree 16). */ + -0x3.d10108c27ebafad4p-4L, + 0x1.cd557caff7d2b202p+0L, + 0x3.819b4856d3995034p+0L, + 0x6.8505cbad03dd3bd8p+0L, + 0xb.c1b2e653aa0b924p+0L, + 0x1.50a53a38f05f72d6p+4L, + 0x2.57ae00cbd06efb34p+4L, + 0x4.2b1563077a577e9p+4L, + 0x7.6989ed790138a7f8p+4L, + 0xd.2dd28417b4f8406p+4L, + 0x1.76e1b71f0710803ap+8L, + 0x2.9a7a096254ac032p+8L, + 0x4.a0e6109e2a039788p+8L, + 0x8.37ea17a93c877b2p+8L, + 0xe.9506a641143612bp+8L, + 0x1.b680ed4ea386d52p+12L, + 0x3.28a2130c8de0ae84p+12L, + /* Interval [-2.75, -2.625] (polynomial degree 15). */ + -0x6.b5d252a56e8a7548p-4L, + 0x1.28d60383da3ac72p+0L, + 0x1.db6513ada8a6703ap+0L, + 0x2.e217118f9d34aa7cp+0L, + 0x4.450112c5cbd6256p+0L, + 0x6.4af99151e972f92p+0L, + 0x9.2db598b5b183cd6p+0L, + 0xd.62bef9c9adcff6ap+0L, + 0x1.379f290d743d9774p+4L, + 0x1.c58271ff823caa26p+4L, + 0x2.93a871b87a06e73p+4L, + 0x3.bf9db66103d7ec98p+4L, + 0x5.73247c111fbf197p+4L, + 0x7.ec8b9973ba27d008p+4L, + 0xb.eca5f9619b39c03p+4L, + 0x1.18f2e46411c78b1cp+8L, + /* Interval [-2.875, -2.75] (polynomial degree 14). */ + -0x8.a41b1e4f36ff88ep-4L, + 0xc.da87d3b69dc0f34p-4L, + 0x1.1474ad5c36158ad2p+0L, + 0x1.761ecb90c5553996p+0L, + 0x1.d279bff9ae234f8p+0L, + 0x2.4e5d0055a16c5414p+0L, + 0x2.d57545a783902f8cp+0L, + 0x3.8514eec263aa9f98p+0L, + 0x4.5235e338245f6fe8p+0L, + 0x5.562b1ef200b256c8p+0L, + 0x6.8ec9782b93bd565p+0L, + 0x8.14baf4836483508p+0L, + 0x9.efaf35dc712ea79p+0L, + 0xc.8431f6a226507a9p+0L, + 0xf.80358289a768401p+0L, + /* Interval [-3, -2.875] (polynomial degree 13). */ + -0xa.046d667e468f3e4p-4L, + 0x9.70b88dcc006c216p-4L, + 0xa.a8a39421c86ce9p-4L, + 0xd.2f4d1363f321e89p-4L, + 0xd.ca9aa1a3ab2f438p-4L, + 0xf.cf09c31f05d02cbp-4L, + 0x1.04b133a195686a38p+0L, + 0x1.22b54799d0072024p+0L, + 0x1.2c5802b869a36ae8p+0L, + 0x1.4aadf23055d7105ep+0L, + 0x1.5794078dd45c55d6p+0L, + 0x1.7759069da18bcf0ap+0L, + 0x1.8e672cefa4623f34p+0L, + 0x1.b2acfa32c17145e6p+0L, + }; + +static const size_t poly_deg[] = + { + 13, + 13, + 14, + 15, + 16, + 15, + 14, + 13, + }; + +static const size_t poly_end[] = + { + 13, + 27, + 42, + 58, + 75, + 91, + 106, + 120, + }; + +/* Compute sin (pi * X) for -0.25 <= X <= 0.5. */ + +static long double +lg_sinpi (long double x) +{ + if (x <= 0.25L) + return __sinl (M_PIl * x); + else + return __cosl (M_PIl * (0.5L - x)); +} + +/* Compute cos (pi * X) for -0.25 <= X <= 0.5. */ + +static long double +lg_cospi (long double x) +{ + if (x <= 0.25L) + return __cosl (M_PIl * x); + else + return __sinl (M_PIl * (0.5L - x)); +} + +/* Compute cot (pi * X) for -0.25 <= X <= 0.5. */ + +static long double +lg_cotpi (long double x) +{ + return lg_cospi (x) / lg_sinpi (x); +} + +/* Compute lgamma of a negative argument -33 < X < -2, setting + *SIGNGAMP accordingly. */ + +long double +__lgamma_negl (long double x, int *signgamp) +{ + /* Determine the half-integer region X lies in, handle exact + integers and determine the sign of the result. */ + int i = __floorl (-2 * x); + if ((i & 1) == 0 && i == -2 * x) + return 1.0L / 0.0L; + long double xn = ((i & 1) == 0 ? -i / 2 : (-i - 1) / 2); + i -= 4; + *signgamp = ((i & 2) == 0 ? -1 : 1); + + SET_RESTORE_ROUNDL (FE_TONEAREST); + + /* Expand around the zero X0 = X0_HI + X0_LO. */ + long double x0_hi = lgamma_zeros[i][0], x0_lo = lgamma_zeros[i][1]; + long double xdiff = x - x0_hi - x0_lo; + + /* For arguments in the range -3 to -2, use polynomial + approximations to an adjusted version of the gamma function. */ + if (i < 2) + { + int j = __floorl (-8 * x) - 16; + long double xm = (-33 - 2 * j) * 0.0625L; + long double x_adj = x - xm; + size_t deg = poly_deg[j]; + size_t end = poly_end[j]; + long double g = poly_coeff[end]; + for (size_t j = 1; j <= deg; j++) + g = g * x_adj + poly_coeff[end - j]; + return __log1pl (g * xdiff / (x - xn)); + } + + /* The result we want is log (sinpi (X0) / sinpi (X)) + + log (gamma (1 - X0) / gamma (1 - X)). */ + long double x_idiff = fabsl (xn - x), x0_idiff = fabsl (xn - x0_hi - x0_lo); + long double log_sinpi_ratio; + if (x0_idiff < x_idiff * 0.5L) + /* Use log not log1p to avoid inaccuracy from log1p of arguments + close to -1. */ + log_sinpi_ratio = __ieee754_logl (lg_sinpi (x0_idiff) + / lg_sinpi (x_idiff)); + else + { + /* Use log1p not log to avoid inaccuracy from log of arguments + close to 1. X0DIFF2 has positive sign if X0 is further from + XN than X is from XN, negative sign otherwise. */ + long double x0diff2 = ((i & 1) == 0 ? xdiff : -xdiff) * 0.5L; + long double sx0d2 = lg_sinpi (x0diff2); + long double cx0d2 = lg_cospi (x0diff2); + log_sinpi_ratio = __log1pl (2 * sx0d2 + * (-sx0d2 + cx0d2 * lg_cotpi (x_idiff))); + } + + long double log_gamma_ratio; + long double y0 = 1 - x0_hi; + long double y0_eps = -x0_hi + (1 - y0) - x0_lo; + long double y = 1 - x; + long double y_eps = -x + (1 - y); + /* We now wish to compute LOG_GAMMA_RATIO + = log (gamma (Y0 + Y0_EPS) / gamma (Y + Y_EPS)). XDIFF + accurately approximates the difference Y0 + Y0_EPS - Y - + Y_EPS. Use Stirling's approximation. First, we may need to + adjust into the range where Stirling's approximation is + sufficiently accurate. */ + long double log_gamma_adj = 0; + if (i < 8) + { + int n_up = (9 - i) / 2; + long double ny0, ny0_eps, ny, ny_eps; + ny0 = y0 + n_up; + ny0_eps = y0 - (ny0 - n_up) + y0_eps; + y0 = ny0; + y0_eps = ny0_eps; + ny = y + n_up; + ny_eps = y - (ny - n_up) + y_eps; + y = ny; + y_eps = ny_eps; + long double prodm1 = __lgamma_productl (xdiff, y - n_up, y_eps, n_up); + log_gamma_adj = -__log1pl (prodm1); + } + long double log_gamma_high + = (xdiff * __log1pl ((y0 - e_hi - e_lo + y0_eps) / e_hi) + + (y - 0.5L + y_eps) * __log1pl (xdiff / y) + log_gamma_adj); + /* Compute the sum of (B_2k / 2k(2k-1))(Y0^-(2k-1) - Y^-(2k-1)). */ + long double y0r = 1 / y0, yr = 1 / y; + long double y0r2 = y0r * y0r, yr2 = yr * yr; + long double rdiff = -xdiff / (y * y0); + long double bterm[NCOEFF]; + long double dlast = rdiff, elast = rdiff * yr * (yr + y0r); + bterm[0] = dlast * lgamma_coeff[0]; + for (size_t j = 1; j < NCOEFF; j++) + { + long double dnext = dlast * y0r2 + elast; + long double enext = elast * yr2; + bterm[j] = dnext * lgamma_coeff[j]; + dlast = dnext; + elast = enext; + } + long double log_gamma_low = 0; + for (size_t j = 0; j < NCOEFF; j++) + log_gamma_low += bterm[NCOEFF - 1 - j]; + log_gamma_ratio = log_gamma_high + log_gamma_low; + + return log_sinpi_ratio + log_gamma_ratio; +} diff --git a/sysdeps/ieee754/ldbl-96/lgamma_product.c b/sysdeps/ieee754/ldbl-96/lgamma_product.c new file mode 100644 index 0000000000..e3ba72d8e4 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/lgamma_product.c @@ -0,0 +1,37 @@ +/* Compute a product of 1 + (T/X), 1 + (T/(X+1)), .... + Copyright (C) 2015-2016 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <float.h> + +/* Compute the product of 1 + (T / (X + X_EPS)), 1 + (T / (X + X_EPS + + 1)), ..., 1 + (T / (X + X_EPS + N - 1)), minus 1. X is such that + all the values X + 1, ..., X + N - 1 are exactly representable, and + X_EPS / X is small enough that factors quadratic in it can be + neglected. */ + +double +__lgamma_product (double t, double x, double x_eps, int n) +{ + long double x_full = (long double) x + (long double) x_eps; + long double ret = 0; + for (int i = 0; i < n; i++) + ret += (t / (x_full + i)) * (1 + ret); + return ret; +} diff --git a/sysdeps/ieee754/ldbl-96/lgamma_productl.c b/sysdeps/ieee754/ldbl-96/lgamma_productl.c new file mode 100644 index 0000000000..de67cbe665 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/lgamma_productl.c @@ -0,0 +1,82 @@ +/* Compute a product of 1 + (T/X), 1 + (T/(X+1)), .... + Copyright (C) 2015-2016 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <float.h> + +/* Calculate X * Y exactly and store the result in *HI + *LO. It is + given that the values are small enough that no overflow occurs and + large enough (or zero) that no underflow occurs. */ + +static void +mul_split (long double *hi, long double *lo, long double x, long double y) +{ +#ifdef __FP_FAST_FMAL + /* Fast built-in fused multiply-add. */ + *hi = x * y; + *lo = __builtin_fmal (x, y, -*hi); +#elif defined FP_FAST_FMAL + /* Fast library fused multiply-add, compiler before GCC 4.6. */ + *hi = x * y; + *lo = __fmal (x, y, -*hi); +#else + /* Apply Dekker's algorithm. */ + *hi = x * y; +# define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1) + long double x1 = x * C; + long double y1 = y * C; +# undef C + x1 = (x - x1) + x1; + y1 = (y - y1) + y1; + long double x2 = x - x1; + long double y2 = y - y1; + *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2; +#endif +} + +/* Compute the product of 1 + (T / (X + X_EPS)), 1 + (T / (X + X_EPS + + 1)), ..., 1 + (T / (X + X_EPS + N - 1)), minus 1. X is such that + all the values X + 1, ..., X + N - 1 are exactly representable, and + X_EPS / X is small enough that factors quadratic in it can be + neglected. */ + +long double +__lgamma_productl (long double t, long double x, long double x_eps, int n) +{ + long double ret = 0, ret_eps = 0; + for (int i = 0; i < n; i++) + { + long double xi = x + i; + long double quot = t / xi; + long double mhi, mlo; + mul_split (&mhi, &mlo, quot, xi); + long double quot_lo = (t - mhi - mlo) / xi - t * x_eps / (xi * xi); + /* We want (1 + RET + RET_EPS) * (1 + QUOT + QUOT_LO) - 1. */ + long double rhi, rlo; + mul_split (&rhi, &rlo, ret, quot); + long double rpq = ret + quot; + long double rpq_eps = (ret - rpq) + quot; + long double nret = rpq + rhi; + long double nret_eps = (rpq - nret) + rhi; + ret_eps += (rpq_eps + nret_eps + rlo + ret_eps * quot + + quot_lo + quot_lo * (ret + ret_eps)); + ret = nret; + } + return ret + ret_eps; +} diff --git a/sysdeps/ieee754/ldbl-96/mpn2ldbl.c b/sysdeps/ieee754/ldbl-96/mpn2ldbl.c index b219de0fca..6159d7dabf 100644 --- a/sysdeps/ieee754/ldbl-96/mpn2ldbl.c +++ b/sysdeps/ieee754/ldbl-96/mpn2ldbl.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1995-2015 Free Software Foundation, Inc. +/* Copyright (C) 1995-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or diff --git a/sysdeps/ieee754/ldbl-96/printf_fphex.c b/sysdeps/ieee754/ldbl-96/printf_fphex.c index 564f062727..65bf538590 100644 --- a/sysdeps/ieee754/ldbl-96/printf_fphex.c +++ b/sysdeps/ieee754/ldbl-96/printf_fphex.c @@ -1,5 +1,5 @@ /* Print floating point number in hexadecimal notation according to ISO C99. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or diff --git a/sysdeps/ieee754/ldbl-96/s_asinhl.c b/sysdeps/ieee754/ldbl-96/s_asinhl.c index 75e47c71ea..da49ea5988 100644 --- a/sysdeps/ieee754/ldbl-96/s_asinhl.c +++ b/sysdeps/ieee754/ldbl-96/s_asinhl.c @@ -45,11 +45,7 @@ long double __asinhl(long double x) GET_LDOUBLE_EXP(hx,x); ix = hx&0x7fff; if(__builtin_expect(ix< 0x3fde, 0)) { /* |x|<2**-34 */ - if (fabsl (x) < LDBL_MIN) - { - long double force_underflow = x * x; - math_force_eval (force_underflow); - } + math_check_force_underflow (x); if(huge+x>one) return x; /* return x inexact except 0 */ } if(__builtin_expect(ix>0x4020,0)) { /* |x| > 2**34 */ diff --git a/sysdeps/ieee754/ldbl-96/s_cbrtl.c b/sysdeps/ieee754/ldbl-96/s_cbrtl.c index bc3c2abfeb..42849ab517 100644 --- a/sysdeps/ieee754/ldbl-96/s_cbrtl.c +++ b/sysdeps/ieee754/ldbl-96/s_cbrtl.c @@ -1,5 +1,5 @@ /* Compute cubic root of double value. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Dirk Alboth <dirka@uni-paderborn.de> and Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/ldbl-96/s_erfl.c b/sysdeps/ieee754/ldbl-96/s_erfl.c index c27de812ca..d00adb1000 100644 --- a/sysdeps/ieee754/ldbl-96/s_erfl.c +++ b/sysdeps/ieee754/ldbl-96/s_erfl.c @@ -274,11 +274,7 @@ __erfl (long double x) { /* Avoid spurious underflow. */ long double ret = 0.0625 * (16.0 * x + (16.0 * efx) * x); - if (fabsl (ret) < LDBL_MIN) - { - long double force_underflow = ret * ret; - math_force_eval (force_underflow); - } + math_check_force_underflow (ret); return ret; } return x + efx * x; diff --git a/sysdeps/ieee754/ldbl-96/s_fma.c b/sysdeps/ieee754/ldbl-96/s_fma.c index 8fd238cfe4..ab45bcfce2 100644 --- a/sysdeps/ieee754/ldbl-96/s_fma.c +++ b/sysdeps/ieee754/ldbl-96/s_fma.c @@ -1,5 +1,5 @@ /* Compute x * y + z as ternary operation. - Copyright (C) 2010-2015 Free Software Foundation, Inc. + Copyright (C) 2010-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jakub@redhat.com>, 2010. @@ -41,7 +41,10 @@ __fma (double x, double y, double z) /* Ensure correct sign of exact 0 + 0. */ if (__glibc_unlikely ((x == 0 || y == 0) && z == 0)) - return x * y + z; + { + x = math_opt_barrier (x); + return x * y + z; + } fenv_t env; feholdexcept (&env); diff --git a/sysdeps/ieee754/ldbl-96/s_fmal.c b/sysdeps/ieee754/ldbl-96/s_fmal.c index eec5a02e0b..f1467fda3d 100644 --- a/sysdeps/ieee754/ldbl-96/s_fmal.c +++ b/sysdeps/ieee754/ldbl-96/s_fmal.c @@ -1,5 +1,5 @@ /* Compute x * y + z as ternary operation. - Copyright (C) 2010-2015 Free Software Foundation, Inc. + Copyright (C) 2010-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jakub@redhat.com>, 2010. @@ -68,7 +68,7 @@ __fmal (long double x, long double y, long double z) if (u.ieee.exponent + v.ieee.exponent > 0x7fff + IEEE854_LONG_DOUBLE_BIAS) return x * y; - /* If x * y is less than 1/4 of LDBL_DENORM_MIN, neither the + /* If x * y is less than 1/4 of LDBL_TRUE_MIN, neither the result nor whether there is underflow depends on its exact value, only on its sign. */ if (u.ieee.exponent + v.ieee.exponent @@ -92,8 +92,8 @@ __fmal (long double x, long double y, long double z) && w.ieee.mantissa1 == 0 && w.ieee.mantissa0 == 0x80000000))) { - volatile long double force_underflow = x * y; - (void) force_underflow; + long double force_underflow = x * y; + math_force_eval (force_underflow); } return v.d * 0x1p-65L; } @@ -119,7 +119,7 @@ __fmal (long double x, long double y, long double z) very small, adjust them up to avoid spurious underflows, rather than down. */ if (u.ieee.exponent + v.ieee.exponent - <= IEEE854_LONG_DOUBLE_BIAS + LDBL_MANT_DIG) + <= IEEE854_LONG_DOUBLE_BIAS + 2 * LDBL_MANT_DIG) { if (u.ieee.exponent > v.ieee.exponent) u.ieee.exponent += 2 * LDBL_MANT_DIG + 2; @@ -177,7 +177,10 @@ __fmal (long double x, long double y, long double z) /* Ensure correct sign of exact 0 + 0. */ if (__glibc_unlikely ((x == 0 || y == 0) && z == 0)) - return x * y + z; + { + x = math_opt_barrier (x); + return x * y + z; + } fenv_t env; feholdexcept (&env); diff --git a/sysdeps/ieee754/ldbl-96/s_isinf_nsl.c b/sysdeps/ieee754/ldbl-96/s_isinf_nsl.c deleted file mode 100644 index 9c7868b490..0000000000 --- a/sysdeps/ieee754/ldbl-96/s_isinf_nsl.c +++ /dev/null @@ -1,18 +0,0 @@ -/* - * Written by Ulrich Drepper <drepper@gmail.com>. - */ - -/* - * __isinf_nsl(x) returns != 0 if x is ±inf, else 0; - */ - -#include <math.h> -#include <math_private.h> - -int -__isinf_nsl (long double x) -{ - int32_t se,hx,lx; - GET_LDOUBLE_WORDS(se,hx,lx,x); - return !(((se & 0x7fff) ^ 0x7fff) | lx | (hx & 0x7fffffff)); -} diff --git a/sysdeps/ieee754/ldbl-96/s_issignalingl.c b/sysdeps/ieee754/ldbl-96/s_issignalingl.c index adbee9823c..73646cac0c 100644 --- a/sysdeps/ieee754/ldbl-96/s_issignalingl.c +++ b/sysdeps/ieee754/ldbl-96/s_issignalingl.c @@ -1,5 +1,5 @@ /* Test for signaling NaN. - Copyright (C) 2013-2015 Free Software Foundation, Inc. + Copyright (C) 2013-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or diff --git a/sysdeps/ieee754/ldbl-96/s_llrintl.c b/sysdeps/ieee754/ldbl-96/s_llrintl.c index 070246cb6a..592d51c607 100644 --- a/sysdeps/ieee754/ldbl-96/s_llrintl.c +++ b/sysdeps/ieee754/ldbl-96/s_llrintl.c @@ -1,6 +1,6 @@ /* Round argument to nearest integral value according to current rounding direction. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -18,6 +18,8 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <fenv.h> +#include <limits.h> #include <math.h> #include <math_private.h> @@ -35,7 +37,7 @@ __llrintl (long double x) int32_t se,j0; u_int32_t i0, i1; long long int result; - volatile long double w; + long double w; long double t; int sx; @@ -50,8 +52,21 @@ __llrintl (long double x) result = (((long long int) i0 << 32) | i1) << (j0 - 63); else { - w = two63[sx] + x; - t = w - two63[sx]; +#if defined FE_INVALID || defined FE_INEXACT + /* X < LLONG_MAX + 1 implied by J0 < 63. */ + if (x > (long double) LLONG_MAX) + { + /* In the event of overflow we must raise the "invalid" + exception, but not "inexact". */ + t = __nearbyintl (x); + feraiseexcept (t == LLONG_MAX ? FE_INEXACT : FE_INVALID); + } + else +#endif + { + w = two63[sx] + x; + t = w - two63[sx]; + } GET_LDOUBLE_WORDS (se, i0, i1, t); j0 = (se & 0x7fff) - 0x3fff; diff --git a/sysdeps/ieee754/ldbl-96/s_llroundl.c b/sysdeps/ieee754/ldbl-96/s_llroundl.c index 8381649269..483199d442 100644 --- a/sysdeps/ieee754/ldbl-96/s_llroundl.c +++ b/sysdeps/ieee754/ldbl-96/s_llroundl.c @@ -1,5 +1,5 @@ /* Round long double value to long long int. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -17,6 +17,8 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <fenv.h> +#include <limits.h> #include <math.h> #include <math_private.h> @@ -64,7 +66,14 @@ __llroundl (long double x) ++result; if (j0 > 31) - result = (result << (j0 - 31)) | (j >> (63 - j0)); + { + result = (result << (j0 - 31)) | (j >> (63 - j0)); +#ifdef FE_INVALID + if (sign == 1 && result == LLONG_MIN) + /* Rounding brought the value out of range. */ + feraiseexcept (FE_INVALID); +#endif + } } } else diff --git a/sysdeps/ieee754/ldbl-96/s_lrintl.c b/sysdeps/ieee754/ldbl-96/s_lrintl.c index 96c6a5ea7d..bd902deb47 100644 --- a/sysdeps/ieee754/ldbl-96/s_lrintl.c +++ b/sysdeps/ieee754/ldbl-96/s_lrintl.c @@ -1,6 +1,6 @@ /* Round argument to nearest integral value according to current rounding direction. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -18,6 +18,8 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <fenv.h> +#include <limits.h> #include <math.h> #include <math_private.h> @@ -35,7 +37,7 @@ __lrintl (long double x) int32_t se,j0; u_int32_t i0, i1; long int result; - volatile long double w; + long double w; long double t; int sx; @@ -46,8 +48,22 @@ __lrintl (long double x) if (j0 < 31) { - w = two63[sx] + x; - t = w - two63[sx]; +#if defined FE_INVALID || defined FE_INEXACT + /* X < LONG_MAX + 1 implied by J0 < 31. */ + if (sizeof (long int) == 4 + && x > (long double) LONG_MAX) + { + /* In the event of overflow we must raise the "invalid" + exception, but not "inexact". */ + t = __nearbyintl (x); + feraiseexcept (t == LONG_MAX ? FE_INEXACT : FE_INVALID); + } + else +#endif + { + w = two63[sx] + x; + t = w - two63[sx]; + } GET_LDOUBLE_WORDS (se, i0, i1, t); j0 = (se & 0x7fff) - 0x3fff; @@ -59,8 +75,22 @@ __lrintl (long double x) result = ((long int) i0 << (j0 - 31)) | (i1 << (j0 - 63)); else { - w = two63[sx] + x; - t = w - two63[sx]; +#if defined FE_INVALID || defined FE_INEXACT + /* X < LONG_MAX + 1 implied by J0 < 63. */ + if (sizeof (long int) == 8 + && x > (long double) LONG_MAX) + { + /* In the event of overflow we must raise the "invalid" + exception, but not "inexact". */ + t = __nearbyintl (x); + feraiseexcept (t == LONG_MAX ? FE_INEXACT : FE_INVALID); + } + else +#endif + { + w = two63[sx] + x; + t = w - two63[sx]; + } GET_LDOUBLE_WORDS (se, i0, i1, t); j0 = (se & 0x7fff) - 0x3fff; @@ -72,8 +102,21 @@ __lrintl (long double x) } else { - /* The number is too large. It is left implementation defined - what happens. */ + /* The number is too large. Unless it rounds to LONG_MIN, + FE_INVALID must be raised and the return value is + unspecified. */ +#if defined FE_INVALID || defined FE_INEXACT + if (sizeof (long int) == 4 + && x < (long double) LONG_MIN + && x > (long double) LONG_MIN - 1.0L) + { + /* If truncation produces LONG_MIN, the cast will not raise + the exception, but may raise "inexact". */ + t = __nearbyintl (x); + feraiseexcept (t == LONG_MIN ? FE_INEXACT : FE_INVALID); + return LONG_MIN; + } +#endif return (long int) x; } diff --git a/sysdeps/ieee754/ldbl-96/s_lroundl.c b/sysdeps/ieee754/ldbl-96/s_lroundl.c index e34edd307d..3b43d77e72 100644 --- a/sysdeps/ieee754/ldbl-96/s_lroundl.c +++ b/sysdeps/ieee754/ldbl-96/s_lroundl.c @@ -1,5 +1,5 @@ /* Round long double value to long int. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -17,6 +17,8 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <fenv.h> +#include <limits.h> #include <math.h> #include <math_private.h> @@ -49,6 +51,13 @@ __lroundl (long double x) } result = j >> (31 - j0); +#ifdef FE_INVALID + if (sizeof (long int) == 4 + && sign == 1 + && result == LONG_MIN) + /* Rounding brought the value out of range. */ + feraiseexcept (FE_INVALID); +#endif } } else if (j0 < (int32_t) (8 * sizeof (long int)) - 1) @@ -58,19 +67,41 @@ __lroundl (long double x) else { u_int32_t j = i1 + (0x80000000 >> (j0 - 31)); + unsigned long int ures = i0; + if (j < i1) - ++i0; + ++ures; if (j0 == 31) - result = (long int) i0; + result = ures; else - result = ((long int) i0 << (j0 - 31)) | (j >> (63 - j0)); + { + result = (ures << (j0 - 31)) | (j >> (63 - j0)); +#ifdef FE_INVALID + if (sizeof (long int) == 8 + && sign == 1 + && result == LONG_MIN) + /* Rounding brought the value out of range. */ + feraiseexcept (FE_INVALID); +#endif + } } } else { - /* The number is too large. It is left implementation defined - what happens. */ + /* The number is too large. Unless it rounds to LONG_MIN, + FE_INVALID must be raised and the return value is + unspecified. */ +#ifdef FE_INVALID + if (sizeof (long int) == 4 + && x <= (long double) LONG_MIN - 0.5L) + { + /* If truncation produces LONG_MIN, the cast will not raise + the exception, but may raise "inexact". */ + feraiseexcept (FE_INVALID); + return LONG_MIN; + } +#endif return (long int) x; } diff --git a/sysdeps/ieee754/ldbl-96/s_nexttoward.c b/sysdeps/ieee754/ldbl-96/s_nexttoward.c index f7a8b2165a..3d0382eac9 100644 --- a/sysdeps/ieee754/ldbl-96/s_nexttoward.c +++ b/sysdeps/ieee754/ldbl-96/s_nexttoward.c @@ -25,6 +25,7 @@ static char rcsid[] = "$NetBSD: $"; * Special cases: */ +#include <errno.h> #include <math.h> #include <math_private.h> #include <float.h> @@ -70,15 +71,14 @@ double __nexttoward(double x, long double y) } hy = hx&0x7ff00000; if(hy>=0x7ff00000) { - x = x+x; /* overflow */ - if (FLT_EVAL_METHOD != 0 && FLT_EVAL_METHOD != 1) - /* Force conversion to double. */ - asm ("" : "+m"(x)); - return x; + double u = x+x; /* overflow */ + math_force_eval (u); + __set_errno (ERANGE); } if(hy<0x00100000) { double u = x*x; /* underflow */ math_force_eval (u); /* raise underflow flag */ + __set_errno (ERANGE); } INSERT_WORDS(x,hx,lx); return x; diff --git a/sysdeps/ieee754/ldbl-96/s_nexttowardf.c b/sysdeps/ieee754/ldbl-96/s_nexttowardf.c index a96f9da2c2..ae7538942f 100644 --- a/sysdeps/ieee754/ldbl-96/s_nexttowardf.c +++ b/sysdeps/ieee754/ldbl-96/s_nexttowardf.c @@ -17,6 +17,7 @@ static char rcsid[] = "$NetBSD: $"; #endif +#include <errno.h> #include <math.h> #include <math_private.h> #include <float.h> @@ -58,15 +59,14 @@ float __nexttowardf(float x, long double y) } hy = hx&0x7f800000; if(hy>=0x7f800000) { - x = x+x; /* overflow */ - if (FLT_EVAL_METHOD != 0) - /* Force conversion to float. */ - asm ("" : "+m"(x)); - return x; + float u = x+x; /* overflow */ + math_force_eval (u); + __set_errno (ERANGE); } if(hy<0x00800000) { float u = x*x; /* underflow */ math_force_eval (u); /* raise underflow flag */ + __set_errno (ERANGE); } SET_FLOAT_WORD(x,hx); return x; diff --git a/sysdeps/ieee754/ldbl-96/s_remquol.c b/sysdeps/ieee754/ldbl-96/s_remquol.c index 230d0fa57c..89b2630d46 100644 --- a/sysdeps/ieee754/ldbl-96/s_remquol.c +++ b/sysdeps/ieee754/ldbl-96/s_remquol.c @@ -1,5 +1,5 @@ /* Compute remainder and a congruent to the quotient. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/ldbl-96/s_roundl.c b/sysdeps/ieee754/ldbl-96/s_roundl.c index 560f7afb99..4f35c4847b 100644 --- a/sysdeps/ieee754/ldbl-96/s_roundl.c +++ b/sysdeps/ieee754/ldbl-96/s_roundl.c @@ -1,5 +1,5 @@ /* Round long double to integer away from zero. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/ldbl-96/s_signbitl.c b/sysdeps/ieee754/ldbl-96/s_signbitl.c index bbe72a62ea..ee5d77e27a 100644 --- a/sysdeps/ieee754/ldbl-96/s_signbitl.c +++ b/sysdeps/ieee754/ldbl-96/s_signbitl.c @@ -1,5 +1,5 @@ /* Return nonzero value if number is negative. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -19,13 +19,8 @@ #include <math.h> -#include <math_private.h> - int __signbitl (long double x) { - int32_t e; - - GET_LDOUBLE_EXP (e, x); - return e & 0x8000; + return __builtin_signbitl (x); } diff --git a/sysdeps/ieee754/ldbl-96/s_sincosl.c b/sysdeps/ieee754/ldbl-96/s_sincosl.c index 37067bdaeb..ab32b73e7d 100644 --- a/sysdeps/ieee754/ldbl-96/s_sincosl.c +++ b/sysdeps/ieee754/ldbl-96/s_sincosl.c @@ -1,5 +1,5 @@ /* Compute sine and cosine of argument. - Copyright (C) 1997-2015 Free Software Foundation, Inc. + Copyright (C) 1997-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -42,7 +42,7 @@ __sincosl (long double x, long double *sinx, long double *cosx) { /* sin(Inf or NaN) is NaN */ *sinx = *cosx = x - x; - if (__isinf_nsl (x)) + if (isinf (x)) __set_errno (EDOM); } else diff --git a/sysdeps/ieee754/ldbl-96/s_tanhl.c b/sysdeps/ieee754/ldbl-96/s_tanhl.c index 7ec6247315..38edf9f75e 100644 --- a/sysdeps/ieee754/ldbl-96/s_tanhl.c +++ b/sysdeps/ieee754/ldbl-96/s_tanhl.c @@ -42,6 +42,7 @@ static char rcsid[] = "$NetBSD: $"; * only tanhl(0)=0 is exact for finite argument. */ +#include <float.h> #include <math.h> #include <math_private.h> @@ -69,7 +70,10 @@ long double __tanhl(long double x) if ((ix|j0|j1) == 0) return x; /* x == +- 0 */ if (ix<0x3fc8) /* |x|<2**-55 */ + { + math_check_force_underflow (x); return x*(one+tiny); /* tanh(small) = small */ + } if (ix>=0x3fff) { /* |x|>=1 */ t = __expm1l(two*fabsl(x)); z = one - two/(t+two); diff --git a/sysdeps/ieee754/ldbl-96/strtod_nan_ldouble.h b/sysdeps/ieee754/ldbl-96/strtod_nan_ldouble.h new file mode 100644 index 0000000000..2694b5ee34 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/strtod_nan_ldouble.h @@ -0,0 +1,30 @@ +/* Convert string for NaN payload to corresponding NaN. For ldbl-96. + Copyright (C) 1997-2016 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#define FLOAT long double +#define SET_MANTISSA(flt, mant) \ + do \ + { \ + union ieee854_long_double u; \ + u.d = (flt); \ + u.ieee_nan.mantissa0 = (mant) >> 32; \ + u.ieee_nan.mantissa1 = (mant); \ + if ((u.ieee.mantissa0 | u.ieee.mantissa1) != 0) \ + (flt) = u.d; \ + } \ + while (0) diff --git a/sysdeps/ieee754/ldbl-96/strtold_l.c b/sysdeps/ieee754/ldbl-96/strtold_l.c index c082e7472e..b51560e18a 100644 --- a/sysdeps/ieee754/ldbl-96/strtold_l.c +++ b/sysdeps/ieee754/ldbl-96/strtold_l.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1999-2015 Free Software Foundation, Inc. +/* Copyright (C) 1999-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -25,19 +25,13 @@ #ifdef USE_WIDE_CHAR # define STRTOF wcstold_l # define __STRTOF __wcstold_l +# define STRTOF_NAN __wcstold_nan #else # define STRTOF strtold_l # define __STRTOF __strtold_l +# define STRTOF_NAN __strtold_nan #endif #define MPN2FLOAT __mpn_construct_long_double #define FLOAT_HUGE_VAL HUGE_VALL -#define SET_MANTISSA(flt, mant) \ - do { union ieee854_long_double u; \ - u.d = (flt); \ - u.ieee_nan.mantissa0 = (mant) >> 32; \ - u.ieee_nan.mantissa1 = (mant); \ - if ((u.ieee.mantissa0 | u.ieee.mantissa1) != 0) \ - (flt) = u.d; \ - } while (0) #include <stdlib/strtod_l.c> diff --git a/sysdeps/ieee754/ldbl-96/t_sincosl.c b/sysdeps/ieee754/ldbl-96/t_sincosl.c index 4efc4d9a02..f00a182061 100644 --- a/sysdeps/ieee754/ldbl-96/t_sincosl.c +++ b/sysdeps/ieee754/ldbl-96/t_sincosl.c @@ -1,5 +1,5 @@ /* Extended-precision floating point sine and cosine tables. - Copyright (C) 1999-2015 Free Software Foundation, Inc. + Copyright (C) 1999-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Based on quad-precision tables by Jakub Jelinek <jj@ultra.linux.cz> diff --git a/sysdeps/ieee754/ldbl-96/w_expl.c b/sysdeps/ieee754/ldbl-96/w_expl.c index 3cc165b39d..5c7a1812bc 100644 --- a/sysdeps/ieee754/ldbl-96/w_expl.c +++ b/sysdeps/ieee754/ldbl-96/w_expl.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2011-2015 Free Software Foundation, Inc. +/* Copyright (C) 2011-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@gmail.com>, 2011. diff --git a/sysdeps/ieee754/ldbl-96/x2y2m1.c b/sysdeps/ieee754/ldbl-96/x2y2m1.c index a6cc82cb29..4119f42acc 100644 --- a/sysdeps/ieee754/ldbl-96/x2y2m1.c +++ b/sysdeps/ieee754/ldbl-96/x2y2m1.c @@ -1,5 +1,5 @@ /* Compute x^2 + y^2 - 1, without large cancellation error. - Copyright (C) 2012-2015 Free Software Foundation, Inc. + Copyright (C) 2012-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -27,8 +27,8 @@ #else /* Return X^2 + Y^2 - 1, computed without large cancellation error. - It is given that 1 > X >= Y >= epsilon / 2, and that either X >= - 0.75 or Y >= 0.5. */ + It is given that 1 > X >= Y >= epsilon / 2, and that X^2 + Y^2 >= + 0.5. */ double __x2y2m1 (double x, double y) diff --git a/sysdeps/ieee754/ldbl-96/x2y2m1l.c b/sysdeps/ieee754/ldbl-96/x2y2m1l.c index 11757c6af8..733742da04 100644 --- a/sysdeps/ieee754/ldbl-96/x2y2m1l.c +++ b/sysdeps/ieee754/ldbl-96/x2y2m1l.c @@ -1,5 +1,5 @@ /* Compute x^2 + y^2 - 1, without large cancellation error. - Copyright (C) 2012-2015 Free Software Foundation, Inc. + Copyright (C) 2012-2016 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -80,32 +80,26 @@ compare (const void *p, const void *q) } /* Return X^2 + Y^2 - 1, computed without large cancellation error. - It is given that 1 > X >= Y >= epsilon / 2, and that either X >= - 0.75 or Y >= 0.5. */ + It is given that 1 > X >= Y >= epsilon / 2, and that X^2 + Y^2 >= + 0.5. */ long double __x2y2m1l (long double x, long double y) { - long double vals[4]; + long double vals[5]; SET_RESTORE_ROUNDL (FE_TONEAREST); mul_split (&vals[1], &vals[0], x, x); mul_split (&vals[3], &vals[2], y, y); - if (x >= 0.75L) - vals[1] -= 1.0L; - else - { - vals[1] -= 0.5L; - vals[3] -= 0.5L; - } - qsort (vals, 4, sizeof (long double), compare); + vals[4] = -1.0L; + qsort (vals, 5, sizeof (long double), compare); /* Add up the values so that each element of VALS has absolute value at most equal to the last set bit of the next nonzero element. */ - for (size_t i = 0; i <= 2; i++) + for (size_t i = 0; i <= 3; i++) { add_split (&vals[i + 1], &vals[i], vals[i + 1], vals[i]); - qsort (vals + i + 1, 3 - i, sizeof (long double), compare); + qsort (vals + i + 1, 4 - i, sizeof (long double), compare); } /* Now any error from this addition will be small. */ - return vals[3] + vals[2] + vals[1] + vals[0]; + return vals[4] + vals[3] + vals[2] + vals[1] + vals[0]; } |