+2012-11-01 Tobias Burnus <burnus@net-b.de>
+
+ * Makefile.am (libquadmath_la_SOURCES): Add new math/* files.
+ * Makefile.in: Regenerated.
+ * math/acoshq.c: Update comment.
+ * math/acosq.c: Ditto.
+ * math/asinhq.c: Ditto.
+ * math/asinq.c: Ditto.
+ * math/atan2q.c: Ditto.
+ * math/atanhq.c: Ditto.
+ * math/ceilq.c: Ditto.
+ * math/copysignq.c: Ditto.
+ * math/cosq.c: Ditto.
+ * math/coshq.c: Ditto.
+ * math/erfq.c: Ditto.
+ * math/fabsq.c: Ditto.
+ * math/finiteq.c: Ditto.
+ * math/floorq.c: Ditto.
+ * math/fmodq.c: Ditto.
+ * math/frexpq.c: Ditto.
+ * math/isnanq.c: Ditto.
+ * math/j0q.c: Ditto.
+ * math/j1q.c: Ditto.
+ * math/ldexpq.c: Ditto.
+ * math/llroundq.c: Ditto.
+ * math/log10q.c: Ditto.
+ * math/log1pq.c: Ditto.
+ * math/log2q.c: Ditto.
+ * math/logq.c: Ditto.
+ * math/lroundq.c: Ditto.
+ * math/modfq.c: Ditto.
+ * math/nextafterq.c: Ditto.
+ * math/powq.c: Ditto.
+ * math/rem_pio2q.c: Ditto.
+ * math/remainderq.c: Ditto.
+ * math/rintq.c: Ditto.
+ * math/roundq.c: Ditto.
+ * math/scalblnq.c: Ditto.
+ * math/scalbnq.c: Ditto.
+ * math/sincosq_kernel.c: Ditto.
+ * math/sinq.c: Ditto.
+ * math/tanq.c: Ditto.
+ * math/expq.c: Ditto.
+ (__expq_table, expq): Renamed local array from __expl_table.
+ * math/cosq_kernel.c (__quadmath_kernel_cosq): Fix sign handling.
+ * math/cacoshq.c: Changes from GLIBC; fix returned sign.
+ * math/casinhq.c: Changes from GLIBC to fix special-case.
+ * math/cbrtq.c: Use modified GLIBC version.
+ * math/complex.c (ccoshd, cexpq, clog10q, clogq, csinhq, csinq,
+ ctanhq, ctanq): Moved to separates files.
+ (mult_c128, div_c128): Removed no longer needed functions.
+ (cexpiq): Call sincosq instead of sinq and cosq.
+ (cosq): Call cosh(-re,im) instead of cosq/sinq/sinh/cosh.
+ * math/ccoshq.c (ccoshq): New file, moved from complex.c and
+ modified based on GLIBC.
+ * math/cexpq.c (cexp): Ditto.
+ * math/clog10q.c (clog10q): Ditto.
+ * math/clogq.c (clogq): Ditto.
+ * math/csinhq.c: Ditto.
+ * math/csinq.c: Ditto.
+ * math/csqrtq.c: Ditto.
+ * math/ctanhq.c: Ditto.
+ * math/ctanq.c: Ditto.
+ * math/fmaq.c (fmaq): Port TININESS_AFTER_ROUNDING handling
+ from GLIBC.
+ * math/ilogbq.c (ilogbq): Add errno = EDOM handling.
+ * math/isinf_nsq.c (__quadmath_isinf_nsq): New file, ported
+ from GLIBC.
+ * math/lgammaq.c (lgammaq): Add signgam handling.
+ * math/sinhq.c (sinhq): Fix sign handling.
+ * math/sinq_kernel.c (__quadmath_kernel_sinq): Ditto.
+ * math/tgammaq.c (tgammaq): Ditto.
+ * math/x2y2m1q.c: New file.
+ * quadmath-imp.h (TININESS_AFTER_ROUNDING): New define.
+ (__quadmath_x2y2m1q, __quadmath_isinf_nsq): New prototypes.
+
2012-10-31 Tobias Burnus <burnus@net-b.de>
Joseph Myers <joseph@codesourcery.com>
David S. Miller <davem@davemloft.net>
libsubincludedir = $(libdir)/gcc/$(target_alias)/$(gcc_version)/include
libquadmath_la_SOURCES = \
- math/acoshq.c math/fmodq.c math/acosq.c math/frexpq.c \
+ math/x2y2m1q.c math/isinf_nsq.c math/acoshq.c math/fmodq.c \
+ math/acosq.c math/frexpq.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/roundq.c math/atanhq.c math/isnanq.c math/scalblnq.c math/atanq.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/llrintq.c math/log2q.c math/lrintq.c math/nearbyintq.c math/remquoq.c \
+ math/ccoshq.c math/cexpq.c math/clog10q.c math/clogq.c math/csinq.c \
+ math/csinhq.c math/csqrtq.c math/ctanq.c math/ctanhq.c \
printf/addmul_1.c printf/add_n.c printf/cmp.c printf/divrem.c \
printf/flt1282mpn.c printf/fpioconst.c printf/lshift.c printf/mul_1.c \
printf/mul_n.c printf/mul.c printf/printf_fphex.c printf/printf_fp.c \
"$(DESTDIR)$(libsubincludedir)"
LTLIBRARIES = $(toolexeclib_LTLIBRARIES)
am__dirstamp = $(am__leading_dot)dirstamp
-@BUILD_LIBQUADMATH_TRUE@am_libquadmath_la_OBJECTS = math/acoshq.lo \
+@BUILD_LIBQUADMATH_TRUE@am_libquadmath_la_OBJECTS = math/x2y2m1q.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/isinf_nsq.lo math/acoshq.lo \
@BUILD_LIBQUADMATH_TRUE@ math/fmodq.lo math/acosq.lo \
@BUILD_LIBQUADMATH_TRUE@ math/frexpq.lo math/rem_pio2q.lo \
@BUILD_LIBQUADMATH_TRUE@ math/asinhq.lo math/hypotq.lo \
@BUILD_LIBQUADMATH_TRUE@ math/ilogbq.lo math/llrintq.lo \
@BUILD_LIBQUADMATH_TRUE@ math/log2q.lo math/lrintq.lo \
@BUILD_LIBQUADMATH_TRUE@ math/nearbyintq.lo math/remquoq.lo \
-@BUILD_LIBQUADMATH_TRUE@ printf/addmul_1.lo printf/add_n.lo \
-@BUILD_LIBQUADMATH_TRUE@ printf/cmp.lo printf/divrem.lo \
-@BUILD_LIBQUADMATH_TRUE@ printf/flt1282mpn.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/ccoshq.lo math/cexpq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/clog10q.lo math/clogq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/csinq.lo math/csinhq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/csqrtq.lo math/ctanq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/ctanhq.lo printf/addmul_1.lo \
+@BUILD_LIBQUADMATH_TRUE@ printf/add_n.lo printf/cmp.lo \
+@BUILD_LIBQUADMATH_TRUE@ printf/divrem.lo printf/flt1282mpn.lo \
@BUILD_LIBQUADMATH_TRUE@ printf/fpioconst.lo printf/lshift.lo \
@BUILD_LIBQUADMATH_TRUE@ printf/mul_1.lo printf/mul_n.lo \
@BUILD_LIBQUADMATH_TRUE@ printf/mul.lo printf/printf_fphex.lo \
@BUILD_LIBQUADMATH_TRUE@nodist_libsubinclude_HEADERS = quadmath.h quadmath_weak.h
@BUILD_LIBQUADMATH_TRUE@libsubincludedir = $(libdir)/gcc/$(target_alias)/$(gcc_version)/include
@BUILD_LIBQUADMATH_TRUE@libquadmath_la_SOURCES = \
-@BUILD_LIBQUADMATH_TRUE@ math/acoshq.c math/fmodq.c math/acosq.c math/frexpq.c \
+@BUILD_LIBQUADMATH_TRUE@ math/x2y2m1q.c math/isinf_nsq.c math/acoshq.c math/fmodq.c \
+@BUILD_LIBQUADMATH_TRUE@ math/acosq.c math/frexpq.c \
@BUILD_LIBQUADMATH_TRUE@ math/rem_pio2q.c math/asinhq.c math/hypotq.c math/remainderq.c \
@BUILD_LIBQUADMATH_TRUE@ math/asinq.c math/rintq.c math/atan2q.c math/isinfq.c \
@BUILD_LIBQUADMATH_TRUE@ math/roundq.c math/atanhq.c math/isnanq.c math/scalblnq.c math/atanq.c \
@BUILD_LIBQUADMATH_TRUE@ math/catanhq.c math/catanq.c math/cimagq.c math/conjq.c math/cprojq.c \
@BUILD_LIBQUADMATH_TRUE@ math/crealq.c math/fdimq.c math/fmaxq.c math/fminq.c math/ilogbq.c \
@BUILD_LIBQUADMATH_TRUE@ math/llrintq.c math/log2q.c math/lrintq.c math/nearbyintq.c math/remquoq.c \
+@BUILD_LIBQUADMATH_TRUE@ math/ccoshq.c math/cexpq.c math/clog10q.c math/clogq.c math/csinq.c \
+@BUILD_LIBQUADMATH_TRUE@ math/csinhq.c math/csqrtq.c math/ctanq.c math/ctanhq.c \
@BUILD_LIBQUADMATH_TRUE@ printf/addmul_1.c printf/add_n.c printf/cmp.c printf/divrem.c \
@BUILD_LIBQUADMATH_TRUE@ printf/flt1282mpn.c printf/fpioconst.c printf/lshift.c printf/mul_1.c \
@BUILD_LIBQUADMATH_TRUE@ printf/mul_n.c printf/mul.c printf/printf_fphex.c printf/printf_fp.c \
math/$(DEPDIR)/$(am__dirstamp):
@$(MKDIR_P) math/$(DEPDIR)
@: > math/$(DEPDIR)/$(am__dirstamp)
+math/x2y2m1q.lo: math/$(am__dirstamp) math/$(DEPDIR)/$(am__dirstamp)
+math/isinf_nsq.lo: math/$(am__dirstamp) math/$(DEPDIR)/$(am__dirstamp)
math/acoshq.lo: math/$(am__dirstamp) math/$(DEPDIR)/$(am__dirstamp)
math/fmodq.lo: math/$(am__dirstamp) math/$(DEPDIR)/$(am__dirstamp)
math/acosq.lo: math/$(am__dirstamp) math/$(DEPDIR)/$(am__dirstamp)
math/nearbyintq.lo: math/$(am__dirstamp) \
math/$(DEPDIR)/$(am__dirstamp)
math/remquoq.lo: math/$(am__dirstamp) math/$(DEPDIR)/$(am__dirstamp)
+math/ccoshq.lo: math/$(am__dirstamp) math/$(DEPDIR)/$(am__dirstamp)
+math/cexpq.lo: math/$(am__dirstamp) math/$(DEPDIR)/$(am__dirstamp)
+math/clog10q.lo: math/$(am__dirstamp) math/$(DEPDIR)/$(am__dirstamp)
+math/clogq.lo: math/$(am__dirstamp) math/$(DEPDIR)/$(am__dirstamp)
+math/csinq.lo: math/$(am__dirstamp) math/$(DEPDIR)/$(am__dirstamp)
+math/csinhq.lo: math/$(am__dirstamp) math/$(DEPDIR)/$(am__dirstamp)
+math/csqrtq.lo: math/$(am__dirstamp) math/$(DEPDIR)/$(am__dirstamp)
+math/ctanq.lo: math/$(am__dirstamp) math/$(DEPDIR)/$(am__dirstamp)
+math/ctanhq.lo: math/$(am__dirstamp) math/$(DEPDIR)/$(am__dirstamp)
printf/$(am__dirstamp):
@$(MKDIR_P) printf
@: > printf/$(am__dirstamp)
-rm -f math/catanq.lo
-rm -f math/cbrtq.$(OBJEXT)
-rm -f math/cbrtq.lo
+ -rm -f math/ccoshq.$(OBJEXT)
+ -rm -f math/ccoshq.lo
-rm -f math/ceilq.$(OBJEXT)
-rm -f math/ceilq.lo
+ -rm -f math/cexpq.$(OBJEXT)
+ -rm -f math/cexpq.lo
-rm -f math/cimagq.$(OBJEXT)
-rm -f math/cimagq.lo
+ -rm -f math/clog10q.$(OBJEXT)
+ -rm -f math/clog10q.lo
+ -rm -f math/clogq.$(OBJEXT)
+ -rm -f math/clogq.lo
-rm -f math/complex.$(OBJEXT)
-rm -f math/complex.lo
-rm -f math/conjq.$(OBJEXT)
-rm -f math/cprojq.lo
-rm -f math/crealq.$(OBJEXT)
-rm -f math/crealq.lo
+ -rm -f math/csinhq.$(OBJEXT)
+ -rm -f math/csinhq.lo
+ -rm -f math/csinq.$(OBJEXT)
+ -rm -f math/csinq.lo
+ -rm -f math/csqrtq.$(OBJEXT)
+ -rm -f math/csqrtq.lo
+ -rm -f math/ctanhq.$(OBJEXT)
+ -rm -f math/ctanhq.lo
+ -rm -f math/ctanq.$(OBJEXT)
+ -rm -f math/ctanq.lo
-rm -f math/erfq.$(OBJEXT)
-rm -f math/erfq.lo
-rm -f math/expm1q.$(OBJEXT)
-rm -f math/hypotq.lo
-rm -f math/ilogbq.$(OBJEXT)
-rm -f math/ilogbq.lo
+ -rm -f math/isinf_nsq.$(OBJEXT)
+ -rm -f math/isinf_nsq.lo
-rm -f math/isinfq.$(OBJEXT)
-rm -f math/isinfq.lo
-rm -f math/isnanq.$(OBJEXT)
-rm -f math/tgammaq.lo
-rm -f math/truncq.$(OBJEXT)
-rm -f math/truncq.lo
+ -rm -f math/x2y2m1q.$(OBJEXT)
+ -rm -f math/x2y2m1q.lo
-rm -f printf/add_n.$(OBJEXT)
-rm -f printf/add_n.lo
-rm -f printf/addmul_1.$(OBJEXT)
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/catanhq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/catanq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/cbrtq.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/ccoshq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/ceilq.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/cexpq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/cimagq.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/clog10q.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/clogq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/complex.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/conjq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/copysignq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/cosq_kernel.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/cprojq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/crealq.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/csinhq.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/csinq.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/csqrtq.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/ctanhq.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/ctanq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/erfq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/expm1q.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/expq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/frexpq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/hypotq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/ilogbq.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/isinf_nsq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/isinfq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/isnanq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/j0q.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/tanq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/tgammaq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/truncq.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/x2y2m1q.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@printf/$(DEPDIR)/add_n.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@printf/$(DEPDIR)/addmul_1.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@printf/$(DEPDIR)/cmp.Plo@am__quote@
-/* e_acoshl.c -- long double version of e_acosh.c.
+/* acoshq.c -- __float128 version of e_acosh.c.
* Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
* ====================================================
*/
-/* __ieee754_acoshl(x)
+/* acoshq(x)
* Method :
* Based on
* acoshl(x) = logl [ x + sqrtl(x*x-1) ]
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
-/* __ieee754_acosl(x)
+/* acosq(x)
* Method :
* acos(x) = pi/2 - asin(x)
* acos(-x) = pi/2 + asin(x)
* if x is NaN, return x itself;
* if |x|>1, return NaN with invalid signal.
*
- * Functions needed: __ieee754_sqrtl.
+ * Functions needed: sqrtq.
*/
#include "quadmath-imp.h"
-/* s_asinhl.c -- long double version of s_asinh.c.
+/* asinhq.c -- __float128 version of s_asinh.c.
* Conversion to long double by Ulrich Drepper,
* Cygnus Support, drepper@cygnus.com.
*/
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
-/* __ieee754_asin(x)
+/* asinq(x)
* Method :
* Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
* we approximate asin(x) on [0,0.5] by
-/* e_atan2l.c -- long double version of e_atan2.c.
+/* atan2q.c -- __float128 version of e_atan2.c.
* Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
* ====================================================
*/
-/* __ieee754_atanhl(x)
+/* atanhq(x)
* Method :
* 1.Reduced x to positive by atanh(-x) = -atanh(x)
* 2.For x>=0.5
__real__ res = 0.0;
__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 (signbit (__real__ res))
+ __real__ res = 0.0Q;
+ }
else
{
__complex128 y;
y = csqrtq (y);
- if (__real__ x < 0.0)
+ if (signbitq (x))
y = -y;
__real__ y += __real__ x;
__imag__ y += __imag__ x;
res = clogq (y);
-
- /* We have to use the positive branch. */
- if (__real__ res < 0.0)
- res = -res;
}
return res;
__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;
+/* cbrtq.c
+ *
+ * Cube root, __float128 precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * __float128 x, y, cbrtq();
+ *
+ * y = cbrtq( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the cube root of the argument, which may be negative.
+ *
+ * Range reduction involves determining the power of 2 of
+ * the argument. A polynomial of degree 2 applied to the
+ * mantissa, and multiplication by the cube root of 1, 2, or 4
+ * approximates the root to within about 0.1%. Then Newton's
+ * iteration is used three times to converge to an accurate
+ * result.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE -8,8 100000 1.3e-34 3.9e-35
+ * IEEE exp(+-707) 100000 1.3e-34 4.3e-35
+ *
+ */
+
+/*
+Cephes Math Library Release 2.2: January, 1991
+Copyright 1984, 1991 by Stephen L. Moshier
+Adapted for glibc October, 2001.
+
+ This 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.
+
+ This 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 this library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
+
#include "quadmath-imp.h"
-#include <math.h>
-#include <float.h>
+
+static const long double CBRT2 = 1.259921049894873164767210607278228350570251Q;
+static const long double CBRT4 = 1.587401051968199474751705639272308260391493Q;
+static const long double CBRT2I = 0.7937005259840997373758528196361541301957467Q;
+static const long double CBRT4I = 0.6299605249474365823836053036391141752851257Q;
+
__float128
-cbrtq (const __float128 x)
+cbrtq ( __float128 x)
{
- __float128 y;
- int exp, i;
+ int e, rem, sign;
+ __float128 z;
+
+ if (!finiteq (x))
+ return x + x;
if (x == 0)
- return x;
-
- if (isnanq (x))
- return x;
-
- if (x <= DBL_MAX && x >= DBL_MIN)
- {
- /* Use double result as starting point. */
- y = cbrt ((double) x);
-
- /* Two Newton iterations. */
- y -= 0.333333333333333333333333333333333333333333333333333Q
- * (y - x / (y * y));
- y -= 0.333333333333333333333333333333333333333333333333333Q
- * (y - x / (y * y));
- return y;
- }
-
-#ifdef HAVE_CBRTL
- if (x <= LDBL_MAX && x >= LDBL_MIN)
- {
- /* Use long double result as starting point. */
- y = cbrtl ((long double) x);
-
- /* One Newton iteration. */
- y -= 0.333333333333333333333333333333333333333333333333333Q
- * (y - x / (y * y));
- return y;
- }
-#endif
-
- /* If we're outside of the range of C types, we have to compute
- the initial guess the hard way. */
- y = frexpq (x, &exp);
-
- i = exp % 3;
- y = (i >= 0 ? i : -i);
- if (i == 1)
- y *= 2, exp--;
- else if (i == 2)
- y *= 4, exp -= 2;
-
- y = cbrt (y);
- y = scalbnq (y, exp / 3);
-
- /* Two Newton iterations. */
- y -= 0.333333333333333333333333333333333333333333333333333Q
- * (y - x / (y * y));
- y -= 0.333333333333333333333333333333333333333333333333333Q
- * (y - x / (y * y));
- return y;
-}
+ return (x);
+
+ if (x > 0)
+ sign = 1;
+ else
+ {
+ sign = -1;
+ x = -x;
+ }
+
+ z = x;
+ /* extract power of 2, leaving mantissa between 0.5 and 1 */
+ x = frexpq (x, &e);
+
+ /* Approximate cube root of number between .5 and 1,
+ peak relative error = 1.2e-6 */
+ x = ((((1.3584464340920900529734e-1L * x
+ - 6.3986917220457538402318e-1L) * x
+ + 1.2875551670318751538055e0L) * x
+ - 1.4897083391357284957891e0L) * x
+ + 1.3304961236013647092521e0L) * x + 3.7568280825958912391243e-1L;
+ /* exponent divided by 3 */
+ if (e >= 0)
+ {
+ rem = e;
+ e /= 3;
+ rem -= 3 * e;
+ if (rem == 1)
+ x *= CBRT2;
+ else if (rem == 2)
+ x *= CBRT4;
+ }
+ else
+ { /* argument less than 1 */
+ e = -e;
+ rem = e;
+ e /= 3;
+ rem -= 3 * e;
+ if (rem == 1)
+ x *= CBRT2I;
+ else if (rem == 2)
+ x *= CBRT4I;
+ e = -e;
+ }
+
+ /* multiply by power of 2 */
+ x = ldexpq (x, e);
+
+ /* Newton iteration */
+ x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
+ x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
+ x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
+
+ if (sign < 0)
+ x = -x;
+ return (x);
+}
--- /dev/null
+/* Complex cosine hyperbole function for complex __float128.
+ Copyright (C) 1997-2012 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+ 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"
+
+#ifdef HAVE_FENV_H
+# include <fenv.h>
+#endif
+
+
+__complex128
+ccoshq (__complex128 x)
+{
+ __complex128 retval;
+ int rcls = fpclassifyq (__real__ x);
+ int icls = fpclassifyq (__imag__ x);
+
+ if (__builtin_expect (rcls >= QUADFP_ZERO, 1))
+ {
+ /* Real part is finite. */
+ if (__builtin_expect (icls >= QUADFP_ZERO, 1))
+ {
+ /* Imaginary part is finite. */
+ const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q);
+ __float128 sinix, cosix;
+
+ if (__builtin_expect (icls != QUADFP_SUBNORMAL, 1))
+ {
+ sincosq (__imag__ x, &sinix, &cosix);
+ }
+ else
+ {
+ sinix = __imag__ x;
+ cosix = 1.0Q;
+ }
+
+ if (fabsq (__real__ x) > t)
+ {
+ __float128 exp_t = expq (t);
+ __float128 rx = fabsq (__real__ x);
+ if (signbitq (__real__ x))
+ sinix = -sinix;
+ rx -= t;
+ sinix *= exp_t / 2.0Q;
+ cosix *= exp_t / 2.0Q;
+ if (rx > t)
+ {
+ rx -= t;
+ sinix *= exp_t;
+ cosix *= exp_t;
+ }
+ if (rx > t)
+ {
+ /* Overflow (original real part of x > 3t). */
+ __real__ retval = FLT128_MAX * cosix;
+ __imag__ retval = FLT128_MAX * sinix;
+ }
+ else
+ {
+ __float128 exp_val = expq (rx);
+ __real__ retval = exp_val * cosix;
+ __imag__ retval = exp_val * sinix;
+ }
+ }
+ else
+ {
+ __real__ retval = coshq (__real__ x) * cosix;
+ __imag__ retval = sinhq (__real__ x) * sinix;
+ }
+ }
+ else
+ {
+ __imag__ retval = __real__ x == 0.0Q ? 0.0Q : nanq ("");
+ __real__ retval = nanq ("") + nanq ("");
+
+#ifdef HAVE_FENV_H
+ if (icls == QUADFP_INFINITE)
+ feraiseexcept (FE_INVALID);
+#endif
+ }
+ }
+ else if (rcls == QUADFP_INFINITE)
+ {
+ /* Real part is infinite. */
+ if (__builtin_expect (icls > QUADFP_ZERO, 1))
+ {
+ /* Imaginary part is finite. */
+ __float128 sinix, cosix;
+
+ if (__builtin_expect (icls != QUADFP_SUBNORMAL, 1))
+ {
+ sincosq (__imag__ x, &sinix, &cosix);
+ }
+ else
+ {
+ sinix = __imag__ x;
+ cosix = 1.0Q;
+ }
+
+ __real__ retval = copysignq (HUGE_VALQ, cosix);
+ __imag__ retval = (copysignq (HUGE_VALQ, sinix)
+ * copysignq (1.0Q, __real__ x));
+ }
+ else if (icls == QUADFP_ZERO)
+ {
+ /* Imaginary part is 0.0. */
+ __real__ retval = HUGE_VALQ;
+ __imag__ retval = __imag__ x * copysignq (1.0Q, __real__ x);
+ }
+ else
+ {
+ /* The addition raises the invalid exception. */
+ __real__ retval = HUGE_VALQ;
+ __imag__ retval = nanq ("") + nanq ("");
+
+#ifdef HAVE_FENV_H
+ if (icls == QUADFP_INFINITE)
+ feraiseexcept (FE_INVALID);
+#endif
+ }
+ }
+ else
+ {
+ __real__ retval = nanq ("");
+ __imag__ retval = __imag__ x == 0.0 ? __imag__ x : nanq ("");
+ }
+
+ return retval;
+}
-/* s_ceill.c -- long double version of s_ceil.c.
+/* ceilq.c -- __float128 version of s_ceil.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
--- /dev/null
+/* Return value of complex exponential function for complex __float128 value.
+ Copyright (C) 1997-2012 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+ 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"
+
+#ifdef HAVE_FENV_H
+# include <fenv.h>
+#endif
+
+
+__complex128
+cexpq (__complex128 x)
+{
+ __complex128 retval;
+ int rcls = fpclassifyq (__real__ x);
+ int icls = fpclassifyq (__imag__ x);
+
+ if (__builtin_expect (rcls >= QUADFP_ZERO, 1))
+ {
+ /* Real part is finite. */
+ if (__builtin_expect (icls >= QUADFP_ZERO, 1))
+ {
+ /* Imaginary part is finite. */
+ const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q);
+ __float128 sinix, cosix;
+
+ if (__builtin_expect (icls != QUADFP_SUBNORMAL, 1))
+ {
+ sincosq (__imag__ x, &sinix, &cosix);
+ }
+ else
+ {
+ sinix = __imag__ x;
+ cosix = 1.0Q;
+ }
+
+ if (__real__ x > t)
+ {
+ __float128 exp_t = expq (t);
+ __real__ x -= t;
+ sinix *= exp_t;
+ cosix *= exp_t;
+ if (__real__ x > t)
+ {
+ __real__ x -= t;
+ sinix *= exp_t;
+ cosix *= exp_t;
+ }
+ }
+ if (__real__ x > t)
+ {
+ /* Overflow (original real part of x > 3t). */
+ __real__ retval = FLT128_MAX * cosix;
+ __imag__ retval = FLT128_MAX * sinix;
+ }
+ else
+ {
+ __float128 exp_val = expq (__real__ x);
+ __real__ retval = exp_val * cosix;
+ __imag__ retval = exp_val * sinix;
+ }
+ }
+ else
+ {
+ /* If the imaginary part is +-inf or NaN and the real part
+ is not +-inf the result is NaN + iNaN. */
+ __real__ retval = nanq ("");
+ __imag__ retval = nanq ("");
+
+#ifdef HAVE_FENV_H
+ feraiseexcept (FE_INVALID);
+#endif
+ }
+ }
+ else if (__builtin_expect (rcls == QUADFP_INFINITE, 1))
+ {
+ /* Real part is infinite. */
+ if (__builtin_expect (icls >= QUADFP_ZERO, 1))
+ {
+ /* Imaginary part is finite. */
+ __float128 value = signbitq (__real__ x) ? 0.0Q : HUGE_VALQ;
+
+ if (icls == QUADFP_ZERO)
+ {
+ /* Imaginary part is 0.0. */
+ __real__ retval = value;
+ __imag__ retval = __imag__ x;
+ }
+ else
+ {
+ __float128 sinix, cosix;
+
+ if (__builtin_expect (icls != QUADFP_SUBNORMAL, 1))
+ {
+ sincosq (__imag__ x, &sinix, &cosix);
+ }
+ else
+ {
+ sinix = __imag__ x;
+ cosix = 1.0Q;
+ }
+
+ __real__ retval = copysignq (value, cosix);
+ __imag__ retval = copysignq (value, sinix);
+ }
+ }
+ else if (signbitq (__real__ x) == 0)
+ {
+ __real__ retval = HUGE_VALQ;
+ __imag__ retval = nanq ("");
+
+#ifdef HAVE_FENV_H
+ if (icls == QUADFP_INFINITE)
+ feraiseexcept (FE_INVALID);
+#endif
+ }
+ else
+ {
+ __real__ retval = 0.0Q;
+ __imag__ retval = copysignq (0.0Q, __imag__ x);
+ }
+ }
+ else
+ {
+ /* If the real part is NaN the result is NaN + iNaN. */
+ __real__ retval = nanq ("");
+ __imag__ retval = nanq ("");
+
+#ifdef HAVE_FENV_H
+ if (rcls != QUADFP_NAN || icls != QUADFP_NAN)
+ feraiseexcept (FE_INVALID);
+#endif
+ }
+
+ return retval;
+}
--- /dev/null
+/* Compute complex base 10 logarithm for complex __float128.
+ Copyright (C) 1997-2012 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+ 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"
+
+
+/* log_10 (2). */
+#define M_LOG10_2q 0.3010299956639811952137388947244930267682Q
+
+
+__complex128
+clog10q (__complex128 x)
+{
+ __complex128 result;
+ int rcls = fpclassifyq (__real__ x);
+ int icls = fpclassifyq (__imag__ x);
+
+ if (__builtin_expect (rcls == QUADFP_ZERO && icls == QUADFP_ZERO, 0))
+ {
+ /* Real and imaginary part are 0.0. */
+ __imag__ result = signbitq (__real__ x) ? M_PIq : 0.0Q;
+ __imag__ result = copysignq (__imag__ result, __imag__ x);
+ /* Yes, the following line raises an exception. */
+ __real__ result = -1.0Q / fabsq (__real__ x);
+ }
+ else if (__builtin_expect (rcls != QUADFP_NAN && icls != QUADFP_NAN, 1))
+ {
+ /* Neither real nor imaginary part is NaN. */
+ __float128 absx = fabsq (__real__ x), absy = fabsq (__imag__ x);
+ int scale = 0;
+
+ if (absx < absy)
+ {
+ __float128 t = absx;
+ absx = absy;
+ absy = t;
+ }
+
+ if (absx > FLT128_MAX / 2.0Q)
+ {
+ scale = -1;
+ absx = scalbnq (absx, scale);
+ absy = (absy >= FLT128_MIN * 2.0Q ? scalbnq (absy, scale) : 0.0Q);
+ }
+ else if (absx < FLT128_MIN && absy < FLT128_MIN)
+ {
+ scale = FLT128_MANT_DIG;
+ absx = scalbnq (absx, scale);
+ absy = scalbnq (absy, scale);
+ }
+
+ if (absx == 1.0Q && scale == 0)
+ {
+ __float128 absy2 = absy * absy;
+ if (absy2 <= FLT128_MIN * 2.0Q * M_LN10q)
+ __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)
+ {
+ __float128 d2m1 = (absx - 1.0Q) * (absx + 1.0Q);
+ if (absy >= FLT128_EPSILON)
+ d2m1 += absy * absy;
+ __real__ result = log1pq (d2m1) * (M_LOG10Eq / 2.0Q);
+ }
+ else if (absx < 1.0Q
+ && absx >= 0.75Q
+ && absy < FLT128_EPSILON / 2.0Q
+ && scale == 0)
+ {
+ __float128 d2m1 = (absx - 1.0Q) * (absx + 1.0Q);
+ __real__ result = log1pq (d2m1) * (M_LOG10Eq / 2.0Q);
+ }
+ else if (absx < 1.0Q && (absx >= 0.75Q || absy >= 0.5Q) && scale == 0)
+ {
+ __float128 d2m1 = __quadmath_x2y2m1q (absx, absy);
+ __real__ result = log1pq (d2m1) * (M_LOG10Eq / 2.0Q);
+ }
+ else
+ {
+ __float128 d = hypotq (absx, absy);
+ __real__ result = log10q (d) - scale * M_LOG10_2q;
+ }
+
+ __imag__ result = M_LOG10Eq * atan2q (__imag__ x, __real__ x);
+ }
+ else
+ {
+ __imag__ result = nanq ("");
+ if (rcls == QUADFP_INFINITE || icls == QUADFP_INFINITE)
+ /* Real or imaginary part is infinite. */
+ __real__ result = HUGE_VALQ;
+ else
+ __real__ result = nanq ("");
+ }
+
+ return result;
+}
--- /dev/null
+/* Compute complex natural logarithm for complex __float128.
+ Copyright (C) 1997-2012 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+ 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"
+
+
+__complex128
+clogq (__complex128 x)
+{
+ __complex128 result;
+ int rcls = fpclassifyq (__real__ x);
+ int icls = fpclassifyq (__imag__ x);
+
+ if (__builtin_expect (rcls == QUADFP_ZERO && icls == QUADFP_ZERO, 0))
+ {
+ /* Real and imaginary part are 0.0. */
+ __imag__ result = signbitq (__real__ x) ? M_PIq : 0.0Q;
+ __imag__ result = copysignq (__imag__ result, __imag__ x);
+ /* Yes, the following line raises an exception. */
+ __real__ result = -1.0Q / fabsq (__real__ x);
+ }
+ else if (__builtin_expect (rcls != QUADFP_NAN && icls != QUADFP_NAN, 1))
+ {
+ /* Neither real nor imaginary part is NaN. */
+ __float128 absx = fabsq (__real__ x), absy = fabsq (__imag__ x);
+ int scale = 0;
+
+ if (absx < absy)
+ {
+ __float128 t = absx;
+ absx = absy;
+ absy = t;
+ }
+
+ if (absx > FLT128_MAX / 2.0)
+ {
+ scale = -1;
+ absx = scalbnq (absx, scale);
+ absy = (absy >= FLT128_MIN * 2.0Q ? scalbnq (absy, scale) : 0.0Q);
+ }
+ else if (absx < FLT128_MIN && absy < FLT128_MIN)
+ {
+ scale = FLT128_MANT_DIG;
+ absx = scalbnq (absx, scale);
+ absy = scalbnq (absy, scale);
+ }
+
+ if (absx == 1.0Q && scale == 0)
+ {
+ __float128 absy2 = absy * absy;
+ if (absy2 <= FLT128_MIN * 2.0Q)
+ __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)
+ {
+ __float128 d2m1 = (absx - 1.0Q) * (absx + 1.0Q);
+ if (absy >= FLT128_EPSILON)
+ d2m1 += absy * absy;
+ __real__ result = log1pq (d2m1) / 2.0Q;
+ }
+ else if (absx < 1.0Q
+ && absx >= 0.75Q
+ && absy < FLT128_EPSILON / 2.0Q
+ && scale == 0)
+ {
+ __float128 d2m1 = (absx - 1.0Q) * (absx + 1.0Q);
+ __real__ result = log1pq (d2m1) / 2.0Q;
+ }
+ else if (absx < 1.0 && (absx >= 0.75Q || absy >= 0.5Q) && scale == 0)
+ {
+ __float128 d2m1 = __quadmath_x2y2m1q (absx, absy);
+ __real__ result = log1pq (d2m1) / 2.0Q;
+ }
+ else
+ {
+ __float128 d = hypotq (absx, absy);
+ __real__ result = logq (d) - scale * M_LN2q;
+ }
+
+ __imag__ result = atan2q (__imag__ x, __real__ x);
+ }
+ else
+ {
+ __imag__ result = nanq ("");
+ if (rcls == QUADFP_INFINITE || icls == QUADFP_INFINITE)
+ /* Real or imaginary part is infinite. */
+ __real__ result = HUGE_VALQ;
+ else
+ __real__ result = nanq ("");
+ }
+
+ return result;
+}
+/* GCC Quad-Precision Math Library
+ Copyright (C) 2010, 2011 Free Software Foundation, Inc.
+ Written by Francois-Xavier Coudert <fxcoudert@gcc.gnu.org>
+
+This file is part of the libquadmath library.
+Libquadmath is free software; you can redistribute it and/or
+modify it under the terms of the GNU Library General Public
+License as published by the Free Software Foundation; either
+version 2 of the License, or (at your option) any later version.
+
+Libquadmath 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
+Library General Public License for more details.
+
+You should have received a copy of the GNU Library General Public
+License along with libquadmath; see the file COPYING.LIB. If
+not, write to the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
+Boston, MA 02110-1301, USA. */
+
#include "quadmath-imp.h"
+#ifdef HAVE_FENV_H
+# include <fenv.h>
+#endif
+
#define REALPART(z) (__real__(z))
#define IMAGPART(z) (__imag__(z))
#define COMPLEX_ASSIGN(z_, r_, i_) {__real__(z_) = (r_); __imag__(z_) = (i_);}
-// Horrible... GCC doesn't know how to multiply or divide these
-// __complex128 things. We have to do it on our own.
-// Protect it around macros so, some day, we can switch it on
-
-#if 0
-
-# define C128_MULT(x,y) ((x)*(y))
-# define C128_DIV(x,y) ((x)/(y))
-
-#else
-
-#define C128_MULT(x,y) mult_c128(x,y)
-#define C128_DIV(x,y) div_c128(x,y)
-
-static inline __complex128 mult_c128 (__complex128 x, __complex128 y)
-{
- __float128 r1 = REALPART(x), i1 = IMAGPART(x);
- __float128 r2 = REALPART(y), i2 = IMAGPART(y);
- __complex128 res;
- COMPLEX_ASSIGN(res, r1*r2 - i1*i2, i2*r1 + i1*r2);
- return res;
-}
-
-
-// Careful: the algorithm for the division sucks. A lot.
-static inline __complex128 div_c128 (__complex128 x, __complex128 y)
-{
- __float128 n = hypotq (REALPART (y), IMAGPART (y));
- __float128 r1 = REALPART(x), i1 = IMAGPART(x);
- __float128 r2 = REALPART(y), i2 = IMAGPART(y);
- __complex128 res;
- COMPLEX_ASSIGN(res, r1*r2 + i1*i2, i1*r2 - i2*r1);
- return res / n;
-}
-
-#endif
-
-
-
__float128
cabsq (__complex128 z)
{
}
-__complex128
-cexpq (__complex128 z)
-{
- __float128 a, b;
- __complex128 v;
-
- a = REALPART (z);
- b = IMAGPART (z);
- COMPLEX_ASSIGN (v, cosq (b), sinq (b));
- return expq (a) * v;
-}
-
-
__complex128
cexpiq (__float128 x)
{
+ __float128 sinix, cosix;
__complex128 v;
- COMPLEX_ASSIGN (v, cosq (x), sinq (x));
+ sincosq (x, &sinix, &cosix);
+ COMPLEX_ASSIGN (v, cosix, sinix);
return v;
}
}
-__complex128
-clogq (__complex128 z)
-{
- __complex128 v;
- COMPLEX_ASSIGN (v, logq (cabsq (z)), cargq (z));
- return v;
-}
-
-
-__complex128
-clog10q (__complex128 z)
-{
- __complex128 v;
- COMPLEX_ASSIGN (v, log10q (cabsq (z)), cargq (z));
- return v;
-}
-
-
__complex128
cpowq (__complex128 base, __complex128 power)
{
- return cexpq (C128_MULT(power, clogq (base)));
+ return cexpq (power * clogq (base));
}
__complex128
-csinq (__complex128 a)
+ccosq (__complex128 x)
{
- __float128 r = REALPART (a), i = IMAGPART (a);
- __complex128 v;
- COMPLEX_ASSIGN (v, sinq (r) * coshq (i), cosq (r) * sinhq (i));
- return v;
-}
-
-
-__complex128
-csinhq (__complex128 a)
-{
- __float128 r = REALPART (a), i = IMAGPART (a);
- __complex128 v;
- COMPLEX_ASSIGN (v, sinhq (r) * cosq (i), coshq (r) * sinq (i));
- return v;
-}
-
-
-__complex128
-ccosq (__complex128 a)
-{
- __float128 r = REALPART (a), i = IMAGPART (a);
- __complex128 v;
- COMPLEX_ASSIGN (v, cosq (r) * coshq (i), - (sinq (r) * sinhq (i)));
- return v;
-}
-
-
-__complex128
-ccoshq (__complex128 a)
-{
- __float128 r = REALPART (a), i = IMAGPART (a);
- __complex128 v;
- COMPLEX_ASSIGN (v, coshq (r) * cosq (i), sinhq (r) * sinq (i));
- return v;
-}
-
-
-__complex128
-ctanq (__complex128 a)
-{
- __float128 rt = tanq (REALPART (a)), it = tanhq (IMAGPART (a));
- __complex128 n, d;
- COMPLEX_ASSIGN (n, rt, it);
- COMPLEX_ASSIGN (d, 1, - (rt * it));
- return C128_DIV(n,d);
-}
-
-
-__complex128
-ctanhq (__complex128 a)
-{
- __float128 rt = tanhq (REALPART (a)), it = tanq (IMAGPART (a));
- __complex128 n, d;
- COMPLEX_ASSIGN (n, rt, it);
- COMPLEX_ASSIGN (d, 1, rt * it);
- return C128_DIV(n,d);
-}
-
-
-/* Square root algorithm from glibc. */
-__complex128
-csqrtq (__complex128 z)
-{
- __float128 re = REALPART(z), im = IMAGPART(z);
- __complex128 v;
+ __complex128 y;
- if (im == 0)
- {
- if (isnanq (re))
- {
- COMPLEX_ASSIGN (v, -re, -re);
- }
- else if (re < 0)
- {
- COMPLEX_ASSIGN (v, 0, copysignq (sqrtq (-re), im));
- }
- else
- {
- COMPLEX_ASSIGN (v, fabsq (sqrtq (re)), copysignq (0, im));
- }
- }
- else if (isinfq (im))
- {
- COMPLEX_ASSIGN (v, fabsq (im), im);
- }
- else if (re == 0)
- {
- __float128 r = sqrtq (0.5 * fabsq (im));
- COMPLEX_ASSIGN (v, r, copysignq (r, im));
- }
- else
- {
- __float128 d = hypotq (re, im);
- __float128 r, s;
-
- /* Use the identity 2 Re res Im res = Im x
- to avoid cancellation error in d +/- Re x. */
- if (re > 0)
- r = sqrtq (0.5 * d + 0.5 * re), s = (0.5 * im) / r;
- else
- s = sqrtq (0.5 * d - 0.5 * re), r = fabsq ((0.5 * im) / s);
-
- COMPLEX_ASSIGN (v, r, copysignq (s, im));
- }
- return v;
+ COMPLEX_ASSIGN (y, -IMAGPART (x), REALPART (x));
+ return ccoshq (y);
}
-
-/* s_copysignl.c -- long double version of s_copysign.c.
+/* copysignq.c -- __float128 version of s_copysign.c.
* Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
-/* __ieee754_coshl(x)
+/* coshq(x)
* Method :
- * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
- * 1. Replace x by |x| (coshl(x) = coshl(-x)).
+ * mathematically coshq(x) if defined to be (exp(x)+exp(-x))/2
+ * 1. Replace x by |x| (coshq(x) = coshq(-x)).
* 2.
* [ exp(x) - 1 ]^2
- * 0 <= x <= ln2/2 : coshl(x) := 1 + -------------------
+ * 0 <= x <= ln2/2 : coshq(x) := 1 + -------------------
* 2*exp(x)
*
* exp(x) + 1/exp(x)
- * ln2/2 <= x <= 22 : coshl(x) := -------------------
+ * ln2/2 <= x <= 22 : coshq(x) := -------------------
* 2
- * 22 <= x <= lnovft : coshl(x) := expl(x)/2
- * lnovft <= x <= ln2ovft: coshl(x) := expl(x/2)/2 * expl(x/2)
- * ln2ovft < x : coshl(x) := huge*huge (overflow)
+ * 22 <= x <= lnovft : coshq(x) := expq(x)/2
+ * lnovft <= x <= ln2ovft: coshq(x) := expq(x/2)/2 * expq(x/2)
+ * ln2ovft < x : coshq(x) := huge*huge (overflow)
*
* Special cases:
- * coshl(x) is |x| if x is +INF, -INF, or NaN.
- * only coshl(0)=1 is exact for finite x.
+ * coshq(x) is |x| if x is +INF, -INF, or NaN.
+ * only coshq(0)=1 is exact for finite x.
*/
#include "quadmath-imp.h"
if (ex >= 0x7fff0000)
return x * x;
- /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
+ /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expq(|x|)) */
if (ex < 0x3ffd62e4) /* 0.3465728759765625 */
{
t = expm1q (u.value);
-/* s_cosl.c -- long double version of s_cos.c.
+/* cosq.c -- __float128 version of s_cos.c.
* Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
* ====================================================
*/
-/* cosl(x)
+/* cosq(x)
* Return cosine function of x.
*
* kernel function:
- * __kernel_sinl ... sine function on [-pi/4,pi/4]
- * __kernel_cosl ... cosine function on [-pi/4,pi/4]
- * __ieee754_rem_pio2l ... argument reduction routine
+ * __quadmath_kernel_sinq ... sine function on [-pi/4,pi/4]
+ * __quadmath_kernel_cosq ... cosine function on [-pi/4,pi/4]
+ * __quadmath_rem_pio2q ... argument reduction routine
*
* Method.
* Let S,C and T denote the sin, cos and tan respectively on
{
/* 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
- possible values for h. We look up cosl(h) and sinl(h) in
- pre-computed tables, compute cosl(l) and sinl(l) using a
+ possible values for h. We look up cosq(h) and sinq(h) in
+ pre-computed tables, compute cosq(l) and sinq(l) using a
Chebyshev polynomial of degree 10(11) and compute
- cosl(h+l) = cosl(h)cosl(l) - sinl(h)sinl(l). */
+ cosq(h+l) = cosq(h)cosq(l) - sinq(h)sinq(l). */
index = 0x3ffe - (tix >> 16);
hix = (tix + (0x200 << index)) & (0xfffffc00 << index);
- x = fabsq (x);
+ if (signbitq (x))
+ {
+ x = -x;
+ y = -y;
+ }
switch (index)
{
case 0: index = ((45 << 10) + hix - 0x3ffe0000) >> 8; break;
--- /dev/null
+/* Complex sine hyperbole function for complex __float128.
+ Copyright (C) 1997-2012 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+ 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"
+
+#ifdef HAVE_FENV_H
+# include <fenv.h>
+#endif
+
+
+__complex128
+csinhq (__complex128 x)
+{
+ __complex128 retval;
+ int negate = signbitq (__real__ x);
+ int rcls = fpclassifyq (__real__ x);
+ int icls = fpclassifyq (__imag__ x);
+
+ __real__ x = fabsq (__real__ x);
+
+ if (__builtin_expect (rcls >= QUADFP_ZERO, 1))
+ {
+ /* Real part is finite. */
+ if (__builtin_expect (icls >= QUADFP_ZERO, 1))
+ {
+ /* Imaginary part is finite. */
+ const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q);
+ __float128 sinix, cosix;
+
+ if (__builtin_expect (icls != QUADFP_SUBNORMAL, 1))
+ {
+ sincosq (__imag__ x, &sinix, &cosix);
+ }
+ else
+ {
+ sinix = __imag__ x;
+ cosix = 1.0Q;
+ }
+
+ if (fabsq (__real__ x) > t)
+ {
+ __float128 exp_t = expq (t);
+ __float128 rx = fabsq (__real__ x);
+ if (signbitq (__real__ x))
+ cosix = -cosix;
+ rx -= t;
+ sinix *= exp_t / 2.0Q;
+ cosix *= exp_t / 2.0Q;
+ if (rx > t)
+ {
+ rx -= t;
+ sinix *= exp_t;
+ cosix *= exp_t;
+ }
+ if (rx > t)
+ {
+ /* Overflow (original real part of x > 3t). */
+ __real__ retval = FLT128_MAX * cosix;
+ __imag__ retval = FLT128_MAX * sinix;
+ }
+ else
+ {
+ __float128 exp_val = expq (rx);
+ __real__ retval = exp_val * cosix;
+ __imag__ retval = exp_val * sinix;
+ }
+ }
+ else
+ {
+ __real__ retval = sinhq (__real__ x) * cosix;
+ __imag__ retval = coshq (__real__ x) * sinix;
+ }
+
+ if (negate)
+ __real__ retval = -__real__ retval;
+ }
+ else
+ {
+ if (rcls == QUADFP_ZERO)
+ {
+ /* Real part is 0.0. */
+ __real__ retval = copysignq (0.0Q, negate ? -1.0Q : 1.0Q);
+ __imag__ retval = nanq ("") + nanq ("");
+
+#ifdef HAVE_FENV_H
+ if (icls == QUADFP_INFINITE)
+ feraiseexcept (FE_INVALID);
+#endif
+ }
+ else
+ {
+ __real__ retval = nanq ("");
+ __imag__ retval = nanq ("");
+
+#ifdef HAVE_FENV_H
+ feraiseexcept (FE_INVALID);
+#endif
+ }
+ }
+ }
+ else if (rcls == QUADFP_INFINITE)
+ {
+ /* Real part is infinite. */
+ if (__builtin_expect (icls > QUADFP_ZERO, 1))
+ {
+ /* Imaginary part is finite. */
+ __float128 sinix, cosix;
+
+ if (__builtin_expect (icls != QUADFP_SUBNORMAL, 1))
+ {
+ sincosq (__imag__ x, &sinix, &cosix);
+ }
+ else
+ {
+ sinix = __imag__ x;
+ cosix = 1.0;
+ }
+
+ __real__ retval = copysignq (HUGE_VALQ, cosix);
+ __imag__ retval = copysignq (HUGE_VALQ, sinix);
+
+ if (negate)
+ __real__ retval = -__real__ retval;
+ }
+ else if (icls == QUADFP_ZERO)
+ {
+ /* Imaginary part is 0.0. */
+ __real__ retval = negate ? -HUGE_VALQ : HUGE_VALQ;
+ __imag__ retval = __imag__ x;
+ }
+ else
+ {
+ /* The addition raises the invalid exception. */
+ __real__ retval = HUGE_VALQ;
+ __imag__ retval = nanq ("") + nanq ("");
+
+#ifdef HAVE_FENV_H
+ if (icls == QUADFP_INFINITE)
+ feraiseexcept (FE_INVALID);
+#endif
+ }
+ }
+ else
+ {
+ __real__ retval = nanq ("");
+ __imag__ retval = __imag__ x == 0.0Q ? __imag__ x : nanq ("");
+ }
+
+ return retval;
+}
--- /dev/null
+/* Complex sine function for complex __float128.
+ Copyright (C) 1997-2012 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+ 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"
+
+#ifdef HAVE_FENV_H
+# include <fenv.h>
+#endif
+
+
+__complex128
+csinq (__complex128 x)
+{
+ __complex128 retval;
+ int negate = signbitq (__real__ x);
+ int rcls = fpclassifyq (__real__ x);
+ int icls = fpclassifyq (__imag__ x);
+
+ __real__ x = fabsq (__real__ x);
+
+ if (__builtin_expect (icls >= QUADFP_ZERO, 1))
+ {
+ /* Imaginary part is finite. */
+ if (__builtin_expect (rcls >= QUADFP_ZERO, 1))
+ {
+ /* Real part is finite. */
+ const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q);
+ __float128 sinix, cosix;
+
+ if (__builtin_expect (rcls != QUADFP_SUBNORMAL, 1))
+ {
+ sincosq (__real__ x, &sinix, &cosix);
+ }
+ else
+ {
+ sinix = __real__ x;
+ cosix = 1.0Q;
+ }
+
+ if (fabsq (__imag__ x) > t)
+ {
+ __float128 exp_t = expq (t);
+ __float128 ix = fabsq (__imag__ x);
+ if (signbitq (__imag__ x))
+ cosix = -cosix;
+ ix -= t;
+ sinix *= exp_t / 2.0Q;
+ cosix *= exp_t / 2.0Q;
+ if (ix > t)
+ {
+ ix -= t;
+ sinix *= exp_t;
+ cosix *= exp_t;
+ }
+ if (ix > t)
+ {
+ /* Overflow (original imaginary part of x > 3t). */
+ __real__ retval = FLT128_MAX * sinix;
+ __imag__ retval = FLT128_MAX * cosix;
+ }
+ else
+ {
+ __float128 exp_val = expq (ix);
+ __real__ retval = exp_val * sinix;
+ __imag__ retval = exp_val * cosix;
+ }
+ }
+ else
+ {
+ __real__ retval = coshq (__imag__ x) * sinix;
+ __imag__ retval = sinhq (__imag__ x) * cosix;
+ }
+
+ if (negate)
+ __real__ retval = -__real__ retval;
+ }
+ else
+ {
+ if (icls == QUADFP_ZERO)
+ {
+ /* Imaginary part is 0.0. */
+ __real__ retval = nanq ("");
+ __imag__ retval = __imag__ x;
+
+#ifdef HAVE_FENV_H
+ if (rcls == QUADFP_INFINITE)
+ feraiseexcept (FE_INVALID);
+#endif
+ }
+ else
+ {
+ __real__ retval = nanq ("");
+ __imag__ retval = nanq ("");
+
+#ifdef HAVE_FENV_H
+ feraiseexcept (FE_INVALID);
+#endif
+ }
+ }
+ }
+ else if (icls == QUADFP_INFINITE)
+ {
+ /* Imaginary part is infinite. */
+ if (rcls == QUADFP_ZERO)
+ {
+ /* Real part is 0.0. */
+ __real__ retval = copysignq (0.0Q, negate ? -1.0Q : 1.0Q);
+ __imag__ retval = __imag__ x;
+ }
+ else if (rcls > QUADFP_ZERO)
+ {
+ /* Real part is finite. */
+ __float128 sinix, cosix;
+
+ if (__builtin_expect (rcls != QUADFP_SUBNORMAL, 1))
+ {
+ sincosq (__real__ x, &sinix, &cosix);
+ }
+ else
+ {
+ sinix = __real__ x;
+ cosix = 1.0;
+ }
+
+ __real__ retval = copysignq (HUGE_VALQ, sinix);
+ __imag__ retval = copysignq (HUGE_VALQ, cosix);
+
+ if (negate)
+ __real__ retval = -__real__ retval;
+ if (signbitq (__imag__ x))
+ __imag__ retval = -__imag__ retval;
+ }
+ else
+ {
+ /* The addition raises the invalid exception. */
+ __real__ retval = nanq ("");
+ __imag__ retval = HUGE_VALQ;
+
+#ifdef HAVE_FENV_H
+ if (rcls == QUADFP_INFINITE)
+ feraiseexcept (FE_INVALID);
+#endif
+ }
+ }
+ else
+ {
+ if (rcls == QUADFP_ZERO)
+ __real__ retval = copysignq (0.0Q, negate ? -1.0Q : 1.0Q);
+ else
+ __real__ retval = nanq ("");
+ __imag__ retval = nanq ("");
+ }
+
+ return retval;
+}
--- /dev/null
+/* Complex square root of __float128 value.
+ Copyright (C) 1997-2012 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Based on an algorithm by Stephen L. Moshier <moshier@world.std.com>.
+ Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+ 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"
+
+#ifdef HAVE_FENV_H
+# include <fenv.h>
+#endif
+
+
+__complex128
+csqrtq (__complex128 x)
+{
+ __complex128 res;
+ int rcls = fpclassifyq (__real__ x);
+ int icls = fpclassifyq (__imag__ x);
+
+ if (__builtin_expect (rcls <= QUADFP_INFINITE || icls <= QUADFP_INFINITE, 0))
+ {
+ if (icls == QUADFP_INFINITE)
+ {
+ __real__ res = HUGE_VALQ;
+ __imag__ res = __imag__ x;
+ }
+ else if (rcls == QUADFP_INFINITE)
+ {
+ if (__real__ x < 0.0Q)
+ {
+ __real__ res = icls == QUADFP_NAN ? nanq ("") : 0;
+ __imag__ res = copysignq (HUGE_VALQ, __imag__ x);
+ }
+ else
+ {
+ __real__ res = __real__ x;
+ __imag__ res = (icls == QUADFP_NAN
+ ? nanq ("") : copysignq (0.0Q, __imag__ x));
+ }
+ }
+ else
+ {
+ __real__ res = nanq ("");
+ __imag__ res = nanq ("");
+ }
+ }
+ else
+ {
+ if (__builtin_expect (icls == QUADFP_ZERO, 0))
+ {
+ if (__real__ x < 0.0Q)
+ {
+ __real__ res = 0.0Q;
+ __imag__ res = copysignq (sqrtq (-__real__ x),
+ __imag__ x);
+ }
+ else
+ {
+ __real__ res = fabsq (sqrtq (__real__ x));
+ __imag__ res = copysignq (0.0Q, __imag__ x);
+ }
+ }
+ else if (__builtin_expect (rcls == QUADFP_ZERO, 0))
+ {
+ __float128 r;
+ if (fabsq (__imag__ x) >= 2.0Q * FLT128_MIN)
+ r = sqrtq (0.5Q * fabsq (__imag__ x));
+ else
+ r = 0.5Q * sqrtq (2.0Q * fabsq (__imag__ x));
+
+ __real__ res = r;
+ __imag__ res = copysignq (r, __imag__ x);
+ }
+ else
+ {
+ __float128 d, r, s;
+ int scale = 0;
+
+ if (fabsq (__real__ x) > FLT128_MAX / 4.0Q)
+ {
+ scale = 1;
+ __real__ x = scalbnq (__real__ x, -2 * scale);
+ __imag__ x = scalbnq (__imag__ x, -2 * scale);
+ }
+ else if (fabsq (__imag__ x) > FLT128_MAX / 4.0Q)
+ {
+ scale = 1;
+ if (fabsq (__real__ x) >= 4.0Q * FLT128_MIN)
+ __real__ x = scalbnq (__real__ x, -2 * scale);
+ else
+ __real__ x = 0.0Q;
+ __imag__ x = scalbnq (__imag__ x, -2 * scale);
+ }
+ else if (fabsq (__real__ x) < FLT128_MIN
+ && fabsq (__imag__ x) < FLT128_MIN)
+ {
+ scale = -(FLT128_MANT_DIG / 2);
+ __real__ x = scalbnq (__real__ x, -2 * scale);
+ __imag__ x = scalbnq (__imag__ x, -2 * scale);
+ }
+
+ d = hypotq (__real__ x, __imag__ x);
+ /* Use the identity 2 Re res Im res = Im x
+ to avoid cancellation error in d +/- Re x. */
+ if (__real__ x > 0)
+ {
+ r = sqrtq (0.5Q * (d + __real__ x));
+ s = 0.5Q * (__imag__ x / r);
+ }
+ else
+ {
+ s = sqrtq (0.5Q * (d - __real__ x));
+ r = fabsq (0.5Q * (__imag__ x / s));
+ }
+
+ if (scale)
+ {
+ r = scalbnq (r, scale);
+ s = scalbnq (s, scale);
+ }
+
+ __real__ res = r;
+ __imag__ res = copysignq (s, __imag__ x);
+ }
+ }
+
+ return res;
+}
--- /dev/null
+/* Complex hyperbole tangent for __float128.
+ Copyright (C) 1997-2012 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+ 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"
+
+#ifdef HAVE_FENV_H
+# include <fenv.h>
+#endif
+
+
+__complex128
+ctanhq (__complex128 x)
+{
+ __complex128 res;
+
+ if (__builtin_expect (!finiteq (__real__ x) || !finiteq (__imag__ x), 0))
+ {
+ if (__quadmath_isinf_nsq (__real__ x))
+ {
+ __real__ res = copysignq (1.0Q, __real__ x);
+ __imag__ res = copysignq (0.0Q, __imag__ x);
+ }
+ else if (__imag__ x == 0.0Q)
+ {
+ res = x;
+ }
+ else
+ {
+ __real__ res = nanq ("");
+ __imag__ res = nanq ("");
+
+#ifdef HAVE_FENV_H
+ if (__quadmath_isinf_nsq (__imag__ x))
+ feraiseexcept (FE_INVALID);
+#endif
+ }
+ }
+ else
+ {
+ __float128 sinix, cosix;
+ __float128 den;
+ 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))
+ = (sinh(x)*cosh(x) + i*sin(y)*cos(y))/(sinh(x)^2 + cos(y)^2). */
+
+ if (__builtin_expect (icls != QUADFP_SUBNORMAL, 1))
+ {
+ sincosq (__imag__ x, &sinix, &cosix);
+ }
+ else
+ {
+ sinix = __imag__ x;
+ cosix = 1.0Q;
+ }
+
+ if (fabsq (__real__ x) > t)
+ {
+ /* Avoid intermediate overflow when the imaginary part of
+ the result may be subnormal. Ignoring negligible terms,
+ the real part is +/- 1, the imaginary part is
+ sin(y)*cos(y)/sinh(x)^2 = 4*sin(y)*cos(y)/exp(2x). */
+ __float128 exp_2t = expq (2 * t);
+
+ __real__ res = copysignq (1.0, __real__ x);
+ __imag__ res = 4 * sinix * cosix;
+ __real__ x = fabsq (__real__ x);
+ __real__ x -= t;
+ __imag__ res /= exp_2t;
+ if (__real__ x > t)
+ {
+ /* Underflow (original real part of x has absolute value
+ > 2t). */
+ __imag__ res /= exp_2t;
+ }
+ else
+ __imag__ res /= expq (2 * __real__ x);
+ }
+ else
+ {
+ __float128 sinhrx, coshrx;
+ if (fabsq (__real__ x) > FLT128_MIN)
+ {
+ sinhrx = sinhq (__real__ x);
+ coshrx = coshq (__real__ x);
+ }
+ else
+ {
+ sinhrx = __real__ x;
+ coshrx = 1.0Q;
+ }
+
+ if (fabsq (sinhrx) > fabsq (cosix) * FLT128_EPSILON)
+ den = sinhrx * sinhrx + cosix * cosix;
+ else
+ den = cosix * cosix;
+ __real__ res = sinhrx * coshrx / den;
+ __imag__ res = sinix * cosix / den;
+ }
+ }
+
+ return res;
+}
--- /dev/null
+/* Complex tangent function for complex __float128.
+ Copyright (C) 1997-2012 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+ 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"
+
+#ifdef HAVE_FENV_H
+# include <fenv.h>
+#endif
+
+
+__complex128
+ctanq (__complex128 x)
+{
+ __complex128 res;
+
+ if (__builtin_expect (!finiteq (__real__ x) || !finiteq (__imag__ x), 0))
+ {
+ if (__quadmath_isinf_nsq (__imag__ x))
+ {
+ __real__ res = copysignq (0.0Q, __real__ x);
+ __imag__ res = copysignq (1.0Q, __imag__ x);
+ }
+ else if (__real__ x == 0.0Q)
+ {
+ res = x;
+ }
+ else
+ {
+ __real__ res = nanq ("");
+ __imag__ res = nanq ("");
+
+#ifdef HAVE_FENV_H
+ if (__quadmath_isinf_nsq (__real__ x))
+ feraiseexcept (FE_INVALID);
+#endif
+ }
+ }
+ else
+ {
+ __float128 sinrx, cosrx;
+ __float128 den;
+ const int t = (int) ((FLT128_MAX_EXP - 1) * M_LN2q / 2.0Q);
+ int rcls = fpclassifyq (__real__ x);
+
+ /* 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). */
+
+ if (__builtin_expect (rcls != QUADFP_SUBNORMAL, 1))
+ {
+ sincosq (__real__ x, &sinrx, &cosrx);
+ }
+ else
+ {
+ sinrx = __real__ x;
+ cosrx = 1.0Q;
+ }
+
+ if (fabsq (__imag__ x) > t)
+ {
+ /* Avoid intermediate overflow when the real part of the
+ result may be subnormal. Ignoring negligible terms, the
+ imaginary part is +/- 1, the real part is
+ sin(x)*cos(x)/sinh(y)^2 = 4*sin(x)*cos(x)/exp(2y). */
+ __float128 exp_2t = expq (2 * t);
+
+ __imag__ res = copysignq (1.0Q, __imag__ x);
+ __real__ res = 4 * sinrx * cosrx;
+ __imag__ x = fabsq (__imag__ x);
+ __imag__ x -= t;
+ __real__ res /= exp_2t;
+ if (__imag__ x > t)
+ {
+ /* Underflow (original imaginary part of x has absolute
+ value > 2t). */
+ __real__ res /= exp_2t;
+ }
+ else
+ __real__ res /= expq (2 * __imag__ x);
+ }
+ else
+ {
+ __float128 sinhix, coshix;
+ if (fabsq (__imag__ x) > FLT128_MIN)
+ {
+ sinhix = sinhq (__imag__ x);
+ coshix = coshq (__imag__ x);
+ }
+ else
+ {
+ sinhix = __imag__ x;
+ coshix = 1.0Q;
+ }
+
+ if (fabsq (sinhix) > fabsq (cosrx) * FLT128_EPSILON)
+ den = cosrx * cosrx + sinhix * sinhix;
+ else
+ den = cosrx * cosrx;
+ __real__ res = sinrx * cosrx / den;
+ __imag__ res = sinhix * coshix / den;
+ }
+ }
+
+ return res;
+}
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
-/* double erf(double x)
- * double erfc(double x)
+/* __float128 erfq(__float128 x)
+ * __float128 erfcq(__float128 x)
* x
* 2 |\
* erf(x) = --------- | exp(-t*t)dt
#endif
-/* __expl_table basically consists of four tables, T_EXPL_ARG{1,2} and
+/* __expq_table basically consists of four tables, T_EXPL_ARG{1,2} and
T_EXPL_RES{1,2}. All tables use positive and negative indexes, the 0 points
are marked by T_EXPL_* defines.
For ARG1 and RES1 tables lets B be 89 and S 256.0, for ARG2 and RES2 B is 65
and S 32768.0.
These table have the property that, for all integers -B <= i <= B
- expl(__expl_table[T_EXPL_ARGN+2*i]+__expl_table[T_EXPL_ARGN+2*i+1]+r) ==
- __expl_table[T_EXPL_RESN+i], __expl_table[T_EXPL_RESN+i] is some exact number
+ expl(__expq_table[T_EXPL_ARGN+2*i]+__expq_table[T_EXPL_ARGN+2*i+1]+r) ==
+ __expq_table[T_EXPL_RESN+i], __expq_table[T_EXPL_RESN+i] is some exact number
with the low 58 bits of the mantissa 0,
- __expl_table[T_EXPL_ARGN+2*i] == i/S+s
+ __expq_table[T_EXPL_ARGN+2*i] == i/S+s
where absl(s) <= 2^-54 and absl(r) <= 2^-212. */
-static const __float128 __expl_table [] = {
+static const __float128 __expq_table [] = {
-3.47656250000000000584188889839535373E-01Q, /* bffd640000000000002b1b04213cf000 */
6.90417668990715641167244540876988960E-32Q, /* 3f97667c3fdb588a6ae1af8748357a17 */
-3.43749999999999981853132895957607418E-01Q, /* bffd5ffffffffffffac4ff5f4050b000 */
/* Compute tval1 = t. */
tval1 = (int) (t * TWO8);
- x -= __expl_table[T_EXPL_ARG1+2*tval1];
- xl -= __expl_table[T_EXPL_ARG1+2*tval1+1];
+ x -= __expq_table[T_EXPL_ARG1+2*tval1];
+ xl -= __expq_table[T_EXPL_ARG1+2*tval1+1];
/* Calculate t/32768. */
t = x + THREEp96;
/* Compute tval2 = t. */
tval2 = (int) (t * TWO15);
- x -= __expl_table[T_EXPL_ARG2+2*tval2];
- xl -= __expl_table[T_EXPL_ARG2+2*tval2+1];
+ x -= __expq_table[T_EXPL_ARG2+2*tval2];
+ xl -= __expq_table[T_EXPL_ARG2+2*tval2+1];
x = x + xl;
/* Compute ex2 = 2^n_0 e^(argtable[tval1]) e^(argtable[tval2]). */
- ex2_u.value = __expl_table[T_EXPL_RES1 + tval1]
- * __expl_table[T_EXPL_RES2 + tval2];
+ ex2_u.value = __expq_table[T_EXPL_RES1 + tval1]
+ * __expq_table[T_EXPL_RES2 + tval2];
n_i = (int)n;
/* 'unsafe' is 1 iff n_1 != 0. */
unsafe = abs(n_i) >= -FLT128_MIN_EXP - 1;
-/* s_fabsl.c -- __float128 version of s_fabs.c.
+/* fabsq.c -- __float128 version of s_fabs.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
-/* s_finitel.c -- long double version of s_finite.c.
+/* finiteq.c -- __float128 version of s_finite.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
-/* s_floorl.c -- long double version of s_floor.c.
+/* floorq.c -- __float128 version of s_floor.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
for proper rounding. */
if (v.ieee.exponent == 226)
{
+ /* If the exponent would be in the normal range when
+ rounding to normal precision with unbounded exponent
+ range, the exact result is known and spurious underflows
+ must be avoided on systems detecting tininess after
+ rounding. */
+ if (TININESS_AFTER_ROUNDING)
+ {
+ w.value = a1 + u.value;
+ if (w.ieee.exponent == 227)
+ return w.value * 0x1p-226L;
+ }
/* v.ieee.mant_low & 2 is LSB bit of the result before rounding,
v.ieee.mant_low & 1 is the round bit and j is our sticky
bit. */
-/* e_fmodl.c -- long double version of e_fmod.c.
+/* fmodq.c -- __float128 version of e_fmod.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
/*
-/* s_frexpl.c -- long double version of s_frexp.c.
+/* frexpq.c -- __float128 version of s_frexp.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
-/* s_ilogbl.c -- long double version of s_ilogb.c.
+/* ilogbq.c -- __float128 version of s_ilogb.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
static char rcsid[] = "$NetBSD: $";
#endif
-/* ilogbl(long double x)
+/* ilogbl(__float128 x)
* return the binary exponent of non-zero x
* ilogbl(0) = FP_ILOGB0
* ilogbl(NaN) = FP_ILOGBNAN (no signal is raised)
#include <limits.h>
#include <math.h>
+#include <errno.h>
#include "quadmath-imp.h"
#ifndef FP_ILOGB0
hx &= 0x7fffffffffffffffLL;
if(hx <= 0x0001000000000000LL) {
if((hx|lx)==0)
+ {
+ errno = EDOM;
+#ifdef USE_FENV_H
+ feraiseexcept (FE_INVALID);
+#endif
return FP_ILOGB0; /* ilogbl(0) = FP_ILOGB0 */
+ }
else /* subnormal x */
if(hx==0) {
for (ix = -16431; lx>0; lx<<=1) ix -=1;
else if (FP_ILOGBNAN != INT_MAX) {
/* ISO C99 requires ilogbl(+-Inf) == INT_MAX. */
if (((hx^0x7fff000000000000LL)|lx) == 0)
+ {
+ errno = EDOM;
+#ifdef USE_FENV_H
+ feraiseexcept (FE_INVALID);
+#endif
return INT_MAX;
+ }
}
+
+ errno = EDOM;
+#ifdef USE_FENV_H
+ feraiseexcept (FE_INVALID);
+#endif
return FP_ILOGBNAN;
}
--- /dev/null
+/*
+ * 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));
+}
+
-/* s_isnanl.c -- long double version of s_isnan.c.
+/* isnanq.c -- __float128 version of s_isnan.c.
* Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
*
* SYNOPSIS:
*
- * long double x, y, j0l();
+ * __float128 x, y, j0l();
*
* y = j0l( x );
*
*
* SYNOPSIS:
*
- * double x, y, y0l();
+ * __float128 x, y, y0l();
*
* y = y0l( x );
*
*
* SYNOPSIS:
*
- * long double x, y, j1l();
+ * __float128 x, y, j1q();
*
- * y = j1l( x );
+ * y = j1q( x );
*
*
*
*
* SYNOPSIS:
*
- * double x, y, y1l();
+ * __float128, y, y1q();
*
- * y = y1l( x );
+ * y = y1q( x );
*
*
*
-/* s_ldexpl.c -- long double version of s_ldexp.c.
+/* ldexpq.c -- __float128 version of s_ldexp.c.
* Conversion to long double by Ulrich Drepper,
* Cygnus Support, drepper@cygnus.com.
*/
* Returns the base e (2.718...) logarithm of the absolute
* value of the gamma function of the argument.
* The sign (+1 or -1) of the gamma function is returned in a
- * global (extern) variable named sgngam.
+ * global (extern) variable named signgam.
*
* The positive domain is partitioned into numerous segments for approximation.
* For x > 10,
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
#include "quadmath-imp.h"
+#include <math.h> /* For extern int signgam. */
static const __float128 PIQ = 3.1415926535897932384626433832795028841972E0Q;
static const __float128 MAXLGM = 1.0485738685148938358098967157129705071571E4928Q;
lgammaq (__float128 x)
{
__float128 p, q, w, z, nx;
- int i, nn, sign;
-
- sign = 1;
+ int i, nn;
+ int sign;
if (! finiteq (x))
return x * x;
if (x == 0.0Q)
{
if (signbitq (x))
- sign = -1;
+ sign = -1;
}
+ signgam = sign;
+
if (x < 0.0Q)
{
q = -x;
z = p - q;
}
z = q * sinq (PIQ * z);
+ signgam = sign;
+
if (z == 0.0Q)
return (sign * huge * huge);
w = lgammaq (q);
{
z = x - 0.75Q;
p = z * neval (z, RN1r75, NRN1r75)
- / deval (z, RD1r75, NRD1r75);
+ / deval (z, RD1r75, NRD1r75);
p += lgam1r75b;
p += lgam1r75a;
}
-/* Round long double value to long long int.
+/* Round __float128 value to long long int.
Copyright (C) 1997, 1999, 2004 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
-/* log10l.c
+/* log10q.c
*
- * Common logarithm, 128-bit long double precision
+ * Common logarithm, 128-bit __float128 precision
*
*
*
* SYNOPSIS:
*
- * long double x, y, log10l();
+ * __float128 x, y, log10l();
*
- * y = log10l( x );
+ * y = log10q( x );
*
*
*
/* log1pl.c
*
* Relative error logarithm
- * Natural logarithm of 1+x, 128-bit long double precision
+ * Natural logarithm of 1+x for __float128 precision
*
*
*
* SYNOPSIS:
*
- * long double x, y, log1pl();
+ * __float128 x, y, log1pl();
*
- * y = log1pl( x );
+ * y = log1pq( x );
*
*
*
-/* log2l.c
- * Base 2 logarithm, 128-bit long double precision
+/* log2q.c
+ * Base 2 logarithm for __float128 precision
*
*
*
* SYNOPSIS:
*
- * long double x, y, log2l();
+ * __float128 x, y, log2q();
*
- * y = log2l( x );
+ * y = log2q( x );
*
*
*
-/* logll.c
+/* logq.c
*
- * Natural logarithm for 128-bit long double precision.
+ * Natural logarithm for __float128 precision.
*
*
*
* SYNOPSIS:
*
- * long double x, y, logl();
+ * __float128 x, y, logq();
*
- * y = logl( x );
+ * y = logq( x );
*
*
*
-/* Round long double value to long int.
+/* Round __float128 value to long int.
Copyright (C) 1997, 1999, 2004 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
-/* s_modfl.c -- long double version of s_modf.c.
+/* modfq.c -- __float128 version of s_modf.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
-/* s_nextafterl.c -- long double version of s_nextafter.c.
+/* nextafterq.c -- __float128 version of s_nextafter.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
-/* __ieee754_powl(x,y) return x**y
+/* powq(x,y) return x**y
*
* n
* Method: Let x = 2 * (1+f)
*/
/*
- * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
+ * __quadmath_kernel_rem_pio2 (x,y,e0,nx,prec,ipio2)
* double x[],y[]; int e0,nx,prec; int ipio2[];
*
- * __kernel_rem_pio2 return the last three digits of N with
+ * __quadmath_kernel_rem_pio2 return the last three digits of N with
* y = x - N*pi/2
* so that |y| < pi/2.
*
-/* e_fmodl.c -- long double version of e_fmod.c.
+/* fmodq.c -- __float128 version of e_fmod.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
/*
-/* s_rintl.c -- long double version of s_rint.c.
+/* rintq.c -- __float128 version of s_rint.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
-/* Round long double to integer away from zero.
+/* Round __float128 to integer away from zero.
Copyright (C) 1997, 1999 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
-/* s_scalblnl.c -- long double version of s_scalbn.c.
+/* scalblnq.c -- __float128 version of s_scalbn.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
-/* s_scalbnl.c -- long double version of s_scalbn.c.
+/* scalbnq.c -- __float128 version of s_scalbn.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
{
/* 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
- possible values for h. We look up cosl(h) and sinl(h) in
- pre-computed tables, compute cosl(l) and sinl(l) using a
+ possible values for h. We look up cosq(h) and sinq(h) in
+ pre-computed tables, compute cosq(l) and sinq(l) using a
Chebyshev polynomial of degree 10(11) and compute
- sinl(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l) and
- cosl(h+l) = cosl(h)cosl(l) - sinl(h)sinl(l). */
+ sinq(h+l) = sinq(h)cosq(l) + cosq(h)sinq(l) and
+ cosq(h+l) = cosq(h)cosq(l) - sinq(h)sinq(l). */
index = 0x3ffe - (tix >> 16);
hix = (tix + (0x200 << index)) & (0xfffffc00 << index);
- x = fabsq (x);
+ if (signbitq (x))
+ {
+ x = -x;
+ y = -y;
+ }
switch (index)
{
case 0: index = ((45 << 10) + hix - 0x3ffe0000) >> 8; break;
-/* e_sinhl.c -- __float128 version of e_sinh.c.
+/* sinhq.c -- __float128 version of e_sinh.c.
* Conversion to __float128 by Ulrich Drepper,
* Cygnus Support, drepper@cygnus.com.
*/
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
-/* __ieee754_sinhl(x)
+/* sinhq(x)
* Method :
* mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
- * 1. Replace x by |x| (sinhl(-x) = -sinhl(x)).
+ * 1. Replace x by |x| (sinhq(-x) = -sinhq(x)).
* 2.
* E + E/(E+1)
- * 0 <= x <= 25 : sinhl(x) := --------------, E=expm1l(x)
+ * 0 <= x <= 25 : sinhq(x) := --------------, E=expm1q(x)
* 2
*
- * 25 <= x <= lnovft : sinhl(x) := expl(x)/2
- * lnovft <= x <= ln2ovft: sinhl(x) := expl(x/2)/2 * expl(x/2)
- * ln2ovft < x : sinhl(x) := x*shuge (overflow)
+ * 25 <= x <= lnovft : sinhq(x) := expq(x)/2
+ * lnovft <= x <= ln2ovft: sinhq(x) := expq(x/2)/2 * expq(x/2)
+ * ln2ovft < x : sinhq(x) := x*shuge (overflow)
*
* Special cases:
- * sinhl(x) is |x| if x is +INF, -INF, or NaN.
- * only sinhl(0)=0 is exact for finite x.
+ * sinhq(x) is |x| if x is +INF, -INF, or NaN.
+ * only sinhq(0)=0 is exact for finite x.
*/
#include "quadmath-imp.h"
return t * w;
}
- /* |x| > overflowthreshold, sinhl(x) overflow */
+ /* |x| > overflowthreshold, sinhq(x) overflow */
return x * shuge;
}
-/* s_sinl.c -- long double version of s_sin.c.
+/* sinq.c -- __float128 version of s_sin.c.
* Conversion to long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
* ====================================================
*/
-/* sinl(x)
+/* sinq(x)
* Return sine function of x.
*
* kernel function:
- * __kernel_sinl ... sine function on [-pi/4,pi/4]
- * __kernel_cosl ... cose function on [-pi/4,pi/4]
- * __ieee754_rem_pio2l ... argument reduction routine
+ * __quadmath_kernel_sinq ... sine function on [-pi/4,pi/4]
+ * __quadmath_kernel_cosq ... cose function on [-pi/4,pi/4]
+ * __quadmath_rem_pio2q ... argument reduction routine
*
* Method.
* Let S,C and T denote the sin, cos and tan respectively on
{
/* 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
- possible values for h. We look up cosl(h) and sinl(h) in
- pre-computed tables, compute cosl(l) and sinl(l) using a
+ possible values for h. We look up cosq(h) and sinq(h) in
+ pre-computed tables, compute cosq(l) and sinq(l) using a
Chebyshev polynomial of degree 10(11) and compute
- sinl(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l). */
+ sinq(h+l) = sinq(h)cosq(l) + cosq(h)sinq(l). */
index = 0x3ffe - (tix >> 16);
hix = (tix + (0x200 << index)) & (0xfffffc00 << index);
x = fabsq (x);
SET_FLT128_WORDS64(h, ((uint64_t)hix) << 32, 0);
if (iy)
- l = y - (h - x);
+ l = (ix < 0 ? -y : y) - (h - x);
else
l = x - h;
z = l * l;
-/* s_tanl.c -- long double version of s_tan.c.
+/* tanq.c -- __float128 version of s_tan.c.
* Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
*/
* Return tangent function of x.
*
* kernel function:
- * __kernel_tanq ... tangent function on [-pi/4,pi/4]
- * __ieee754_rem_pio2q ... argument reduction routine
+ * __quadmath_kernel_tanq ... tangent function on [-pi/4,pi/4]
+ * __quadmath_rem_pio2q ... argument reduction routine
*
* Method.
* Let S,C and T denote the sin, cos and tan respectively on
conditions we must check some values separately. */
int64_t hx;
uint64_t lx;
+ __float128 res;
+ int sign;
GET_FLT128_WORDS64 (hx, lx, x);
return x - x;
/* XXX FIXME. */
- return expq (lgammaq (x));
+ res = expq (lgammaq (x));
+ return signbitq (x) ? -res : res;
}
--- /dev/null
+/* Compute x^2 + y^2 - 1, without large cancellation error.
+ Copyright (C) 2012 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"
+#include <stdlib.h>
+
+/* Calculate X + Y exactly and store the result in *HI + *LO. It is
+ given that |X| >= |Y| and the values are small enough that no
+ overflow occurs. */
+
+static inline void
+add_split (__float128 *hi, __float128 *lo, __float128 x, __float128 y)
+{
+ /* Apply Dekker's algorithm. */
+ *hi = x + y;
+ *lo = (x - *hi) + y;
+}
+
+/* Calculate X * Y exactly and store the result in *HI + *LO. It is
+ given that the values are small enough that no overflow occurs and
+ large enough (or zero) that no underflow occurs. */
+
+static inline void
+mul_split (__float128 *hi, __float128 *lo, __float128 x, __float128 y)
+{
+ /* Fast built-in fused multiply-add. */
+ *hi = x * y;
+ *lo = fmaq (x, y, -*hi);
+}
+
+/* Compare absolute values of floating-point values pointed to by P
+ and Q for qsort. */
+
+static int
+compare (const void *p, const void *q)
+{
+ __float128 pld = fabsq (*(const __float128 *) p);
+ __float128 qld = fabsq (*(const __float128 *) q);
+ if (pld < qld)
+ return -1;
+ else if (pld == qld)
+ return 0;
+ else
+ return 1;
+}
+
+/* Return X^2 + Y^2 - 1, computed without large cancellation error.
+ It is given that 1 > X >= Y >= epsilon / 2, and that either X >=
+ 0.75 or Y >= 0.5. */
+
+__float128
+__quadmath_x2y2m1q (__float128 x, __float128 y)
+{
+ __float128 vals[4];
+ size_t i;
+
+ /* FIXME: SET_RESTORE_ROUNDL (FE_TONEAREST); */
+ mul_split (&vals[1], &vals[0], x, x);
+ mul_split (&vals[3], &vals[2], y, y);
+ if (x >= 0.75Q)
+ vals[1] -= 1.0Q;
+ else
+ {
+ vals[1] -= 0.5Q;
+ vals[3] -= 0.5Q;
+ }
+ qsort (vals, 4, sizeof (__float128), compare);
+ /* Add up the values so that each element of VALS has absolute value
+ at most equal to the last set bit of the next nonzero
+ element. */
+ for (i = 0; i <= 2; i++)
+ {
+ add_split (&vals[i + 1], &vals[i], vals[i + 1], vals[i]);
+ qsort (vals + i + 1, 3 - i, sizeof (__float128), compare);
+ }
+ /* Now any error from this addition will be small. */
+ return vals[3] + vals[2] + vals[1] + vals[0];
+}
#include "config.h"
+/* Under IEEE 754, an architecture may determine tininess of
+ floating-point results either "before rounding" or "after
+ rounding", but must do so in the same way for all operations
+ returning binary results. Define TININESS_AFTER_ROUNDING to 1 for
+ "after rounding" architectures, 0 for "before rounding"
+ architectures. */
+
+#define TININESS_AFTER_ROUNDING 1
+
+
/* Prototypes for internal functions. */
extern int32_t __quadmath_rem_pio2q (__float128, __float128 *);
extern void __quadmath_kernel_sincosq (__float128, __float128, __float128 *,
__float128 *, int);
extern __float128 __quadmath_kernel_sinq (__float128, __float128, int);
extern __float128 __quadmath_kernel_cosq (__float128, __float128);
+extern __float128 __quadmath_x2y2m1q (__float128 x, __float128 y);
+extern int __quadmath_isinf_nsq (__float128 x);
+
+
projects / gcc.git / commitdiff