git
/
gcc
mirror of git://gcc.gnu.org/git/gcc.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Makefile.am (libquadmath_la_SOURCES): Add new math/* files. 

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.

From-SVN: r193063
This commit is contained in:
Tobias Burnus 2012-11-01 17:14:42 +01:00 committed by Tobias Burnus
parent 621cf8af1d
commit f029f4be17
65 changed files with 1837 additions and 359 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
libquadmath/ChangeLog
View File

@ -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>

5
libquadmath/Makefile.am
View File

@ -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 \

62
libquadmath/Makefile.in
View File

@ -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@

4
libquadmath/math/acoshq.c
View File

@ -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) ]

4
libquadmath/math/acosq.c
View File

@ -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"

2
libquadmath/math/asinhq.c
View File

@ -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.
*/

2
libquadmath/math/asinq.c
View File

@ -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

2
libquadmath/math/atan2q.c
View File

@ -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.
*/

2
libquadmath/math/atanhq.c
View File

@ -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

16
libquadmath/math/cacoshq.c
View File

@ -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;

5
libquadmath/math/casinhq.c
View File

@ -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;

168
libquadmath/math/cbrtq.c
View File

@ -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;
return (x);
if (isnanq (x))
return x;
if (x > 0)
sign = 1;
else
{
sign = -1;
x = -x;
}
if (x <= DBL_MAX && x >= DBL_MIN)
{
/* Use double result as starting point. */
y = cbrt ((double) x);
z = x;
/* extract power of 2, leaving mantissa between 0.5 and 1 */
x = frexpq (x, &e);
/* Two Newton iterations. */
y -= 0.333333333333333333333333333333333333333333333333333Q
* (y - x / (y * y));
y -= 0.333333333333333333333333333333333333333333333333333Q
* (y - x / (y * y));
return y;
}
/* 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;
#ifdef HAVE_CBRTL
if (x <= LDBL_MAX && x >= LDBL_MIN)
{
/* Use long double result as starting point. */
y = cbrtl ((long double) x);
/* 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;
}
/* One Newton iteration. */
y -= 0.333333333333333333333333333333333333333333333333333Q
* (y - x / (y * y));
return y;
}
#endif
/* multiply by power of 2 */
x = ldexpq (x, e);
/* If we're outside of the range of C types, we have to compute
the initial guess the hard way. */
y = frexpq (x, &exp);
/* Newton iteration */
x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
x -= (x - (z / (x * x))) * 0.3333333333333333333333333333333333333333L;
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;
if (sign < 0)
x = -x;
return (x);
}

145
libquadmath/math/ccoshq.c Normal file
View File

@ -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;
}

2
libquadmath/math/ceilq.c
View File

@ -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.
*/

152
libquadmath/math/cexpq.c Normal file
View File

@ -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;
}

116
libquadmath/math/clog10q.c Normal file
View File

@ -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;
}

111
libquadmath/math/clogq.c Normal file
View File

@ -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;
}

212
libquadmath/math/complex.c
View File

@ -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)
{
@ -52,24 +37,13 @@ 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;
}
@ -81,138 +55,18 @@ 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 y;
COMPLEX_ASSIGN (y, -IMAGPART (x), REALPART (x));
return ccoshq (y);
}
__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;
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;
}

2
libquadmath/math/copysignq.c
View File

@ -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.
*/

22
libquadmath/math/coshq.c
View File

@ -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);

10
libquadmath/math/cosq.c
View File

@ -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

12
libquadmath/math/cosq_kernel.c
View File

@ -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;

166
libquadmath/math/csinhq.c Normal file
View File

@ -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;
}

171
libquadmath/math/csinq.c Normal file
View File

@ -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;
}

143
libquadmath/math/csqrtq.c Normal file
View File

@ -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;
}

120
libquadmath/math/ctanhq.c Normal file
View File

@ -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;
}

120
libquadmath/math/ctanq.c Normal file
View File

@ -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;
}

4
libquadmath/math/erfq.c
View File

@ -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

22
libquadmath/math/expq.c
View File

@ -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;

2
libquadmath/math/fabsq.c
View File

@ -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.
*/

2
libquadmath/math/finiteq.c
View File

@ -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.
*/

2
libquadmath/math/floorq.c
View File

@ -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.
*/

11
libquadmath/math/fmaq.c
View File

@ -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. */

2
libquadmath/math/fmodq.c
View File

@ -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.
*/
/*

2
libquadmath/math/frexpq.c
View File

@ -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.
*/

22
libquadmath/math/ilogbq.c
View File

@ -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;
}

19
libquadmath/math/isinf_nsq.c Normal file
View File

@ -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));
}

2
libquadmath/math/isnanq.c
View File

@ -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.
*/

4
libquadmath/math/j0q.c
View File

@ -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 );
*

8
libquadmath/math/j1q.c
View File

@ -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 );
*
*
*

2
libquadmath/math/ldexpq.c
View File

@ -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.
*/

16
libquadmath/math/lgammaq.c
View File

@ -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;
}

2
libquadmath/math/llroundq.c
View File

@ -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

8
libquadmath/math/log10q.c
View File

@ -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 );
*
*
*

6
libquadmath/math/log1pq.c
View File

@ -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 );
*
*
*

8
libquadmath/math/log2q.c
View File

@ -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 );
*
*
*

8
libquadmath/math/logq.c
View File

@ -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 );
*
*
*

2
libquadmath/math/lroundq.c
View File

@ -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

2
libquadmath/math/modfq.c
View File

@ -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.
*/

2
libquadmath/math/nextafterq.c
View File

@ -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.
*/

2
libquadmath/math/powq.c
View File

@ -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)

4
libquadmath/math/rem_pio2q.c
View File

@ -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.
*

2
libquadmath/math/remainderq.c
View File

@ -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.
*/
/*

2
libquadmath/math/rintq.c
View File

@ -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.
*/

2
libquadmath/math/roundq.c
View File

@ -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

2
libquadmath/math/scalblnq.c
View File

@ -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.
*/

2
libquadmath/math/scalbnq.c
View File

@ -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.
*/

14
libquadmath/math/sincosq_kernel.c
View File

@ -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;

20
libquadmath/math/sinhq.c
View File

@ -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;
}

10
libquadmath/math/sinq.c
View File

@ -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

8
libquadmath/math/sinq_kernel.c
View File

@ -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;

6
libquadmath/math/tanq.c
View File

@ -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

5
libquadmath/math/tgammaq.c
View File

@ -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;
}

93
libquadmath/math/x2y2m1q.c Normal file
View File

@ -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];
}

14
libquadmath/quadmath-imp.h
View File

@ -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);