diff options
Diffstat (limited to 'libquadmath')
65 files changed, 1843 insertions, 365 deletions
diff --git a/libquadmath/ChangeLog b/libquadmath/ChangeLog index 37cc5a25028..efb95d179d2 100644 --- a/libquadmath/ChangeLog +++ b/libquadmath/ChangeLog @@ -1,3 +1,79 @@ +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> diff --git a/libquadmath/Makefile.am b/libquadmath/Makefile.am index c7be3e55414..6c97ee81c5c 100644 --- a/libquadmath/Makefile.am +++ b/libquadmath/Makefile.am @@ -43,7 +43,8 @@ nodist_libsubinclude_HEADERS = quadmath.h quadmath_weak.h 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 \ @@ -60,6 +61,8 @@ libquadmath_la_SOURCES = \ 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 \ diff --git a/libquadmath/Makefile.in b/libquadmath/Makefile.in index 6e389cf6af1..92c5d256d5a 100644 --- a/libquadmath/Makefile.in +++ b/libquadmath/Makefile.in @@ -87,7 +87,8 @@ am__installdirs = "$(DESTDIR)$(toolexeclibdir)" "$(DESTDIR)$(infodir)" \ "$(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 \ @@ -126,9 +127,13 @@ am__dirstamp = $(am__leading_dot)dirstamp @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 \ @@ -321,7 +326,8 @@ AUTOMAKE_OPTIONS = 1.8 foreign @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 \ @@ -338,6 +344,8 @@ AUTOMAKE_OPTIONS = 1.8 foreign @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 \ @@ -504,6 +512,8 @@ math/$(am__dirstamp): 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) @@ -588,6 +598,15 @@ math/lrintq.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) @@ -669,10 +688,18 @@ mostlyclean-compile: -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) @@ -689,6 +716,16 @@ mostlyclean-compile: -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) @@ -717,6 +754,8 @@ mostlyclean-compile: -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) @@ -795,6 +834,8 @@ mostlyclean-compile: -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) @@ -851,8 +892,12 @@ distclean-compile: @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@ @@ -861,6 +906,11 @@ distclean-compile: @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@ @@ -875,6 +925,7 @@ distclean-compile: @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@ @@ -914,6 +965,7 @@ distclean-compile: @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@ diff --git a/libquadmath/math/acoshq.c b/libquadmath/math/acoshq.c index b0a1e4e5496..9845a8e364c 100644 --- a/libquadmath/math/acoshq.c +++ b/libquadmath/math/acoshq.c @@ -1,4 +1,4 @@ -/* 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. */ @@ -13,7 +13,7 @@ * ==================================================== */ -/* __ieee754_acoshl(x) +/* acoshq(x) * Method : * Based on * acoshl(x) = logl [ x + sqrtl(x*x-1) ] diff --git a/libquadmath/math/acosq.c b/libquadmath/math/acosq.c index a8a361d23bb..7ef79474651 100644 --- a/libquadmath/math/acosq.c +++ b/libquadmath/math/acosq.c @@ -31,7 +31,7 @@ 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) @@ -51,7 +51,7 @@ * 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" diff --git a/libquadmath/math/asinhq.c b/libquadmath/math/asinhq.c index be044dcd87f..9be0aa1f053 100644 --- a/libquadmath/math/asinhq.c +++ b/libquadmath/math/asinhq.c @@ -1,4 +1,4 @@ -/* 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. */ diff --git a/libquadmath/math/asinq.c b/libquadmath/math/asinq.c index d4321a58e0a..7bd4d768c97 100644 --- a/libquadmath/math/asinq.c +++ b/libquadmath/math/asinq.c @@ -31,7 +31,7 @@ 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 diff --git a/libquadmath/math/atan2q.c b/libquadmath/math/atan2q.c index fbe64d62b95..daa303efba3 100644 --- a/libquadmath/math/atan2q.c +++ b/libquadmath/math/atan2q.c @@ -1,4 +1,4 @@ -/* 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. */ diff --git a/libquadmath/math/atanhq.c b/libquadmath/math/atanhq.c index 73db957d3b0..b34036715d7 100644 --- a/libquadmath/math/atanhq.c +++ b/libquadmath/math/atanhq.c @@ -14,7 +14,7 @@ * ==================================================== */ -/* __ieee754_atanhl(x) +/* atanhq(x) * Method : * 1.Reduced x to positive by atanh(-x) = -atanh(x) * 2.For x>=0.5 diff --git a/libquadmath/math/cacoshq.c b/libquadmath/math/cacoshq.c index 8acc570de76..263e03d0c11 100644 --- a/libquadmath/math/cacoshq.c +++ b/libquadmath/math/cacoshq.c @@ -63,6 +63,16 @@ cacoshq (__complex128 x) __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; @@ -72,17 +82,13 @@ cacoshq (__complex128 x) 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; diff --git a/libquadmath/math/casinhq.c b/libquadmath/math/casinhq.c index ffa45fa81d9..11487b967fc 100644 --- a/libquadmath/math/casinhq.c +++ b/libquadmath/math/casinhq.c @@ -72,6 +72,11 @@ casinhq (__complex128 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; diff --git a/libquadmath/math/cbrtq.c b/libquadmath/math/cbrtq.c index f61f32513ee..f1f05cac789 100644 --- a/libquadmath/math/cbrtq.c +++ b/libquadmath/math/cbrtq.c @@ -1,64 +1,132 @@ +/* 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); +} diff --git a/libquadmath/math/ccoshq.c b/libquadmath/math/ccoshq.c new file mode 100644 index 00000000000..8d55ad3a99d --- /dev/null +++ b/libquadmath/math/ccoshq.c @@ -0,0 +1,145 @@ +/* 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; +} diff --git a/libquadmath/math/ceilq.c b/libquadmath/math/ceilq.c index 577d8cd9895..0d9bb8b87f3 100644 --- a/libquadmath/math/ceilq.c +++ b/libquadmath/math/ceilq.c @@ -1,4 +1,4 @@ -/* 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. */ diff --git a/libquadmath/math/cexpq.c b/libquadmath/math/cexpq.c new file mode 100644 index 00000000000..bd4be1ebe71 --- /dev/null +++ b/libquadmath/math/cexpq.c @@ -0,0 +1,152 @@ +/* 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; +} diff --git a/libquadmath/math/clog10q.c b/libquadmath/math/clog10q.c new file mode 100644 index 00000000000..c379bec72b6 --- /dev/null +++ b/libquadmath/math/clog10q.c @@ -0,0 +1,116 @@ +/* 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; +} diff --git a/libquadmath/math/clogq.c b/libquadmath/math/clogq.c new file mode 100644 index 00000000000..1a772cd434d --- /dev/null +++ b/libquadmath/math/clogq.c @@ -0,0 +1,111 @@ +/* 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; +} diff --git a/libquadmath/math/complex.c b/libquadmath/math/complex.c index 9953f52ca6b..8cf9e9808cc 100644 --- a/libquadmath/math/complex.c +++ b/libquadmath/math/complex.c @@ -1,50 +1,35 @@ +/* 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) { @@ -53,23 +38,12 @@ 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; } @@ -82,137 +56,17 @@ cargq (__complex128 z) __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); } - diff --git a/libquadmath/math/copysignq.c b/libquadmath/math/copysignq.c index b59fcc5492d..054de2d2eb3 100644 --- a/libquadmath/math/copysignq.c +++ b/libquadmath/math/copysignq.c @@ -1,4 +1,4 @@ -/* 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. */ diff --git a/libquadmath/math/coshq.c b/libquadmath/math/coshq.c index a6e0eb5f193..77f4b98338b 100644 --- a/libquadmath/math/coshq.c +++ b/libquadmath/math/coshq.c @@ -30,25 +30,25 @@ 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" @@ -73,7 +73,7 @@ coshq (__float128 x) 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); diff --git a/libquadmath/math/cosq.c b/libquadmath/math/cosq.c index dc321a27d17..28630b93a6f 100644 --- a/libquadmath/math/cosq.c +++ b/libquadmath/math/cosq.c @@ -1,4 +1,4 @@ -/* 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. */ @@ -13,13 +13,13 @@ * ==================================================== */ -/* 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 diff --git a/libquadmath/math/cosq_kernel.c b/libquadmath/math/cosq_kernel.c index 86f39551c32..42d29adfd8c 100644 --- a/libquadmath/math/cosq_kernel.c +++ b/libquadmath/math/cosq_kernel.c @@ -99,13 +99,17 @@ __quadmath_kernel_cosq (__float128 x, __float128 y) { /* 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; diff --git a/libquadmath/math/csinhq.c b/libquadmath/math/csinhq.c new file mode 100644 index 00000000000..c16d576f4da --- /dev/null +++ b/libquadmath/math/csinhq.c @@ -0,0 +1,166 @@ +/* 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; +} diff --git a/libquadmath/math/csinq.c b/libquadmath/math/csinq.c new file mode 100644 index 00000000000..c837e50b87f --- /dev/null +++ b/libquadmath/math/csinq.c @@ -0,0 +1,171 @@ +/* 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; +} diff --git a/libquadmath/math/csqrtq.c b/libquadmath/math/csqrtq.c new file mode 100644 index 00000000000..5279e4378a2 --- /dev/null +++ b/libquadmath/math/csqrtq.c @@ -0,0 +1,143 @@ +/* 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; +} diff --git a/libquadmath/math/ctanhq.c b/libquadmath/math/ctanhq.c new file mode 100644 index 00000000000..8934cfad59f --- /dev/null +++ b/libquadmath/math/ctanhq.c @@ -0,0 +1,120 @@ +/* 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; +} diff --git a/libquadmath/math/ctanq.c b/libquadmath/math/ctanq.c new file mode 100644 index 00000000000..d390511ca20 --- /dev/null +++ b/libquadmath/math/ctanq.c @@ -0,0 +1,120 @@ +/* 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; +} diff --git a/libquadmath/math/erfq.c b/libquadmath/math/erfq.c index 50db88ae821..8d383e9ca70 100644 --- a/libquadmath/math/erfq.c +++ b/libquadmath/math/erfq.c @@ -30,8 +30,8 @@ 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 diff --git a/libquadmath/math/expq.c b/libquadmath/math/expq.c index 2740b4e2cc9..70c638d7a08 100644 --- a/libquadmath/math/expq.c +++ b/libquadmath/math/expq.c @@ -30,19 +30,19 @@ #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 */ @@ -1122,8 +1122,8 @@ expq (__float128 x) /* 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; @@ -1132,14 +1132,14 @@ expq (__float128 x) /* 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; diff --git a/libquadmath/math/fabsq.c b/libquadmath/math/fabsq.c index 7ec2850feb7..a103f840f38 100644 --- a/libquadmath/math/fabsq.c +++ b/libquadmath/math/fabsq.c @@ -1,4 +1,4 @@ -/* 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. */ diff --git a/libquadmath/math/finiteq.c b/libquadmath/math/finiteq.c index 663c9289263..c5326d2194b 100644 --- a/libquadmath/math/finiteq.c +++ b/libquadmath/math/finiteq.c @@ -1,4 +1,4 @@ -/* 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. */ diff --git a/libquadmath/math/floorq.c b/libquadmath/math/floorq.c index 781cf40e3cb..0d74f94e27d 100644 --- a/libquadmath/math/floorq.c +++ b/libquadmath/math/floorq.c @@ -1,4 +1,4 @@ -/* 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. */ diff --git a/libquadmath/math/fmaq.c b/libquadmath/math/fmaq.c index 23e3188669e..5616c1a2781 100644 --- a/libquadmath/math/fmaq.c +++ b/libquadmath/math/fmaq.c @@ -228,6 +228,17 @@ fmaq (__float128 x, __float128 y, __float128 z) 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. */ diff --git a/libquadmath/math/fmodq.c b/libquadmath/math/fmodq.c index 036e3dd531c..55eb18ffb0a 100644 --- a/libquadmath/math/fmodq.c +++ b/libquadmath/math/fmodq.c @@ -1,4 +1,4 @@ -/* 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. */ /* diff --git a/libquadmath/math/frexpq.c b/libquadmath/math/frexpq.c index fa3d7836d06..b0b305af574 100644 --- a/libquadmath/math/frexpq.c +++ b/libquadmath/math/frexpq.c @@ -1,4 +1,4 @@ -/* 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. */ diff --git a/libquadmath/math/ilogbq.c b/libquadmath/math/ilogbq.c index 47986f5b311..7f95e9c2240 100644 --- a/libquadmath/math/ilogbq.c +++ b/libquadmath/math/ilogbq.c @@ -1,4 +1,4 @@ -/* 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. */ @@ -17,7 +17,7 @@ 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) @@ -26,6 +26,7 @@ static char rcsid[] = "$NetBSD: $"; #include <limits.h> #include <math.h> +#include <errno.h> #include "quadmath-imp.h" #ifndef FP_ILOGB0 @@ -45,7 +46,13 @@ ilogbq (__float128 x) 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; @@ -58,7 +65,18 @@ ilogbq (__float128 x) 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; } diff --git a/libquadmath/math/isinf_nsq.c b/libquadmath/math/isinf_nsq.c new file mode 100644 index 00000000000..2f0834361c5 --- /dev/null +++ b/libquadmath/math/isinf_nsq.c @@ -0,0 +1,19 @@ +/* + * 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)); +} + diff --git a/libquadmath/math/isnanq.c b/libquadmath/math/isnanq.c index ab9df1658c8..4462d0139e0 100644 --- a/libquadmath/math/isnanq.c +++ b/libquadmath/math/isnanq.c @@ -1,4 +1,4 @@ -/* 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. */ diff --git a/libquadmath/math/j0q.c b/libquadmath/math/j0q.c index fecbe627746..3ef9201a8a4 100644 --- a/libquadmath/math/j0q.c +++ b/libquadmath/math/j0q.c @@ -6,7 +6,7 @@ * * SYNOPSIS: * - * long double x, y, j0l(); + * __float128 x, y, j0l(); * * y = j0l( x ); * @@ -52,7 +52,7 @@ * * SYNOPSIS: * - * double x, y, y0l(); + * __float128 x, y, y0l(); * * y = y0l( x ); * diff --git a/libquadmath/math/j1q.c b/libquadmath/math/j1q.c index f6bf2a256ce..8a18e2779b2 100644 --- a/libquadmath/math/j1q.c +++ b/libquadmath/math/j1q.c @@ -6,9 +6,9 @@ * * SYNOPSIS: * - * long double x, y, j1l(); + * __float128 x, y, j1q(); * - * y = j1l( x ); + * y = j1q( x ); * * * @@ -52,9 +52,9 @@ * * SYNOPSIS: * - * double x, y, y1l(); + * __float128, y, y1q(); * - * y = y1l( x ); + * y = y1q( x ); * * * diff --git a/libquadmath/math/ldexpq.c b/libquadmath/math/ldexpq.c index 394d4590cca..c18968b03bc 100644 --- a/libquadmath/math/ldexpq.c +++ b/libquadmath/math/ldexpq.c @@ -1,4 +1,4 @@ -/* 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. */ diff --git a/libquadmath/math/lgammaq.c b/libquadmath/math/lgammaq.c index 6e7697ac6a1..361f7037bc3 100644 --- a/libquadmath/math/lgammaq.c +++ b/libquadmath/math/lgammaq.c @@ -18,7 +18,7 @@ * 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, @@ -69,6 +69,7 @@ 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; @@ -757,9 +758,8 @@ __float128 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; @@ -767,9 +767,11 @@ lgammaq (__float128 x) if (x == 0.0Q) { if (signbitq (x)) - sign = -1; + sign = -1; } + signgam = sign; + if (x < 0.0Q) { q = -x; @@ -788,6 +790,8 @@ lgammaq (__float128 x) z = p - q; } z = q * sinq (PIQ * z); + signgam = sign; + if (z == 0.0Q) return (sign * huge * huge); w = lgammaq (q); @@ -855,7 +859,7 @@ lgammaq (__float128 x) { z = x - 0.75Q; p = z * neval (z, RN1r75, NRN1r75) - / deval (z, RD1r75, NRD1r75); + / deval (z, RD1r75, NRD1r75); p += lgam1r75b; p += lgam1r75a; } diff --git a/libquadmath/math/llroundq.c b/libquadmath/math/llroundq.c index c108e7ad18a..d22180d6bba 100644 --- a/libquadmath/math/llroundq.c +++ b/libquadmath/math/llroundq.c @@ -1,4 +1,4 @@ -/* 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 diff --git a/libquadmath/math/log10q.c b/libquadmath/math/log10q.c index 50caf18b7fd..9eeb9ae3fc4 100644 --- a/libquadmath/math/log10q.c +++ b/libquadmath/math/log10q.c @@ -1,14 +1,14 @@ -/* 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 ); * * * diff --git a/libquadmath/math/log1pq.c b/libquadmath/math/log1pq.c index a466dc8921e..2de13fe2adc 100644 --- a/libquadmath/math/log1pq.c +++ b/libquadmath/math/log1pq.c @@ -1,15 +1,15 @@ /* 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 ); * * * diff --git a/libquadmath/math/log2q.c b/libquadmath/math/log2q.c index 963b38c8483..f8275369b37 100644 --- a/libquadmath/math/log2q.c +++ b/libquadmath/math/log2q.c @@ -1,13 +1,13 @@ -/* 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 ); * * * diff --git a/libquadmath/math/logq.c b/libquadmath/math/logq.c index cd1a48631fa..7aae9b101ad 100644 --- a/libquadmath/math/logq.c +++ b/libquadmath/math/logq.c @@ -1,14 +1,14 @@ -/* 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 ); * * * diff --git a/libquadmath/math/lroundq.c b/libquadmath/math/lroundq.c index 4b12d245267..59c883a1464 100644 --- a/libquadmath/math/lroundq.c +++ b/libquadmath/math/lroundq.c @@ -1,4 +1,4 @@ -/* 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 diff --git a/libquadmath/math/modfq.c b/libquadmath/math/modfq.c index dc564fe3555..8c5db54bb76 100644 --- a/libquadmath/math/modfq.c +++ b/libquadmath/math/modfq.c @@ -1,4 +1,4 @@ -/* 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. */ diff --git a/libquadmath/math/nextafterq.c b/libquadmath/math/nextafterq.c index 01bfa657901..04d63deb8e2 100644 --- a/libquadmath/math/nextafterq.c +++ b/libquadmath/math/nextafterq.c @@ -1,4 +1,4 @@ -/* 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. */ diff --git a/libquadmath/math/powq.c b/libquadmath/math/powq.c index d3863243825..12b87d536d7 100644 --- a/libquadmath/math/powq.c +++ b/libquadmath/math/powq.c @@ -30,7 +30,7 @@ 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) diff --git a/libquadmath/math/rem_pio2q.c b/libquadmath/math/rem_pio2q.c index 47ee8ef2051..60653c8d1d3 100644 --- a/libquadmath/math/rem_pio2q.c +++ b/libquadmath/math/rem_pio2q.c @@ -15,10 +15,10 @@ */ /* - * __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. * diff --git a/libquadmath/math/remainderq.c b/libquadmath/math/remainderq.c index 421f728ffe3..c3f5641293f 100644 --- a/libquadmath/math/remainderq.c +++ b/libquadmath/math/remainderq.c @@ -1,4 +1,4 @@ -/* 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. */ /* diff --git a/libquadmath/math/rintq.c b/libquadmath/math/rintq.c index 8a93fdbfa78..4a50503721c 100644 --- a/libquadmath/math/rintq.c +++ b/libquadmath/math/rintq.c @@ -1,4 +1,4 @@ -/* 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. */ diff --git a/libquadmath/math/roundq.c b/libquadmath/math/roundq.c index d7f577bd941..7c9d640e933 100644 --- a/libquadmath/math/roundq.c +++ b/libquadmath/math/roundq.c @@ -1,4 +1,4 @@ -/* 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 diff --git a/libquadmath/math/scalblnq.c b/libquadmath/math/scalblnq.c index b2332250c61..a414045b51f 100644 --- a/libquadmath/math/scalblnq.c +++ b/libquadmath/math/scalblnq.c @@ -1,4 +1,4 @@ -/* 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. */ diff --git a/libquadmath/math/scalbnq.c b/libquadmath/math/scalbnq.c index f0852ee038d..9975a47154c 100644 --- a/libquadmath/math/scalbnq.c +++ b/libquadmath/math/scalbnq.c @@ -1,4 +1,4 @@ -/* 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. */ diff --git a/libquadmath/math/sincosq_kernel.c b/libquadmath/math/sincosq_kernel.c index 578d1828f75..f6341a4d948 100644 --- a/libquadmath/math/sincosq_kernel.c +++ b/libquadmath/math/sincosq_kernel.c @@ -126,14 +126,18 @@ __quadmath_kernel_sincosq(__float128 x, __float128 y, __float128 *sinx, { /* 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; diff --git a/libquadmath/math/sinhq.c b/libquadmath/math/sinhq.c index 5492180a2a2..eff8149b515 100644 --- a/libquadmath/math/sinhq.c +++ b/libquadmath/math/sinhq.c @@ -1,4 +1,4 @@ -/* 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. */ @@ -35,22 +35,22 @@ 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" @@ -106,6 +106,6 @@ sinhq (__float128 x) return t * w; } - /* |x| > overflowthreshold, sinhl(x) overflow */ + /* |x| > overflowthreshold, sinhq(x) overflow */ return x * shuge; } diff --git a/libquadmath/math/sinq.c b/libquadmath/math/sinq.c index 76254a37302..989f679d6d6 100644 --- a/libquadmath/math/sinq.c +++ b/libquadmath/math/sinq.c @@ -1,4 +1,4 @@ -/* 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. */ @@ -13,13 +13,13 @@ * ==================================================== */ -/* 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 diff --git a/libquadmath/math/sinq_kernel.c b/libquadmath/math/sinq_kernel.c index 395fcaba9cb..86034551d43 100644 --- a/libquadmath/math/sinq_kernel.c +++ b/libquadmath/math/sinq_kernel.c @@ -99,10 +99,10 @@ __quadmath_kernel_sinq (__float128 x, __float128 y, int iy) { /* 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); @@ -116,7 +116,7 @@ __quadmath_kernel_sinq (__float128 x, __float128 y, int iy) 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; diff --git a/libquadmath/math/tanq.c b/libquadmath/math/tanq.c index e1ec6aae86c..690d94b782c 100644 --- a/libquadmath/math/tanq.c +++ b/libquadmath/math/tanq.c @@ -160,7 +160,7 @@ __quadmath_kernel_tanq (__float128 x, __float128 y, int iy) -/* 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. */ @@ -180,8 +180,8 @@ __quadmath_kernel_tanq (__float128 x, __float128 y, int iy) * 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 diff --git a/libquadmath/math/tgammaq.c b/libquadmath/math/tgammaq.c index 3305b6484cd..a07d5831de0 100644 --- a/libquadmath/math/tgammaq.c +++ b/libquadmath/math/tgammaq.c @@ -30,6 +30,8 @@ tgammaq (__float128 x) conditions we must check some values separately. */ int64_t hx; uint64_t lx; + __float128 res; + int sign; GET_FLT128_WORDS64 (hx, lx, x); @@ -46,5 +48,6 @@ tgammaq (__float128 x) return x - x; /* XXX FIXME. */ - return expq (lgammaq (x)); + res = expq (lgammaq (x)); + return signbitq (x) ? -res : res; } diff --git a/libquadmath/math/x2y2m1q.c b/libquadmath/math/x2y2m1q.c new file mode 100644 index 00000000000..90bbc2f605d --- /dev/null +++ b/libquadmath/math/x2y2m1q.c @@ -0,0 +1,93 @@ +/* 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]; +} diff --git a/libquadmath/quadmath-imp.h b/libquadmath/quadmath-imp.h index bac714d1c8b..40b346b6ff2 100644 --- a/libquadmath/quadmath-imp.h +++ b/libquadmath/quadmath-imp.h @@ -27,12 +27,26 @@ Boston, MA 02110-1301, USA. */ #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); + + |