Commit 4239f144 by Joseph Myers Committed by Joseph Myers

Update most of libquadmath/math/ from glibc, automate update (PR libquadmath/68686).

libquadmath sources are mostly based on glibc sources at present, but
derived from them by a manual editing / substitution process and with
subsequent manual merges.  The manual effort involved in merges means
they are sometimes incomplete and long-delayed.

Since libquadmath was first created, glibc's support for this format
has undergone significant changes so that it can also be used in glibc
to provide *f128 functions for the _Float128 type from TS 18661-3.
This makes it significantly easier to use it for libquadmath in a more
automated fashion, since glibc has a float128_private.h header that
redefines many identifiers as macros as needed for building *f128
functions.

Simply using float128_private.h directly in libquadmath, with
unmodified glibc sources except for changing function names in that
one header to be *q instead of *f128, would be tricky, given its
dependence on lots of other glibc-internal headers (whereas
libquadmath supports non-glibc systems), and also given how some libm
functions in glibc are built from type-generic templates using a
further set of macros rather than from separate function
implementations for each type.

So instead this patch adds a script update-quadmath.py to convert
glibc sources into libquadmath ones, and the script reads
float128_private.h to identify many of the substitutions it should
make.  quadmath-imp.h is updated with various new internal
definitions, taken from glibc as needed; this is the main place
expected to need updating manually when subsequent merges from glibc
are done using the script.  No attempt is made to make the script
output match the details of existing formatting, although the
differences are of a size that makes a rough comparison (ignoring
whitespace) possible.

Two new public interfaces are added to libquadmath, exp2q and
issignalingq, at a new QUADMATH_1.2 symbol version, since those
interfaces are used internally by some of the glibc sources being
merged into libquadmath; although there is a new symbol version, no
change however is made to the libtool version in the libtool-version
file.  Although there are various other interfaces now in glibc libm
but not in libquadmath, this patch does nothing to add such interfaces
(although adding many of them would in fact be easy to do, given the
script).

One internal file (not providing any public interfaces),
math/isinf_nsq.c, is removed, as no longer used by anything in
libquadmath after the merge.

Conditionals in individual source files on <fenv.h> availability or
features are moved into quadmath-imp.h (providing trivial macro
versions of the functions if real implementations aren't available),
to simplify the substitutions in individual source files.  Note
however that I haven't tested for any configurations lacking <fenv.h>,
so further changes could well be needed there.

Two files in libquadmath/math/ are based on glibc sources but not
updated in this patch: fmaq.c and rem_pio2q.c.  Both could be updated
after further changes to the script (and quadmath-imp.h as needed); in
the case of rem_pio2q.c, based on two separate glibc source files,
those separate files would naturally be split out into separate
libquadmath source files in the process (as done in this patch with
expq_table.h and tanq_kernel.c, where previously two glibc source
files had been merged into one libquadmath source file).  complex.c,
nanq.c and sqrtq.c are not based on glibc sources (though four of the
(trivial) functions in complex.c could readily be replaced by instead
using the four corresponding files from glibc, if desired).

libquadmath also has printf/ and strtod/ sources based on glibc, also
mostly not updated for a long time.  Again the script could no doubt
be made to generate those automatically, although that would be a
larger change (effectively some completely separate logic in the
script, not sharing much if anything with the existing code).

Bootstrapped with no regressions on x86_64-pc-linux-gnu.

	PR libquadmath/68686
	* Makefile.am: (libquadmath_la_SOURCES): Remove math/isinf_nsq.c.
	Add math/exp2q.c math/issignalingq.c math/lgammaq_neg.c
	math/lgammaq_product.c math/tanq_kernel.c math/tgammaq_product.c
	math/casinhq_kernel.c.
	* Makefile.in: Regenerate.
	* libquadmath.texi (exp2q, issignalingq): Document.
	* quadmath-imp.h: Include <errno.h>, <limits.h>, <stdbool.h> and
	<fenv.h>.
	(HIGH_ORDER_BIT_IS_SET_FOR_SNAN, FIX_FLT128_LONG_CONVERT_OVERFLOW)
	(FIX_FLT128_LLONG_CONVERT_OVERFLOW, __quadmath_kernel_tanq)
	(__quadmath_gamma_productq, __quadmath_gammaq_r)
	(__quadmath_lgamma_negq, __quadmath_lgamma_productq)
	(__quadmath_lgammaq_r, __quadmath_kernel_casinhq, mul_splitq)
	(math_check_force_underflow_complex, __glibc_likely)
	(__glibc_unlikely, struct rm_ctx, SET_RESTORE_ROUNDF128)
	(libc_feholdsetround_ctx, libc_feresetround_ctx): New.
	(feraiseexcept, fenv_t, feholdexcept, fesetround, feupdateenv)
	(fesetenv, fetestexcept, feclearexcept): Define if not supported
	through <fenv.h>.
	(__quadmath_isinf_nsq): Remove.
	* quadmath.h (exp2q, issignalingq): New.
	* quadmath.map (QUADMATH_1.2): New.
	* quadmath_weak.h (exp2q, issignalingq): New.
	* update-quadmath.py: New file.
	* math/isinf_nsq.c: Remove file.
	* math/casinhq_kernel.c, math/exp2q.c, math/expq_table.h,
	math/issignalingq.c, math/lgammaq_neg.c, math/lgammaq_product.c,
	math/tanq_kernel.c, math/tgammaq_product.c: New files.  Generated
	from glibc sources with update-quadmath.py.
	* math/acoshq.c, math/acosq.c, math/asinhq.c, math/asinq.c,
	math/atan2q.c, math/atanhq.c, math/atanq.c, math/cacoshq.c,
	math/cacosq.c, math/casinhq.c, math/casinq.c, math/catanhq.c,
	math/catanq.c, math/cbrtq.c, math/ccoshq.c, math/ceilq.c,
	math/cexpq.c, math/cimagq.c, math/clog10q.c, math/clogq.c,
	math/conjq.c, math/copysignq.c, math/coshq.c, math/cosq.c,
	math/cosq_kernel.c, math/cprojq.c, math/crealq.c, math/csinhq.c,
	math/csinq.c, math/csqrtq.c, math/ctanhq.c, math/ctanq.c,
	math/erfq.c, math/expm1q.c, math/expq.c, math/fabsq.c,
	math/fdimq.c, math/finiteq.c, math/floorq.c, math/fmaxq.c,
	math/fminq.c, math/fmodq.c, math/frexpq.c, math/hypotq.c,
	math/ilogbq.c, math/isinfq.c, math/isnanq.c, math/j0q.c,
	math/j1q.c, math/jnq.c, math/ldexpq.c, math/lgammaq.c,
	math/llrintq.c, math/llroundq.c, math/log10q.c, math/log1pq.c,
	math/log2q.c, math/logbq.c, math/logq.c, math/lrintq.c,
	math/lroundq.c, math/modfq.c, math/nearbyintq.c,
	math/nextafterq.c, math/powq.c, math/remainderq.c, math/remquoq.c,
	math/rintq.c, math/roundq.c, math/scalblnq.c, math/scalbnq.c,
	math/signbitq.c, math/sincos_table.c, math/sincosq.c,
	math/sincosq_kernel.c, math/sinhq.c, math/sinq.c,
	math/sinq_kernel.c, math/tanhq.c, math/tanq.c, math/tgammaq.c,
	math/truncq.c, math/x2y2m1q.c: Regenerate from glibc sources with
	update-quadmath.py.

From-SVN: r265822
parent 0f31a750
2018-11-05 Joseph Myers <joseph@codesourcery.com>
PR libquadmath/68686
* Makefile.am: (libquadmath_la_SOURCES): Remove math/isinf_nsq.c.
Add math/exp2q.c math/issignalingq.c math/lgammaq_neg.c
math/lgammaq_product.c math/tanq_kernel.c math/tgammaq_product.c
math/casinhq_kernel.c.
* Makefile.in: Regenerate.
* libquadmath.texi (exp2q, issignalingq): Document.
* quadmath-imp.h: Include <errno.h>, <limits.h>, <stdbool.h> and
<fenv.h>.
(HIGH_ORDER_BIT_IS_SET_FOR_SNAN, FIX_FLT128_LONG_CONVERT_OVERFLOW)
(FIX_FLT128_LLONG_CONVERT_OVERFLOW, __quadmath_kernel_tanq)
(__quadmath_gamma_productq, __quadmath_gammaq_r)
(__quadmath_lgamma_negq, __quadmath_lgamma_productq)
(__quadmath_lgammaq_r, __quadmath_kernel_casinhq, mul_splitq)
(math_check_force_underflow_complex, __glibc_likely)
(__glibc_unlikely, struct rm_ctx, SET_RESTORE_ROUNDF128)
(libc_feholdsetround_ctx, libc_feresetround_ctx): New.
(feraiseexcept, fenv_t, feholdexcept, fesetround, feupdateenv)
(fesetenv, fetestexcept, feclearexcept): Define if not supported
through <fenv.h>.
(__quadmath_isinf_nsq): Remove.
* quadmath.h (exp2q, issignalingq): New.
* quadmath.map (QUADMATH_1.2): New.
* quadmath_weak.h (exp2q, issignalingq): New.
* update-quadmath.py: New file.
* math/isinf_nsq.c: Remove file.
* math/casinhq_kernel.c, math/exp2q.c, math/expq_table.h,
math/issignalingq.c, math/lgammaq_neg.c, math/lgammaq_product.c,
math/tanq_kernel.c, math/tgammaq_product.c: New files. Generated
from glibc sources with update-quadmath.py.
* math/acoshq.c, math/acosq.c, math/asinhq.c, math/asinq.c,
math/atan2q.c, math/atanhq.c, math/atanq.c, math/cacoshq.c,
math/cacosq.c, math/casinhq.c, math/casinq.c, math/catanhq.c,
math/catanq.c, math/cbrtq.c, math/ccoshq.c, math/ceilq.c,
math/cexpq.c, math/cimagq.c, math/clog10q.c, math/clogq.c,
math/conjq.c, math/copysignq.c, math/coshq.c, math/cosq.c,
math/cosq_kernel.c, math/cprojq.c, math/crealq.c, math/csinhq.c,
math/csinq.c, math/csqrtq.c, math/ctanhq.c, math/ctanq.c,
math/erfq.c, math/expm1q.c, math/expq.c, math/fabsq.c,
math/fdimq.c, math/finiteq.c, math/floorq.c, math/fmaxq.c,
math/fminq.c, math/fmodq.c, math/frexpq.c, math/hypotq.c,
math/ilogbq.c, math/isinfq.c, math/isnanq.c, math/j0q.c,
math/j1q.c, math/jnq.c, math/ldexpq.c, math/lgammaq.c,
math/llrintq.c, math/llroundq.c, math/log10q.c, math/log1pq.c,
math/log2q.c, math/logbq.c, math/logq.c, math/lrintq.c,
math/lroundq.c, math/modfq.c, math/nearbyintq.c,
math/nextafterq.c, math/powq.c, math/remainderq.c, math/remquoq.c,
math/rintq.c, math/roundq.c, math/scalblnq.c, math/scalbnq.c,
math/signbitq.c, math/sincos_table.c, math/sincosq.c,
math/sincosq_kernel.c, math/sinhq.c, math/sinq.c,
math/sinq_kernel.c, math/tanhq.c, math/tanq.c, math/tgammaq.c,
math/truncq.c, math/x2y2m1q.c: Regenerate from glibc sources with
update-quadmath.py.
2018-10-31 Joseph Myers <joseph@codesourcery.com> 2018-10-31 Joseph Myers <joseph@codesourcery.com>
PR bootstrap/82856 PR bootstrap/82856
......
...@@ -44,7 +44,7 @@ nodist_libsubinclude_HEADERS = quadmath.h quadmath_weak.h ...@@ -44,7 +44,7 @@ nodist_libsubinclude_HEADERS = quadmath.h quadmath_weak.h
libsubincludedir = $(libdir)/gcc/$(target_alias)/$(gcc_version)/include libsubincludedir = $(libdir)/gcc/$(target_alias)/$(gcc_version)/include
libquadmath_la_SOURCES = \ libquadmath_la_SOURCES = \
math/x2y2m1q.c math/isinf_nsq.c math/acoshq.c math/fmodq.c \ math/x2y2m1q.c math/acoshq.c math/fmodq.c \
math/acosq.c math/frexpq.c \ math/acosq.c math/frexpq.c \
math/rem_pio2q.c math/asinhq.c math/hypotq.c math/remainderq.c \ math/rem_pio2q.c math/asinhq.c math/hypotq.c math/remainderq.c \
math/asinq.c math/rintq.c math/atan2q.c math/isinfq.c \ math/asinq.c math/rintq.c math/atan2q.c math/isinfq.c \
...@@ -58,6 +58,8 @@ libquadmath_la_SOURCES = \ ...@@ -58,6 +58,8 @@ libquadmath_la_SOURCES = \
math/tanhq.c math/expq.c math/modfq.c math/tanq.c math/fabsq.c \ math/tanhq.c math/expq.c math/modfq.c math/tanq.c math/fabsq.c \
math/nanq.c math/tgammaq.c math/finiteq.c math/nextafterq.c \ math/nanq.c math/tgammaq.c math/finiteq.c math/nextafterq.c \
math/truncq.c math/floorq.c math/powq.c math/fmaq.c math/logbq.c \ math/truncq.c math/floorq.c math/powq.c math/fmaq.c math/logbq.c \
math/exp2q.c math/issignalingq.c math/lgammaq_neg.c math/lgammaq_product.c \
math/tanq_kernel.c math/tgammaq_product.c math/casinhq_kernel.c \
math/cacoshq.c math/cacosq.c math/casinhq.c math/casinq.c \ math/cacoshq.c math/cacosq.c math/casinhq.c math/casinq.c \
math/catanhq.c math/catanq.c math/cimagq.c math/conjq.c math/cprojq.c \ math/catanhq.c math/catanq.c math/cimagq.c math/conjq.c math/cprojq.c \
math/crealq.c math/fdimq.c math/fmaxq.c math/fminq.c math/ilogbq.c \ math/crealq.c math/fdimq.c math/fmaxq.c math/fminq.c math/ilogbq.c \
......
...@@ -157,6 +157,7 @@ The following mathematical functions are available: ...@@ -157,6 +157,7 @@ The following mathematical functions are available:
@item @code{cosq}: cosine function @item @code{cosq}: cosine function
@item @code{erfq}: error function @item @code{erfq}: error function
@item @code{erfcq}: complementary error function @item @code{erfcq}: complementary error function
@item @code{exp2q}: base 2 exponential function
@item @code{expq}: exponential function @item @code{expq}: exponential function
@item @code{expm1q}: exponential minus 1 function @item @code{expm1q}: exponential minus 1 function
@need 800 @need 800
...@@ -173,6 +174,7 @@ The following mathematical functions are available: ...@@ -173,6 +174,7 @@ The following mathematical functions are available:
@item @code{ilogbq}: get exponent of the value @item @code{ilogbq}: get exponent of the value
@item @code{isinfq}: check for infinity @item @code{isinfq}: check for infinity
@item @code{isnanq}: check for not a number @item @code{isnanq}: check for not a number
@item @code{issignalingq}: check for signaling not a number
@item @code{j0q}: Bessel function of the first kind, first order @item @code{j0q}: Bessel function of the first kind, first order
@item @code{j1q}: Bessel function of the first kind, second order @item @code{j1q}: Bessel function of the first kind, second order
@item @code{jnq}: Bessel function of the first kind, @var{n}-th order @item @code{jnq}: Bessel function of the first kind, @var{n}-th order
......
/* acoshq.c -- __float128 version of e_acosh.c. /* e_acoshl.c -- long double version of e_acosh.c.
* Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz. * Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz.
*/ */
...@@ -16,11 +16,11 @@ ...@@ -16,11 +16,11 @@
/* acoshq(x) /* acoshq(x)
* Method : * Method :
* Based on * Based on
* acoshl(x) = logl [ x + sqrtl(x*x-1) ] * acoshl(x) = logq [ x + sqrtq(x*x-1) ]
* we have * we have
* acoshl(x) := logl(x)+ln2, if x is large; else * acoshl(x) := logq(x)+ln2, if x is large; else
* acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else * acoshl(x) := logq(2x-1/(sqrtq(x*x-1)+x)) if x>2; else
* acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1. * acoshl(x) := log1pq(t+sqrtq(2.0*t+t*t)); where t=x-1.
* *
* Special cases: * Special cases:
* acoshl(x) is NaN with signal if x<1. * acoshl(x) is NaN with signal if x<1.
...@@ -30,11 +30,11 @@ ...@@ -30,11 +30,11 @@
#include "quadmath-imp.h" #include "quadmath-imp.h"
static const __float128 static const __float128
one = 1.0Q, one = 1.0,
ln2 = 0.6931471805599453094172321214581766Q; ln2 = 0.6931471805599453094172321214581766Q;
__float128 __float128
acoshq (__float128 x) acoshq(__float128 x)
{ {
__float128 t; __float128 t;
uint64_t lx; uint64_t lx;
...@@ -44,16 +44,16 @@ acoshq (__float128 x) ...@@ -44,16 +44,16 @@ acoshq (__float128 x)
return (x-x)/(x-x); return (x-x)/(x-x);
} else if(hx >=0x4035000000000000LL) { /* x > 2**54 */ } else if(hx >=0x4035000000000000LL) { /* x > 2**54 */
if(hx >=0x7fff000000000000LL) { /* x is inf of NaN */ if(hx >=0x7fff000000000000LL) { /* x is inf of NaN */
return x+x; return x+x;
} else } else
return logq(x)+ln2; /* acoshl(huge)=logl(2x) */ return logq(x)+ln2; /* acoshl(huge)=logq(2x) */
} else if(((hx-0x3fff000000000000LL)|lx)==0) { } else if(((hx-0x3fff000000000000LL)|lx)==0) {
return 0.0Q; /* acosh(1) = 0 */ return 0; /* acosh(1) = 0 */
} else if (hx > 0x4000000000000000LL) { /* 2**28 > x > 2 */ } else if (hx > 0x4000000000000000LL) { /* 2**28 > x > 2 */
t=x*x; t=x*x;
return logq(2.0Q*x-one/(x+sqrtq(t-one))); return logq(2*x-one/(x+sqrtq(t-one)));
} else { /* 1<x<2 */ } else { /* 1<x<2 */
t = x-one; t = x-one;
return log1pq(t+sqrtq(2.0Q*t+t*t)); return log1pq(t+sqrtq(2*t+t*t));
} }
} }
...@@ -10,11 +10,11 @@ ...@@ -10,11 +10,11 @@
*/ */
/* /*
__float128 expansions are Long double expansions are
Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov> Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
and are incorporated herein by permission of the author. The author and are incorporated herein by permission of the author. The author
reserves the right to distribute this material elsewhere under different reserves the right to distribute this material elsewhere under different
copying permissions. These modifications are distributed here under copying permissions. These modifications are distributed here under
the following terms: the following terms:
This library is free software; you can redistribute it and/or This library is free software; you can redistribute it and/or
...@@ -28,8 +28,8 @@ ...@@ -28,8 +28,8 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software License along with this library; if not, see
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ <http://www.gnu.org/licenses/>. */
/* acosq(x) /* acosq(x)
* Method : * Method :
...@@ -57,7 +57,7 @@ ...@@ -57,7 +57,7 @@
#include "quadmath-imp.h" #include "quadmath-imp.h"
static const __float128 static const __float128
one = 1.0Q, one = 1,
pio2_hi = 1.5707963267948966192313216916397514420986Q, pio2_hi = 1.5707963267948966192313216916397514420986Q,
pio2_lo = 4.3359050650618905123985220130216759843812E-35Q, pio2_lo = 4.3359050650618905123985220130216759843812E-35Q,
...@@ -272,11 +272,11 @@ acosq (__float128 x) ...@@ -272,11 +272,11 @@ acosq (__float128 x)
s = sqrtq (z); s = sqrtq (z);
/* Compute an extended precision square root from /* Compute an extended precision square root from
the Newton iteration s -> 0.5 * (s + z / s). the Newton iteration s -> 0.5 * (s + z / s).
The change w from s to the improved value is The change w from s to the improved value is
w = 0.5 * (s + z / s) - s = (s^2 + z)/2s - s = (z - s^2)/2s. w = 0.5 * (s + z / s) - s = (s^2 + z)/2s - s = (z - s^2)/2s.
Express s = f1 + f2 where f1 * f1 is exactly representable. Express s = f1 + f2 where f1 * f1 is exactly representable.
w = (z - s^2)/2s = (z - f1^2 - 2 f1 f2 - f2^2)/2s . w = (z - s^2)/2s = (z - f1^2 - 2 f1 f2 - f2^2)/2s .
s + w has extended precision. */ s + w has extended precision. */
u.value = s; u.value = s;
u.words32.w2 = 0; u.words32.w2 = 0;
u.words32.w3 = 0; u.words32.w3 = 0;
......
/* asinhq.c -- __float128 version of s_asinh.c. /* s_asinhl.c -- long double version of s_asinh.c.
* Conversion to long double by Ulrich Drepper, * Conversion to long double by Ulrich Drepper,
* Cygnus Support, drepper@cygnus.com. * Cygnus Support, drepper@cygnus.com.
*/ */
...@@ -14,21 +14,25 @@ ...@@ -14,21 +14,25 @@
* ==================================================== * ====================================================
*/ */
/* asinhl(x) #if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/* asinhq(x)
* Method : * Method :
* Based on * Based on
* asinhl(x) = signl(x) * logl [ |x| + sqrtl(x*x+1) ] * asinhq(x) = signl(x) * logq [ |x| + sqrtq(x*x+1) ]
* we have * we have
* asinhl(x) := x if 1+x*x=1, * asinhq(x) := x if 1+x*x=1,
* := signl(x)*(logl(x)+ln2)) for large |x|, else * := signl(x)*(logq(x)+ln2)) for large |x|, else
* := signl(x)*logl(2|x|+1/(|x|+sqrtl(x*x+1))) if|x|>2, else * := signl(x)*logq(2|x|+1/(|x|+sqrtq(x*x+1))) if|x|>2, else
* := signl(x)*log1pl(|x| + x^2/(1 + sqrtl(1+x^2))) * := signl(x)*log1pq(|x| + x^2/(1 + sqrtq(1+x^2)))
*/ */
#include "quadmath-imp.h" #include "quadmath-imp.h"
static const __float128 static const __float128
one = 1.0Q, one = 1,
ln2 = 6.931471805599453094172321214581765681e-1Q, ln2 = 6.931471805599453094172321214581765681e-1Q,
huge = 1.0e+4900Q; huge = 1.0e+4900Q;
......
...@@ -10,11 +10,11 @@ ...@@ -10,11 +10,11 @@
*/ */
/* /*
__float128 expansions are Long double expansions are
Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov> Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
and are incorporated herein by permission of the author. The author and are incorporated herein by permission of the author. The author
reserves the right to distribute this material elsewhere under different reserves the right to distribute this material elsewhere under different
copying permissions. These modifications are distributed here under the copying permissions. These modifications are distributed here under the
following terms: following terms:
This library is free software; you can redistribute it and/or This library is free software; you can redistribute it and/or
...@@ -28,10 +28,10 @@ ...@@ -28,10 +28,10 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software License along with this library; if not, see
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ <http://www.gnu.org/licenses/>. */
/* asinq(x) /* __ieee754_asin(x)
* Method : * Method :
* Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ... * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
* we approximate asin(x) on [0,0.5] by * we approximate asin(x) on [0,0.5] by
...@@ -58,11 +58,10 @@ ...@@ -58,11 +58,10 @@
* *
*/ */
#include "quadmath-imp.h" #include "quadmath-imp.h"
static const __float128 static const __float128
one = 1.0Q, one = 1,
huge = 1.0e+4932Q, huge = 1.0e+4932Q,
pio2_hi = 1.5707963267948966192313216916397514420986Q, pio2_hi = 1.5707963267948966192313216916397514420986Q,
pio2_lo = 4.3359050650618905123985220130216759843812E-35Q, pio2_lo = 4.3359050650618905123985220130216759843812E-35Q,
...@@ -129,8 +128,7 @@ static const __float128 ...@@ -129,8 +128,7 @@ static const __float128
__float128 __float128
asinq (__float128 x) asinq (__float128 x)
{ {
__float128 t = 0; __float128 t, w, p, q, c, r, s;
__float128 w, p, q, c, r, s;
int32_t ix, sign, flag; int32_t ix, sign, flag;
ieee854_float128 u; ieee854_float128 u;
......
/* atan2q.c -- __float128 version of e_atan2.c. /* e_atan2l.c -- long double version of e_atan2.c.
* Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz. * Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz.
*/ */
...@@ -15,9 +15,9 @@ ...@@ -15,9 +15,9 @@
/* atan2q(y,x) /* atan2q(y,x)
* Method : * Method :
* 1. Reduce y to positive by atan2q(y,x)=-atan2q(-y,x). * 1. Reduce y to positive by atan2l(y,x)=-atan2l(-y,x).
* 2. Reduce x to positive by (if x and y are unexceptional): * 2. Reduce x to positive by (if x and y are unexceptional):
* ARG (x+iy) = arctan(y/x) ... if x > 0, * ARG (x+iy) = arctan(y/x) ... if x > 0,
* ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
* *
* Special cases: * Special cases:
...@@ -51,7 +51,7 @@ pi = 3.14159265358979323846264338327950280e+00Q, /* 4000921fb54442d18469898 ...@@ -51,7 +51,7 @@ pi = 3.14159265358979323846264338327950280e+00Q, /* 4000921fb54442d18469898
pi_lo = 8.67181013012378102479704402604335225e-35Q; /* 3f8dcd129024e088a67cc74020bbea64 */ pi_lo = 8.67181013012378102479704402604335225e-35Q; /* 3f8dcd129024e088a67cc74020bbea64 */
__float128 __float128
atan2q (__float128 y, __float128 x) atan2q(__float128 y, __float128 x)
{ {
__float128 z; __float128 z;
int64_t k,m,hx,hy,ix,iy; int64_t k,m,hx,hy,ix,iy;
...@@ -64,14 +64,14 @@ atan2q (__float128 y, __float128 x) ...@@ -64,14 +64,14 @@ atan2q (__float128 y, __float128 x)
if(((ix|((lx|-lx)>>63))>0x7fff000000000000LL)|| if(((ix|((lx|-lx)>>63))>0x7fff000000000000LL)||
((iy|((ly|-ly)>>63))>0x7fff000000000000LL)) /* x or y is NaN */ ((iy|((ly|-ly)>>63))>0x7fff000000000000LL)) /* x or y is NaN */
return x+y; return x+y;
if(((hx-0x3fff000000000000LL)|lx)==0) return atanq(y); /* x=1.0Q */ if(((hx-0x3fff000000000000LL)|lx)==0) return atanq(y); /* x=1.0L */
m = ((hy>>63)&1)|((hx>>62)&2); /* 2*sign(x)+sign(y) */ m = ((hy>>63)&1)|((hx>>62)&2); /* 2*sign(x)+sign(y) */
/* when y = 0 */ /* when y = 0 */
if((iy|ly)==0) { if((iy|ly)==0) {
switch(m) { switch(m) {
case 0: case 0:
case 1: return y; /* atan(+-0,+anything)=+-0 */ case 1: return y; /* atan(+-0,+anything)=+-0 */
case 2: return pi+tiny;/* atan(+0,-anything) = pi */ case 2: return pi+tiny;/* atan(+0,-anything) = pi */
case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */ case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
} }
...@@ -85,8 +85,8 @@ atan2q (__float128 y, __float128 x) ...@@ -85,8 +85,8 @@ atan2q (__float128 y, __float128 x)
switch(m) { switch(m) {
case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */ case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */
case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */ case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
case 2: return 3.0Q*pi_o_4+tiny;/*atan(+INF,-INF)*/ case 2: return 3*pi_o_4+tiny;/*atan(+INF,-INF)*/
case 3: return -3.0Q*pi_o_4-tiny;/*atan(-INF,-INF)*/ case 3: return -3*pi_o_4-tiny;/*atan(-INF,-INF)*/
} }
} else { } else {
switch(m) { switch(m) {
...@@ -102,19 +102,19 @@ atan2q (__float128 y, __float128 x) ...@@ -102,19 +102,19 @@ atan2q (__float128 y, __float128 x)
/* compute y/x */ /* compute y/x */
k = (iy-ix)>>48; k = (iy-ix)>>48;
if(k > 120) z=pi_o_2+0.5Q*pi_lo; /* |y/x| > 2**120 */ if(k > 120) z=pi_o_2+0.5Q*pi_lo; /* |y/x| > 2**120 */
else if(hx<0&&k<-120) z=0.0Q; /* |y|/x < -2**120 */ else if(hx<0&&k<-120) z=0; /* |y|/x < -2**120 */
else z=atanq(fabsq(y/x)); /* safe to do y/x */ else z=atanq(fabsq(y/x)); /* safe to do y/x */
switch (m) { switch (m) {
case 0: return z ; /* atan(+,+) */ case 0: return z ; /* atan(+,+) */
case 1: { case 1: {
uint64_t zh; uint64_t zh;
GET_FLT128_MSW64(zh,z); GET_FLT128_MSW64(zh,z);
SET_FLT128_MSW64(z,zh ^ 0x8000000000000000ULL); SET_FLT128_MSW64(z,zh ^ 0x8000000000000000ULL);
} }
return z ; /* atan(-,+) */ return z ; /* atan(-,+) */
case 2: return pi-(z-pi_lo);/* atan(+,-) */ case 2: return pi-(z-pi_lo);/* atan(+,-) */
default: /* case 3 */ default: /* case 3 */
return (z-pi_lo)-pi;/* atan(-,-) */ return (z-pi_lo)-pi;/* atan(-,-) */
} }
} }
/* s_atanhl.c -- __float128 version of s_atan.c. /* s_atanhl.c -- long double version of s_atan.c.
* Conversion to __float128 by Ulrich Drepper, * Conversion to long double by Ulrich Drepper,
* Cygnus Support, drepper@cygnus.com. * Cygnus Support, drepper@cygnus.com.
*/ */
...@@ -22,8 +22,8 @@ ...@@ -22,8 +22,8 @@
* atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) * atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
* 2 1 - x 1 - x * 2 1 - x 1 - x
* *
* For x<0.5 * For x<0.5
* atanhl(x) = 0.5*log1pl(2x+2x*x/(1-x)) * atanhl(x) = 0.5*log1pq(2x+2x*x/(1-x))
* *
* Special cases: * Special cases:
* atanhl(x) is NaN if |x| > 1 with signal; * atanhl(x) is NaN if |x| > 1 with signal;
...@@ -34,11 +34,12 @@ ...@@ -34,11 +34,12 @@
#include "quadmath-imp.h" #include "quadmath-imp.h"
static const __float128 one = 1.0Q, huge = 1e4900Q; static const __float128 one = 1, huge = 1e4900Q;
static const __float128 zero = 0.0Q;
static const __float128 zero = 0;
__float128 __float128
atanhq (__float128 x) atanhq(__float128 x)
{ {
__float128 t; __float128 t;
uint32_t jx, ix; uint32_t jx, ix;
......
/* s_atanl.c /* s_atanl.c
* *
* Inverse circular tangent for 128-bit __float128 precision * Inverse circular tangent for 128-bit long double precision
* (arctangent) * (arctangent)
* *
* *
* *
* SYNOPSIS: * SYNOPSIS:
* *
* __float128 x, y, atanl(); * long double x, y, atanq();
* *
* y = atanl( x ); * y = atanq( x );
* *
* *
* *
...@@ -55,9 +55,8 @@ ...@@ -55,9 +55,8 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software License along with this library; if not, see
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ <http://www.gnu.org/licenses/>. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
...@@ -228,7 +227,7 @@ atanq (__float128 x) ...@@ -228,7 +227,7 @@ atanq (__float128 x)
/* Index of nearest table element. /* Index of nearest table element.
Roundoff to integer is asymmetrical to avoid cancellation when t < 0 Roundoff to integer is asymmetrical to avoid cancellation when t < 0
(cf. fdlibm). */ (cf. fdlibm). */
k = 8.0Q * x + 0.25Q; k = 8.0 * x + 0.25;
u = 0.125Q * k; u = 0.125Q * k;
/* Small arctan argument. */ /* Small arctan argument. */
t = (x - u) / (1.0 + x * u); t = (x - u) / (1.0 + x * u);
......
/* Return arc hyperbole cosine for __float128 value. /* Return arc hyperbolic cosine for a complex type.
Copyright (C) 1997, 1998, 2006 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -14,13 +14,11 @@ ...@@ -14,13 +14,11 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
__complex128 __complex128
cacoshq (__complex128 x) cacoshq (__complex128 x)
{ {
...@@ -38,57 +36,54 @@ cacoshq (__complex128 x) ...@@ -38,57 +36,54 @@ cacoshq (__complex128 x)
__imag__ res = nanq (""); __imag__ res = nanq ("");
else else
__imag__ res = copysignq ((rcls == QUADFP_INFINITE __imag__ res = copysignq ((rcls == QUADFP_INFINITE
? (__real__ x < 0.0 ? (__real__ x < 0
? M_PIq - M_PI_4q : M_PI_4q) ? M_PIq - M_PI_4q
: M_PI_2q), __imag__ x); : M_PI_4q)
: M_PI_2q), __imag__ x);
} }
else if (rcls == QUADFP_INFINITE) else if (rcls == QUADFP_INFINITE)
{ {
__real__ res = HUGE_VALQ; __real__ res = HUGE_VALQ;
if (icls >= QUADFP_ZERO) if (icls >= QUADFP_ZERO)
__imag__ res = copysignq (signbitq (__real__ x) ? M_PIq : 0.0, __imag__ res = copysignq (signbitq (__real__ x)
__imag__ x); ? M_PIq : 0, __imag__ x);
else else
__imag__ res = nanq (""); __imag__ res = nanq ("");
} }
else else
{ {
__real__ res = nanq (""); __real__ res = nanq ("");
__imag__ res = nanq (""); if (rcls == QUADFP_ZERO)
__imag__ res = M_PI_2q;
else
__imag__ res = nanq ("");
} }
} }
else if (rcls == QUADFP_ZERO && icls == QUADFP_ZERO) else if (rcls == QUADFP_ZERO && icls == QUADFP_ZERO)
{ {
__real__ res = 0.0; __real__ res = 0;
__imag__ res = copysignq (M_PI_2q, __imag__ x); __imag__ res = copysignq (M_PI_2q, __imag__ x);
} }
/* The factor 16 is just a guess. */
else if (16.0Q * fabsq (__imag__ x) < fabsq (__real__ x))
{
/* Kahan's formula which avoid cancellation through subtraction in
some cases. */
res = 2.0Q * clogq (csqrtq ((x + 1.0Q) / 2.0Q)
+ csqrtq ((x - 1.0Q) / 2.0Q));
if (signbitq (__real__ res))
__real__ res = 0.0Q;
}
else else
{ {
__complex128 y; __complex128 y;
__real__ y = (__real__ x - __imag__ x) * (__real__ x + __imag__ x) - 1.0; __real__ y = -__imag__ x;
__imag__ y = 2.0 * __real__ x * __imag__ x; __imag__ y = __real__ x;
y = csqrtq (y);
if (signbitq (x)) y = __quadmath_kernel_casinhq (y, 1);
y = -y;
__real__ y += __real__ x; if (signbitq (__imag__ x))
__imag__ y += __imag__ x; {
__real__ res = __real__ y;
res = clogq (y); __imag__ res = -__imag__ y;
}
else
{
__real__ res = -__real__ y;
__imag__ res = __imag__ y;
}
} }
return res; return res;
......
/* Return cosine of complex __float128 value. /* Return cosine of a complex type.
Copyright (C) 1997, 1998 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -14,9 +14,8 @@ ...@@ -14,9 +14,8 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
...@@ -25,11 +24,29 @@ cacosq (__complex128 x) ...@@ -25,11 +24,29 @@ cacosq (__complex128 x)
{ {
__complex128 y; __complex128 y;
__complex128 res; __complex128 res;
int rcls = fpclassifyq (__real__ x);
y = casinq (x); int icls = fpclassifyq (__imag__ x);
__real__ res = M_PI_2q - __real__ y; if (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE
__imag__ res = -__imag__ y; || (rcls == QUADFP_ZERO && icls == QUADFP_ZERO))
{
y = casinq (x);
__real__ res = (__float128) M_PI_2q - __real__ y;
if (__real__ res == 0)
__real__ res = 0;
__imag__ res = -__imag__ y;
}
else
{
__real__ y = -__imag__ x;
__imag__ y = __real__ x;
y = __quadmath_kernel_casinhq (y, 1);
__real__ res = __imag__ y;
__imag__ res = __real__ y;
}
return res; return res;
} }
/* Return arc hyperbole sine for __float128 value. /* Return arc hyperbolic sine for a complex float type.
Copyright (C) 1997, 1998 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -14,13 +14,11 @@ ...@@ -14,13 +14,11 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
__complex128 __complex128
casinhq (__complex128 x) casinhq (__complex128 x)
{ {
...@@ -37,15 +35,16 @@ casinhq (__complex128 x) ...@@ -37,15 +35,16 @@ casinhq (__complex128 x)
if (rcls == QUADFP_NAN) if (rcls == QUADFP_NAN)
__imag__ res = nanq (""); __imag__ res = nanq ("");
else else
__imag__ res = copysignq (rcls >= QUADFP_ZERO ? M_PI_2q : M_PI_4q, __imag__ res = copysignq ((rcls >= QUADFP_ZERO
__imag__ x); ? M_PI_2q : M_PI_4q),
__imag__ x);
} }
else if (rcls <= QUADFP_INFINITE) else if (rcls <= QUADFP_INFINITE)
{ {
__real__ res = __real__ x; __real__ res = __real__ x;
if ((rcls == QUADFP_INFINITE && icls >= QUADFP_ZERO) if ((rcls == QUADFP_INFINITE && icls >= QUADFP_ZERO)
|| (rcls == QUADFP_NAN && icls == QUADFP_ZERO)) || (rcls == QUADFP_NAN && icls == QUADFP_ZERO))
__imag__ res = copysignq (0.0, __imag__ x); __imag__ res = copysignq (0, __imag__ x);
else else
__imag__ res = nanq (""); __imag__ res = nanq ("");
} }
...@@ -61,22 +60,7 @@ casinhq (__complex128 x) ...@@ -61,22 +60,7 @@ casinhq (__complex128 x)
} }
else else
{ {
__complex128 y; res = __quadmath_kernel_casinhq (x, 0);
__real__ y = (__real__ x - __imag__ x) * (__real__ x + __imag__ x) + 1.0;
__imag__ y = 2.0 * __real__ x * __imag__ x;
y = csqrtq (y);
__real__ y += __real__ x;
__imag__ y += __imag__ x;
res = clogq (y);
/* Ensure zeros have correct sign and results are correct if
very close to branch cuts. */
__real__ res = copysignq (__real__ res, __real__ x);
__imag__ res = copysignq (__imag__ res, __imag__ x);
} }
return res; return res;
......
/* Return arc hyperbolic sine for a complex float type, with the
imaginary part of the result possibly adjusted for use in
computing other functions.
Copyright (C) 1997-2018 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 "quadmath-imp.h"
/* Return the complex inverse hyperbolic sine of finite nonzero Z,
with the imaginary part of the result subtracted from pi/2 if ADJ
is nonzero. */
__complex128
__quadmath_kernel_casinhq (__complex128 x, int adj)
{
__complex128 res;
__float128 rx, ix;
__complex128 y;
/* Avoid cancellation by reducing to the first quadrant. */
rx = fabsq (__real__ x);
ix = fabsq (__imag__ x);
if (rx >= 1 / FLT128_EPSILON || ix >= 1 / FLT128_EPSILON)
{
/* For large x in the first quadrant, x + csqrt (1 + x * x)
is sufficiently close to 2 * x to make no significant
difference to the result; avoid possible overflow from
the squaring and addition. */
__real__ y = rx;
__imag__ y = ix;
if (adj)
{
__float128 t = __real__ y;
__real__ y = copysignq (__imag__ y, __imag__ x);
__imag__ y = t;
}
res = clogq (y);
__real__ res += (__float128) M_LN2q;
}
else if (rx >= 0.5Q && ix < FLT128_EPSILON / 8)
{
__float128 s = hypotq (1, rx);
__real__ res = logq (rx + s);
if (adj)
__imag__ res = atan2q (s, __imag__ x);
else
__imag__ res = atan2q (ix, s);
}
else if (rx < FLT128_EPSILON / 8 && ix >= 1.5Q)
{
__float128 s = sqrtq ((ix + 1) * (ix - 1));
__real__ res = logq (ix + s);
if (adj)
__imag__ res = atan2q (rx, copysignq (s, __imag__ x));
else
__imag__ res = atan2q (s, rx);
}
else if (ix > 1 && ix < 1.5Q && rx < 0.5Q)
{
if (rx < FLT128_EPSILON * FLT128_EPSILON)
{
__float128 ix2m1 = (ix + 1) * (ix - 1);
__float128 s = sqrtq (ix2m1);
__real__ res = log1pq (2 * (ix2m1 + ix * s)) / 2;
if (adj)
__imag__ res = atan2q (rx, copysignq (s, __imag__ x));
else
__imag__ res = atan2q (s, rx);
}
else
{
__float128 ix2m1 = (ix + 1) * (ix - 1);
__float128 rx2 = rx * rx;
__float128 f = rx2 * (2 + rx2 + 2 * ix * ix);
__float128 d = sqrtq (ix2m1 * ix2m1 + f);
__float128 dp = d + ix2m1;
__float128 dm = f / dp;
__float128 r1 = sqrtq ((dm + rx2) / 2);
__float128 r2 = rx * ix / r1;
__real__ res = log1pq (rx2 + dp + 2 * (rx * r1 + ix * r2)) / 2;
if (adj)
__imag__ res = atan2q (rx + r1, copysignq (ix + r2, __imag__ x));
else
__imag__ res = atan2q (ix + r2, rx + r1);
}
}
else if (ix == 1 && rx < 0.5Q)
{
if (rx < FLT128_EPSILON / 8)
{
__real__ res = log1pq (2 * (rx + sqrtq (rx))) / 2;
if (adj)
__imag__ res = atan2q (sqrtq (rx), copysignq (1, __imag__ x));
else
__imag__ res = atan2q (1, sqrtq (rx));
}
else
{
__float128 d = rx * sqrtq (4 + rx * rx);
__float128 s1 = sqrtq ((d + rx * rx) / 2);
__float128 s2 = sqrtq ((d - rx * rx) / 2);
__real__ res = log1pq (rx * rx + d + 2 * (rx * s1 + s2)) / 2;
if (adj)
__imag__ res = atan2q (rx + s1, copysignq (1 + s2, __imag__ x));
else
__imag__ res = atan2q (1 + s2, rx + s1);
}
}
else if (ix < 1 && rx < 0.5Q)
{
if (ix >= FLT128_EPSILON)
{
if (rx < FLT128_EPSILON * FLT128_EPSILON)
{
__float128 onemix2 = (1 + ix) * (1 - ix);
__float128 s = sqrtq (onemix2);
__real__ res = log1pq (2 * rx / s) / 2;
if (adj)
__imag__ res = atan2q (s, __imag__ x);
else
__imag__ res = atan2q (ix, s);
}
else
{
__float128 onemix2 = (1 + ix) * (1 - ix);
__float128 rx2 = rx * rx;
__float128 f = rx2 * (2 + rx2 + 2 * ix * ix);
__float128 d = sqrtq (onemix2 * onemix2 + f);
__float128 dp = d + onemix2;
__float128 dm = f / dp;
__float128 r1 = sqrtq ((dp + rx2) / 2);
__float128 r2 = rx * ix / r1;
__real__ res = log1pq (rx2 + dm + 2 * (rx * r1 + ix * r2)) / 2;
if (adj)
__imag__ res = atan2q (rx + r1, copysignq (ix + r2,
__imag__ x));
else
__imag__ res = atan2q (ix + r2, rx + r1);
}
}
else
{
__float128 s = hypotq (1, rx);
__real__ res = log1pq (2 * rx * (rx + s)) / 2;
if (adj)
__imag__ res = atan2q (s, __imag__ x);
else
__imag__ res = atan2q (ix, s);
}
math_check_force_underflow_nonneg (__real__ res);
}
else
{
__real__ y = (rx - ix) * (rx + ix) + 1;
__imag__ y = 2 * rx * ix;
y = csqrtq (y);
__real__ y += rx;
__imag__ y += ix;
if (adj)
{
__float128 t = __real__ y;
__real__ y = copysignq (__imag__ y, __imag__ x);
__imag__ y = t;
}
res = clogq (y);
}
/* Give results the correct sign for the original argument. */
__real__ res = copysignq (__real__ res, __real__ x);
__imag__ res = copysignq (__imag__ res, (adj ? 1 : __imag__ x));
return res;
}
/* Return arc sine of complex __float128 value. /* Return arc sine of a complex float type.
Copyright (C) 1997 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -14,13 +14,11 @@ ...@@ -14,13 +14,11 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
__complex128 __complex128
casinq (__complex128 x) casinq (__complex128 x)
{ {
...@@ -28,7 +26,7 @@ casinq (__complex128 x) ...@@ -28,7 +26,7 @@ casinq (__complex128 x)
if (isnanq (__real__ x) || isnanq (__imag__ x)) if (isnanq (__real__ x) || isnanq (__imag__ x))
{ {
if (__real__ x == 0.0) if (__real__ x == 0)
{ {
res = x; res = x;
} }
......
/* Return arc hyperbole tangent for __float128 value. /* Return arc hyperbolic tangent for a complex float type.
Copyright (C) 1997, 1998 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -14,13 +14,11 @@ ...@@ -14,13 +14,11 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
__complex128 __complex128
catanhq (__complex128 x) catanhq (__complex128 x)
{ {
...@@ -28,16 +26,16 @@ catanhq (__complex128 x) ...@@ -28,16 +26,16 @@ catanhq (__complex128 x)
int rcls = fpclassifyq (__real__ x); int rcls = fpclassifyq (__real__ x);
int icls = fpclassifyq (__imag__ x); int icls = fpclassifyq (__imag__ x);
if (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE) if (__glibc_unlikely (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE))
{ {
if (icls == QUADFP_INFINITE) if (icls == QUADFP_INFINITE)
{ {
__real__ res = copysignq (0.0, __real__ x); __real__ res = copysignq (0, __real__ x);
__imag__ res = copysignq (M_PI_2q, __imag__ x); __imag__ res = copysignq (M_PI_2q, __imag__ x);
} }
else if (rcls == QUADFP_INFINITE || rcls == QUADFP_ZERO) else if (rcls == QUADFP_INFINITE || rcls == QUADFP_ZERO)
{ {
__real__ res = copysignq (0.0, __real__ x); __real__ res = copysignq (0, __real__ x);
if (icls >= QUADFP_ZERO) if (icls >= QUADFP_ZERO)
__imag__ res = copysignq (M_PI_2q, __imag__ x); __imag__ res = copysignq (M_PI_2q, __imag__ x);
else else
...@@ -49,27 +47,83 @@ catanhq (__complex128 x) ...@@ -49,27 +47,83 @@ catanhq (__complex128 x)
__imag__ res = nanq (""); __imag__ res = nanq ("");
} }
} }
else if (rcls == QUADFP_ZERO && icls == QUADFP_ZERO) else if (__glibc_unlikely (rcls == QUADFP_ZERO && icls == QUADFP_ZERO))
{ {
res = x; res = x;
} }
else else
{ {
__float128 i2, num, den; if (fabsq (__real__ x) >= 16 / FLT128_EPSILON
|| fabsq (__imag__ x) >= 16 / FLT128_EPSILON)
{
__imag__ res = copysignq (M_PI_2q, __imag__ x);
if (fabsq (__imag__ x) <= 1)
__real__ res = 1 / __real__ x;
else if (fabsq (__real__ x) <= 1)
__real__ res = __real__ x / __imag__ x / __imag__ x;
else
{
__float128 h = hypotq (__real__ x / 2, __imag__ x / 2);
__real__ res = __real__ x / h / h / 4;
}
}
else
{
if (fabsq (__real__ x) == 1
&& fabsq (__imag__ x) < FLT128_EPSILON * FLT128_EPSILON)
__real__ res = (copysignq (0.5Q, __real__ x)
* ((__float128) M_LN2q
- logq (fabsq (__imag__ x))));
else
{
__float128 i2 = 0;
if (fabsq (__imag__ x) >= FLT128_EPSILON * FLT128_EPSILON)
i2 = __imag__ x * __imag__ x;
__float128 num = 1 + __real__ x;
num = i2 + num * num;
i2 = __imag__ x * __imag__ x; __float128 den = 1 - __real__ x;
den = i2 + den * den;
num = 1.0 + __real__ x; __float128 f = num / den;
num = i2 + num * num; if (f < 0.5Q)
__real__ res = 0.25Q * logq (f);
else
{
num = 4 * __real__ x;
__real__ res = 0.25Q * log1pq (num / den);
}
}
den = 1.0 - __real__ x; __float128 absx, absy, den;
den = i2 + den * den;
__real__ res = 0.25 * (logq (num) - logq (den)); absx = fabsq (__real__ x);
absy = fabsq (__imag__ x);
if (absx < absy)
{
__float128 t = absx;
absx = absy;
absy = t;
}
den = 1 - __real__ x * __real__ x - i2; if (absy < FLT128_EPSILON / 2)
{
den = (1 - absx) * (1 + absx);
if (den == 0)
den = 0;
}
else if (absx >= 1)
den = (1 - absx) * (1 + absx) - absy * absy;
else if (absx >= 0.75Q || absy >= 0.5Q)
den = -__quadmath_x2y2m1q (absx, absy);
else
den = (1 - absx) * (1 + absx) - absy * absy;
__imag__ res = 0.5Q * atan2q (2 * __imag__ x, den);
}
__imag__ res = 0.5 * atan2q (2.0 * __imag__ x, den); math_check_force_underflow_complex (res);
} }
return res; return res;
......
/* Return arc tangent of complex __float128 value. /* Return arc tangent of complex float type.
Copyright (C) 1997, 1998 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -14,13 +14,11 @@ ...@@ -14,13 +14,11 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
__complex128 __complex128
catanq (__complex128 x) catanq (__complex128 x)
{ {
...@@ -28,12 +26,12 @@ catanq (__complex128 x) ...@@ -28,12 +26,12 @@ catanq (__complex128 x)
int rcls = fpclassifyq (__real__ x); int rcls = fpclassifyq (__real__ x);
int icls = fpclassifyq (__imag__ x); int icls = fpclassifyq (__imag__ x);
if (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE) if (__glibc_unlikely (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE))
{ {
if (rcls == QUADFP_INFINITE) if (rcls == QUADFP_INFINITE)
{ {
__real__ res = copysignq (M_PI_2q, __real__ x); __real__ res = copysignq (M_PI_2q, __real__ x);
__imag__ res = copysignq (0.0, __imag__ x); __imag__ res = copysignq (0, __imag__ x);
} }
else if (icls == QUADFP_INFINITE) else if (icls == QUADFP_INFINITE)
{ {
...@@ -41,12 +39,12 @@ catanq (__complex128 x) ...@@ -41,12 +39,12 @@ catanq (__complex128 x)
__real__ res = copysignq (M_PI_2q, __real__ x); __real__ res = copysignq (M_PI_2q, __real__ x);
else else
__real__ res = nanq (""); __real__ res = nanq ("");
__imag__ res = copysignq (0.0, __imag__ x); __imag__ res = copysignq (0, __imag__ x);
} }
else if (icls == QUADFP_ZERO || icls == QUADFP_INFINITE) else if (icls == QUADFP_ZERO || icls == QUADFP_INFINITE)
{ {
__real__ res = nanq (""); __real__ res = nanq ("");
__imag__ res = copysignq (0.0, __imag__ x); __imag__ res = copysignq (0, __imag__ x);
} }
else else
{ {
...@@ -54,27 +52,84 @@ catanq (__complex128 x) ...@@ -54,27 +52,84 @@ catanq (__complex128 x)
__imag__ res = nanq (""); __imag__ res = nanq ("");
} }
} }
else if (rcls == QUADFP_ZERO && icls == QUADFP_ZERO) else if (__glibc_unlikely (rcls == QUADFP_ZERO && icls == QUADFP_ZERO))
{ {
res = x; res = x;
} }
else else
{ {
__float128 r2, num, den; if (fabsq (__real__ x) >= 16 / FLT128_EPSILON
|| fabsq (__imag__ x) >= 16 / FLT128_EPSILON)
{
__real__ res = copysignq (M_PI_2q, __real__ x);
if (fabsq (__real__ x) <= 1)
__imag__ res = 1 / __imag__ x;
else if (fabsq (__imag__ x) <= 1)
__imag__ res = __imag__ x / __real__ x / __real__ x;
else
{
__float128 h = hypotq (__real__ x / 2, __imag__ x / 2);
__imag__ res = __imag__ x / h / h / 4;
}
}
else
{
__float128 den, absx, absy;
r2 = __real__ x * __real__ x; absx = fabsq (__real__ x);
absy = fabsq (__imag__ x);
if (absx < absy)
{
__float128 t = absx;
absx = absy;
absy = t;
}
den = 1 - r2 - __imag__ x * __imag__ x; if (absy < FLT128_EPSILON / 2)
{
den = (1 - absx) * (1 + absx);
if (den == 0)
den = 0;
}
else if (absx >= 1)
den = (1 - absx) * (1 + absx) - absy * absy;
else if (absx >= 0.75Q || absy >= 0.5Q)
den = -__quadmath_x2y2m1q (absx, absy);
else
den = (1 - absx) * (1 + absx) - absy * absy;
__real__ res = 0.5Q * atan2q (2 * __real__ x, den);
if (fabsq (__imag__ x) == 1
&& fabsq (__real__ x) < FLT128_EPSILON * FLT128_EPSILON)
__imag__ res = (copysignq (0.5Q, __imag__ x)
* ((__float128) M_LN2q
- logq (fabsq (__real__ x))));
else
{
__float128 r2 = 0, num, f;
if (fabsq (__real__ x) >= FLT128_EPSILON * FLT128_EPSILON)
r2 = __real__ x * __real__ x;
__real__ res = 0.5 * atan2q (2.0 * __real__ x, den); num = __imag__ x + 1;
num = r2 + num * num;
num = __imag__ x + 1.0; den = __imag__ x - 1;
num = r2 + num * num; den = r2 + den * den;
den = __imag__ x - 1.0; f = num / den;
den = r2 + den * den; if (f < 0.5Q)
__imag__ res = 0.25Q * logq (f);
else
{
num = 4 * __imag__ x;
__imag__ res = 0.25Q * log1pq (num / den);
}
}
}
__imag__ res = 0.25 * logq (num / den); math_check_force_underflow_complex (res);
} }
return res; return res;
......
/* cbrtq.c /* cbrtq.c
* *
* Cube root, __float128 precision * Cube root, long double precision
* *
* *
* *
* SYNOPSIS: * SYNOPSIS:
* *
* __float128 x, y, cbrtq(); * long double x, y, cbrtq();
* *
* y = cbrtq( x ); * y = cbrtq( x );
* *
...@@ -53,7 +53,6 @@ Adapted for glibc October, 2001. ...@@ -53,7 +53,6 @@ Adapted for glibc October, 2001.
License along with this library; if not, see License along with this library; if not, see
<http://www.gnu.org/licenses/>. */ <http://www.gnu.org/licenses/>. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
static const __float128 CBRT2 = 1.259921049894873164767210607278228350570251Q; static const __float128 CBRT2 = 1.259921049894873164767210607278228350570251Q;
...@@ -63,7 +62,7 @@ static const __float128 CBRT4I = 0.6299605249474365823836053036391141752851257Q; ...@@ -63,7 +62,7 @@ static const __float128 CBRT4I = 0.6299605249474365823836053036391141752851257Q;
__float128 __float128
cbrtq ( __float128 x) cbrtq (__float128 x)
{ {
int e, rem, sign; int e, rem, sign;
__float128 z; __float128 z;
......
/* Complex cosine hyperbole function for complex __float128. /* Complex cosine hyperbolic function for float types.
Copyright (C) 1997-2012 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -19,11 +19,6 @@ ...@@ -19,11 +19,6 @@
#include "quadmath-imp.h" #include "quadmath-imp.h"
#ifdef HAVE_FENV_H
# include <fenv.h>
#endif
__complex128 __complex128
ccoshq (__complex128 x) ccoshq (__complex128 x)
{ {
...@@ -31,23 +26,23 @@ ccoshq (__complex128 x) ...@@ -31,23 +26,23 @@ ccoshq (__complex128 x)
int rcls = fpclassifyq (__real__ x); int rcls = fpclassifyq (__real__ x);
int icls = fpclassifyq (__imag__ x); int icls = fpclassifyq (__imag__ x);
if (__builtin_expect (rcls >= QUADFP_ZERO, 1)) if (__glibc_likely (rcls >= QUADFP_ZERO))
{ {
/* Real part is finite. */ /* Real part is finite. */
if (__builtin_expect (icls >= QUADFP_ZERO, 1)) if (__glibc_likely (icls >= QUADFP_ZERO))
{ {
/* Imaginary part is finite. */ /* Imaginary part is finite. */
const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q); const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q);
__float128 sinix, cosix; __float128 sinix, cosix;
if (__builtin_expect (icls != QUADFP_SUBNORMAL, 1)) if (__glibc_likely (fabsq (__imag__ x) > FLT128_MIN))
{ {
sincosq (__imag__ x, &sinix, &cosix); sincosq (__imag__ x, &sinix, &cosix);
} }
else else
{ {
sinix = __imag__ x; sinix = __imag__ x;
cosix = 1.0Q; cosix = 1;
} }
if (fabsq (__real__ x) > t) if (fabsq (__real__ x) > t)
...@@ -57,8 +52,8 @@ ccoshq (__complex128 x) ...@@ -57,8 +52,8 @@ ccoshq (__complex128 x)
if (signbitq (__real__ x)) if (signbitq (__real__ x))
sinix = -sinix; sinix = -sinix;
rx -= t; rx -= t;
sinix *= exp_t / 2.0Q; sinix *= exp_t / 2;
cosix *= exp_t / 2.0Q; cosix *= exp_t / 2;
if (rx > t) if (rx > t)
{ {
rx -= t; rx -= t;
...@@ -83,62 +78,53 @@ ccoshq (__complex128 x) ...@@ -83,62 +78,53 @@ ccoshq (__complex128 x)
__real__ retval = coshq (__real__ x) * cosix; __real__ retval = coshq (__real__ x) * cosix;
__imag__ retval = sinhq (__real__ x) * sinix; __imag__ retval = sinhq (__real__ x) * sinix;
} }
math_check_force_underflow_complex (retval);
} }
else else
{ {
__imag__ retval = __real__ x == 0.0Q ? 0.0Q : nanq (""); __imag__ retval = __real__ x == 0 ? 0 : nanq ("");
__real__ retval = nanq ("") + nanq (""); __real__ retval = __imag__ x - __imag__ x;
}
#ifdef HAVE_FENV_H
if (icls == QUADFP_INFINITE)
feraiseexcept (FE_INVALID);
#endif
}
} }
else if (rcls == QUADFP_INFINITE) else if (rcls == QUADFP_INFINITE)
{ {
/* Real part is infinite. */ /* Real part is infinite. */
if (__builtin_expect (icls > QUADFP_ZERO, 1)) if (__glibc_likely (icls > QUADFP_ZERO))
{ {
/* Imaginary part is finite. */ /* Imaginary part is finite. */
__float128 sinix, cosix; __float128 sinix, cosix;
if (__builtin_expect (icls != QUADFP_SUBNORMAL, 1)) if (__glibc_likely (fabsq (__imag__ x) > FLT128_MIN))
{ {
sincosq (__imag__ x, &sinix, &cosix); sincosq (__imag__ x, &sinix, &cosix);
} }
else else
{ {
sinix = __imag__ x; sinix = __imag__ x;
cosix = 1.0Q; cosix = 1;
} }
__real__ retval = copysignq (HUGE_VALQ, cosix); __real__ retval = copysignq (HUGE_VALQ, cosix);
__imag__ retval = (copysignq (HUGE_VALQ, sinix) __imag__ retval = (copysignq (HUGE_VALQ, sinix)
* copysignq (1.0Q, __real__ x)); * copysignq (1, __real__ x));
} }
else if (icls == QUADFP_ZERO) else if (icls == QUADFP_ZERO)
{ {
/* Imaginary part is 0.0. */ /* Imaginary part is 0.0. */
__real__ retval = HUGE_VALQ; __real__ retval = HUGE_VALQ;
__imag__ retval = __imag__ x * copysignq (1.0Q, __real__ x); __imag__ retval = __imag__ x * copysignq (1, __real__ x);
} }
else else
{ {
/* The addition raises the invalid exception. */
__real__ retval = HUGE_VALQ; __real__ retval = HUGE_VALQ;
__imag__ retval = nanq ("") + nanq (""); __imag__ retval = __imag__ x - __imag__ x;
}
#ifdef HAVE_FENV_H
if (icls == QUADFP_INFINITE)
feraiseexcept (FE_INVALID);
#endif
}
} }
else else
{ {
__real__ retval = nanq (""); __real__ retval = nanq ("");
__imag__ retval = __imag__ x == 0.0 ? __imag__ x : nanq (""); __imag__ retval = __imag__ x == 0 ? __imag__ x : nanq ("");
} }
return retval; return retval;
......
/* ceilq.c -- __float128 version of s_ceil.c. /* s_ceill.c -- long double version of s_ceil.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz. * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/ */
...@@ -13,10 +13,22 @@ ...@@ -13,10 +13,22 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/*
* ceilq(x)
* Return x rounded toward -inf to integral value
* Method:
* Bit twiddling.
*/
#define NO_MATH_REDIRECT
#include "quadmath-imp.h" #include "quadmath-imp.h"
__float128 __float128 ceilq(__float128 x)
ceilq (__float128 x)
{ {
int64_t i0,i1,j0; int64_t i0,i1,j0;
uint64_t i,j; uint64_t i,j;
......
/* Return value of complex exponential function for complex __float128 value. /* Return value of complex exponential function for a float type.
Copyright (C) 1997-2012 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -19,11 +19,6 @@ ...@@ -19,11 +19,6 @@
#include "quadmath-imp.h" #include "quadmath-imp.h"
#ifdef HAVE_FENV_H
# include <fenv.h>
#endif
__complex128 __complex128
cexpq (__complex128 x) cexpq (__complex128 x)
{ {
...@@ -31,23 +26,23 @@ cexpq (__complex128 x) ...@@ -31,23 +26,23 @@ cexpq (__complex128 x)
int rcls = fpclassifyq (__real__ x); int rcls = fpclassifyq (__real__ x);
int icls = fpclassifyq (__imag__ x); int icls = fpclassifyq (__imag__ x);
if (__builtin_expect (rcls >= QUADFP_ZERO, 1)) if (__glibc_likely (rcls >= QUADFP_ZERO))
{ {
/* Real part is finite. */ /* Real part is finite. */
if (__builtin_expect (icls >= QUADFP_ZERO, 1)) if (__glibc_likely (icls >= QUADFP_ZERO))
{ {
/* Imaginary part is finite. */ /* Imaginary part is finite. */
const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q); const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q);
__float128 sinix, cosix; __float128 sinix, cosix;
if (__builtin_expect (icls != QUADFP_SUBNORMAL, 1)) if (__glibc_likely (fabsq (__imag__ x) > FLT128_MIN))
{ {
sincosq (__imag__ x, &sinix, &cosix); sincosq (__imag__ x, &sinix, &cosix);
} }
else else
{ {
sinix = __imag__ x; sinix = __imag__ x;
cosix = 1.0Q; cosix = 1;
} }
if (__real__ x > t) if (__real__ x > t)
...@@ -75,6 +70,7 @@ cexpq (__complex128 x) ...@@ -75,6 +70,7 @@ cexpq (__complex128 x)
__real__ retval = exp_val * cosix; __real__ retval = exp_val * cosix;
__imag__ retval = exp_val * sinix; __imag__ retval = exp_val * sinix;
} }
math_check_force_underflow_complex (retval);
} }
else else
{ {
...@@ -83,18 +79,16 @@ cexpq (__complex128 x) ...@@ -83,18 +79,16 @@ cexpq (__complex128 x)
__real__ retval = nanq (""); __real__ retval = nanq ("");
__imag__ retval = nanq (""); __imag__ retval = nanq ("");
#ifdef HAVE_FENV_H
feraiseexcept (FE_INVALID); feraiseexcept (FE_INVALID);
#endif
} }
} }
else if (__builtin_expect (rcls == QUADFP_INFINITE, 1)) else if (__glibc_likely (rcls == QUADFP_INFINITE))
{ {
/* Real part is infinite. */ /* Real part is infinite. */
if (__builtin_expect (icls >= QUADFP_ZERO, 1)) if (__glibc_likely (icls >= QUADFP_ZERO))
{ {
/* Imaginary part is finite. */ /* Imaginary part is finite. */
__float128 value = signbitq (__real__ x) ? 0.0Q : HUGE_VALQ; __float128 value = signbitq (__real__ x) ? 0 : HUGE_VALQ;
if (icls == QUADFP_ZERO) if (icls == QUADFP_ZERO)
{ {
...@@ -106,14 +100,14 @@ cexpq (__complex128 x) ...@@ -106,14 +100,14 @@ cexpq (__complex128 x)
{ {
__float128 sinix, cosix; __float128 sinix, cosix;
if (__builtin_expect (icls != QUADFP_SUBNORMAL, 1)) if (__glibc_likely (fabsq (__imag__ x) > FLT128_MIN))
{ {
sincosq (__imag__ x, &sinix, &cosix); sincosq (__imag__ x, &sinix, &cosix);
} }
else else
{ {
sinix = __imag__ x; sinix = __imag__ x;
cosix = 1.0Q; cosix = 1;
} }
__real__ retval = copysignq (value, cosix); __real__ retval = copysignq (value, cosix);
...@@ -123,29 +117,28 @@ cexpq (__complex128 x) ...@@ -123,29 +117,28 @@ cexpq (__complex128 x)
else if (signbitq (__real__ x) == 0) else if (signbitq (__real__ x) == 0)
{ {
__real__ retval = HUGE_VALQ; __real__ retval = HUGE_VALQ;
__imag__ retval = nanq (""); __imag__ retval = __imag__ x - __imag__ x;
#ifdef HAVE_FENV_H
if (icls == QUADFP_INFINITE)
feraiseexcept (FE_INVALID);
#endif
} }
else else
{ {
__real__ retval = 0.0Q; __real__ retval = 0;
__imag__ retval = copysignq (0.0Q, __imag__ x); __imag__ retval = copysignq (0, __imag__ x);
} }
} }
else else
{ {
/* If the real part is NaN the result is NaN + iNaN. */ /* If the real part is NaN the result is NaN + iNaN unless the
imaginary part is zero. */
__real__ retval = nanq (""); __real__ retval = nanq ("");
__imag__ retval = nanq (""); if (icls == QUADFP_ZERO)
__imag__ retval = __imag__ x;
else
{
__imag__ retval = nanq ("");
#ifdef HAVE_FENV_H if (rcls != QUADFP_NAN || icls != QUADFP_NAN)
if (rcls != QUADFP_NAN || icls != QUADFP_NAN) feraiseexcept (FE_INVALID);
feraiseexcept (FE_INVALID); }
#endif
} }
return retval; return retval;
......
/* Return imaginary part of complex __float128 value. /* Return imaginary part of complex float type.
Copyright (C) 1997, 1998 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -14,9 +14,8 @@ ...@@ -14,9 +14,8 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
......
/* Compute complex base 10 logarithm for complex __float128. /* Compute complex base 10 logarithm.
Copyright (C) 1997-2012 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -19,10 +19,11 @@ ...@@ -19,10 +19,11 @@
#include "quadmath-imp.h" #include "quadmath-imp.h"
/* log_10 (2). */ /* log_10 (2). */
#define M_LOG10_2q 0.3010299956639811952137388947244930267682Q #define LOG10_2 0.3010299956639811952137388947244930267682Q
/* pi * log10 (e). */
#define PI_LOG10E 1.364376353841841347485783625431355770210Q
__complex128 __complex128
clog10q (__complex128 x) clog10q (__complex128 x)
...@@ -31,15 +32,15 @@ clog10q (__complex128 x) ...@@ -31,15 +32,15 @@ clog10q (__complex128 x)
int rcls = fpclassifyq (__real__ x); int rcls = fpclassifyq (__real__ x);
int icls = fpclassifyq (__imag__ x); int icls = fpclassifyq (__imag__ x);
if (__builtin_expect (rcls == QUADFP_ZERO && icls == QUADFP_ZERO, 0)) if (__glibc_unlikely (rcls == QUADFP_ZERO && icls == QUADFP_ZERO))
{ {
/* Real and imaginary part are 0.0. */ /* Real and imaginary part are 0.0. */
__imag__ result = signbitq (__real__ x) ? M_PIq : 0.0Q; __imag__ result = signbitq (__real__ x) ? PI_LOG10E : 0;
__imag__ result = copysignq (__imag__ result, __imag__ x); __imag__ result = copysignq (__imag__ result, __imag__ x);
/* Yes, the following line raises an exception. */ /* Yes, the following line raises an exception. */
__real__ result = -1.0Q / fabsq (__real__ x); __real__ result = -1 / fabsq (__real__ x);
} }
else if (__builtin_expect (rcls != QUADFP_NAN && icls != QUADFP_NAN, 1)) else if (__glibc_likely (rcls != QUADFP_NAN && icls != QUADFP_NAN))
{ {
/* Neither real nor imaginary part is NaN. */ /* Neither real nor imaginary part is NaN. */
__float128 absx = fabsq (__real__ x), absy = fabsq (__imag__ x); __float128 absx = fabsq (__real__ x), absy = fabsq (__imag__ x);
...@@ -52,11 +53,11 @@ clog10q (__complex128 x) ...@@ -52,11 +53,11 @@ clog10q (__complex128 x)
absy = t; absy = t;
} }
if (absx > FLT128_MAX / 2.0Q) if (absx > FLT128_MAX / 2)
{ {
scale = -1; scale = -1;
absx = scalbnq (absx, scale); absx = scalbnq (absx, scale);
absy = (absy >= FLT128_MIN * 2.0Q ? scalbnq (absy, scale) : 0.0Q); absy = (absy >= FLT128_MIN * 2 ? scalbnq (absy, scale) : 0);
} }
else if (absx < FLT128_MIN && absy < FLT128_MIN) else if (absx < FLT128_MIN && absy < FLT128_MIN)
{ {
...@@ -65,39 +66,39 @@ clog10q (__complex128 x) ...@@ -65,39 +66,39 @@ clog10q (__complex128 x)
absy = scalbnq (absy, scale); absy = scalbnq (absy, scale);
} }
if (absx == 1.0Q && scale == 0) if (absx == 1 && scale == 0)
{ {
__float128 absy2 = absy * absy; __real__ result = (log1pq (absy * absy)
if (absy2 <= FLT128_MIN * 2.0Q * M_LN10q) * ((__float128) M_LOG10Eq / 2));
__real__ result math_check_force_underflow_nonneg (__real__ result);
= (absy2 / 2.0Q - absy2 * absy2 / 4.0Q) * M_LOG10Eq;
else
__real__ result = log1pq (absy2) * (M_LOG10Eq / 2.0Q);
} }
else if (absx > 1.0Q && absx < 2.0Q && absy < 1.0Q && scale == 0) else if (absx > 1 && absx < 2 && absy < 1 && scale == 0)
{ {
__float128 d2m1 = (absx - 1.0Q) * (absx + 1.0Q); __float128 d2m1 = (absx - 1) * (absx + 1);
if (absy >= FLT128_EPSILON) if (absy >= FLT128_EPSILON)
d2m1 += absy * absy; d2m1 += absy * absy;
__real__ result = log1pq (d2m1) * (M_LOG10Eq / 2.0Q); __real__ result = log1pq (d2m1) * ((__float128) M_LOG10Eq / 2);
} }
else if (absx < 1.0Q else if (absx < 1
&& absx >= 0.75Q && absx >= 0.5Q
&& absy < FLT128_EPSILON / 2.0Q && absy < FLT128_EPSILON / 2
&& scale == 0) && scale == 0)
{ {
__float128 d2m1 = (absx - 1.0Q) * (absx + 1.0Q); __float128 d2m1 = (absx - 1) * (absx + 1);
__real__ result = log1pq (d2m1) * (M_LOG10Eq / 2.0Q); __real__ result = log1pq (d2m1) * ((__float128) M_LOG10Eq / 2);
} }
else if (absx < 1.0Q && (absx >= 0.75Q || absy >= 0.5Q) && scale == 0) else if (absx < 1
&& absx >= 0.5Q
&& scale == 0
&& absx * absx + absy * absy >= 0.5Q)
{ {
__float128 d2m1 = __quadmath_x2y2m1q (absx, absy); __float128 d2m1 = __quadmath_x2y2m1q (absx, absy);
__real__ result = log1pq (d2m1) * (M_LOG10Eq / 2.0Q); __real__ result = log1pq (d2m1) * ((__float128) M_LOG10Eq / 2);
} }
else else
{ {
__float128 d = hypotq (absx, absy); __float128 d = hypotq (absx, absy);
__real__ result = log10q (d) - scale * M_LOG10_2q; __real__ result = log10q (d) - scale * LOG10_2;
} }
__imag__ result = M_LOG10Eq * atan2q (__imag__ x, __real__ x); __imag__ result = M_LOG10Eq * atan2q (__imag__ x, __real__ x);
......
/* Compute complex natural logarithm for complex __float128. /* Compute complex natural logarithm.
Copyright (C) 1997-2012 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -19,7 +19,6 @@ ...@@ -19,7 +19,6 @@
#include "quadmath-imp.h" #include "quadmath-imp.h"
__complex128 __complex128
clogq (__complex128 x) clogq (__complex128 x)
{ {
...@@ -27,15 +26,15 @@ clogq (__complex128 x) ...@@ -27,15 +26,15 @@ clogq (__complex128 x)
int rcls = fpclassifyq (__real__ x); int rcls = fpclassifyq (__real__ x);
int icls = fpclassifyq (__imag__ x); int icls = fpclassifyq (__imag__ x);
if (__builtin_expect (rcls == QUADFP_ZERO && icls == QUADFP_ZERO, 0)) if (__glibc_unlikely (rcls == QUADFP_ZERO && icls == QUADFP_ZERO))
{ {
/* Real and imaginary part are 0.0. */ /* Real and imaginary part are 0.0. */
__imag__ result = signbitq (__real__ x) ? M_PIq : 0.0Q; __imag__ result = signbitq (__real__ x) ? (__float128) M_PIq : 0;
__imag__ result = copysignq (__imag__ result, __imag__ x); __imag__ result = copysignq (__imag__ result, __imag__ x);
/* Yes, the following line raises an exception. */ /* Yes, the following line raises an exception. */
__real__ result = -1.0Q / fabsq (__real__ x); __real__ result = -1 / fabsq (__real__ x);
} }
else if (__builtin_expect (rcls != QUADFP_NAN && icls != QUADFP_NAN, 1)) else if (__glibc_likely (rcls != QUADFP_NAN && icls != QUADFP_NAN))
{ {
/* Neither real nor imaginary part is NaN. */ /* Neither real nor imaginary part is NaN. */
__float128 absx = fabsq (__real__ x), absy = fabsq (__imag__ x); __float128 absx = fabsq (__real__ x), absy = fabsq (__imag__ x);
...@@ -48,11 +47,11 @@ clogq (__complex128 x) ...@@ -48,11 +47,11 @@ clogq (__complex128 x)
absy = t; absy = t;
} }
if (absx > FLT128_MAX / 2.0) if (absx > FLT128_MAX / 2)
{ {
scale = -1; scale = -1;
absx = scalbnq (absx, scale); absx = scalbnq (absx, scale);
absy = (absy >= FLT128_MIN * 2.0Q ? scalbnq (absy, scale) : 0.0Q); absy = (absy >= FLT128_MIN * 2 ? scalbnq (absy, scale) : 0);
} }
else if (absx < FLT128_MIN && absy < FLT128_MIN) else if (absx < FLT128_MIN && absy < FLT128_MIN)
{ {
...@@ -61,38 +60,38 @@ clogq (__complex128 x) ...@@ -61,38 +60,38 @@ clogq (__complex128 x)
absy = scalbnq (absy, scale); absy = scalbnq (absy, scale);
} }
if (absx == 1.0Q && scale == 0) if (absx == 1 && scale == 0)
{ {
__float128 absy2 = absy * absy; __real__ result = log1pq (absy * absy) / 2;
if (absy2 <= FLT128_MIN * 2.0Q) math_check_force_underflow_nonneg (__real__ result);
__real__ result = absy2 / 2.0Q - absy2 * absy2 / 4.0Q;
else
__real__ result = log1pq (absy2) / 2.0Q;
} }
else if (absx > 1.0Q && absx < 2.0Q && absy < 1.0Q && scale == 0) else if (absx > 1 && absx < 2 && absy < 1 && scale == 0)
{ {
__float128 d2m1 = (absx - 1.0Q) * (absx + 1.0Q); __float128 d2m1 = (absx - 1) * (absx + 1);
if (absy >= FLT128_EPSILON) if (absy >= FLT128_EPSILON)
d2m1 += absy * absy; d2m1 += absy * absy;
__real__ result = log1pq (d2m1) / 2.0Q; __real__ result = log1pq (d2m1) / 2;
} }
else if (absx < 1.0Q else if (absx < 1
&& absx >= 0.75Q && absx >= 0.5Q
&& absy < FLT128_EPSILON / 2.0Q && absy < FLT128_EPSILON / 2
&& scale == 0) && scale == 0)
{ {
__float128 d2m1 = (absx - 1.0Q) * (absx + 1.0Q); __float128 d2m1 = (absx - 1) * (absx + 1);
__real__ result = log1pq (d2m1) / 2.0Q; __real__ result = log1pq (d2m1) / 2;
} }
else if (absx < 1.0 && (absx >= 0.75Q || absy >= 0.5Q) && scale == 0) else if (absx < 1
&& absx >= 0.5Q
&& scale == 0
&& absx * absx + absy * absy >= 0.5Q)
{ {
__float128 d2m1 = __quadmath_x2y2m1q (absx, absy); __float128 d2m1 = __quadmath_x2y2m1q (absx, absy);
__real__ result = log1pq (d2m1) / 2.0Q; __real__ result = log1pq (d2m1) / 2;
} }
else else
{ {
__float128 d = hypotq (absx, absy); __float128 d = hypotq (absx, absy);
__real__ result = logq (d) - scale * M_LN2q; __real__ result = logq (d) - scale * (__float128) M_LN2q;
} }
__imag__ result = atan2q (__imag__ x, __real__ x); __imag__ result = atan2q (__imag__ x, __real__ x);
......
/* Return complex conjugate of complex __float128 value. /* Return complex conjugate of complex float type.
Copyright (C) 1997, 1998 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -14,9 +14,8 @@ ...@@ -14,9 +14,8 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
......
/* copysignq.c -- __float128 version of s_copysign.c. /* s_copysignl.c -- long double version of s_copysign.c.
* Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz. * Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz.
*/ */
...@@ -13,14 +13,26 @@ ...@@ -13,14 +13,26 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/*
* copysignq(long double x, long double y)
* copysignq(x,y) returns a value with the magnitude of x and
* with the sign bit of y.
*/
#define NO_MATH_REDIRECT
#include "quadmath-imp.h" #include "quadmath-imp.h"
__float128 __float128 copysignq(__float128 x, __float128 y)
copysignq (__float128 x, __float128 y)
{ {
uint64_t hx,hy; uint64_t hx,hy;
GET_FLT128_MSW64(hx,x); GET_FLT128_MSW64(hx,x);
GET_FLT128_MSW64(hy,y); GET_FLT128_MSW64(hy,y);
SET_FLT128_MSW64(x,(hx&0x7fffffffffffffffULL)|(hy&0x8000000000000000ULL)); SET_FLT128_MSW64(x,(hx&0x7fffffffffffffffULL)
return x; |(hy&0x8000000000000000ULL));
return x;
} }
...@@ -9,11 +9,11 @@ ...@@ -9,11 +9,11 @@
* ==================================================== * ====================================================
*/ */
/* Changes for 128-bit __float128 are /* Changes for 128-bit long double are
Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov> Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
and are incorporated herein by permission of the author. The author and are incorporated herein by permission of the author. The author
reserves the right to distribute this material elsewhere under different reserves the right to distribute this material elsewhere under different
copying permissions. These modifications are distributed here under copying permissions. These modifications are distributed here under
the following terms: the following terms:
This library is free software; you can redistribute it and/or This library is free software; you can redistribute it and/or
...@@ -27,34 +27,34 @@ ...@@ -27,34 +27,34 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software License along with this library; if not, see
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ <http://www.gnu.org/licenses/>. */
/* coshq(x) /* coshq(x)
* Method : * Method :
* mathematically coshq(x) if defined to be (exp(x)+exp(-x))/2 * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
* 1. Replace x by |x| (coshq(x) = coshq(-x)). * 1. Replace x by |x| (coshl(x) = coshl(-x)).
* 2. * 2.
* [ exp(x) - 1 ]^2 * [ exp(x) - 1 ]^2
* 0 <= x <= ln2/2 : coshq(x) := 1 + ------------------- * 0 <= x <= ln2/2 : coshl(x) := 1 + -------------------
* 2*exp(x) * 2*exp(x)
* *
* exp(x) + 1/exp(x) * exp(x) + 1/exp(x)
* ln2/2 <= x <= 22 : coshq(x) := ------------------- * ln2/2 <= x <= 22 : coshl(x) := -------------------
* 2 * 2
* 22 <= x <= lnovft : coshq(x) := expq(x)/2 * 22 <= x <= lnovft : coshl(x) := expq(x)/2
* lnovft <= x <= ln2ovft: coshq(x) := expq(x/2)/2 * expq(x/2) * lnovft <= x <= ln2ovft: coshl(x) := expq(x/2)/2 * expq(x/2)
* ln2ovft < x : coshq(x) := huge*huge (overflow) * ln2ovft < x : coshl(x) := huge*huge (overflow)
* *
* Special cases: * Special cases:
* coshq(x) is |x| if x is +INF, -INF, or NaN. * coshl(x) is |x| if x is +INF, -INF, or NaN.
* only coshq(0)=1 is exact for finite x. * only coshl(0)=1 is exact for finite x.
*/ */
#include "quadmath-imp.h" #include "quadmath-imp.h"
static const __float128 one = 1.0Q, half = 0.5Q, huge = 1.0e4900Q, static const __float128 one = 1.0, half = 0.5, huge = 1.0e4900Q,
ovf_thresh = 1.1357216553474703894801348310092223067821E4Q; ovf_thresh = 1.1357216553474703894801348310092223067821E4Q;
__float128 __float128
coshq (__float128 x) coshq (__float128 x)
...@@ -73,7 +73,7 @@ coshq (__float128 x) ...@@ -73,7 +73,7 @@ coshq (__float128 x)
if (ex >= 0x7fff0000) if (ex >= 0x7fff0000)
return x * x; return x * x;
/* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expq(|x|)) */ /* |x| in [0,0.5*ln2], return 1+expm1q(|x|)^2/(2*expq(|x|)) */
if (ex < 0x3ffd62e4) /* 0.3465728759765625 */ if (ex < 0x3ffd62e4) /* 0.3465728759765625 */
{ {
if (ex < 0x3fb80000) /* |x| < 2^-116 */ if (ex < 0x3fb80000) /* |x| < 2^-116 */
......
/* cosq.c -- __float128 version of s_cos.c. /* s_cosl.c -- long double version of s_cos.c.
* Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz. * Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz.
*/ */
...@@ -17,8 +17,8 @@ ...@@ -17,8 +17,8 @@
* Return cosine function of x. * Return cosine function of x.
* *
* kernel function: * kernel function:
* __quadmath_kernel_sinq ... sine function on [-pi/4,pi/4] * __quadmath_kernel_sinq ... sine function on [-pi/4,pi/4]
* __quadmath_kernel_cosq ... cosine function on [-pi/4,pi/4] * __quadmath_kernel_cosq ... cosine function on [-pi/4,pi/4]
* __quadmath_rem_pio2q ... argument reduction routine * __quadmath_rem_pio2q ... argument reduction routine
* *
* Method. * Method.
...@@ -46,10 +46,9 @@ ...@@ -46,10 +46,9 @@
#include "quadmath-imp.h" #include "quadmath-imp.h"
__float128 __float128 cosq(__float128 x)
cosq (__float128 x)
{ {
__float128 y[2],z=0.0Q; __float128 y[2],z=0;
int64_t n, ix; int64_t n, ix;
/* High word of x. */ /* High word of x. */
...@@ -64,6 +63,8 @@ cosq (__float128 x) ...@@ -64,6 +63,8 @@ cosq (__float128 x)
else if (ix>=0x7fff000000000000LL) { else if (ix>=0x7fff000000000000LL) {
if (ix == 0x7fff000000000000LL) { if (ix == 0x7fff000000000000LL) {
GET_FLT128_LSW64(n,x); GET_FLT128_LSW64(n,x);
if (n == 0)
errno = EDOM;
} }
return x-x; return x-x;
} }
......
/* Quad-precision floating point cosine on <-pi/4,pi/4>. /* Quad-precision floating point cosine on <-pi/4,pi/4>.
Copyright (C) 1999 Free Software Foundation, Inc. Copyright (C) 1999-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Jakub Jelinek <jj@ultra.linux.cz> Contributed by Jakub Jelinek <jj@ultra.linux.cz>
...@@ -14,9 +14,8 @@ ...@@ -14,9 +14,8 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
...@@ -70,14 +69,14 @@ static const __float128 c[] = { ...@@ -70,14 +69,14 @@ static const __float128 c[] = {
-2.50521016467996193495359189395805639E-08Q, /* bfe5ae644ee90c47dc71839de75b2787 */ -2.50521016467996193495359189395805639E-08Q, /* bfe5ae644ee90c47dc71839de75b2787 */
}; };
#define SINCOSQ_COS_HI 0 #define SINCOSL_COS_HI 0
#define SINCOSQ_COS_LO 1 #define SINCOSL_COS_LO 1
#define SINCOSQ_SIN_HI 2 #define SINCOSL_SIN_HI 2
#define SINCOSQ_SIN_LO 3 #define SINCOSL_SIN_LO 3
extern const __float128 __sincosq_table[]; extern const __float128 __sincosq_table[];
__float128 __float128
__quadmath_kernel_cosq (__float128 x, __float128 y) __quadmath_kernel_cosq(__float128 x, __float128 y)
{ {
__float128 h, l, z, sin_l, cos_l_m1; __float128 h, l, z, sin_l, cos_l_m1;
int64_t ix; int64_t ix;
...@@ -98,7 +97,7 @@ __quadmath_kernel_cosq (__float128 x, __float128 y) ...@@ -98,7 +97,7 @@ __quadmath_kernel_cosq (__float128 x, __float128 y)
else else
{ {
/* So that we don't have to use too large polynomial, we find /* So that we don't have to use too large polynomial, we find
l and h such that x = l + h, where fabsl(l) <= 1.0/256 with 83 l and h such that x = l + h, where fabsq(l) <= 1.0/256 with 83
possible values for h. We look up cosq(h) and sinq(h) in possible values for h. We look up cosq(h) and sinq(h) in
pre-computed tables, compute cosq(l) and sinq(l) using a pre-computed tables, compute cosq(l) and sinq(l) using a
Chebyshev polynomial of degree 10(11) and compute Chebyshev polynomial of degree 10(11) and compute
...@@ -106,10 +105,10 @@ __quadmath_kernel_cosq (__float128 x, __float128 y) ...@@ -106,10 +105,10 @@ __quadmath_kernel_cosq (__float128 x, __float128 y)
index = 0x3ffe - (tix >> 16); index = 0x3ffe - (tix >> 16);
hix = (tix + (0x200 << index)) & (0xfffffc00 << index); hix = (tix + (0x200 << index)) & (0xfffffc00 << index);
if (signbitq (x)) if (signbitq (x))
{ {
x = -x; x = -x;
y = -y; y = -y;
} }
switch (index) switch (index)
{ {
case 0: index = ((45 << 10) + hix - 0x3ffe0000) >> 8; break; case 0: index = ((45 << 10) + hix - 0x3ffe0000) >> 8; break;
...@@ -123,9 +122,9 @@ __quadmath_kernel_cosq (__float128 x, __float128 y) ...@@ -123,9 +122,9 @@ __quadmath_kernel_cosq (__float128 x, __float128 y)
z = l * l; z = l * l;
sin_l = l*(ONE+z*(SSIN1+z*(SSIN2+z*(SSIN3+z*(SSIN4+z*SSIN5))))); sin_l = l*(ONE+z*(SSIN1+z*(SSIN2+z*(SSIN3+z*(SSIN4+z*SSIN5)))));
cos_l_m1 = z*(SCOS1+z*(SCOS2+z*(SCOS3+z*(SCOS4+z*SCOS5)))); cos_l_m1 = z*(SCOS1+z*(SCOS2+z*(SCOS3+z*(SCOS4+z*SCOS5))));
return __sincosq_table [index + SINCOSQ_COS_HI] return __sincosq_table [index + SINCOSL_COS_HI]
+ (__sincosq_table [index + SINCOSQ_COS_LO] + (__sincosq_table [index + SINCOSL_COS_LO]
- (__sincosq_table [index + SINCOSQ_SIN_HI] * sin_l - (__sincosq_table [index + SINCOSL_SIN_HI] * sin_l
- __sincosq_table [index + SINCOSQ_COS_HI] * cos_l_m1)); - __sincosq_table [index + SINCOSL_COS_HI] * cos_l_m1));
} }
} }
/* Compute projection of complex __float128 value to Riemann sphere. /* Compute projection of complex float type value to Riemann sphere.
Copyright (C) 1997, 1999, 2010 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -14,24 +14,20 @@ ...@@ -14,24 +14,20 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
__complex128 __complex128
cprojq (__complex128 x) cprojq (__complex128 x)
{ {
if (isnanq (__real__ x) && isnanq (__imag__ x)) if (isinfq (__real__ x) || isinfq (__imag__ x))
return x;
else if (!finiteq (__real__ x) || !finiteq (__imag__ x))
{ {
__complex128 res; __complex128 res;
__real__ res = __builtin_inf (); __real__ res = __builtin_inf ();
__imag__ res = copysignq (0.0, __imag__ x); __imag__ res = copysignq (0, __imag__ x);
return res; return res;
} }
......
/* Return real part of complex __float128 value. /* Return real part of complex float type.
Copyright (C) 1997, 1998 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -14,9 +14,8 @@ ...@@ -14,9 +14,8 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
......
/* Complex sine hyperbole function for complex __float128. /* Complex sine hyperbole function for float types.
Copyright (C) 1997-2012 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -19,11 +19,6 @@ ...@@ -19,11 +19,6 @@
#include "quadmath-imp.h" #include "quadmath-imp.h"
#ifdef HAVE_FENV_H
# include <fenv.h>
#endif
__complex128 __complex128
csinhq (__complex128 x) csinhq (__complex128 x)
{ {
...@@ -34,25 +29,28 @@ csinhq (__complex128 x) ...@@ -34,25 +29,28 @@ csinhq (__complex128 x)
__real__ x = fabsq (__real__ x); __real__ x = fabsq (__real__ x);
if (__builtin_expect (rcls >= QUADFP_ZERO, 1)) if (__glibc_likely (rcls >= QUADFP_ZERO))
{ {
/* Real part is finite. */ /* Real part is finite. */
if (__builtin_expect (icls >= QUADFP_ZERO, 1)) if (__glibc_likely (icls >= QUADFP_ZERO))
{ {
/* Imaginary part is finite. */ /* Imaginary part is finite. */
const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q); const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q);
__float128 sinix, cosix; __float128 sinix, cosix;
if (__builtin_expect (icls != QUADFP_SUBNORMAL, 1)) if (__glibc_likely (fabsq (__imag__ x) > FLT128_MIN))
{ {
sincosq (__imag__ x, &sinix, &cosix); sincosq (__imag__ x, &sinix, &cosix);
} }
else else
{ {
sinix = __imag__ x; sinix = __imag__ x;
cosix = 1.0Q; cosix = 1;
} }
if (negate)
cosix = -cosix;
if (fabsq (__real__ x) > t) if (fabsq (__real__ x) > t)
{ {
__float128 exp_t = expq (t); __float128 exp_t = expq (t);
...@@ -60,8 +58,8 @@ csinhq (__complex128 x) ...@@ -60,8 +58,8 @@ csinhq (__complex128 x)
if (signbitq (__real__ x)) if (signbitq (__real__ x))
cosix = -cosix; cosix = -cosix;
rx -= t; rx -= t;
sinix *= exp_t / 2.0Q; sinix *= exp_t / 2;
cosix *= exp_t / 2.0Q; cosix *= exp_t / 2;
if (rx > t) if (rx > t)
{ {
rx -= t; rx -= t;
...@@ -87,49 +85,41 @@ csinhq (__complex128 x) ...@@ -87,49 +85,41 @@ csinhq (__complex128 x)
__imag__ retval = coshq (__real__ x) * sinix; __imag__ retval = coshq (__real__ x) * sinix;
} }
if (negate) math_check_force_underflow_complex (retval);
__real__ retval = -__real__ retval;
} }
else else
{ {
if (rcls == QUADFP_ZERO) if (rcls == QUADFP_ZERO)
{ {
/* Real part is 0.0. */ /* Real part is 0.0. */
__real__ retval = copysignq (0.0Q, negate ? -1.0Q : 1.0Q); __real__ retval = copysignq (0, negate ? -1 : 1);
__imag__ retval = nanq ("") + nanq (""); __imag__ retval = __imag__ x - __imag__ x;
#ifdef HAVE_FENV_H
if (icls == QUADFP_INFINITE)
feraiseexcept (FE_INVALID);
#endif
} }
else else
{ {
__real__ retval = nanq (""); __real__ retval = nanq ("");
__imag__ retval = nanq (""); __imag__ retval = nanq ("");
#ifdef HAVE_FENV_H
feraiseexcept (FE_INVALID); feraiseexcept (FE_INVALID);
#endif
} }
} }
} }
else if (rcls == QUADFP_INFINITE) else if (rcls == QUADFP_INFINITE)
{ {
/* Real part is infinite. */ /* Real part is infinite. */
if (__builtin_expect (icls > QUADFP_ZERO, 1)) if (__glibc_likely (icls > QUADFP_ZERO))
{ {
/* Imaginary part is finite. */ /* Imaginary part is finite. */
__float128 sinix, cosix; __float128 sinix, cosix;
if (__builtin_expect (icls != QUADFP_SUBNORMAL, 1)) if (__glibc_likely (fabsq (__imag__ x) > FLT128_MIN))
{ {
sincosq (__imag__ x, &sinix, &cosix); sincosq (__imag__ x, &sinix, &cosix);
} }
else else
{ {
sinix = __imag__ x; sinix = __imag__ x;
cosix = 1.0; cosix = 1;
} }
__real__ retval = copysignq (HUGE_VALQ, cosix); __real__ retval = copysignq (HUGE_VALQ, cosix);
...@@ -146,20 +136,14 @@ csinhq (__complex128 x) ...@@ -146,20 +136,14 @@ csinhq (__complex128 x)
} }
else else
{ {
/* The addition raises the invalid exception. */
__real__ retval = HUGE_VALQ; __real__ retval = HUGE_VALQ;
__imag__ retval = nanq ("") + nanq (""); __imag__ retval = __imag__ x - __imag__ x;
#ifdef HAVE_FENV_H
if (icls == QUADFP_INFINITE)
feraiseexcept (FE_INVALID);
#endif
} }
} }
else else
{ {
__real__ retval = nanq (""); __real__ retval = nanq ("");
__imag__ retval = __imag__ x == 0.0Q ? __imag__ x : nanq (""); __imag__ retval = __imag__ x == 0 ? __imag__ x : nanq ("");
} }
return retval; return retval;
......
/* Complex sine function for complex __float128. /* Complex sine function for float types.
Copyright (C) 1997-2012 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -19,11 +19,6 @@ ...@@ -19,11 +19,6 @@
#include "quadmath-imp.h" #include "quadmath-imp.h"
#ifdef HAVE_FENV_H
# include <fenv.h>
#endif
__complex128 __complex128
csinq (__complex128 x) csinq (__complex128 x)
{ {
...@@ -34,25 +29,28 @@ csinq (__complex128 x) ...@@ -34,25 +29,28 @@ csinq (__complex128 x)
__real__ x = fabsq (__real__ x); __real__ x = fabsq (__real__ x);
if (__builtin_expect (icls >= QUADFP_ZERO, 1)) if (__glibc_likely (icls >= QUADFP_ZERO))
{ {
/* Imaginary part is finite. */ /* Imaginary part is finite. */
if (__builtin_expect (rcls >= QUADFP_ZERO, 1)) if (__glibc_likely (rcls >= QUADFP_ZERO))
{ {
/* Real part is finite. */ /* Real part is finite. */
const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q); const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q);
__float128 sinix, cosix; __float128 sinix, cosix;
if (__builtin_expect (rcls != QUADFP_SUBNORMAL, 1)) if (__glibc_likely (__real__ x > FLT128_MIN))
{ {
sincosq (__real__ x, &sinix, &cosix); sincosq (__real__ x, &sinix, &cosix);
} }
else else
{ {
sinix = __real__ x; sinix = __real__ x;
cosix = 1.0Q; cosix = 1;
} }
if (negate)
sinix = -sinix;
if (fabsq (__imag__ x) > t) if (fabsq (__imag__ x) > t)
{ {
__float128 exp_t = expq (t); __float128 exp_t = expq (t);
...@@ -60,8 +58,8 @@ csinq (__complex128 x) ...@@ -60,8 +58,8 @@ csinq (__complex128 x)
if (signbitq (__imag__ x)) if (signbitq (__imag__ x))
cosix = -cosix; cosix = -cosix;
ix -= t; ix -= t;
sinix *= exp_t / 2.0Q; sinix *= exp_t / 2;
cosix *= exp_t / 2.0Q; cosix *= exp_t / 2;
if (ix > t) if (ix > t)
{ {
ix -= t; ix -= t;
...@@ -87,30 +85,22 @@ csinq (__complex128 x) ...@@ -87,30 +85,22 @@ csinq (__complex128 x)
__imag__ retval = sinhq (__imag__ x) * cosix; __imag__ retval = sinhq (__imag__ x) * cosix;
} }
if (negate) math_check_force_underflow_complex (retval);
__real__ retval = -__real__ retval;
} }
else else
{ {
if (icls == QUADFP_ZERO) if (icls == QUADFP_ZERO)
{ {
/* Imaginary part is 0.0. */ /* Imaginary part is 0.0. */
__real__ retval = nanq (""); __real__ retval = __real__ x - __real__ x;
__imag__ retval = __imag__ x; __imag__ retval = __imag__ x;
#ifdef HAVE_FENV_H
if (rcls == QUADFP_INFINITE)
feraiseexcept (FE_INVALID);
#endif
} }
else else
{ {
__real__ retval = nanq (""); __real__ retval = nanq ("");
__imag__ retval = nanq (""); __imag__ retval = nanq ("");
#ifdef HAVE_FENV_H
feraiseexcept (FE_INVALID); feraiseexcept (FE_INVALID);
#endif
} }
} }
} }
...@@ -120,7 +110,7 @@ csinq (__complex128 x) ...@@ -120,7 +110,7 @@ csinq (__complex128 x)
if (rcls == QUADFP_ZERO) if (rcls == QUADFP_ZERO)
{ {
/* Real part is 0.0. */ /* Real part is 0.0. */
__real__ retval = copysignq (0.0Q, negate ? -1.0Q : 1.0Q); __real__ retval = copysignq (0, negate ? -1 : 1);
__imag__ retval = __imag__ x; __imag__ retval = __imag__ x;
} }
else if (rcls > QUADFP_ZERO) else if (rcls > QUADFP_ZERO)
...@@ -128,14 +118,14 @@ csinq (__complex128 x) ...@@ -128,14 +118,14 @@ csinq (__complex128 x)
/* Real part is finite. */ /* Real part is finite. */
__float128 sinix, cosix; __float128 sinix, cosix;
if (__builtin_expect (rcls != QUADFP_SUBNORMAL, 1)) if (__glibc_likely (__real__ x > FLT128_MIN))
{ {
sincosq (__real__ x, &sinix, &cosix); sincosq (__real__ x, &sinix, &cosix);
} }
else else
{ {
sinix = __real__ x; sinix = __real__ x;
cosix = 1.0; cosix = 1;
} }
__real__ retval = copysignq (HUGE_VALQ, sinix); __real__ retval = copysignq (HUGE_VALQ, sinix);
...@@ -148,20 +138,14 @@ csinq (__complex128 x) ...@@ -148,20 +138,14 @@ csinq (__complex128 x)
} }
else else
{ {
/* The addition raises the invalid exception. */ __real__ retval = __real__ x - __real__ x;
__real__ retval = nanq ("");
__imag__ retval = HUGE_VALQ; __imag__ retval = HUGE_VALQ;
#ifdef HAVE_FENV_H
if (rcls == QUADFP_INFINITE)
feraiseexcept (FE_INVALID);
#endif
} }
} }
else else
{ {
if (rcls == QUADFP_ZERO) if (rcls == QUADFP_ZERO)
__real__ retval = copysignq (0.0Q, negate ? -1.0Q : 1.0Q); __real__ retval = copysignq (0, negate ? -1 : 1);
else else
__real__ retval = nanq (""); __real__ retval = nanq ("");
__imag__ retval = nanq (""); __imag__ retval = nanq ("");
......
/* Complex square root of __float128 value. /* Complex square root of a float type.
Copyright (C) 1997-2012 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Based on an algorithm by Stephen L. Moshier <moshier@world.std.com>. Based on an algorithm by Stephen L. Moshier <moshier@world.std.com>.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -20,11 +20,6 @@ ...@@ -20,11 +20,6 @@
#include "quadmath-imp.h" #include "quadmath-imp.h"
#ifdef HAVE_FENV_H
# include <fenv.h>
#endif
__complex128 __complex128
csqrtq (__complex128 x) csqrtq (__complex128 x)
{ {
...@@ -32,7 +27,7 @@ csqrtq (__complex128 x) ...@@ -32,7 +27,7 @@ csqrtq (__complex128 x)
int rcls = fpclassifyq (__real__ x); int rcls = fpclassifyq (__real__ x);
int icls = fpclassifyq (__imag__ x); int icls = fpclassifyq (__imag__ x);
if (__builtin_expect (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE, 0)) if (__glibc_unlikely (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE))
{ {
if (icls == QUADFP_INFINITE) if (icls == QUADFP_INFINITE)
{ {
...@@ -41,7 +36,7 @@ csqrtq (__complex128 x) ...@@ -41,7 +36,7 @@ csqrtq (__complex128 x)
} }
else if (rcls == QUADFP_INFINITE) else if (rcls == QUADFP_INFINITE)
{ {
if (__real__ x < 0.0Q) if (__real__ x < 0)
{ {
__real__ res = icls == QUADFP_NAN ? nanq ("") : 0; __real__ res = icls == QUADFP_NAN ? nanq ("") : 0;
__imag__ res = copysignq (HUGE_VALQ, __imag__ x); __imag__ res = copysignq (HUGE_VALQ, __imag__ x);
...@@ -50,7 +45,7 @@ csqrtq (__complex128 x) ...@@ -50,7 +45,7 @@ csqrtq (__complex128 x)
{ {
__real__ res = __real__ x; __real__ res = __real__ x;
__imag__ res = (icls == QUADFP_NAN __imag__ res = (icls == QUADFP_NAN
? nanq ("") : copysignq (0.0Q, __imag__ x)); ? nanq ("") : copysignq (0, __imag__ x));
} }
} }
else else
...@@ -61,27 +56,26 @@ csqrtq (__complex128 x) ...@@ -61,27 +56,26 @@ csqrtq (__complex128 x)
} }
else else
{ {
if (__builtin_expect (icls == QUADFP_ZERO, 0)) if (__glibc_unlikely (icls == QUADFP_ZERO))
{ {
if (__real__ x < 0.0Q) if (__real__ x < 0)
{ {
__real__ res = 0.0Q; __real__ res = 0;
__imag__ res = copysignq (sqrtq (-__real__ x), __imag__ res = copysignq (sqrtq (-__real__ x), __imag__ x);
__imag__ x);
} }
else else
{ {
__real__ res = fabsq (sqrtq (__real__ x)); __real__ res = fabsq (sqrtq (__real__ x));
__imag__ res = copysignq (0.0Q, __imag__ x); __imag__ res = copysignq (0, __imag__ x);
} }
} }
else if (__builtin_expect (rcls == QUADFP_ZERO, 0)) else if (__glibc_unlikely (rcls == QUADFP_ZERO))
{ {
__float128 r; __float128 r;
if (fabsq (__imag__ x) >= 2.0Q * FLT128_MIN) if (fabsq (__imag__ x) >= 2 * FLT128_MIN)
r = sqrtq (0.5Q * fabsq (__imag__ x)); r = sqrtq (0.5Q * fabsq (__imag__ x));
else else
r = 0.5Q * sqrtq (2.0Q * fabsq (__imag__ x)); r = 0.5Q * sqrtq (2 * fabsq (__imag__ x));
__real__ res = r; __real__ res = r;
__imag__ res = copysignq (r, __imag__ x); __imag__ res = copysignq (r, __imag__ x);
...@@ -91,25 +85,25 @@ csqrtq (__complex128 x) ...@@ -91,25 +85,25 @@ csqrtq (__complex128 x)
__float128 d, r, s; __float128 d, r, s;
int scale = 0; int scale = 0;
if (fabsq (__real__ x) > FLT128_MAX / 4.0Q) if (fabsq (__real__ x) > FLT128_MAX / 4)
{ {
scale = 1; scale = 1;
__real__ x = scalbnq (__real__ x, -2 * scale); __real__ x = scalbnq (__real__ x, -2 * scale);
__imag__ x = scalbnq (__imag__ x, -2 * scale); __imag__ x = scalbnq (__imag__ x, -2 * scale);
} }
else if (fabsq (__imag__ x) > FLT128_MAX / 4.0Q) else if (fabsq (__imag__ x) > FLT128_MAX / 4)
{ {
scale = 1; scale = 1;
if (fabsq (__real__ x) >= 4.0Q * FLT128_MIN) if (fabsq (__real__ x) >= 4 * FLT128_MIN)
__real__ x = scalbnq (__real__ x, -2 * scale); __real__ x = scalbnq (__real__ x, -2 * scale);
else else
__real__ x = 0.0Q; __real__ x = 0;
__imag__ x = scalbnq (__imag__ x, -2 * scale); __imag__ x = scalbnq (__imag__ x, -2 * scale);
} }
else if (fabsq (__real__ x) < FLT128_MIN else if (fabsq (__real__ x) < 2 * FLT128_MIN
&& fabsq (__imag__ x) < FLT128_MIN) && fabsq (__imag__ x) < 2 * FLT128_MIN)
{ {
scale = -(FLT128_MANT_DIG / 2); scale = -((FLT128_MANT_DIG + 1) / 2);
__real__ x = scalbnq (__real__ x, -2 * scale); __real__ x = scalbnq (__real__ x, -2 * scale);
__imag__ x = scalbnq (__imag__ x, -2 * scale); __imag__ x = scalbnq (__imag__ x, -2 * scale);
} }
...@@ -120,12 +114,28 @@ csqrtq (__complex128 x) ...@@ -120,12 +114,28 @@ csqrtq (__complex128 x)
if (__real__ x > 0) if (__real__ x > 0)
{ {
r = sqrtq (0.5Q * (d + __real__ x)); r = sqrtq (0.5Q * (d + __real__ x));
s = 0.5Q * (__imag__ x / r); if (scale == 1 && fabsq (__imag__ x) < 1)
{
/* Avoid possible intermediate underflow. */
s = __imag__ x / r;
r = scalbnq (r, scale);
scale = 0;
}
else
s = 0.5Q * (__imag__ x / r);
} }
else else
{ {
s = sqrtq (0.5Q * (d - __real__ x)); s = sqrtq (0.5Q * (d - __real__ x));
r = fabsq (0.5Q * (__imag__ x / s)); if (scale == 1 && fabsq (__imag__ x) < 1)
{
/* Avoid possible intermediate underflow. */
r = fabsq (__imag__ x / s);
s = scalbnq (s, scale);
scale = 0;
}
else
r = fabsq (0.5Q * (__imag__ x / s));
} }
if (scale) if (scale)
...@@ -134,6 +144,9 @@ csqrtq (__complex128 x) ...@@ -134,6 +144,9 @@ csqrtq (__complex128 x)
s = scalbnq (s, scale); s = scalbnq (s, scale);
} }
math_check_force_underflow (r);
math_check_force_underflow (s);
__real__ res = r; __real__ res = r;
__imag__ res = copysignq (s, __imag__ x); __imag__ res = copysignq (s, __imag__ x);
} }
......
/* Complex hyperbole tangent for __float128. /* Complex hyperbolic tangent for float types.
Copyright (C) 1997-2012 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -19,36 +19,39 @@ ...@@ -19,36 +19,39 @@
#include "quadmath-imp.h" #include "quadmath-imp.h"
#ifdef HAVE_FENV_H
# include <fenv.h>
#endif
__complex128 __complex128
ctanhq (__complex128 x) ctanhq (__complex128 x)
{ {
__complex128 res; __complex128 res;
if (__builtin_expect (!finiteq (__real__ x) || !finiteq (__imag__ x), 0)) if (__glibc_unlikely (!finiteq (__real__ x) || !finiteq (__imag__ x)))
{ {
if (__quadmath_isinf_nsq (__real__ x)) if (isinfq (__real__ x))
{ {
__real__ res = copysignq (1.0Q, __real__ x); __real__ res = copysignq (1, __real__ x);
__imag__ res = copysignq (0.0Q, __imag__ x); if (finiteq (__imag__ x) && fabsq (__imag__ x) > 1)
{
__float128 sinix, cosix;
sincosq (__imag__ x, &sinix, &cosix);
__imag__ res = copysignq (0, sinix * cosix);
}
else
__imag__ res = copysignq (0, __imag__ x);
} }
else if (__imag__ x == 0.0Q) else if (__imag__ x == 0)
{ {
res = x; res = x;
} }
else else
{ {
__real__ res = nanq (""); if (__real__ x == 0)
__real__ res = __real__ x;
else
__real__ res = nanq ("");
__imag__ res = nanq (""); __imag__ res = nanq ("");
#ifdef HAVE_FENV_H if (isinfq (__imag__ x))
if (__quadmath_isinf_nsq (__imag__ x))
feraiseexcept (FE_INVALID); feraiseexcept (FE_INVALID);
#endif
} }
} }
else else
...@@ -56,19 +59,18 @@ ctanhq (__complex128 x) ...@@ -56,19 +59,18 @@ ctanhq (__complex128 x)
__float128 sinix, cosix; __float128 sinix, cosix;
__float128 den; __float128 den;
const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q / 2); const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q / 2);
int icls = fpclassifyq (__imag__ x);
/* tanh(x+iy) = (sinh(2x) + i*sin(2y))/(cosh(2x) + cos(2y)) /* tanh(x+iy) = (sinh(2x) + i*sin(2y))/(cosh(2x) + cos(2y))
= (sinh(x)*cosh(x) + i*sin(y)*cos(y))/(sinh(x)^2 + cos(y)^2). */ = (sinh(x)*cosh(x) + i*sin(y)*cos(y))/(sinh(x)^2 + cos(y)^2). */
if (__builtin_expect (icls != QUADFP_SUBNORMAL, 1)) if (__glibc_likely (fabsq (__imag__ x) > FLT128_MIN))
{ {
sincosq (__imag__ x, &sinix, &cosix); sincosq (__imag__ x, &sinix, &cosix);
} }
else else
{ {
sinix = __imag__ x; sinix = __imag__ x;
cosix = 1.0Q; cosix = 1;
} }
if (fabsq (__real__ x) > t) if (fabsq (__real__ x) > t)
...@@ -79,7 +81,7 @@ ctanhq (__complex128 x) ...@@ -79,7 +81,7 @@ ctanhq (__complex128 x)
sin(y)*cos(y)/sinh(x)^2 = 4*sin(y)*cos(y)/exp(2x). */ sin(y)*cos(y)/sinh(x)^2 = 4*sin(y)*cos(y)/exp(2x). */
__float128 exp_2t = expq (2 * t); __float128 exp_2t = expq (2 * t);
__real__ res = copysignq (1.0, __real__ x); __real__ res = copysignq (1, __real__ x);
__imag__ res = 4 * sinix * cosix; __imag__ res = 4 * sinix * cosix;
__real__ x = fabsq (__real__ x); __real__ x = fabsq (__real__ x);
__real__ x -= t; __real__ x -= t;
...@@ -104,7 +106,7 @@ ctanhq (__complex128 x) ...@@ -104,7 +106,7 @@ ctanhq (__complex128 x)
else else
{ {
sinhrx = __real__ x; sinhrx = __real__ x;
coshrx = 1.0Q; coshrx = 1;
} }
if (fabsq (sinhrx) > fabsq (cosix) * FLT128_EPSILON) if (fabsq (sinhrx) > fabsq (cosix) * FLT128_EPSILON)
...@@ -114,6 +116,7 @@ ctanhq (__complex128 x) ...@@ -114,6 +116,7 @@ ctanhq (__complex128 x)
__real__ res = sinhrx * coshrx / den; __real__ res = sinhrx * coshrx / den;
__imag__ res = sinix * cosix / den; __imag__ res = sinix * cosix / den;
} }
math_check_force_underflow_complex (res);
} }
return res; return res;
......
/* Complex tangent function for complex __float128. /* Complex tangent function for a complex float type.
Copyright (C) 1997-2012 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -19,56 +19,58 @@ ...@@ -19,56 +19,58 @@
#include "quadmath-imp.h" #include "quadmath-imp.h"
#ifdef HAVE_FENV_H
# include <fenv.h>
#endif
__complex128 __complex128
ctanq (__complex128 x) ctanq (__complex128 x)
{ {
__complex128 res; __complex128 res;
if (__builtin_expect (!finiteq (__real__ x) || !finiteq (__imag__ x), 0)) if (__glibc_unlikely (!finiteq (__real__ x) || !finiteq (__imag__ x)))
{ {
if (__quadmath_isinf_nsq (__imag__ x)) if (isinfq (__imag__ x))
{ {
__real__ res = copysignq (0.0Q, __real__ x); if (finiteq (__real__ x) && fabsq (__real__ x) > 1)
__imag__ res = copysignq (1.0Q, __imag__ x); {
__float128 sinrx, cosrx;
sincosq (__real__ x, &sinrx, &cosrx);
__real__ res = copysignq (0, sinrx * cosrx);
}
else
__real__ res = copysignq (0, __real__ x);
__imag__ res = copysignq (1, __imag__ x);
} }
else if (__real__ x == 0.0Q) else if (__real__ x == 0)
{ {
res = x; res = x;
} }
else else
{ {
__real__ res = nanq (""); __real__ res = nanq ("");
__imag__ res = nanq (""); if (__imag__ x == 0)
__imag__ res = __imag__ x;
else
__imag__ res = nanq ("");
#ifdef HAVE_FENV_H if (isinfq (__real__ x))
if (__quadmath_isinf_nsq (__real__ x))
feraiseexcept (FE_INVALID); feraiseexcept (FE_INVALID);
#endif
} }
} }
else else
{ {
__float128 sinrx, cosrx; __float128 sinrx, cosrx;
__float128 den; __float128 den;
const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q / 2.0Q); const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q / 2);
int rcls = fpclassifyq (__real__ x);
/* tan(x+iy) = (sin(2x) + i*sinh(2y))/(cos(2x) + cosh(2y)) /* tan(x+iy) = (sin(2x) + i*sinh(2y))/(cos(2x) + cosh(2y))
= (sin(x)*cos(x) + i*sinh(y)*cosh(y)/(cos(x)^2 + sinh(y)^2). */ = (sin(x)*cos(x) + i*sinh(y)*cosh(y)/(cos(x)^2 + sinh(y)^2). */
if (__builtin_expect (rcls != QUADFP_SUBNORMAL, 1)) if (__glibc_likely (fabsq (__real__ x) > FLT128_MIN))
{ {
sincosq (__real__ x, &sinrx, &cosrx); sincosq (__real__ x, &sinrx, &cosrx);
} }
else else
{ {
sinrx = __real__ x; sinrx = __real__ x;
cosrx = 1.0Q; cosrx = 1;
} }
if (fabsq (__imag__ x) > t) if (fabsq (__imag__ x) > t)
...@@ -79,7 +81,7 @@ ctanq (__complex128 x) ...@@ -79,7 +81,7 @@ ctanq (__complex128 x)
sin(x)*cos(x)/sinh(y)^2 = 4*sin(x)*cos(x)/exp(2y). */ sin(x)*cos(x)/sinh(y)^2 = 4*sin(x)*cos(x)/exp(2y). */
__float128 exp_2t = expq (2 * t); __float128 exp_2t = expq (2 * t);
__imag__ res = copysignq (1.0Q, __imag__ x); __imag__ res = copysignq (1, __imag__ x);
__real__ res = 4 * sinrx * cosrx; __real__ res = 4 * sinrx * cosrx;
__imag__ x = fabsq (__imag__ x); __imag__ x = fabsq (__imag__ x);
__imag__ x -= t; __imag__ x -= t;
...@@ -104,7 +106,7 @@ ctanq (__complex128 x) ...@@ -104,7 +106,7 @@ ctanq (__complex128 x)
else else
{ {
sinhix = __imag__ x; sinhix = __imag__ x;
coshix = 1.0Q; coshix = 1;
} }
if (fabsq (sinhix) > fabsq (cosrx) * FLT128_EPSILON) if (fabsq (sinhix) > fabsq (cosrx) * FLT128_EPSILON)
...@@ -114,6 +116,7 @@ ctanq (__complex128 x) ...@@ -114,6 +116,7 @@ ctanq (__complex128 x)
__real__ res = sinrx * cosrx / den; __real__ res = sinrx * cosrx / den;
__imag__ res = sinhix * coshix / den; __imag__ res = sinhix * coshix / den;
} }
math_check_force_underflow_complex (res);
} }
return res; return res;
......
...@@ -27,11 +27,11 @@ ...@@ -27,11 +27,11 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software License along with this library; if not, see
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ <http://www.gnu.org/licenses/>. */
/* __float128 erfq(__float128 x) /* double erf(double x)
* __float128 erfcq(__float128 x) * double erfc(double x)
* x * x
* 2 |\ * 2 |\
* erf(x) = --------- | exp(-t*t)dt * erf(x) = --------- | exp(-t*t)dt
...@@ -96,14 +96,8 @@ ...@@ -96,14 +96,8 @@
* erfc/erf(NaN) is NaN * erfc/erf(NaN) is NaN
*/ */
#include <errno.h>
#include "quadmath-imp.h" #include "quadmath-imp.h"
__float128 erfcq (__float128);
/* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */ /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
static __float128 static __float128
...@@ -143,8 +137,8 @@ deval (__float128 x, const __float128 *p, int n) ...@@ -143,8 +137,8 @@ deval (__float128 x, const __float128 *p, int n)
static const __float128 static const __float128
tiny = 1e-4931Q, tiny = 1e-4931Q,
one = 1.0Q, one = 1,
two = 2.0Q, two = 2,
/* 2/sqrt(pi) - 1 */ /* 2/sqrt(pi) - 1 */
efx = 1.2837916709551257389615890312154517168810E-1Q; efx = 1.2837916709551257389615890312154517168810E-1Q;
...@@ -810,7 +804,7 @@ erfq (__float128 x) ...@@ -810,7 +804,7 @@ erfq (__float128 x)
__float128 __float128
erfcq (__float128 x) erfcq (__float128 x)
{ {
__float128 y = 0.0Q, z, p, r; __float128 y, z, p, r;
int32_t i, ix, sign; int32_t i, ix, sign;
ieee854_float128 u; ieee854_float128 u;
...@@ -868,7 +862,7 @@ erfcq (__float128 x) ...@@ -868,7 +862,7 @@ erfcq (__float128 x)
y += C18a; y += C18a;
break; break;
case 8: case 8:
z = x - 1.0Q; z = x - 1;
y = C19b + z * neval (z, RNr19, NRNr19) / deval (z, RDr19, NRDr19); y = C19b + z * neval (z, RNr19, NRNr19) / deval (z, RDr19, NRDr19);
y += C19a; y += C19a;
break; break;
...@@ -879,7 +873,7 @@ erfcq (__float128 x) ...@@ -879,7 +873,7 @@ erfcq (__float128 x)
break; break;
} }
if (sign & 0x80000000) if (sign & 0x80000000)
y = 2.0Q - y; y = 2 - y;
return y; return y;
} }
/* 1.25 < |x| < 107 */ /* 1.25 < |x| < 107 */
...@@ -924,7 +918,8 @@ erfcq (__float128 x) ...@@ -924,7 +918,8 @@ erfcq (__float128 x)
u.words32.w3 = 0; u.words32.w3 = 0;
u.words32.w2 &= 0xfe000000; u.words32.w2 &= 0xfe000000;
z = u.value; z = u.value;
r = expq (-z * z - 0.5625) * expq ((z - x) * (z + x) + p); r = expq (-z * z - 0.5625) *
expq ((z - x) * (z + x) + p);
if ((sign & 0x80000000) == 0) if ((sign & 0x80000000) == 0)
{ {
__float128 ret = r / x; __float128 ret = r / x;
......
/* Compute 2^x.
Copyright (C) 2012-2018 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 "quadmath-imp.h"
__float128
exp2q (__float128 x)
{
if (__glibc_likely (__builtin_isless (x, (__float128) FLT128_MAX_EXP)))
{
if (__builtin_expect (__builtin_isgreaterequal (x, (__float128) (FLT128_MIN_EXP - FLT128_MANT_DIG
- 1)), 1))
{
int intx = (int) x;
__float128 fractx = x - intx;
__float128 result;
if (fabsq (fractx) < FLT128_EPSILON / 4)
result = scalbnq (1 + fractx, intx);
else
result = scalbnq (expq (M_LN2q * fractx), intx);
math_check_force_underflow_nonneg (result);
return result;
}
else
{
/* Underflow or exact zero. */
if (isinfq (x))
return 0;
else
return FLT128_MIN * FLT128_MIN;
}
}
else
/* Infinity, NaN or overflow. */
return FLT128_MAX * x;
}
/* expm1l.c /* expm1q.c
* *
* Exponential function, minus 1 * Exponential function, minus 1
* 128-bit __float128 precision * 128-bit long double precision
* *
* *
* *
* SYNOPSIS: * SYNOPSIS:
* *
* __float128 x, y, expm1l(); * long double x, y, expm1q();
* *
* y = expm1l( x ); * y = expm1q( x );
* *
* *
* *
...@@ -48,12 +48,9 @@ ...@@ -48,12 +48,9 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software License along with this library; if not, see
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ <http://www.gnu.org/licenses/>. */
#include <errno.h>
#include "quadmath-imp.h" #include "quadmath-imp.h"
/* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x) /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
...@@ -83,7 +80,7 @@ static const __float128 ...@@ -83,7 +80,7 @@ static const __float128
C1 = 6.93145751953125E-1Q, C1 = 6.93145751953125E-1Q,
C2 = 1.428606820309417232121458176568075500134E-6Q, C2 = 1.428606820309417232121458176568075500134E-6Q,
/* ln 2^-114 */ /* ln 2^-114 */
minarg = -7.9018778583833765273564461846232128760607E1Q; minarg = -7.9018778583833765273564461846232128760607E1Q, big = 1e4932Q;
__float128 __float128
...@@ -108,7 +105,7 @@ expm1q (__float128 x) ...@@ -108,7 +105,7 @@ expm1q (__float128 x)
{ {
/* Infinity (which must be negative infinity). */ /* Infinity (which must be negative infinity). */
if (((ix & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3) == 0) if (((ix & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3) == 0)
return -1.0Q; return -1;
/* NaN. Invalid exception if signaling. */ /* NaN. Invalid exception if signaling. */
return x + x; return x + x;
} }
...@@ -119,7 +116,7 @@ expm1q (__float128 x) ...@@ -119,7 +116,7 @@ expm1q (__float128 x)
/* Minimum value. */ /* Minimum value. */
if (x < minarg) if (x < minarg)
return (4.0/HUGE_VALQ - 1.0Q); return (4.0/big - 1);
/* Avoid internal underflow when result does not underflow, while /* Avoid internal underflow when result does not underflow, while
ensuring underflow (without returning a zero of the wrong sign) ensuring underflow (without returning a zero of the wrong sign)
...@@ -156,7 +153,7 @@ expm1q (__float128 x) ...@@ -156,7 +153,7 @@ expm1q (__float128 x)
exp(x) - 1 = 2^k (qx + 1) - 1 exp(x) - 1 = 2^k (qx + 1) - 1
= 2^k qx + 2^k - 1. */ = 2^k qx + 2^k - 1. */
px = ldexpq (1.0Q, k); px = ldexpq (1, k);
x = px * qx + (px - 1.0); x = px * qx + (px - 1.0);
return x; return x;
} }
/* fabsq.c -- __float128 version of s_fabs.c. /* s_fabsl.c -- long double version of s_fabs.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz. * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/ */
...@@ -13,13 +13,20 @@ ...@@ -13,13 +13,20 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/*
* fabsq(x) returns the absolute value of x.
*/
#include "quadmath-imp.h" #include "quadmath-imp.h"
__float128 __float128 fabsq(__float128 x)
fabsq (__float128 x)
{ {
uint64_t hx; uint64_t hx;
GET_FLT128_MSW64(hx,x); GET_FLT128_MSW64(hx,x);
SET_FLT128_MSW64(x,hx&0x7fffffffffffffffLL); SET_FLT128_MSW64(x,hx&0x7fffffffffffffffLL);
return x; return x;
} }
/* Return positive difference between arguments. /* Return positive difference between arguments.
Copyright (C) 1997, 2004, 2009 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -14,29 +14,19 @@ ...@@ -14,29 +14,19 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include <errno.h>
#include "quadmath-imp.h" #include "quadmath-imp.h"
__float128 __float128
fdimq (__float128 x, __float128 y) fdimq (__float128 x, __float128 y)
{ {
int clsx = fpclassifyq (x); if (__builtin_islessequal (x, y))
int clsy = fpclassifyq (y); return 0;
if (clsx == QUADFP_NAN || clsy == QUADFP_NAN __float128 r = math_narrow_eval (x - y);
|| (y < 0 && clsx == QUADFP_INFINITE && clsy == QUADFP_INFINITE)) if (isinfq (r) && !isinfq (x) && !isinfq (y))
/* Raise invalid flag. */
return x - y;
if (x <= y)
return 0.0Q;
__float128 r = x - y;
if (isinfq (r))
errno = ERANGE; errno = ERANGE;
return r; return r;
......
/* finiteq.c -- __float128 version of s_finite.c. /* s_finitel.c -- long double version of s_finite.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz. * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/ */
...@@ -13,18 +13,21 @@ ...@@ -13,18 +13,21 @@
* ==================================================== * ====================================================
*/ */
#include "quadmath-imp.h" #if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/* /*
* finiteq(x) returns 1 is x is finite, else 0; * finiteq(x) returns 1 is x is finite, else 0;
* no branching! * no branching!
*/ */
int #include "quadmath-imp.h"
finiteq (const __float128 x)
int finiteq(__float128 x)
{ {
int64_t hx; int64_t hx;
GET_FLT128_MSW64(hx,x); GET_FLT128_MSW64(hx,x);
return (int)((uint64_t)((hx&0x7fff000000000000LL) return (int)((uint64_t)((hx&0x7fff000000000000LL)
-0x7fff000000000000LL)>>63); -0x7fff000000000000LL)>>63);
} }
/* floorq.c -- __float128 version of s_floor.c. /* s_floorl.c -- long double version of s_floor.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz. * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/ */
...@@ -13,10 +13,22 @@ ...@@ -13,10 +13,22 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/*
* floorq(x)
* Return x rounded toward -inf to integral value
* Method:
* Bit twiddling.
*/
#define NO_MATH_REDIRECT
#include "quadmath-imp.h" #include "quadmath-imp.h"
__float128 __float128 floorq(__float128 x)
floorq (__float128 x)
{ {
int64_t i0,i1,j0; int64_t i0,i1,j0;
uint64_t i,j; uint64_t i,j;
......
/* Return maximum numeric value of X and Y. /* Return maximum numeric value of X and Y.
Copyright (C) 1997 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -14,15 +14,20 @@ ...@@ -14,15 +14,20 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
__float128 __float128
fmaxq (__float128 x, __float128 y) fmaxq (__float128 x, __float128 y)
{ {
return (__builtin_isgreaterequal (x, y) || isnanq (y)) ? x : y; if (__builtin_isgreaterequal (x, y))
return x;
else if (__builtin_isless (x, y))
return y;
else if (issignalingq (x) || issignalingq (y))
return x + y;
else
return isnanq (y) ? x : y;
} }
/* Return minimum numeric value of X and Y. /* Return minimum numeric value of X and Y.
Copyright (C) 1997 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
...@@ -14,15 +14,20 @@ ...@@ -14,15 +14,20 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
__float128 __float128
fminq (__float128 x, __float128 y) fminq (__float128 x, __float128 y)
{ {
return (__builtin_islessequal (x, y) || isnanq (y)) ? x : y; if (__builtin_islessequal (x, y))
return x;
else if (__builtin_isgreater (x, y))
return y;
else if (issignalingq (x) || issignalingq (y))
return x + y;
else
return isnanq (y) ? x : y;
} }
/* fmodq.c -- __float128 version of e_fmod.c. /* e_fmodl.c -- long double version of e_fmod.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz. * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/ */
/* /*
* ==================================================== * ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (C) 1993, 2011 by Sun Microsystems, Inc. All rights reserved.
* *
* Developed at SunPro, a Sun Microsystems, Inc. business. * Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this * Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice * software is freely granted, provided that this notice
* is preserved. * is preserved.
* ==================================================== * ====================================================
*/ */
/* /*
* fmodq(x,y) * fmodq(x,y)
* Return x mod y in exact arithmetic * Return x mod y in exact arithmetic
* Method: shift and subtract * Method: shift and subtract
...@@ -25,105 +25,105 @@ static const __float128 one = 1.0, Zero[] = {0.0, -0.0,}; ...@@ -25,105 +25,105 @@ static const __float128 one = 1.0, Zero[] = {0.0, -0.0,};
__float128 __float128
fmodq (__float128 x, __float128 y) fmodq (__float128 x, __float128 y)
{ {
int64_t n,hx,hy,hz,ix,iy,sx,i; int64_t n,hx,hy,hz,ix,iy,sx,i;
uint64_t lx,ly,lz; uint64_t lx,ly,lz;
GET_FLT128_WORDS64(hx,lx,x); GET_FLT128_WORDS64(hx,lx,x);
GET_FLT128_WORDS64(hy,ly,y); GET_FLT128_WORDS64(hy,ly,y);
sx = hx&0x8000000000000000ULL; /* sign of x */ sx = hx&0x8000000000000000ULL; /* sign of x */
hx ^=sx; /* |x| */ hx ^=sx; /* |x| */
hy &= 0x7fffffffffffffffLL; /* |y| */ hy &= 0x7fffffffffffffffLL; /* |y| */
/* purge off exception values */ /* purge off exception values */
if((hy|ly)==0||(hx>=0x7fff000000000000LL)|| /* y=0,or x not finite */ if((hy|ly)==0||(hx>=0x7fff000000000000LL)|| /* y=0,or x not finite */
((hy|((ly|-ly)>>63))>0x7fff000000000000LL)) /* or y is NaN */ ((hy|((ly|-ly)>>63))>0x7fff000000000000LL)) /* or y is NaN */
return (x*y)/(x*y); return (x*y)/(x*y);
if(hx<=hy) { if(hx<=hy) {
if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */ if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
if(lx==ly) if(lx==ly)
return Zero[(uint64_t)sx>>63]; /* |x|=|y| return x*0*/ return Zero[(uint64_t)sx>>63]; /* |x|=|y| return x*0*/
} }
/* determine ix = ilogb(x) */ /* determine ix = ilogb(x) */
if(hx<0x0001000000000000LL) { /* subnormal x */ if(hx<0x0001000000000000LL) { /* subnormal x */
if(hx==0) { if(hx==0) {
for (ix = -16431, i=lx; i>0; i<<=1) ix -=1; for (ix = -16431, i=lx; i>0; i<<=1) ix -=1;
} else { } else {
for (ix = -16382, i=hx<<15; i>0; i<<=1) ix -=1; for (ix = -16382, i=hx<<15; i>0; i<<=1) ix -=1;
} }
} else ix = (hx>>48)-0x3fff; } else ix = (hx>>48)-0x3fff;
/* determine iy = ilogb(y) */ /* determine iy = ilogb(y) */
if(hy<0x0001000000000000LL) { /* subnormal y */ if(hy<0x0001000000000000LL) { /* subnormal y */
if(hy==0) { if(hy==0) {
for (iy = -16431, i=ly; i>0; i<<=1) iy -=1; for (iy = -16431, i=ly; i>0; i<<=1) iy -=1;
} else { } else {
for (iy = -16382, i=hy<<15; i>0; i<<=1) iy -=1; for (iy = -16382, i=hy<<15; i>0; i<<=1) iy -=1;
} }
} else iy = (hy>>48)-0x3fff; } else iy = (hy>>48)-0x3fff;
/* set up {hx,lx}, {hy,ly} and align y to x */ /* set up {hx,lx}, {hy,ly} and align y to x */
if(ix >= -16382) if(ix >= -16382)
hx = 0x0001000000000000LL|(0x0000ffffffffffffLL&hx); hx = 0x0001000000000000LL|(0x0000ffffffffffffLL&hx);
else { /* subnormal x, shift x to normal */ else { /* subnormal x, shift x to normal */
n = -16382-ix; n = -16382-ix;
if(n<=63) { if(n<=63) {
hx = (hx<<n)|(lx>>(64-n)); hx = (hx<<n)|(lx>>(64-n));
lx <<= n; lx <<= n;
} else { } else {
hx = lx<<(n-64); hx = lx<<(n-64);
lx = 0; lx = 0;
} }
} }
if(iy >= -16382) if(iy >= -16382)
hy = 0x0001000000000000LL|(0x0000ffffffffffffLL&hy); hy = 0x0001000000000000LL|(0x0000ffffffffffffLL&hy);
else { /* subnormal y, shift y to normal */ else { /* subnormal y, shift y to normal */
n = -16382-iy; n = -16382-iy;
if(n<=63) { if(n<=63) {
hy = (hy<<n)|(ly>>(64-n)); hy = (hy<<n)|(ly>>(64-n));
ly <<= n; ly <<= n;
} else { } else {
hy = ly<<(n-64); hy = ly<<(n-64);
ly = 0; ly = 0;
} }
} }
/* fix point fmod */ /* fix point fmod */
n = ix - iy; n = ix - iy;
while(n--) { while(n--) {
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;} if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;}
else { else {
if((hz|lz)==0) /* return sign(x)*0 */ if((hz|lz)==0) /* return sign(x)*0 */
return Zero[(uint64_t)sx>>63]; return Zero[(uint64_t)sx>>63];
hx = hz+hz+(lz>>63); lx = lz+lz; hx = hz+hz+(lz>>63); lx = lz+lz;
} }
} }
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz>=0) {hx=hz;lx=lz;} if(hz>=0) {hx=hz;lx=lz;}
/* convert back to floating value and restore the sign */ /* convert back to floating value and restore the sign */
if((hx|lx)==0) /* return sign(x)*0 */ if((hx|lx)==0) /* return sign(x)*0 */
return Zero[(uint64_t)sx>>63]; return Zero[(uint64_t)sx>>63];
while(hx<0x0001000000000000LL) { /* normalize x */ while(hx<0x0001000000000000LL) { /* normalize x */
hx = hx+hx+(lx>>63); lx = lx+lx; hx = hx+hx+(lx>>63); lx = lx+lx;
iy -= 1; iy -= 1;
} }
if(iy>= -16382) { /* normalize output */ if(iy>= -16382) { /* normalize output */
hx = ((hx-0x0001000000000000LL)|((iy+16383)<<48)); hx = ((hx-0x0001000000000000LL)|((iy+16383)<<48));
SET_FLT128_WORDS64(x,hx|sx,lx); SET_FLT128_WORDS64(x,hx|sx,lx);
} else { /* subnormal output */ } else { /* subnormal output */
n = -16382 - iy; n = -16382 - iy;
if(n<=48) { if(n<=48) {
lx = (lx>>n)|((uint64_t)hx<<(64-n)); lx = (lx>>n)|((uint64_t)hx<<(64-n));
hx >>= n; hx >>= n;
} else if (n<=63) { } else if (n<=63) {
lx = (hx<<(64-n))|(lx>>n); hx = sx; lx = (hx<<(64-n))|(lx>>n); hx = sx;
} else { } else {
lx = hx>>(n-64); hx = sx; lx = hx>>(n-64); hx = sx;
} }
SET_FLT128_WORDS64(x,hx|sx,lx); SET_FLT128_WORDS64(x,hx|sx,lx);
x *= one; /* create necessary signal */ x *= one; /* create necessary signal */
} }
return x; /* exact output */ return x; /* exact output */
} }
/* frexpq.c -- __float128 version of s_frexp.c. /* s_frexpl.c -- long double version of s_frexp.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz. * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/ */
...@@ -13,10 +13,14 @@ ...@@ -13,10 +13,14 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/* /*
* for non-zero x * for non-zero x
* x = frexpq(arg,&exp); * x = frexpq(arg,&exp);
* return a __float128 fp quantity x such that 0.5 <= |x| <1.0 * return a long double fp quantity x such that 0.5 <= |x| <1.0
* and the corresponding binary exponent "exp". That is * and the corresponding binary exponent "exp". That is
* arg = x*2^exp. * arg = x*2^exp.
* If arg is inf, 0.0, or NaN, then frexpq(arg,&exp) returns arg * If arg is inf, 0.0, or NaN, then frexpq(arg,&exp) returns arg
...@@ -28,22 +32,21 @@ ...@@ -28,22 +32,21 @@
static const __float128 static const __float128
two114 = 2.0769187434139310514121985316880384E+34Q; /* 0x4071000000000000, 0 */ two114 = 2.0769187434139310514121985316880384E+34Q; /* 0x4071000000000000, 0 */
__float128 __float128 frexpq(__float128 x, int *eptr)
frexpq (__float128 x, int *eptr)
{ {
uint64_t hx, lx, ix; uint64_t hx, lx, ix;
GET_FLT128_WORDS64(hx,lx,x); GET_FLT128_WORDS64(hx,lx,x);
ix = 0x7fffffffffffffffULL&hx; ix = 0x7fffffffffffffffULL&hx;
*eptr = 0; *eptr = 0;
if(ix>=0x7fff000000000000ULL||((ix|lx)==0)) return x + x;/* 0,inf,nan */ if(ix>=0x7fff000000000000ULL||((ix|lx)==0)) return x + x;/* 0,inf,nan */
if (ix<0x0001000000000000ULL) { /* subnormal */ if (ix<0x0001000000000000ULL) { /* subnormal */
x *= two114; x *= two114;
GET_FLT128_MSW64(hx,x); GET_FLT128_MSW64(hx,x);
ix = hx&0x7fffffffffffffffULL; ix = hx&0x7fffffffffffffffULL;
*eptr = -114; *eptr = -114;
} }
*eptr += (ix>>48)-16382; *eptr += (ix>>48)-16382;
hx = (hx&0x8000ffffffffffffULL) | 0x3ffe000000000000ULL; hx = (hx&0x8000ffffffffffffULL) | 0x3ffe000000000000ULL;
SET_FLT128_MSW64(x,hx); SET_FLT128_MSW64(x,hx);
return x; return x;
} }
/* e_hypotl.c -- long double version of e_hypot.c.
* Conversion to long double by Jakub Jelinek, jakub@redhat.com.
*/
/* /*
* ==================================================== * ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
...@@ -9,19 +13,14 @@ ...@@ -9,19 +13,14 @@
* ==================================================== * ====================================================
*/ */
/* From e_hypotl.c -- long double version of e_hypot.c.
* Conversion to long double by Jakub Jelinek, jakub@redhat.com.
* Conversion to __float128 by FX Coudert, fxcoudert@gcc.gnu.org.
*/
/* hypotq(x,y) /* hypotq(x,y)
* *
* Method : * Method :
* If (assume round-to-nearest) z=x*x+y*y * If (assume round-to-nearest) z=x*x+y*y
* has error less than sqrtl(2)/2 ulp, than * has error less than sqrtq(2)/2 ulp, than
* sqrtl(z) has error less than 1 ulp (exercise). * sqrtq(z) has error less than 1 ulp (exercise).
* *
* So, compute sqrtl(x*x+y*y) with some care as * So, compute sqrtq(x*x+y*y) with some care as
* follows to get the error below 1 ulp: * follows to get the error below 1 ulp:
* *
* Assume x>y>0; * Assume x>y>0;
...@@ -38,100 +37,102 @@ ...@@ -38,100 +37,102 @@
* large or too tiny * large or too tiny
* *
* Special cases: * Special cases:
* hypotq(x,y) is INF if x or y is +INF or -INF; else * hypotl(x,y) is INF if x or y is +INF or -INF; else
* hypotq(x,y) is NAN if x or y is NAN. * hypotl(x,y) is NAN if x or y is NAN.
* *
* Accuracy: * Accuracy:
* hypotq(x,y) returns sqrtl(x^2+y^2) with error less * hypotl(x,y) returns sqrtq(x^2+y^2) with error less
* than 1 ulps (units in the last place) * than 1 ulps (units in the last place)
*/ */
#include "quadmath-imp.h" #include "quadmath-imp.h"
__float128 __float128
hypotq (__float128 x, __float128 y) hypotq(__float128 x, __float128 y)
{ {
__float128 a, b, t1, t2, y1, y2, w; __float128 a,b,t1,t2,y1,y2,w;
int64_t j, k, ha, hb; int64_t j,k,ha,hb;
GET_FLT128_MSW64(ha,x); GET_FLT128_MSW64(ha,x);
ha &= 0x7fffffffffffffffLL; ha &= 0x7fffffffffffffffLL;
GET_FLT128_MSW64(hb,y); GET_FLT128_MSW64(hb,y);
hb &= 0x7fffffffffffffffLL; hb &= 0x7fffffffffffffffLL;
if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;} if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
SET_FLT128_MSW64(a,ha); /* a <- |a| */ SET_FLT128_MSW64(a,ha); /* a <- |a| */
SET_FLT128_MSW64(b,hb); /* b <- |b| */ SET_FLT128_MSW64(b,hb); /* b <- |b| */
if((ha-hb)>0x78000000000000LL) {return a+b;} /* x/y > 2**120 */ if((ha-hb)>0x78000000000000LL) {return a+b;} /* x/y > 2**120 */
k=0; k=0;
if(ha > 0x5f3f000000000000LL) { /* a>2**8000 */ if(ha > 0x5f3f000000000000LL) { /* a>2**8000 */
if(ha >= 0x7fff000000000000LL) { /* Inf or NaN */ if(ha >= 0x7fff000000000000LL) { /* Inf or NaN */
uint64_t low; uint64_t low;
w = a+b; /* for sNaN */ w = a+b; /* for sNaN */
GET_FLT128_LSW64(low,a); if (issignalingq (a) || issignalingq (b))
if(((ha&0xffffffffffffLL)|low)==0) w = a; return w;
GET_FLT128_LSW64(low,b); GET_FLT128_LSW64(low,a);
if(((hb^0x7fff000000000000LL)|low)==0) w = b; if(((ha&0xffffffffffffLL)|low)==0) w = a;
return w; GET_FLT128_LSW64(low,b);
} if(((hb^0x7fff000000000000LL)|low)==0) w = b;
/* scale a and b by 2**-9600 */ return w;
ha -= 0x2580000000000000LL; }
hb -= 0x2580000000000000LL; k += 9600; /* scale a and b by 2**-9600 */
SET_FLT128_MSW64(a,ha); ha -= 0x2580000000000000LL;
SET_FLT128_MSW64(b,hb); hb -= 0x2580000000000000LL; k += 9600;
} SET_FLT128_MSW64(a,ha);
if(hb < 0x20bf000000000000LL) { /* b < 2**-8000 */ SET_FLT128_MSW64(b,hb);
if(hb <= 0x0000ffffffffffffLL) { /* subnormal b or 0 */ }
uint64_t low; if(hb < 0x20bf000000000000LL) { /* b < 2**-8000 */
GET_FLT128_LSW64(low,b); if(hb <= 0x0000ffffffffffffLL) { /* subnormal b or 0 */
if((hb|low)==0) return a; uint64_t low;
t1=0; GET_FLT128_LSW64(low,b);
SET_FLT128_MSW64(t1,0x7ffd000000000000LL); /* t1=2^16382 */ if((hb|low)==0) return a;
b *= t1; t1=0;
a *= t1; SET_FLT128_MSW64(t1,0x7ffd000000000000LL); /* t1=2^16382 */
k -= 16382; b *= t1;
GET_FLT128_MSW64 (ha, a); a *= t1;
GET_FLT128_MSW64 (hb, b); k -= 16382;
if (hb > ha) GET_FLT128_MSW64 (ha, a);
{ GET_FLT128_MSW64 (hb, b);
t1 = a; if (hb > ha)
a = b; {
b = t1; t1 = a;
j = ha; a = b;
ha = hb; b = t1;
hb = j; j = ha;
} ha = hb;
} else { /* scale a and b by 2^9600 */ hb = j;
ha += 0x2580000000000000LL; /* a *= 2^9600 */ }
hb += 0x2580000000000000LL; /* b *= 2^9600 */ } else { /* scale a and b by 2^9600 */
k -= 9600; ha += 0x2580000000000000LL; /* a *= 2^9600 */
SET_FLT128_MSW64(a,ha); hb += 0x2580000000000000LL; /* b *= 2^9600 */
SET_FLT128_MSW64(b,hb); k -= 9600;
} SET_FLT128_MSW64(a,ha);
} SET_FLT128_MSW64(b,hb);
}
}
/* medium size a and b */ /* medium size a and b */
w = a-b; w = a-b;
if (w>b) { if (w>b) {
t1 = 0; t1 = 0;
SET_FLT128_MSW64(t1,ha); SET_FLT128_MSW64(t1,ha);
t2 = a-t1; t2 = a-t1;
w = sqrtq(t1*t1-(b*(-b)-t2*(a+t1))); w = sqrtq(t1*t1-(b*(-b)-t2*(a+t1)));
} else { } else {
a = a+a; a = a+a;
y1 = 0; y1 = 0;
SET_FLT128_MSW64(y1,hb); SET_FLT128_MSW64(y1,hb);
y2 = b - y1; y2 = b - y1;
t1 = 0; t1 = 0;
SET_FLT128_MSW64(t1,ha+0x0001000000000000LL); SET_FLT128_MSW64(t1,ha+0x0001000000000000LL);
t2 = a - t1; t2 = a - t1;
w = sqrtq(t1*y1-(w*(-w)-(t1*y2+t2*b))); w = sqrtq(t1*y1-(w*(-w)-(t1*y2+t2*b)));
} }
if(k!=0) { if(k!=0) {
uint64_t high; uint64_t high;
t1 = 1.0Q; t1 = 1;
GET_FLT128_MSW64(high,t1); GET_FLT128_MSW64(high,t1);
SET_FLT128_MSW64(t1,high+(k<<48)); SET_FLT128_MSW64(t1,high+(k<<48));
w *= t1; w *= t1;
math_check_force_underflow_nonneg (w); math_check_force_underflow_nonneg (w);
return w; return w;
} else return w; } else return w;
} }
/* ilogbq.c -- __float128 version of s_ilogb.c. /* s_ilogbl.c -- long double version of s_ilogb.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz. * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/ */
...@@ -17,18 +17,15 @@ ...@@ -17,18 +17,15 @@
static char rcsid[] = "$NetBSD: $"; static char rcsid[] = "$NetBSD: $";
#endif #endif
/* ilogbl(__float128 x) /* ilogbl(long double x)
* return the binary exponent of non-zero x * return the binary exponent of non-zero x
* ilogbl(0) = FP_ILOGB0 * ilogbl(0) = FP_ILOGB0
* ilogbl(NaN) = FP_ILOGBNAN (no signal is raised) * ilogbl(NaN) = FP_ILOGBNAN (no signal is raised)
* ilogbl(+-Inf) = INT_MAX (no signal is raised) * ilogbl(+-Inf) = INT_MAX (no signal is raised)
*/ */
#include <limits.h>
#include <math.h> #include <math.h>
#include <errno.h>
#include "quadmath-imp.h" #include "quadmath-imp.h"
#ifndef FP_ILOGB0 #ifndef FP_ILOGB0
# define FP_ILOGB0 INT_MIN # define FP_ILOGB0 INT_MIN
#endif #endif
...@@ -36,8 +33,7 @@ static char rcsid[] = "$NetBSD: $"; ...@@ -36,8 +33,7 @@ static char rcsid[] = "$NetBSD: $";
# define FP_ILOGBNAN INT_MAX # define FP_ILOGBNAN INT_MAX
#endif #endif
int int ilogbq (__float128 x)
ilogbq (__float128 x)
{ {
int64_t hx,lx; int64_t hx,lx;
int ix; int ix;
...@@ -46,13 +42,7 @@ ilogbq (__float128 x) ...@@ -46,13 +42,7 @@ ilogbq (__float128 x)
hx &= 0x7fffffffffffffffLL; hx &= 0x7fffffffffffffffLL;
if(hx <= 0x0001000000000000LL) { if(hx <= 0x0001000000000000LL) {
if((hx|lx)==0) if((hx|lx)==0)
{ { errno = EDOM; feraiseexcept (FE_INVALID); return FP_ILOGB0; } /* ilogbl(0) = FP_ILOGB0 */
errno = EDOM;
#ifdef USE_FENV_H
feraiseexcept (FE_INVALID);
#endif
return FP_ILOGB0; /* ilogbl(0) = FP_ILOGB0 */
}
else /* subnormal x */ else /* subnormal x */
if(hx==0) { if(hx==0) {
for (ix = -16431; lx>0; lx<<=1) ix -=1; for (ix = -16431; lx>0; lx<<=1) ix -=1;
...@@ -65,18 +55,7 @@ ilogbq (__float128 x) ...@@ -65,18 +55,7 @@ ilogbq (__float128 x)
else if (FP_ILOGBNAN != INT_MAX) { else if (FP_ILOGBNAN != INT_MAX) {
/* ISO C99 requires ilogbl(+-Inf) == INT_MAX. */ /* ISO C99 requires ilogbl(+-Inf) == INT_MAX. */
if (((hx^0x7fff000000000000LL)|lx) == 0) if (((hx^0x7fff000000000000LL)|lx) == 0)
{ { errno = EDOM; feraiseexcept (FE_INVALID); return INT_MAX; }
errno = EDOM;
#ifdef USE_FENV_H
feraiseexcept (FE_INVALID);
#endif
return INT_MAX;
}
} }
{ errno = EDOM; feraiseexcept (FE_INVALID); return FP_ILOGBNAN; }
errno = EDOM;
#ifdef USE_FENV_H
feraiseexcept (FE_INVALID);
#endif
return FP_ILOGBNAN;
} }
/*
* Written by Ulrich Drepper <drepper@gmail.com>
*/
/*
* __quadmath_isinf_nsq (x) returns != 0 if x is ±inf, else 0;
* no branching!
*/
#include "quadmath-imp.h"
int
__quadmath_isinf_nsq (__float128 x)
{
int64_t hx,lx;
GET_FLT128_WORDS64(hx,lx,x);
return !(lx | ((hx & 0x7fffffffffffffffLL) ^ 0x7fff000000000000LL));
}
...@@ -4,14 +4,23 @@ ...@@ -4,14 +4,23 @@
* Public domain. * Public domain.
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/*
* isinfq(x) returns 1 if x is inf, -1 if x is -inf, else 0;
* no branching!
*/
#include "quadmath-imp.h" #include "quadmath-imp.h"
int int
isinfq (__float128 x) isinfq (__float128 x)
{ {
int64_t hx,lx; int64_t hx,lx;
GET_FLT128_WORDS64(hx,lx,x); GET_FLT128_WORDS64(hx,lx,x);
lx |= (hx & 0x7fffffffffffffffLL) ^ 0x7fff000000000000LL; lx |= (hx & 0x7fffffffffffffffLL) ^ 0x7fff000000000000LL;
lx |= -lx; lx |= -lx;
return ~(lx >> 63) & (hx >> 62); return ~(lx >> 63) & (hx >> 62);
} }
/* isnanq.c -- __float128 version of s_isnan.c. /* s_isnanl.c -- long double version of s_isnan.c.
* Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz. * Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz.
*/ */
...@@ -13,15 +13,23 @@ ...@@ -13,15 +13,23 @@
* ==================================================== * ====================================================
*/ */
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: $";
#endif
/*
* isnanq(x) returns 1 is x is nan, else 0;
* no branching!
*/
#include "quadmath-imp.h" #include "quadmath-imp.h"
int int isnanq(__float128 x)
isnanq (const __float128 x)
{ {
int64_t hx,lx; int64_t hx,lx;
GET_FLT128_WORDS64(hx,lx,x); GET_FLT128_WORDS64(hx,lx,x);
hx &= 0x7fffffffffffffffLL; hx &= 0x7fffffffffffffffLL;
hx |= (uint64_t)(lx|(-lx))>>63; hx |= (uint64_t)(lx|(-lx))>>63;
hx = 0x7fff000000000000LL - hx; hx = 0x7fff000000000000LL - hx;
return (int)((uint64_t)hx>>63); return (int)((uint64_t)hx>>63);
} }
/* Test for signaling NaN.
Copyright (C) 2013-2018 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 "quadmath-imp.h"
int
issignalingq (__float128 x)
{
uint64_t hxi, lxi __attribute__ ((unused));
GET_FLT128_WORDS64 (hxi, lxi, x);
#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN
/* We only have to care about the high-order bit of x's significand, because
having it set (sNaN) already makes the significand different from that
used to designate infinity. */
return ((hxi & UINT64_C (0x7fff800000000000))
== UINT64_C (0x7fff800000000000));
#else
/* To keep the following comparison simple, toggle the quiet/signaling bit,
so that it is set for sNaNs. This is inverse to IEEE 754-2008 (as well as
common practice for IEEE 754-1985). */
hxi ^= UINT64_C (0x0000800000000000);
/* If lxi != 0, then set any suitable bit of the significand in hxi. */
hxi |= (lxi | -lxi) >> 63;
/* We have to compare for greater (instead of greater or equal), because x's
significand being all-zero designates infinity not NaN. */
return (hxi & UINT64_C (0x7fffffffffffffff)) > UINT64_C (0x7fff800000000000);
#endif
}
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* *
* SYNOPSIS: * SYNOPSIS:
* *
* __float128 x, y, j0l(); * long double x, y, j0l();
* *
* y = j0l( x ); * y = j0l( x );
* *
...@@ -52,7 +52,7 @@ ...@@ -52,7 +52,7 @@
* *
* SYNOPSIS: * SYNOPSIS:
* *
* __float128 x, y, y0l(); * double x, y, y0l();
* *
* y = y0l( x ); * y = y0l( x );
* *
...@@ -88,8 +88,8 @@ ...@@ -88,8 +88,8 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software License along with this library; if not, see
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ <http://www.gnu.org/licenses/>. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
...@@ -97,7 +97,7 @@ ...@@ -97,7 +97,7 @@
static const __float128 ONEOSQPI = 5.6418958354775628694807945156077258584405E-1Q; static const __float128 ONEOSQPI = 5.6418958354775628694807945156077258584405E-1Q;
/* 2 / pi */ /* 2 / pi */
static const __float128 TWOOPI = 6.3661977236758134307553505349005744813784E-1Q; static const __float128 TWOOPI = 6.3661977236758134307553505349005744813784E-1Q;
static const __float128 zero = 0.0Q; static const __float128 zero = 0;
/* J0(x) = 1 - x^2/4 + x^2 x^2 R(x^2) /* J0(x) = 1 - x^2/4 + x^2 x^2 R(x^2)
Peak relative error 3.4e-37 Peak relative error 3.4e-37
...@@ -683,21 +683,21 @@ j0q (__float128 x) ...@@ -683,21 +683,21 @@ j0q (__float128 x)
if (x != x) if (x != x)
return x + x; return x + x;
else else
return 0.0Q; return 0;
} }
if (x == 0.0Q) if (x == 0)
return 1.0Q; return 1;
xx = fabsq (x); xx = fabsq (x);
if (xx <= 2.0Q) if (xx <= 2)
{ {
if (xx < 0x1p-57Q) if (xx < 0x1p-57Q)
return 1.0Q; return 1;
/* 0 <= x <= 2 */ /* 0 <= x <= 2 */
z = xx * xx; z = xx * xx;
p = z * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D); p = z * z * neval (z, J0_2N, NJ0_2N) / deval (z, J0_2D, NJ0_2D);
p -= 0.25Q * z; p -= 0.25Q * z;
p += 1.0Q; p += 1;
return p; return p;
} }
...@@ -711,7 +711,7 @@ j0q (__float128 x) ...@@ -711,7 +711,7 @@ j0q (__float128 x)
sincosq (xx, &s, &c); sincosq (xx, &s, &c);
ss = s - c; ss = s - c;
cc = s + c; cc = s + c;
if (xx <= FLT128_MAX / 2.0Q) if (xx <= FLT128_MAX / 2)
{ {
z = -cosq (xx + xx); z = -cosq (xx + xx);
if ((s * c) < 0) if ((s * c) < 0)
...@@ -723,7 +723,7 @@ j0q (__float128 x) ...@@ -723,7 +723,7 @@ j0q (__float128 x)
if (xx > 0x1p256Q) if (xx > 0x1p256Q)
return ONEOSQPI * cc / sqrtq (xx); return ONEOSQPI * cc / sqrtq (xx);
xinv = 1.0Q / xx; xinv = 1 / xx;
z = xinv * xinv; z = xinv * xinv;
if (xinv <= 0.25) if (xinv <= 0.25)
{ {
...@@ -781,7 +781,7 @@ j0q (__float128 x) ...@@ -781,7 +781,7 @@ j0q (__float128 x)
q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D); q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
} }
} }
p = 1.0Q + z * p; p = 1 + z * p;
q = z * xinv * q; q = z * xinv * q;
q = q - 0.125Q * xinv; q = q - 0.125Q * xinv;
z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx); z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx);
...@@ -789,11 +789,12 @@ j0q (__float128 x) ...@@ -789,11 +789,12 @@ j0q (__float128 x)
} }
/* Y0(x) = 2/pi * log(x) * J0(x) + R(x^2) /* Y0(x) = 2/pi * log(x) * J0(x) + R(x^2)
Peak absolute error 1.7e-36 (relative where Y0 > 1) Peak absolute error 1.7e-36 (relative where Y0 > 1)
0 <= x <= 2 */ 0 <= x <= 2 */
#define NY0_2N 7 #define NY0_2N 7
static __float128 Y0_2N[NY0_2N + 1] = { static const __float128 Y0_2N[NY0_2N + 1] = {
-1.062023609591350692692296993537002558155E19Q, -1.062023609591350692692296993537002558155E19Q,
2.542000883190248639104127452714966858866E19Q, 2.542000883190248639104127452714966858866E19Q,
-1.984190771278515324281415820316054696545E18Q, -1.984190771278515324281415820316054696545E18Q,
...@@ -804,7 +805,7 @@ static __float128 Y0_2N[NY0_2N + 1] = { ...@@ -804,7 +805,7 @@ static __float128 Y0_2N[NY0_2N + 1] = {
8.230845651379566339707130644134372793322E6Q, 8.230845651379566339707130644134372793322E6Q,
}; };
#define NY0_2D 7 #define NY0_2D 7
static __float128 Y0_2D[NY0_2D + 1] = { static const __float128 Y0_2D[NY0_2D + 1] = {
1.438972634353286978700329883122253752192E20Q, 1.438972634353286978700329883122253752192E20Q,
1.856409101981569254247700169486907405500E18Q, 1.856409101981569254247700169486907405500E18Q,
1.219693352678218589553725579802986255614E16Q, 1.219693352678218589553725579802986255614E16Q,
...@@ -821,22 +822,22 @@ static const __float128 U0 = -7.3804295108687225274343927948483016310862e-02Q; ...@@ -821,22 +822,22 @@ static const __float128 U0 = -7.3804295108687225274343927948483016310862e-02Q;
/* Bessel function of the second kind, order zero. */ /* Bessel function of the second kind, order zero. */
__float128 __float128
y0q (__float128 x) y0q(__float128 x)
{ {
__float128 xx, xinv, z, p, q, c, s, cc, ss; __float128 xx, xinv, z, p, q, c, s, cc, ss;
if (! finiteq (x)) if (! finiteq (x))
return 1 / (x + x * x); return 1 / (x + x * x);
if (x <= 0.0Q) if (x <= 0)
{ {
if (x < 0.0Q) if (x < 0)
return (zero / (zero * x)); return (zero / (zero * x));
return -1 / zero; /* -inf and divide by zero exception. */ return -1 / zero; /* -inf and divide by zero exception. */
} }
xx = fabsq (x); xx = fabsq (x);
if (xx <= 0x1p-57) if (xx <= 0x1p-57)
return U0 + TWOOPI * logq (x); return U0 + TWOOPI * logq (x);
if (xx <= 2.0Q) if (xx <= 2)
{ {
/* 0 <= x <= 2 */ /* 0 <= x <= 2 */
z = xx * xx; z = xx * xx;
...@@ -855,7 +856,7 @@ y0q (__float128 x) ...@@ -855,7 +856,7 @@ y0q (__float128 x)
sincosq (x, &s, &c); sincosq (x, &s, &c);
ss = s - c; ss = s - c;
cc = s + c; cc = s + c;
if (xx <= FLT128_MAX / 2.0Q) if (xx <= FLT128_MAX / 2)
{ {
z = -cosq (x + x); z = -cosq (x + x);
if ((s * c) < 0) if ((s * c) < 0)
...@@ -867,7 +868,7 @@ y0q (__float128 x) ...@@ -867,7 +868,7 @@ y0q (__float128 x)
if (xx > 0x1p256Q) if (xx > 0x1p256Q)
return ONEOSQPI * ss / sqrtq (x); return ONEOSQPI * ss / sqrtq (x);
xinv = 1.0Q / xx; xinv = 1 / xx;
z = xinv * xinv; z = xinv * xinv;
if (xinv <= 0.25) if (xinv <= 0.25)
{ {
...@@ -925,7 +926,7 @@ y0q (__float128 x) ...@@ -925,7 +926,7 @@ y0q (__float128 x)
q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D); q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
} }
} }
p = 1.0Q + z * p; p = 1 + z * p;
q = z * xinv * q; q = z * xinv * q;
q = q - 0.125Q * xinv; q = q - 0.125Q * xinv;
z = ONEOSQPI * (p * ss + q * cc) / sqrtq (x); z = ONEOSQPI * (p * ss + q * cc) / sqrtq (x);
......
...@@ -6,9 +6,9 @@ ...@@ -6,9 +6,9 @@
* *
* SYNOPSIS: * SYNOPSIS:
* *
* __float128 x, y, j1q(); * long double x, y, j1l();
* *
* y = j1q( x ); * y = j1l( x );
* *
* *
* *
...@@ -52,9 +52,9 @@ ...@@ -52,9 +52,9 @@
* *
* SYNOPSIS: * SYNOPSIS:
* *
* __float128, y, y1q(); * double x, y, y1l();
* *
* y = y1q( x ); * y = y1l( x );
* *
* *
* *
...@@ -92,17 +92,16 @@ ...@@ -92,17 +92,16 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software License along with this library; if not, see
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ <http://www.gnu.org/licenses/>. */
#include <errno.h>
#include "quadmath-imp.h" #include "quadmath-imp.h"
/* 1 / sqrt(pi) */ /* 1 / sqrt(pi) */
static const __float128 ONEOSQPI = 5.6418958354775628694807945156077258584405E-1Q; static const __float128 ONEOSQPI = 5.6418958354775628694807945156077258584405E-1Q;
/* 2 / pi */ /* 2 / pi */
static const __float128 TWOOPI = 6.3661977236758134307553505349005744813784E-1Q; static const __float128 TWOOPI = 6.3661977236758134307553505349005744813784E-1Q;
static const __float128 zero = 0.0Q; static const __float128 zero = 0;
/* J1(x) = .5x + x x^2 R(x^2) /* J1(x) = .5x + x x^2 R(x^2)
Peak relative error 1.9e-35 Peak relative error 1.9e-35
...@@ -126,7 +125,7 @@ static const __float128 J0_2D[NJ0_2D + 1] = { ...@@ -126,7 +125,7 @@ static const __float128 J0_2D[NJ0_2D + 1] = {
5.673775894803172808323058205986256928794E8Q, 5.673775894803172808323058205986256928794E8Q,
1.080329960080981204840966206372671147224E6Q, 1.080329960080981204840966206372671147224E6Q,
1.411951256636576283942477881535283304912E3Q, 1.411951256636576283942477881535283304912E3Q,
/* 1.000000000000000000000000000000000000000E0Q */ /* 1.000000000000000000000000000000000000000E0L */
}; };
/* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2), /* J1(x)cosX + Y1(x)sinX = sqrt( 2/(pi x)) P1(x), P1(x) = 1 + 1/x^2 R(1/x^2),
...@@ -690,9 +689,9 @@ j1q (__float128 x) ...@@ -690,9 +689,9 @@ j1q (__float128 x)
if (x != x) if (x != x)
return x + x; return x + x;
else else
return 0.0Q; return 0;
} }
if (x == 0.0Q) if (x == 0)
return x; return x;
xx = fabsq (x); xx = fabsq (x);
if (xx <= 0x1p-58Q) if (xx <= 0x1p-58Q)
...@@ -703,7 +702,7 @@ j1q (__float128 x) ...@@ -703,7 +702,7 @@ j1q (__float128 x)
errno = ERANGE; errno = ERANGE;
return ret; return ret;
} }
if (xx <= 2.0Q) if (xx <= 2)
{ {
/* 0 <= x <= 2 */ /* 0 <= x <= 2 */
z = xx * xx; z = xx * xx;
...@@ -723,7 +722,7 @@ j1q (__float128 x) ...@@ -723,7 +722,7 @@ j1q (__float128 x)
sincosq (xx, &s, &c); sincosq (xx, &s, &c);
ss = -s - c; ss = -s - c;
cc = s - c; cc = s - c;
if (xx <= FLT128_MAX / 2.0Q) if (xx <= FLT128_MAX / 2)
{ {
z = cosq (xx + xx); z = cosq (xx + xx);
if ((s * c) > 0) if ((s * c) > 0)
...@@ -740,7 +739,7 @@ j1q (__float128 x) ...@@ -740,7 +739,7 @@ j1q (__float128 x)
return z; return z;
} }
xinv = 1.0Q / xx; xinv = 1 / xx;
z = xinv * xinv; z = xinv * xinv;
if (xinv <= 0.25) if (xinv <= 0.25)
{ {
...@@ -798,7 +797,7 @@ j1q (__float128 x) ...@@ -798,7 +797,7 @@ j1q (__float128 x)
q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D); q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
} }
} }
p = 1.0Q + z * p; p = 1 + z * p;
q = z * q; q = z * q;
q = q * xinv + 0.375Q * xinv; q = q * xinv + 0.375Q * xinv;
z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx); z = ONEOSQPI * (p * cc - q * ss) / sqrtq (xx);
...@@ -808,11 +807,12 @@ j1q (__float128 x) ...@@ -808,11 +807,12 @@ j1q (__float128 x)
} }
/* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2) /* Y1(x) = 2/pi * (log(x) * J1(x) - 1/x) + x R(x^2)
Peak relative error 6.2e-38 Peak relative error 6.2e-38
0 <= x <= 2 */ 0 <= x <= 2 */
#define NY0_2N 7 #define NY0_2N 7
static __float128 Y0_2N[NY0_2N + 1] = { static const __float128 Y0_2N[NY0_2N + 1] = {
-6.804415404830253804408698161694720833249E19Q, -6.804415404830253804408698161694720833249E19Q,
1.805450517967019908027153056150465849237E19Q, 1.805450517967019908027153056150465849237E19Q,
-8.065747497063694098810419456383006737312E17Q, -8.065747497063694098810419456383006737312E17Q,
...@@ -823,7 +823,7 @@ static __float128 Y0_2N[NY0_2N + 1] = { ...@@ -823,7 +823,7 @@ static __float128 Y0_2N[NY0_2N + 1] = {
9.541172044989995856117187515882879304461E5Q, 9.541172044989995856117187515882879304461E5Q,
}; };
#define NY0_2D 7 #define NY0_2D 7
static __float128 Y0_2D[NY0_2D + 1] = { static const __float128 Y0_2D[NY0_2D + 1] = {
3.470629591820267059538637461549677594549E20Q, 3.470629591820267059538637461549677594549E20Q,
4.120796439009916326855848107545425217219E18Q, 4.120796439009916326855848107545425217219E18Q,
2.477653371652018249749350657387030814542E16Q, 2.477653371652018249749350657387030814542E16Q,
...@@ -845,9 +845,9 @@ y1q (__float128 x) ...@@ -845,9 +845,9 @@ y1q (__float128 x)
if (! finiteq (x)) if (! finiteq (x))
return 1 / (x + x * x); return 1 / (x + x * x);
if (x <= 0.0Q) if (x <= 0)
{ {
if (x < 0.0Q) if (x < 0)
return (zero / (zero * x)); return (zero / (zero * x));
return -1 / zero; /* -inf and divide by zero exception. */ return -1 / zero; /* -inf and divide by zero exception. */
} }
...@@ -859,10 +859,10 @@ y1q (__float128 x) ...@@ -859,10 +859,10 @@ y1q (__float128 x)
errno = ERANGE; errno = ERANGE;
return z; return z;
} }
if (xx <= 2.0Q) if (xx <= 2)
{ {
/* 0 <= x <= 2 */ /* 0 <= x <= 2 */
/* FIXME: SET_RESTORE_ROUNDL (FE_TONEAREST); */ SET_RESTORE_ROUNDF128 (FE_TONEAREST);
z = xx * xx; z = xx * xx;
p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D); p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
p = -TWOOPI / xx + p; p = -TWOOPI / xx + p;
...@@ -879,7 +879,7 @@ y1q (__float128 x) ...@@ -879,7 +879,7 @@ y1q (__float128 x)
sincosq (xx, &s, &c); sincosq (xx, &s, &c);
ss = -s - c; ss = -s - c;
cc = s - c; cc = s - c;
if (xx <= FLT128_MAX / 2.0Q) if (xx <= FLT128_MAX / 2)
{ {
z = cosq (xx + xx); z = cosq (xx + xx);
if ((s * c) > 0) if ((s * c) > 0)
...@@ -891,7 +891,7 @@ y1q (__float128 x) ...@@ -891,7 +891,7 @@ y1q (__float128 x)
if (xx > 0x1p256Q) if (xx > 0x1p256Q)
return ONEOSQPI * ss / sqrtq (xx); return ONEOSQPI * ss / sqrtq (xx);
xinv = 1.0Q / xx; xinv = 1 / xx;
z = xinv * xinv; z = xinv * xinv;
if (xinv <= 0.25) if (xinv <= 0.25)
{ {
...@@ -949,7 +949,7 @@ y1q (__float128 x) ...@@ -949,7 +949,7 @@ y1q (__float128 x)
q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D); q = neval (z, Q2_2r3N, NQ2_2r3N) / deval (z, Q2_2r3D, NQ2_2r3D);
} }
} }
p = 1.0Q + z * p; p = 1 + z * p;
q = z * q; q = z * q;
q = q * xinv + 0.375Q * xinv; q = q * xinv + 0.375Q * xinv;
z = ONEOSQPI * (p * ss + q * cc) / sqrtq (xx); z = ONEOSQPI * (p * ss + q * cc) / sqrtq (xx);
......
/* ldexpq.c -- __float128 version of s_ldexp.c. /* @(#)s_ldexp.c 5.1 93/09/24 */
* Conversion to long double by Ulrich Drepper,
* Cygnus Support, drepper@cygnus.com.
*/
/* /*
* ==================================================== * ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
...@@ -14,14 +10,23 @@ ...@@ -14,14 +10,23 @@
* ==================================================== * ====================================================
*/ */
#include <errno.h> #if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: s_ldexp.c,v 1.6 1995/05/10 20:47:40 jtc Exp $";
#endif
#include "quadmath-imp.h" #include "quadmath-imp.h"
__float128 __float128
ldexpq (__float128 value, int exp) ldexpq (__float128 value, int exp)
{ {
if(!finiteq(value)||value==0.0Q) return value; if(!finiteq(value)||value==0) return value + value;
value = scalbnq(value,exp); value = scalbnq(value,exp);
if(!finiteq(value)||value==0.0Q) errno = ERANGE; if(!finiteq(value)||value==0) errno = ERANGE;
return value; return value;
} }
/* Note, versioning issues are punted to ldbl-opt in this case. */
/* Compute a product of 1 + (T/X), 1 + (T/(X+1)), ....
Copyright (C) 2015-2018 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 "quadmath-imp.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. */
__float128
__quadmath_lgamma_productq (__float128 t, __float128 x, __float128 x_eps, int n)
{
__float128 ret = 0, ret_eps = 0;
for (int i = 0; i < n; i++)
{
__float128 xi = x + i;
__float128 quot = t / xi;
__float128 mhi, mlo;
mul_splitq (&mhi, &mlo, quot, xi);
__float128 quot_lo = (t - mhi - mlo) / xi - t * x_eps / (xi * xi);
/* We want (1 + RET + RET_EPS) * (1 + QUOT + QUOT_LO) - 1. */
__float128 rhi, rlo;
mul_splitq (&rhi, &rlo, ret, quot);
__float128 rpq = ret + quot;
__float128 rpq_eps = (ret - rpq) + quot;
__float128 nret = rpq + rhi;
__float128 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;
}
/* Round argument to nearest integral value according to current rounding /* Round argument to nearest integral value according to current rounding
direction. direction.
Copyright (C) 1997-2017 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
Jakub Jelinek <jj@ultra.linux.cz>, 1999. Jakub Jelinek <jj@ultra.linux.cz>, 1999.
...@@ -16,9 +16,8 @@ ...@@ -16,9 +16,8 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
...@@ -53,9 +52,7 @@ llrintq (__float128 x) ...@@ -53,9 +52,7 @@ llrintq (__float128 x)
/* In the event of overflow we must raise the "invalid" /* In the event of overflow we must raise the "invalid"
exception, but not "inexact". */ exception, but not "inexact". */
t = nearbyintq (x); t = nearbyintq (x);
#ifdef USE_FENV_H
feraiseexcept (t == LLONG_MAX ? FE_INEXACT : FE_INVALID); feraiseexcept (t == LLONG_MAX ? FE_INEXACT : FE_INVALID);
#endif
} }
else else
#endif #endif
...@@ -82,21 +79,21 @@ llrintq (__float128 x) ...@@ -82,21 +79,21 @@ llrintq (__float128 x)
unspecified. */ unspecified. */
#if defined FE_INVALID || defined FE_INEXACT #if defined FE_INVALID || defined FE_INEXACT
if (x < (__float128) LLONG_MIN if (x < (__float128) LLONG_MIN
&& x > (__float128) LLONG_MIN - 1.0Q) && x > (__float128) LLONG_MIN - 1)
{ {
/* If truncation produces LLONG_MIN, the cast will not raise /* If truncation produces LLONG_MIN, the cast will not raise
the exception, but may raise "inexact". */ the exception, but may raise "inexact". */
t = nearbyintq (x); t = nearbyintq (x);
#ifdef USE_FENV_H
feraiseexcept (t == LLONG_MIN ? FE_INEXACT : FE_INVALID); feraiseexcept (t == LLONG_MIN ? FE_INEXACT : FE_INVALID);
#endif
return LLONG_MIN; return LLONG_MIN;
} }
else if (FIX_FLT128_LLONG_CONVERT_OVERFLOW && x != (__float128) LLONG_MIN)
{
feraiseexcept (FE_INVALID);
return sx == 0 ? LLONG_MAX : LLONG_MIN;
}
#endif #endif
/* The number is too large. It is left implementation defined
what happens. */
return (long long int) x; return (long long int) x;
} }
......
/* Round __float128 value to long long int. /* Round long double value to long long int.
Copyright (C) 1997-2017 Free Software Foundation, Inc. Copyright (C) 1997-2018 Free Software Foundation, Inc.
This file is part of the GNU C Library. This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
Jakub Jelinek <jj@ultra.linux.cz>, 1999. Jakub Jelinek <jj@ultra.linux.cz>, 1999.
...@@ -15,13 +15,11 @@ ...@@ -15,13 +15,11 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free License along with the GNU C Library; if not, see
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA <http://www.gnu.org/licenses/>. */
02111-1307 USA. */
#include "quadmath-imp.h" #include "quadmath-imp.h"
long long int long long int
llroundq (__float128 x) llroundq (__float128 x)
{ {
...@@ -61,7 +59,7 @@ llroundq (__float128 x) ...@@ -61,7 +59,7 @@ llroundq (__float128 x)
else else
{ {
result = ((long long int) i0 << (j0 - 48)) | (j >> (112 - j0)); result = ((long long int) i0 << (j0 - 48)) | (j >> (112 - j0));
#if defined FE_INVALID && defined USE_FENV_H #ifdef FE_INVALID
if (sign == 1 && result == LLONG_MIN) if (sign == 1 && result == LLONG_MIN)
/* Rounding brought the value out of range. */ /* Rounding brought the value out of range. */
feraiseexcept (FE_INVALID); feraiseexcept (FE_INVALID);
...@@ -75,13 +73,18 @@ llroundq (__float128 x) ...@@ -75,13 +73,18 @@ llroundq (__float128 x)
FE_INVALID must be raised and the return value is FE_INVALID must be raised and the return value is
unspecified. */ unspecified. */
#ifdef FE_INVALID #ifdef FE_INVALID
if (x <= (__float128) LLONG_MIN - 0.5Q) if (FIX_FLT128_LLONG_CONVERT_OVERFLOW
&& !(sign == -1 && x > (__float128) LLONG_MIN - 0.5Q))
{
feraiseexcept (FE_INVALID);
return sign == 1 ? LLONG_MAX : LLONG_MIN;
}
else if (!FIX_FLT128_LLONG_CONVERT_OVERFLOW
&& x <= (__float128) LLONG_MIN - 0.5Q)
{ {
/* If truncation produces LLONG_MIN, the cast will not raise /* If truncation produces LLONG_MIN, the cast will not raise
the exception, but may raise "inexact". */ the exception, but may raise "inexact". */
#ifdef USE_FENV_H
feraiseexcept (FE_INVALID); feraiseexcept (FE_INVALID);
#endif
return LLONG_MIN; return LLONG_MIN;
} }
#endif #endif
......
/* log10q.c /* log10l.c
* *
* Common logarithm, 128-bit __float128 precision * Common logarithm, 128-bit long double precision
* *
* *
* *
* SYNOPSIS: * SYNOPSIS:
* *
* __float128 x, y, log10l(); * long double x, y, log10l();
* *
* y = log10q( x ); * y = log10l( x );
* *
* *
* *
...@@ -57,9 +57,7 @@ ...@@ -57,9 +57,7 @@
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software License along with this library; if not, see <http://www.gnu.org/licenses/>.
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/ */
#include "quadmath-imp.h" #include "quadmath-imp.h"
...@@ -99,7 +97,7 @@ static const __float128 Q[12] = ...@@ -99,7 +97,7 @@ static const __float128 Q[12] =
9.147150349299596453976674231612674085381E3Q, 9.147150349299596453976674231612674085381E3Q,
9.104928120962988414618126155557301584078E2Q, 9.104928120962988414618126155557301584078E2Q,
4.839208193348159620282142911143429644326E1Q 4.839208193348159620282142911143429644326E1Q
/* 1.000000000000000000000000000000000000000E0Q, */ /* 1.000000000000000000000000000000000000000E0L, */
}; };
/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2), /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
...@@ -125,7 +123,7 @@ static const __float128 S[6] = ...@@ -125,7 +123,7 @@ static const __float128 S[6] =
-5.748542087379434595104154610899551484314E4Q, -5.748542087379434595104154610899551484314E4Q,
3.998526750980007367835804959888064681098E3Q, 3.998526750980007367835804959888064681098E3Q,
-1.186359407982897997337150403816839480438E2Q -1.186359407982897997337150403816839480438E2Q
/* 1.000000000000000000000000000000000000000E0Q, */ /* 1.000000000000000000000000000000000000000E0L, */
}; };
static const __float128 static const __float128
...@@ -188,14 +186,14 @@ log10q (__float128 x) ...@@ -188,14 +186,14 @@ log10q (__float128 x)
/* Test for domain */ /* Test for domain */
GET_FLT128_WORDS64 (hx, lx, x); GET_FLT128_WORDS64 (hx, lx, x);
if (((hx & 0x7fffffffffffffffLL) | lx) == 0) if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
return (-1.0Q / fabsq (x)); /* log10l(+-0)=-inf */ return (-1 / fabsq (x)); /* log10l(+-0)=-inf */
if (hx < 0) if (hx < 0)
return (x - x) / (x - x); return (x - x) / (x - x);
if (hx >= 0x7fff000000000000LL) if (hx >= 0x7fff000000000000LL)
return (x + x); return (x + x);
if (x == 1.0Q) if (x == 1)
return 0.0Q; return 0;
/* separate mantissa from exponent */ /* separate mantissa from exponent */
...@@ -234,11 +232,11 @@ log10q (__float128 x) ...@@ -234,11 +232,11 @@ log10q (__float128 x)
if (x < SQRTH) if (x < SQRTH)
{ {
e -= 1; e -= 1;
x = 2.0 * x - 1.0Q; /* 2x - 1 */ x = 2.0 * x - 1; /* 2x - 1 */
} }
else else
{ {
x = x - 1.0Q; x = x - 1;
} }
z = x * x; z = x * x;
y = x * (z * neval (x, P, 12) / deval (x, Q, 11)); y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment