diff options
| author | joseph <joseph@7b3dc134-2b1b-0410-93df-9e9f96275f8d> | 2011-10-25 16:50:31 +0000 |
|---|---|---|
| committer | joseph <joseph@7b3dc134-2b1b-0410-93df-9e9f96275f8d> | 2011-10-25 16:50:31 +0000 |
| commit | c86ab84d63b20aff7cf391414009a38477fe7137 (patch) | |
| tree | eb92a78a0fd85b83317300f94a1dbed8ca414fb9 /libc | |
| parent | 4bbe4e2185c5484328182720ff7b3bb4f9593bff (diff) | |
Merge changes between r15532 and r15557 from /fsf/trunk.
git-svn-id: svn://svn.eglibc.org/trunk@15558 7b3dc134-2b1b-0410-93df-9e9f96275f8d
Diffstat (limited to 'libc')
99 files changed, 1826 insertions, 792 deletions
diff --git a/libc/ChangeLog b/libc/ChangeLog index c41dfd840..ef639398b 100644 --- a/libc/ChangeLog +++ b/libc/ChangeLog @@ -1,5 +1,179 @@ +2011-10-25 Andreas Schwab <schwab@redhat.com> + + * wcsmbs/wcscmp.c (WCSCMP): Compare as wchar_t, not wint_t. + * wcsmbs/wmemcmp.c (WMEMCMP): Likewise. + + * string/test-strchr.c (do_test): Don't generate NUL bytes. + +2011-10-25 Ulrich Drepper <drepper@gmail.com> + + * sysdeps/ieee754/dbl-64/e_atanh.c: Use math_force_eval instead of a + useful if() expression. + * sysdeps/ieee754/dbl-64/e_j0.c: Likewise. + * sysdeps/ieee754/dbl-64/s_ceil.c: Likewise. + * sysdeps/ieee754/dbl-64/s_expm1.c: Likewise. + * sysdeps/ieee754/dbl-64/s_floor.c: Likewise. + * sysdeps/ieee754/dbl-64/s_log1p.c: Likewise. + * sysdeps/ieee754/dbl-64/s_round.c: Likewise. + * sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c: Likewise. + * sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c: Likewise. + * sysdeps/ieee754/dbl-64/wordsize-64/s_round.c: Likewise. + * sysdeps/ieee754/flt-32/e_atanhf.c: Likewise. + * sysdeps/ieee754/flt-32/e_j0f.c: Likewise. + * sysdeps/ieee754/flt-32/s_ceilf.c: Likewise. + * sysdeps/ieee754/flt-32/s_expm1f.c: Likewise. + * sysdeps/ieee754/flt-32/s_floorf.c: Likewise. + * sysdeps/ieee754/flt-32/s_log1pf.c: Likewise. + * sysdeps/ieee754/flt-32/s_roundf.c: Likewise. + * sysdeps/ieee754/ldbl-96/e_atanhl.c: Likewise. + * sysdeps/ieee754/ldbl-96/e_j0l.c: Likewise. + * sysdeps/ieee754/ldbl-96/s_roundl.c: Likewise. + + * sysdeps/x86_64/fpu/math_private.h: Use VEX encoding when possible. + +2011-10-25 Andreas Schwab <schwab@redhat.com> + + * elf/dl-deps.c (_dl_map_object_deps): Remove always true + condition. + * elf/dl-fini.c (_dl_sort_fini): Likewise. + +2011-10-25 Ulrich Drepper <drepper@gmail.com> + + * sysdeps/ieee754/dbl-64/branred.c: Move FMA4 code into separate + .text section. Avoid duplicate constants. + * sysdeps/ieee754/dbl-64/doasin.c: Likewise. + * sysdeps/ieee754/dbl-64/dosincos.c: Likewise. + * sysdeps/ieee754/dbl-64/e_asin.c: Likewise. + * sysdeps/ieee754/dbl-64/e_atan2.c: Likewise. + * sysdeps/ieee754/dbl-64/e_exp.c: Likewise. + * sysdeps/ieee754/dbl-64/e_log.c: Likewise. + * sysdeps/ieee754/dbl-64/e_pow.c: Likewise. + * sysdeps/ieee754/dbl-64/halfulp.c: Likewise. + * sysdeps/ieee754/dbl-64/mpa.c: Likewise. + * sysdeps/ieee754/dbl-64/mpa.h: Likewise. + * sysdeps/ieee754/dbl-64/mpatan.c: Likewise. + * sysdeps/ieee754/dbl-64/mpatan.h: Likewise. + * sysdeps/ieee754/dbl-64/mpatan2.c: Likewise. + * sysdeps/ieee754/dbl-64/mpexp.c: Likewise. + * sysdeps/ieee754/dbl-64/mpexp.h: Likewise. + * sysdeps/ieee754/dbl-64/mpsqrt.c: Likewise. + * sysdeps/ieee754/dbl-64/mpsqrt.h: Likewise. + * sysdeps/ieee754/dbl-64/mptan.c: Likewise. + * sysdeps/ieee754/dbl-64/s_sin.c: Likewise. + * sysdeps/ieee754/dbl-64/s_tan.c: Likewise. + * sysdeps/ieee754/dbl-64/sincos32.c: Likewise. + * sysdeps/ieee754/dbl-64/slowexp.c: Likewise. + * sysdeps/ieee754/dbl-64/slowpow.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/brandred-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/doasin-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/dosincos-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/e_asin-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/e_atan2-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/e_exp-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/e_log-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/e_pow-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/halfulp-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/mpa-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/mpatan-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/mpatan2-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/mpexp-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/mplog-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/mpsqrt-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/mptan-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/s_atan-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/s_sin-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/s_tan-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/sincos32-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/slowexp-fma4.c: Likewise. + * sysdeps/x86_64/fpu/multiarch/slowpow-fma4.c: Likewise. + +2011-10-24 Ulrich Drepper <drepper@gmail.com> + + * sysdeps/x86_64/dla.h: Move to ... + * sysdeps/x86_64/fpu/dla.h: ...here. + (DLA_FMS): Some compilers fail to inline __builtin_fma in some + situations. Use __builtin_fma only for gcc 4.6 and up. + + * config.make.in: Add have-mfma4 entry. + * configure.in: Substitute libc_cv_cc_fma4. + * math/Makefile (dbl-only-routines): Add sincostab. + * sysdeps/ieee754/dbl-64/dosincos.c: Don't include sincos.tbl. + Use __sincostab not sincos. + * sysdeps/ieee754/dbl-64/e_asin.c: Don't define aliases when function + name is a macro. + * sysdeps/ieee754/dbl-64/e_exp.c: Likewise. + * sysdeps/ieee754/dbl-64/e_log.c: Likewise. + * sysdeps/ieee754/dbl-64/e_pow.c: Likewise. + * sysdeps/ieee754/dbl-64/e_atan2.c: Likewise. Define singArctan2 + using __copysign. + * sysdeps/ieee754/dbl-64/mpa.c: Don't export __acr. Don't define + __cr and __cpymn. Define __cpy unless NO___CPY is defined. Define + norm, denorm, and __mp_dbl unless NO___MP_DBL is defined. + * sysdeps/ieee754/dbl-64/mpa.h: Don't declare __acr, __cr, __cpymn, + and __inv. + * sysdeps/ieee754/dbl-64/mpsqrt.c: Make fastiroot static. + * sysdeps/ieee754/dbl-64/s_atan.c: Define __signArctan using + __copysign. + * sysdeps/ieee754/dbl-64/s_sin.c: Use __sincostab not sincos. Don't + define aliases when function name is a macro. + * sysdeps/ieee754/dbl-64/sincostab.c: Renamed from + sysdeps/ieee754/dbl-64/sincos.tbl. + * sysdeps/x86_64/fpu/multiarch/Makefile: Add entries to build + fma4-enabled routines. + * sysdeps/x86_64/fpu/multiarch/brandred-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/doasin-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/dosincos-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/e_asin-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/e_asin.c: New file. + * sysdeps/x86_64/fpu/multiarch/e_atan2-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/e_atan2.c: New file. + * sysdeps/x86_64/fpu/multiarch/e_exp-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/e_exp.c: New file. + * sysdeps/x86_64/fpu/multiarch/e_log-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/e_log.c: New file. + * sysdeps/x86_64/fpu/multiarch/e_pow-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/e_pow.c: New file. + * sysdeps/x86_64/fpu/multiarch/halfulp-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/mpa-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/mpatan-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/mpatan2-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/mpexp-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/mplog-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/mpsqrt-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/mptan-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/s_atan-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/s_atan.c: New file. + * sysdeps/x86_64/fpu/multiarch/s_sin-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/s_sin.c: New file. + * sysdeps/x86_64/fpu/multiarch/s_tan-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/s_tan.c: New file. + * sysdeps/x86_64/fpu/multiarch/sincos32-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/slowexp-fma4.c: New file. + * sysdeps/x86_64/fpu/multiarch/slowpow-fma4.c: New file. + + * sysdeps/ieee754/dbl-64/doasin.c: Adjust for DLA_FMA -> DLA_FMS + rename. + * sysdeps/ieee754/dbl-64/dosincos.c: Likewise. + * sysdeps/ieee754/dbl-64/dosincos.c: Likewise. + * sysdeps/ieee754/dbl-64/e_atan2.c: Likewise. + * sysdeps/ieee754/dbl-64/e_log.c: Likewise. + * sysdeps/ieee754/dbl-64/e_pow.c: Likewise. + * sysdeps/ieee754/dbl-64/halfulp.c: Likewise. + * sysdeps/ieee754/dbl-64/s_atan.c: Likewise. + * sysdeps/ieee754/dbl-64/s_tan.c: Likewise. + +2011-10-24 Andreas Schwab <schwab@redhat.com> + + * wcsmbs/wcslen.c: Don't define WCSLEN, reverse logic. + 2011-10-23 Ulrich Drepper <drepper@gmail.com> + * sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c: New file. + + * sysdeps/ieee754/dbl-64/e_fmod.c (__ieee754_fmod): Add some branch + prediction. + * sysdeps/ieee754/dbl-64/wordsize-64/e_fmod.c: New file. + * string/strnlen.c: Don't define STRNLEN, reverse logic. Remove unused variable magic_bits. * sysdeps/i386/i686/multiarch/rtld-strnlen.c: New file. diff --git a/libc/config.make.in b/libc/config.make.in index 2b7ae58c6..2c3bbe82d 100644 --- a/libc/config.make.in +++ b/libc/config.make.in @@ -59,6 +59,7 @@ have-cpp-asm-debuginfo = @libc_cv_cpp_asm_debuginfo@ enable-check-abi = @enable_check_abi@ have-forced-unwind = @libc_cv_forced_unwind@ have-fpie = @libc_cv_fpie@ +have-mfma4 = @libc_cv_cc_fma4@ gnu89-inline-CFLAGS = @gnu89_inline@ have-ssp = @libc_cv_ssp@ have-selinux = @have_selinux@ diff --git a/libc/configure b/libc/configure index ce163c961..fdbb52386 100755 --- a/libc/configure +++ b/libc/configure @@ -658,6 +658,7 @@ elf ldd_rewrite_script use_ldconfig libc_cv_as_i686 +libc_cv_cc_fma4 libc_cv_cc_novzeroupper libc_cv_cc_avx libc_cv_cc_sse4 @@ -8941,6 +8942,7 @@ fi + if test $elf = yes; then cat >>confdefs.h <<\_ACEOF #define HAVE_ELF 1 diff --git a/libc/configure.in b/libc/configure.in index 3877f53e3..614ed18ef 100644 --- a/libc/configure.in +++ b/libc/configure.in @@ -2368,6 +2368,7 @@ AC_SUBST(libc_cv_cpp_asm_debuginfo) AC_SUBST(libc_cv_cc_sse4) AC_SUBST(libc_cv_cc_avx) AC_SUBST(libc_cv_cc_novzeroupper) +AC_SUBST(libc_cv_cc_fma4) AC_SUBST(libc_cv_as_i686) AC_SUBST(use_ldconfig) diff --git a/libc/elf/dl-deps.c b/libc/elf/dl-deps.c index 326fc16a5..fb180ac99 100644 --- a/libc/elf/dl-deps.c +++ b/libc/elf/dl-deps.c @@ -627,7 +627,7 @@ Filters not supported with LD_TRACE_PRELINKING")); while (1) { /* Keep track of which object we looked at this round. */ - seen[i] += seen[i] < 2; + ++seen[i]; struct link_map *thisp = l_initfini[i]; /* Find the last object in the list for which the current one is diff --git a/libc/elf/dl-fini.c b/libc/elf/dl-fini.c index 953849d92..d10d89f40 100644 --- a/libc/elf/dl-fini.c +++ b/libc/elf/dl-fini.c @@ -44,7 +44,7 @@ _dl_sort_fini (struct link_map **maps, size_t nmaps, char *used, Lmid_t ns) while (1) { /* Keep track of which object we looked at this round. */ - seen[i] += seen[i] < 2; + ++seen[i]; struct link_map *thisp = maps[i]; /* Do not handle ld.so in secondary namespaces and object which diff --git a/libc/math/Makefile b/libc/math/Makefile index 95fd6b066..10f14be8c 100644 --- a/libc/math/Makefile +++ b/libc/math/Makefile @@ -68,7 +68,7 @@ include ../Makeconfig dbl-only-routines := branred doasin dosincos halfulp mpa mpatan2 \ mpatan mpexp mplog mpsqrt mptan sincos32 slowexp \ - slowpow + slowpow sincostab libm-routines = $(strip $(libm-support) $(libm-calls) \ $(patsubst %_rf,%f_r,$(libm-calls:=f)) \ $(long-m-$(long-double-fcts))) \ diff --git a/libc/nptl/ChangeLog b/libc/nptl/ChangeLog index 7b0676381..f4b2f04c0 100644 --- a/libc/nptl/ChangeLog +++ b/libc/nptl/ChangeLog @@ -1,3 +1,10 @@ +2011-10-24 Ulrich Drepper <drepper@gmail.com> + + * tst-cancel7.c: Avoid warning. + * tst-mutex6.c: Likewise. + * tst-mutex9.c: Likewise. + * tst-mutexpi6.c: Likewise. + 2011-10-23 Ulrich Drepper <drepper@gmail.com> * sysdeps/i386/tls.h: Remove #include <list.h>. diff --git a/libc/nptl/tst-cancel7.c b/libc/nptl/tst-cancel7.c index be9b1c606..af0d18f1e 100644 --- a/libc/nptl/tst-cancel7.c +++ b/libc/nptl/tst-cancel7.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2002, 2003 Free Software Foundation, Inc. +/* Copyright (C) 2002, 2003, 2011 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jakub@redhat.com>, 2002. @@ -204,7 +204,7 @@ do_cleanup (void) case OPT_PIDFILE: \ pidfile = optarg; \ break; -// #define CLEANUP_HANDLER do_cleanup () +#define CLEANUP_HANDLER do_cleanup () #define PREPARE(argc, argv) do_prepare (argc, argv) #define TEST_FUNCTION do_test () #define TIMEOUT 5 diff --git a/libc/nptl/tst-mutex6.c b/libc/nptl/tst-mutex6.c index de64bdb43..19611ee94 100644 --- a/libc/nptl/tst-mutex6.c +++ b/libc/nptl/tst-mutex6.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2002, 2006 Free Software Foundation, Inc. +/* Copyright (C) 2002, 2006, 2011 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@redhat.com>, 2002. @@ -25,7 +25,8 @@ #ifndef ATTR -# define ATTR NULL +pthread_mutexattr_t *attr; +# define ATTR attr #endif diff --git a/libc/nptl/tst-mutex9.c b/libc/nptl/tst-mutex9.c index f9d379343..bdf1dc842 100644 --- a/libc/nptl/tst-mutex9.c +++ b/libc/nptl/tst-mutex9.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2003, 2004, 2006 Free Software Foundation, Inc. +/* Copyright (C) 2003, 2004, 2006, 2011 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@redhat.com>, 2003. @@ -37,7 +37,6 @@ do_test (void) pthread_mutex_t *m; pthread_mutexattr_t a; pid_t pid; - char *p; fd = mkstemp (tmpfname); if (fd == -1) @@ -68,7 +67,6 @@ do_test (void) m = (pthread_mutex_t *) (((uintptr_t) mem + __alignof (pthread_mutex_t)) & ~(__alignof (pthread_mutex_t) - 1)); - p = (char *) (m + 1); if (pthread_mutexattr_init (&a) != 0) { diff --git a/libc/nptl/tst-mutexpi6.c b/libc/nptl/tst-mutexpi6.c index 42cda377d..8881a1d2c 100644 --- a/libc/nptl/tst-mutexpi6.c +++ b/libc/nptl/tst-mutexpi6.c @@ -3,11 +3,13 @@ #include <stdlib.h> -static pthread_mutexattr_t a; +pthread_mutexattr_t a; +pthread_mutexattr_t *attr; static void prepare (void) { + attr = &a; if (pthread_mutexattr_init (&a) != 0) { puts ("mutexattr_init failed"); @@ -23,5 +25,5 @@ prepare (void) #define PREPARE(argc, argv) prepare () -#define ATTR &a +#define ATTR attr #include "tst-mutex6.c" diff --git a/libc/string/test-strchr.c b/libc/string/test-strchr.c index 518a4dce3..3bbc2f570 100644 --- a/libc/string/test-strchr.c +++ b/libc/string/test-strchr.c @@ -63,8 +63,8 @@ stupid_STRCHR (const CHAR *s, int c) return NULL; } -IMPL (stupid_STRCHR, 1) -IMPL (simple_STRCHR, 1) +IMPL (stupid_STRCHR, 0) +IMPL (simple_STRCHR, 0) IMPL (STRCHR, 1) static void @@ -100,15 +100,15 @@ do_one_test (impl_t *impl, const CHAR *s, int c, const CHAR *exp_res) static void do_test (size_t align, size_t pos, size_t len, int seek_char, int max_char) -/* for wcschr: align here means align not in bytes, - * but in wchar_ts, in bytes it will equal to align * (sizeof (wchar_t)) - * len for wcschr here isn't in bytes but it's number of wchar_t symbols */ +/* For wcschr: align here means align not in bytes, + but in wchar_ts, in bytes it will equal to align * (sizeof (wchar_t)) + len for wcschr here isn't in bytes but it's number of wchar_t symbols. */ { size_t i; CHAR *result; - CHAR *buf = (CHAR *) (buf1); + CHAR *buf = (CHAR *) buf1; align &= 15; - if ((align + len) * sizeof(CHAR) >= page_size) + if ((align + len) * sizeof (CHAR) >= page_size) return; for (i = 0; i < len; ++i) @@ -116,6 +116,8 @@ do_test (size_t align, size_t pos, size_t len, int seek_char, int max_char) buf[align + i] = 32 + 23 * i % max_char; if (buf[align + i] == seek_char) buf[align + i] = seek_char + 1; + else if (buf[align + i] == 0) + buf[align + i] = 1; } buf[align + len] = 0; @@ -130,7 +132,8 @@ do_test (size_t align, size_t pos, size_t len, int seek_char, int max_char) result = NULL; if (HP_TIMING_AVAIL) - printf ("Length %4zd, alignment in bytes %2zd:", pos, align * sizeof (CHAR)); + printf ("Length %4zd, alignment in bytes %2zd:", + pos, align * sizeof (CHAR)); FOR_EACH_IMPL (impl, 0) do_one_test (impl, buf + align, seek_char, result); @@ -149,14 +152,15 @@ do_random_tests (void) for (n = 0; n < ITERATIONS; n++) { -/* for wcschr: align here means align not in bytes, but in wchar_ts, - * in bytes it will equal to align * (sizeof (wchar_t)) */ + /* For wcschr: align here means align not in bytes, but in wchar_ts, + in bytes it will equal to align * (sizeof (wchar_t)). */ align = random () & 15; pos = random () & 511; seek_char = random () & 255; if (pos + align >= 511) pos = 510 - align - (random () & 7); -/* len for wcschr here isn't in bytes but it's number of wchar_t symbols */ + /* len for wcschr here isn't in bytes but it's number of wchar_t + symbols. */ len = random () & 511; if ((pos == len && seek_char) || (pos > len && (random () & 1))) diff --git a/libc/sysdeps/ieee754/dbl-64/branred.c b/libc/sysdeps/ieee754/dbl-64/branred.c index 76015f0c5..c8483034a 100644 --- a/libc/sysdeps/ieee754/dbl-64/branred.c +++ b/libc/sysdeps/ieee754/dbl-64/branred.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation, Inc. + * Copyright (C) 2001, 2011 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -38,6 +38,10 @@ #include "branred.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + /*******************************************************************/ /* Routine branred() performs range reduction of a double number */ @@ -45,7 +49,9 @@ /* x=n*pi/2+(a+aa), abs(a+aa)<pi/4, n=0,+-1,+-2,.... */ /* Routine return integer (n mod 4) */ /*******************************************************************/ -int __branred(double x, double *a, double *aa) +int +SECTION +__branred(double x, double *a, double *aa) { int i,k; #if 0 diff --git a/libc/sysdeps/ieee754/dbl-64/doasin.c b/libc/sysdeps/ieee754/dbl-64/doasin.c index c21d4b7df..14958b5ca 100644 --- a/libc/sysdeps/ieee754/dbl-64/doasin.c +++ b/libc/sysdeps/ieee754/dbl-64/doasin.c @@ -34,11 +34,17 @@ #include <dla.h> #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + /********************************************************************/ /* Compute arcsin(x,dx,v) of double-length number (x+dx) the result */ /* stored in v where v= v[0]+v[1] =arcsin(x+dx) */ /********************************************************************/ -void __doasin(double x, double dx, double v[]) { +void +SECTION +__doasin(double x, double dx, double v[]) { #include "doasin.h" @@ -53,7 +59,7 @@ void __doasin(double x, double dx, double v[]) { double xx,p,pp,u,uu,r,s; double tc,tcc; -#ifndef DLA_FMA +#ifndef DLA_FMS double hx,tx,hy,ty,tp,tq; #endif diff --git a/libc/sysdeps/ieee754/dbl-64/dosincos.c b/libc/sysdeps/ieee754/dbl-64/dosincos.c index 4ae88c31c..e8890ff8d 100644 --- a/libc/sysdeps/ieee754/dbl-64/dosincos.c +++ b/libc/sysdeps/ieee754/dbl-64/dosincos.c @@ -35,11 +35,20 @@ #include "endian.h" #include "mydefs.h" -#include "sincos.tbl" #include <dla.h> #include "dosincos.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + +extern const union +{ + int4 i[880]; + double x[440]; +} __sincostab attribute_hidden; + /***********************************************************************/ /* Routine receive Double-Length number (x+dx) and computing sin(x+dx) */ /* as Double-Length number and store it at array v .It computes it by */ @@ -47,10 +56,12 @@ /*(x+dx) between 0 and PI/4 */ /***********************************************************************/ -void __dubsin(double x, double dx, double v[]) { +void +SECTION +__dubsin(double x, double dx, double v[]) { double r,s,c,cc,d,dd,d2,dd2,e,ee, sn,ssn,cs,ccs,ds,dss,dc,dcc; -#ifndef DLA_FMA +#ifndef DLA_FMS double p,hx,tx,hy,ty,q; #endif #if 0 @@ -66,10 +77,10 @@ void __dubsin(double x, double dx, double v[]) { dd=(x-d)+dx; /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */ MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc); - sn=sincos.x[k]; /* */ - ssn=sincos.x[k+1]; /* sin(Xi) and cos(Xi) */ - cs=sincos.x[k+2]; /* */ - ccs=sincos.x[k+3]; /* */ + sn=__sincostab.x[k]; /* */ + ssn=__sincostab.x[k+1]; /* sin(Xi) and cos(Xi) */ + cs=__sincostab.x[k+2]; /* */ + ccs=__sincostab.x[k+3]; /* */ MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* Taylor */ ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s); MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); /* series */ @@ -101,10 +112,12 @@ void __dubsin(double x, double dx, double v[]) { /*(x+dx) between 0 and PI/4 */ /**********************************************************************/ -void __dubcos(double x, double dx, double v[]) { +void +SECTION +__dubcos(double x, double dx, double v[]) { double r,s,c,cc,d,dd,d2,dd2,e,ee, sn,ssn,cs,ccs,ds,dss,dc,dcc; -#ifndef DLA_FMA +#ifndef DLA_FMS double p,hx,tx,hy,ty,q; #endif #if 0 @@ -118,10 +131,10 @@ void __dubcos(double x, double dx, double v[]) { d=x+dx; dd=(x-d)+dx; /* cos(x+dx)=cos(Xi+t)=cos(Xi)cos(t) - sin(Xi)sin(t) */ MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc); - sn=sincos.x[k]; /* */ - ssn=sincos.x[k+1]; /* sin(Xi) and cos(Xi) */ - cs=sincos.x[k+2]; /* */ - ccs=sincos.x[k+3]; /* */ + sn=__sincostab.x[k]; /* */ + ssn=__sincostab.x[k+1]; /* sin(Xi) and cos(Xi) */ + cs=__sincostab.x[k+2]; /* */ + ccs=__sincostab.x[k+3]; /* */ MUL2(d2,dd2,s7.x,ss7.x,ds,dss,p,hx,tx,hy,ty,q,c,cc); ADD2(ds,dss,s5.x,ss5.x,ds,dss,r,s); MUL2(d2,dd2,ds,dss,ds,dss,p,hx,tx,hy,ty,q,c,cc); @@ -167,7 +180,9 @@ void __dubcos(double x, double dx, double v[]) { /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */ /* as Double-Length number and store it in array v */ /**********************************************************************/ -void __docos(double x, double dx, double v[]) { +void +SECTION +__docos(double x, double dx, double v[]) { double y,yy,p,w[2]; if (x>0) {y=x; yy=dx;} else {y=-x; yy=-dx;} diff --git a/libc/sysdeps/ieee754/dbl-64/e_asin.c b/libc/sysdeps/ieee754/dbl-64/e_asin.c index 02efb7ad2..65319c0b5 100644 --- a/libc/sysdeps/ieee754/dbl-64/e_asin.c +++ b/libc/sysdeps/ieee754/dbl-64/e_asin.c @@ -42,6 +42,10 @@ #include "uasncs.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + void __doasin(double x, double dx, double w[]); void __dubsin(double x, double dx, double v[]); void __dubcos(double x, double dx, double v[]); @@ -53,7 +57,9 @@ double __cos32(double x, double res, double res1); /* An ultimate asin routine. Given an IEEE double machine number x */ /* it computes the correctly rounded (to nearest) value of arcsin(x) */ /***************************************************************************/ -double __ieee754_asin(double x){ +double +SECTION +__ieee754_asin(double x){ double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2]; mynumber u,v; int4 k,m,n; @@ -324,7 +330,9 @@ double __ieee754_asin(double x){ return u.x/v.x; /* NaN */ } } +#ifndef __ieee754_asin strong_alias (__ieee754_asin, __asin_finite) +#endif /*******************************************************************/ /* */ @@ -332,7 +340,9 @@ strong_alias (__ieee754_asin, __asin_finite) /* */ /*******************************************************************/ -double __ieee754_acos(double x) +double +SECTION +__ieee754_acos(double x) { double x1,x2,xx,s1,s2,res1,p,t,res,r,cor,cc,y,c,z,w[2],eps; #if 0 @@ -636,4 +646,6 @@ double __ieee754_acos(double x) return u.x/v.x; } } +#ifndef __ieee754_acos strong_alias (__ieee754_acos, __acos_finite) +#endif diff --git a/libc/sysdeps/ieee754/dbl-64/e_atan2.c b/libc/sysdeps/ieee754/dbl-64/e_atan2.c index f8f678bc5..64dae3e8d 100644 --- a/libc/sysdeps/ieee754/dbl-64/e_atan2.c +++ b/libc/sysdeps/ieee754/dbl-64/e_atan2.c @@ -44,6 +44,10 @@ #include "atnat2.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + /************************************************************************/ /* An ultimate atan2 routine. Given two IEEE double machine numbers y,x */ /* it computes the correctly rounded (to nearest) value of atan2(y,x). */ @@ -51,11 +55,17 @@ /* round to nearest mode of IEEE 754 standard. */ /************************************************************************/ static double atan2Mp(double ,double ,const int[]); -static double signArctan2(double ,double); + /* Fix the sign and return after stage 1 or stage 2 */ +static double signArctan2(double y,double z) +{ + return __copysign(z, y); +} static double normalized(double ,double,double ,double); void __mpatan2(mp_no *,mp_no *,mp_no *,int); -double __ieee754_atan2(double y,double x) { +double +SECTION +__ieee754_atan2(double y,double x) { int i,de,ux,dx,uy,dy; #if 0 @@ -64,7 +74,7 @@ double __ieee754_atan2(double y,double x) { static const int pr[MM]={6,8,10,20,32}; double ax,ay,u,du,u9,ua,v,vv,dv,t1,t2,t3,t7,t8, z,zz,cor,s1,ss1,s2,ss2; -#ifndef DLA_FMA +#ifndef DLA_FMS double t4,t5,t6; #endif #if 0 @@ -375,10 +385,14 @@ double __ieee754_atan2(double y,double x) { } } } +#ifndef __ieee754_atan2 strong_alias (__ieee754_atan2, __atan2_finite) +#endif /* Treat the Denormalized case */ -static double normalized(double ax,double ay,double y, double z) +static double +SECTION +normalized(double ax,double ay,double y, double z) { int p; mp_no mpx,mpy,mpz,mperr,mpz2,mpt1; p=6; @@ -387,13 +401,10 @@ static double normalized(double ax,double ay,double y, double z) __sub(&mpz,&mperr,&mpz2,p); __mp_dbl(&mpz2,&z,p); return signArctan2(y,z); } - /* Fix the sign and return after stage 1 or stage 2 */ -static double signArctan2(double y,double z) -{ - return ((y<ZERO) ? -z : z); -} /* Stage 3: Perform a multi-Precision computation */ -static double atan2Mp(double x,double y,const int pr[]) +static double +SECTION +atan2Mp(double x,double y,const int pr[]) { double z1,z2; int i,p; diff --git a/libc/sysdeps/ieee754/dbl-64/e_atanh.c b/libc/sysdeps/ieee754/dbl-64/e_atanh.c index de3bc8f14..1f83e3198 100644 --- a/libc/sysdeps/ieee754/dbl-64/e_atanh.c +++ b/libc/sysdeps/ieee754/dbl-64/e_atanh.c @@ -49,8 +49,11 @@ __ieee754_atanh (double x) double t; if (xa < 0.5) { - if (__builtin_expect (xa < 0x1.0p-28, 0) && (huge + x) > 0.0) - return x; + if (__builtin_expect (xa < 0x1.0p-28, 0)) + { + math_force_eval (huge + x); + return x; + } t = xa + xa; t = 0.5 * __log1p (t + t * xa / (1.0 - xa)); diff --git a/libc/sysdeps/ieee754/dbl-64/e_exp.c b/libc/sysdeps/ieee754/dbl-64/e_exp.c index f4b34a636..e7a839d42 100644 --- a/libc/sysdeps/ieee754/dbl-64/e_exp.c +++ b/libc/sysdeps/ieee754/dbl-64/e_exp.c @@ -40,13 +40,19 @@ #include "uexp.tbl" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + double __slowexp(double); /***************************************************************************/ /* An ultimate exp routine. Given an IEEE double machine number x */ /* it computes the correctly rounded (to nearest) value of e^x */ /***************************************************************************/ -double __ieee754_exp(double x) { +double +SECTION +__ieee754_exp(double x) { double bexp, t, eps, del, base, y, al, bet, res, rem, cor; mynumber junk1, junk2, binexp = {{0,0}}; #if 0 @@ -145,7 +151,9 @@ double __ieee754_exp(double x) { else return __slowexp(x); } } +#ifndef __ieee754_exp strong_alias (__ieee754_exp, __exp_finite) +#endif /************************************************************************/ /* Compute e^(x+xx)(Double-Length number) .The routine also receive */ @@ -154,7 +162,9 @@ strong_alias (__ieee754_exp, __exp_finite) /*else return e^(x + xx) (always positive ) */ /************************************************************************/ -double __exp1(double x, double xx, double error) { +double +SECTION +__exp1(double x, double xx, double error) { double bexp, t, eps, del, base, y, al, bet, res, rem, cor; mynumber junk1, junk2, binexp = {{0,0}}; #if 0 diff --git a/libc/sysdeps/ieee754/dbl-64/e_fmod.c b/libc/sysdeps/ieee754/dbl-64/e_fmod.c index a575f616b..0328b011d 100644 --- a/libc/sysdeps/ieee754/dbl-64/e_fmod.c +++ b/libc/sysdeps/ieee754/dbl-64/e_fmod.c @@ -1,4 +1,3 @@ -/* @(#)e_fmod.c 5.1 93/09/24 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. @@ -44,7 +43,7 @@ __ieee754_fmod (double x, double y) } /* determine ix = ilogb(x) */ - if(hx<0x00100000) { /* subnormal x */ + if(__builtin_expect(hx<0x00100000, 0)) { /* subnormal x */ if(hx==0) { for (ix = -1043, i=lx; i>0; i<<=1) ix -=1; } else { @@ -53,7 +52,7 @@ __ieee754_fmod (double x, double y) } else ix = (hx>>20)-1023; /* determine iy = ilogb(y) */ - if(hy<0x00100000) { /* subnormal y */ + if(__builtin_expect(hy<0x00100000, 0)) { /* subnormal y */ if(hy==0) { for (iy = -1043, i=ly; i>0; i<<=1) iy -=1; } else { @@ -62,7 +61,7 @@ __ieee754_fmod (double x, double y) } else iy = (hy>>20)-1023; /* set up {hx,lx}, {hy,ly} and align y to x */ - if(ix >= -1022) + if(__builtin_expect(ix >= -1022, 1)) hx = 0x00100000|(0x000fffff&hx); else { /* subnormal x, shift x to normal */ n = -1022-ix; @@ -74,7 +73,7 @@ __ieee754_fmod (double x, double y) lx = 0; } } - if(iy >= -1022) + if(__builtin_expect(iy >= -1022, 1)) hy = 0x00100000|(0x000fffff&hy); else { /* subnormal y, shift y to normal */ n = -1022-iy; @@ -108,7 +107,7 @@ __ieee754_fmod (double x, double y) hx = hx+hx+(lx>>31); lx = lx+lx; iy -= 1; } - if(iy>= -1022) { /* normalize output */ + if(__builtin_expect(iy>= -1022, 1)) { /* normalize output */ hx = ((hx-0x00100000)|((iy+1023)<<20)); INSERT_WORDS(x,hx|sx,lx); } else { /* subnormal output */ diff --git a/libc/sysdeps/ieee754/dbl-64/e_j0.c b/libc/sysdeps/ieee754/dbl-64/e_j0.c index 5ebf2056b..48584a60b 100644 --- a/libc/sysdeps/ieee754/dbl-64/e_j0.c +++ b/libc/sysdeps/ieee754/dbl-64/e_j0.c @@ -111,10 +111,9 @@ __ieee754_j0(double x) return z; } if(ix<0x3f200000) { /* |x| < 2**-13 */ - if(huge+x>one) { /* raise inexact if x != 0 */ - if(ix<0x3e400000) return one; /* |x|<2**-27 */ - else return one - 0.25*x*x; - } + math_force_eval(huge+x); /* raise inexact if x != 0 */ + if(ix<0x3e400000) return one; /* |x|<2**-27 */ + else return one - 0.25*x*x; } z = x*x; #ifdef DO_NOT_USE_THIS diff --git a/libc/sysdeps/ieee754/dbl-64/e_log.c b/libc/sysdeps/ieee754/dbl-64/e_log.c index 14851638a..e45520eba 100644 --- a/libc/sysdeps/ieee754/dbl-64/e_log.c +++ b/libc/sysdeps/ieee754/dbl-64/e_log.c @@ -41,13 +41,19 @@ #include "MathLib.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + void __mplog(mp_no *, mp_no *, int); /*********************************************************************/ /* An ultimate log routine. Given an IEEE double machine number x */ /* it computes the correctly rounded (to nearest) value of log(x). */ /*********************************************************************/ -double __ieee754_log(double x) { +double +SECTION +__ieee754_log(double x) { #define M 4 static const int pr[M]={8,10,18,32}; int i,j,n,ux,dx,p; @@ -58,7 +64,7 @@ double __ieee754_log(double x) { sij,ssij,ttij,A,B,B0,y,y1,y2,polI,polII,sa,sb, t1,t2,t7,t8,t,ra,rb,ww, a0,aa0,s1,s2,ss2,s3,ss3,a1,aa1,a,aa,b,bb,c; -#ifndef DLA_FMA +#ifndef DLA_FMS double t3,t4,t5,t6; #endif number num; @@ -207,4 +213,6 @@ double __ieee754_log(double x) { } return y1; } +#ifndef __ieee754_log strong_alias (__ieee754_log, __log_finite) +#endif diff --git a/libc/sysdeps/ieee754/dbl-64/e_pow.c b/libc/sysdeps/ieee754/dbl-64/e_pow.c index 789054015..350e93986 100644 --- a/libc/sysdeps/ieee754/dbl-64/e_pow.c +++ b/libc/sysdeps/ieee754/dbl-64/e_pow.c @@ -43,6 +43,10 @@ #include "upow.tbl" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + double __exp1(double x, double xx, double error); static double log1(double x, double *delta, double *error); @@ -55,7 +59,9 @@ static int checkint(double x); /* An ultimate power routine. Given two IEEE double machine numbers y,x */ /* it computes the correctly rounded (to nearest) value of X^y. */ /***************************************************************************/ -double __ieee754_pow(double x, double y) { +double +SECTION +__ieee754_pow(double x, double y) { double z,a,aa,error, t,a1,a2,y1,y2; #if 0 double gor=1.0; @@ -153,12 +159,16 @@ double __ieee754_pow(double x, double y) { if (y<0) return (x<1.0)?INF.x:0; return 0; /* unreachable, to make the compiler happy */ } +#ifndef __ieee754_pow strong_alias (__ieee754_pow, __pow_finite) +#endif /**************************************************************************/ /* Computing x^y using more accurate but more slow log routine */ /**************************************************************************/ -static double power1(double x, double y) { +static double +SECTION +power1(double x, double y) { double z,a,aa,error, t,a1,a2,y1,y2; z = my_log2(x,&aa,&error); t = y*134217729.0; @@ -181,7 +191,9 @@ static double power1(double x, double y) { /* + the parameter delta. */ /* The result is bounded by error (rightmost argument) */ /****************************************************************************/ -static double log1(double x, double *delta, double *error) { +static double +SECTION +log1(double x, double *delta, double *error) { int i,j,m; #if 0 int n; @@ -273,7 +285,9 @@ static double log1(double x, double *delta, double *error) { /* Computing log(x)(x is left argument).The result is return double + delta.*/ /* The result is bounded by error (right argument) */ /****************************************************************************/ -static double my_log2(double x, double *delta, double *error) { +static double +SECTION +my_log2(double x, double *delta, double *error) { int i,j,m; #if 0 int n; @@ -284,7 +298,7 @@ static double my_log2(double x, double *delta, double *error) { #endif double ou1,ou2,lu1,lu2,ov,lv1,lv2,a,a1,a2; double y,yy,z,zz,j1,j2,j7,j8; -#ifndef DLA_FMA +#ifndef DLA_FMS double j3,j4,j5,j6; #endif mynumber u,v; @@ -367,7 +381,9 @@ static double my_log2(double x, double *delta, double *error) { /* Routine receives a double x and checks if it is an integer. If not */ /* it returns 0, else it returns 1 if even or -1 if odd. */ /**********************************************************************/ -static int checkint(double x) { +static int +SECTION +checkint(double x) { union {int4 i[2]; double x;} u; int k,m,n; #if 0 diff --git a/libc/sysdeps/ieee754/dbl-64/halfulp.c b/libc/sysdeps/ieee754/dbl-64/halfulp.c index 373d40522..601830942 100644 --- a/libc/sysdeps/ieee754/dbl-64/halfulp.c +++ b/libc/sysdeps/ieee754/dbl-64/halfulp.c @@ -40,6 +40,10 @@ #include <dla.h> #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + static const int4 tab54[32] = { 262143, 11585, 1782, 511, 210, 107, 63, 42, 30, 22, 17, 14, 12, 10, 9, 7, @@ -47,11 +51,13 @@ static const int4 tab54[32] = { 3, 3, 3, 3, 3, 3, 3, 3 }; -double __halfulp(double x, double y) +double +SECTION +__halfulp(double x, double y) { mynumber v; double z,u,uu; -#ifndef DLA_FMA +#ifndef DLA_FMS double j1,j2,j3,j4,j5; #endif int4 k,l,m,n; diff --git a/libc/sysdeps/ieee754/dbl-64/mpa.c b/libc/sysdeps/ieee754/dbl-64/mpa.c index 68647ba33..39c640882 100644 --- a/libc/sysdeps/ieee754/dbl-64/mpa.c +++ b/libc/sysdeps/ieee754/dbl-64/mpa.c @@ -1,8 +1,7 @@ - /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -48,11 +47,18 @@ #include "mpa.h" #include "mpa2.h" #include <sys/param.h> /* For MIN() */ + +#ifndef SECTION +# define SECTION +#endif + +#ifndef NO___ACR /* mcr() compares the sizes of the mantissas of two multiple precision */ /* numbers. Mantissas are compared regardless of the signs of the */ /* numbers, even if x->d[0] or y->d[0] are zero. Exponents are also */ /* disregarded. */ -static int mcr(const mp_no *x, const mp_no *y, int p) { +static int +mcr(const mp_no *x, const mp_no *y, int p) { int i; for (i=1; i<=p; i++) { if (X[i] == Y[i]) continue; @@ -62,9 +68,9 @@ static int mcr(const mp_no *x, const mp_no *y, int p) { } - /* acr() compares the absolute values of two multiple precision numbers */ -int __acr(const mp_no *x, const mp_no *y, int p) { +int +__acr(const mp_no *x, const mp_no *y, int p) { int i; if (X[0] == ZERO) { @@ -80,10 +86,12 @@ int __acr(const mp_no *x, const mp_no *y, int p) { return i; } +#endif -/* cr90 compares the values of two multiple precision numbers */ -int __cr(const mp_no *x, const mp_no *y, int p) { +#if 0 +/* cr() compares the values of two multiple precision numbers */ +static int __cr(const mp_no *x, const mp_no *y, int p) { int i; if (X[0] > Y[0]) i= 1; @@ -93,36 +101,37 @@ int __cr(const mp_no *x, const mp_no *y, int p) { return i; } +#endif +#ifndef NO___CPY /* Copy a multiple precision number. Set *y=*x. x=y is permissible. */ void __cpy(const mp_no *x, mp_no *y, int p) { - int i; - EY = EX; - for (i=0; i <= p; i++) Y[i] = X[i]; - - return; + for (int i=0; i <= p; i++) Y[i] = X[i]; } +#endif +#if 0 /* Copy a multiple precision number x of precision m into a */ /* multiple precision number y of precision n. In case n>m, */ /* the digits of y beyond the m'th are set to zero. In case */ /* n<m, the digits of x beyond the n'th are ignored. */ /* x=y is permissible. */ -void __cpymn(const mp_no *x, int m, mp_no *y, int n) { +static void __cpymn(const mp_no *x, int m, mp_no *y, int n) { int i,k; EY = EX; k=MIN(m,n); for (i=0; i <= k; i++) Y[i] = X[i]; for ( ; i <= n; i++) Y[i] = ZERO; - - return; } +#endif + +#ifndef NO___MP_DBL /* Convert a multiple precision number *x into a double precision */ /* number *y, normalized case (|x| >= 2**(-1022))) */ static void norm(const mp_no *x, double *y, int p) @@ -141,7 +150,7 @@ static void norm(const mp_no *x, double *y, int p) } else { for (a=ONE, z[1]=X[1]; z[1] < TWO23; ) - {a *= TWO; z[1] *= TWO; } + {a *= TWO; z[1] *= TWO; } for (i=2; i<5; i++) { z[i] = X[i]*a; @@ -157,10 +166,10 @@ static void norm(const mp_no *x, double *y, int p) if (v == TWO18) { if (z[4] == ZERO) { - for (i=5; i <= p; i++) { - if (X[i] == ZERO) continue; - else {z[3] += ONE; break; } - } + for (i=5; i <= p; i++) { + if (X[i] == ZERO) continue; + else {z[3] += ONE; break; } + } } else z[3] += ONE; } @@ -174,7 +183,6 @@ static void norm(const mp_no *x, double *y, int p) for (i=1; i>EX; i--) c *= RADIXI; *y = c; - return; #undef R } @@ -222,8 +230,6 @@ static void denorm(const mp_no *x, double *y, int p) c = X[0]*((z[1] + R*(z[2] + R*z[3])) - TWO10); *y = c*TWOM1032; - return; - #undef R } @@ -242,13 +248,16 @@ void __mp_dbl(const mp_no *x, double *y, int p) { else if (EX==-42 && X[1]>=TWO10) norm(x,y,p); else denorm(x,y,p); } +#endif /* dbl_mp() converts a double precision number x into a multiple precision */ /* number *y. If the precision p is too small the result is truncated. x is */ /* left unchanged. */ -void __dbl_mp(double x, mp_no *y, int p) { +void +SECTION +__dbl_mp(double x, mp_no *y, int p) { int i,n; double u; @@ -269,7 +278,6 @@ void __dbl_mp(double x, mp_no *y, int p) { if (u>x) u -= ONE; Y[i] = u; x -= u; x *= RADIX; } for ( ; i<=p; i++) Y[i] = ZERO; - return; } @@ -279,7 +287,9 @@ void __dbl_mp(double x, mp_no *y, int p) { /* No guard digit is used. The result equals the exact sum, truncated. */ /* *x & *y are left unchanged. */ -static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { +static void +SECTION +add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { int i,j,k; @@ -321,7 +331,9 @@ static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { /* or y&z. One guard digit is used. The error is less than one ulp. */ /* *x & *y are left unchanged. */ -static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { +static void +SECTION +sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { int i,j,k; @@ -336,11 +348,11 @@ static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { else { i=p; j=p+1-j; k=p; if (Y[j] > ZERO) { - Z[k+1] = RADIX - Y[j--]; - Z[k] = MONE; } + Z[k+1] = RADIX - Y[j--]; + Z[k] = MONE; } else { - Z[k+1] = ZERO; - Z[k] = ZERO; j--;} + Z[k+1] = ZERO; + Z[k] = ZERO; j--;} } } @@ -368,8 +380,6 @@ static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { Z[k++] = Z[i++]; for (; k <= p; ) Z[k++] = ZERO; - - return; } @@ -377,7 +387,9 @@ static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { /* but not x&z or y&z. One guard digit is used. The error is less than */ /* one ulp. *x & *y are left unchanged. */ -void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) { +void +SECTION +__add(const mp_no *x, const mp_no *y, mp_no *z, int p) { int n; @@ -393,7 +405,6 @@ void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) { else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = Y[0]; } else Z[0] = ZERO; } - return; } @@ -401,7 +412,9 @@ void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) { /* overlap but not x&z or y&z. One guard digit is used. The error is */ /* less than one ulp. *x & *y are left unchanged. */ -void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) { +void +SECTION +__sub(const mp_no *x, const mp_no *y, mp_no *z, int p) { int n; @@ -417,7 +430,6 @@ void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) { else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = -Y[0]; } else Z[0] = ZERO; } - return; } @@ -426,16 +438,18 @@ void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) { /* truncated to p digits. In case p>3 the error is bounded by 1.001 ulp. */ /* *x & *y are left unchanged. */ -void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { +void +SECTION +__mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { int i, i1, i2, j, k, k2; double u; - /* Is z=0? */ + /* Is z=0? */ if (X[0]*Y[0]==ZERO) { Z[0]=ZERO; return; } - /* Multiply, add and carry */ + /* Multiply, add and carry */ k2 = (p<3) ? p+p : p+3; Z[k2]=ZERO; for (k=k2; k>1; ) { @@ -449,7 +463,7 @@ void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { Z[--k] = u*RADIXI; } - /* Is there a carry beyond the most significant digit? */ + /* Is there a carry beyond the most significant digit? */ if (Z[1] == ZERO) { for (i=1; i<=p; i++) Z[i]=Z[i+1]; EZ = EX + EY - 1; } @@ -457,7 +471,6 @@ void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { EZ = EX + EY; Z[0] = X[0] * Y[0]; - return; } @@ -466,6 +479,8 @@ void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { /* 2.001*r**(1-p) for p>3. */ /* *x=0 is not permissible. *x is left unchanged. */ +static +SECTION void __inv(const mp_no *x, mp_no *y, int p) { int i; #if 0 @@ -474,11 +489,11 @@ void __inv(const mp_no *x, mp_no *y, int p) { double t; mp_no z,w; static const int np1[] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,3, - 4,4,4,4,4,4,4,4,4,4,4,4,4,4,4}; + 4,4,4,4,4,4,4,4,4,4,4,4,4,4,4}; const mp_no mptwo = {1,{1.0,2.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; __cpy(x,&z,p); z.e=0; __mp_dbl(&z,&t,p); t=ONE/t; __dbl_mp(t,y,p); EY -= EX; @@ -489,7 +504,6 @@ void __inv(const mp_no *x, mp_no *y, int p) { __sub(&mptwo,y,&z,p); __mul(&w,&z,y,p); } - return; } @@ -498,11 +512,12 @@ void __inv(const mp_no *x, mp_no *y, int p) { /* Relative error bound = 2.001*r**(1-p) for p=2, 2.063*r**(1-p) for p=3 */ /* and 3.001*r**(1-p) for p>3. *y=0 is not permissible. */ -void __dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) { +void +SECTION +__dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) { mp_no w; if (X[0] == ZERO) Z[0] = ZERO; else {__inv(y,&w,p); __mul(x,&w,z,p);} - return; } diff --git a/libc/sysdeps/ieee754/dbl-64/mpa.h b/libc/sysdeps/ieee754/dbl-64/mpa.h index 4aec48e90..5647ab7b4 100644 --- a/libc/sysdeps/ieee754/dbl-64/mpa.h +++ b/libc/sysdeps/ieee754/dbl-64/mpa.h @@ -1,8 +1,7 @@ - /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation, Inc. + * Copyright (C) 2001, 2011 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -45,14 +44,14 @@ typedef struct {/* This structure holds the details of a multi-precision */ int e; /* floating point number, x: d[0] holds its sign (-1,0 or 1) */ double d[40]; /* e holds its exponent (...,-2,-1,0,1,2,...) and */ } mp_no; /* d[1]...d[p] hold its mantissa digits. The value of x is, */ - /* x = d[1]*r**(e-1) + d[2]*r**(e-2) + ... + d[p]*r**(e-p). */ - /* Here r = 2**24, 0 <= d[i] < r and 1 <= p <= 32. */ - /* p is a global variable. A multi-precision number is */ - /* always normalized. Namely, d[1] > 0. An exception is */ - /* a zero which is characterized by d[0] = 0. The terms */ - /* d[p+1], d[p+2], ... of a none zero number have no */ - /* significance and so are the terms e, d[1],d[2],... */ - /* of a zero. */ + /* x = d[1]*r**(e-1) + d[2]*r**(e-2) + ... + d[p]*r**(e-p). */ + /* Here r = 2**24, 0 <= d[i] < r and 1 <= p <= 32. */ + /* p is a global variable. A multi-precision number is */ + /* always normalized. Namely, d[1] > 0. An exception is */ + /* a zero which is characterized by d[0] = 0. The terms */ + /* d[p+1], d[p+2], ... of a none zero number have no */ + /* significance and so are the terms e, d[1],d[2],... */ + /* of a zero. */ typedef union { int i[2]; double d; } number; @@ -66,15 +65,15 @@ typedef union { int i[2]; double d; } number; #define ABS(x) ((x) < 0 ? -(x) : (x)) int __acr(const mp_no *, const mp_no *, int); -int __cr(const mp_no *, const mp_no *, int); +// int __cr(const mp_no *, const mp_no *, int); void __cpy(const mp_no *, mp_no *, int); -void __cpymn(const mp_no *, int, mp_no *, int); +// void __cpymn(const mp_no *, int, mp_no *, int); void __mp_dbl(const mp_no *, double *, int); void __dbl_mp(double, mp_no *, int); void __add(const mp_no *, const mp_no *, mp_no *, int); void __sub(const mp_no *, const mp_no *, mp_no *, int); void __mul(const mp_no *, const mp_no *, mp_no *, int); -void __inv(const mp_no *, mp_no *, int); +// void __inv(const mp_no *, mp_no *, int); void __dvd(const mp_no *, const mp_no *, mp_no *, int); extern void __mpatan (mp_no *, mp_no *, int); diff --git a/libc/sysdeps/ieee754/dbl-64/mpatan.c b/libc/sysdeps/ieee754/dbl-64/mpatan.c index ee21c2513..f40873ea5 100644 --- a/libc/sysdeps/ieee754/dbl-64/mpatan.c +++ b/libc/sysdeps/ieee754/dbl-64/mpatan.c @@ -1,8 +1,7 @@ - /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -34,11 +33,19 @@ #include "endian.h" #include "mpa.h" -void __mpsqrt(mp_no *, mp_no *, int); -void __mpatan(mp_no *x, mp_no *y, int p) { +#ifndef SECTION +# define SECTION +#endif + #include "mpatan.h" +void __mpsqrt(mp_no *, mp_no *, int); + +void +SECTION +__mpatan(mp_no *x, mp_no *y, int p) { + int i,m,n; double dx; mp_no @@ -54,19 +61,19 @@ void __mpatan(mp_no *x, mp_no *y, int p) { mp_no mps,mpsm,mpt,mpt1,mpt2,mpt3; - /* Choose m and initiate mpone, mptwo & mptwoim1 */ + /* Choose m and initiate mpone, mptwo & mptwoim1 */ if (EX>0) m=7; else if (EX<0) m=0; else { __mp_dbl(x,&dx,p); dx=ABS(dx); for (m=6; m>0; m--) - {if (dx>xm[m].d) break;} + {if (dx>__atan_xm[m].d) break;} } mpone.e = mptwo.e = mptwoim1.e = 1; mpone.d[0] = mpone.d[1] = mptwo.d[0] = mptwoim1.d[0] = ONE; mptwo.d[1] = TWO; - /* Reduce x m times */ + /* Reduce x m times */ __mul(x,x,&mpsm,p); if (m==0) __cpy(x,&mps,p); else { @@ -82,8 +89,8 @@ void __mpatan(mp_no *x, mp_no *y, int p) { __mpsqrt(&mpsm,&mps,p); mps.d[0] = X[0]; } - /* Evaluate a truncated power series for Atan(s) */ - n=np[p]; mptwoim1.d[1] = twonm1[p].d; + /* Evaluate a truncated power series for Atan(s) */ + n=__atan_np[p]; mptwoim1.d[1] = __atan_twonm1[p].d; __dvd(&mpsm,&mptwoim1,&mpt,p); for (i=n-1; i>1; i--) { mptwoim1.d[1] -= TWO; @@ -94,8 +101,8 @@ void __mpatan(mp_no *x, mp_no *y, int p) { __mul(&mps,&mpt,&mpt1,p); __sub(&mps,&mpt1,&mpt,p); - /* Compute Atan(x) */ - mptwoim1.d[1] = twom[m].d; + /* Compute Atan(x) */ + mptwoim1.d[1] = __atan_twom[m].d; __mul(&mptwoim1,&mpt,y,p); return; diff --git a/libc/sysdeps/ieee754/dbl-64/mpatan.h b/libc/sysdeps/ieee754/dbl-64/mpatan.h index d420ff340..003b06c69 100644 --- a/libc/sysdeps/ieee754/dbl-64/mpatan.h +++ b/libc/sysdeps/ieee754/dbl-64/mpatan.h @@ -1,8 +1,7 @@ - /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation, Inc. + * Copyright (C) 2001, 2011 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -29,9 +28,18 @@ #ifndef MPATAN_H #define MPATAN_H +extern const number __atan_xm[8] attribute_hidden; +extern const number __atan_twonm1[33] attribute_hidden; +extern const number __atan_twom[8] attribute_hidden; +extern const number __atan_one attribute_hidden; +extern const number __atan_two attribute_hidden; +extern const int __atan_np[33] attribute_hidden; + + +#ifndef AVOID_MPATAN_H #ifdef BIG_ENDI - static const number - xm[8] = { /* x[m] */ + const number + __atan_xm[8] = { /* x[m] */ /**/ {{0x00000000, 0x00000000} }, /* 0.0 */ /**/ {{0x3f8930be, 0x00000000} }, /* 0.0123 */ /**/ {{0x3f991687, 0x00000000} }, /* 0.0245 */ @@ -40,9 +48,9 @@ /**/ {{0x3fc95810, 0x00000000} }, /* 0.198 */ /**/ {{0x3fda7ef9, 0x00000000} }, /* 0.414 */ /**/ {{0x3ff00000, 0x00000000} }, /* 1.0 */ - }; - static const number - twonm1[33] = { /* 2n-1 */ + }; + const number + __atan_twonm1[33] = { /* 2n-1 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ @@ -76,10 +84,10 @@ /**/ {{0x405b4000, 0x00000000} }, /* 109 */ /**/ {{0x405c4000, 0x00000000} }, /* 113 */ /**/ {{0x405d4000, 0x00000000} }, /* 117 */ - }; + }; - static const number - twom[8] = { /* 2**m */ + const number + __atan_twom[8] = { /* 2**m */ /**/ {{0x3ff00000, 0x00000000} }, /* 1.0 */ /**/ {{0x40000000, 0x00000000} }, /* 2.0 */ /**/ {{0x40100000, 0x00000000} }, /* 4.0 */ @@ -88,17 +96,17 @@ /**/ {{0x40400000, 0x00000000} }, /* 32.0 */ /**/ {{0x40500000, 0x00000000} }, /* 64.0 */ /**/ {{0x40600000, 0x00000000} }, /* 128.0 */ - }; + }; - static const number -/**/ one = {{0x3ff00000, 0x00000000} }, /* 1 */ -/**/ two = {{0x40000000, 0x00000000} }; /* 2 */ + const number +/**/ __atan_one = {{0x3ff00000, 0x00000000} }, /* 1 */ +/**/ __atan_two = {{0x40000000, 0x00000000} }; /* 2 */ #else #ifdef LITTLE_ENDI - static const number - xm[8] = { /* x[m] */ + const number + __atan_xm[8] = { /* x[m] */ /**/ {{0x00000000, 0x00000000} }, /* 0.0 */ /**/ {{0x00000000, 0x3f8930be} }, /* 0.0123 */ /**/ {{0x00000000, 0x3f991687} }, /* 0.0245 */ @@ -107,9 +115,9 @@ /**/ {{0x00000000, 0x3fc95810} }, /* 0.198 */ /**/ {{0x00000000, 0x3fda7ef9} }, /* 0.414 */ /**/ {{0x00000000, 0x3ff00000} }, /* 1.0 */ - }; - static const number - twonm1[33] = { /* 2n-1 */ + }; + const number +__atan_twonm1[33] = { /* 2n-1 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ @@ -143,10 +151,10 @@ /**/ {{0x00000000, 0x405b4000} }, /* 109 */ /**/ {{0x00000000, 0x405c4000} }, /* 113 */ /**/ {{0x00000000, 0x405d4000} }, /* 117 */ - }; + }; - static const number - twom[8] = { /* 2**m */ + const number + __atan_twom[8] = { /* 2**m */ /**/ {{0x00000000, 0x3ff00000} }, /* 1.0 */ /**/ {{0x00000000, 0x40000000} }, /* 2.0 */ /**/ {{0x00000000, 0x40100000} }, /* 4.0 */ @@ -155,20 +163,21 @@ /**/ {{0x00000000, 0x40400000} }, /* 32.0 */ /**/ {{0x00000000, 0x40500000} }, /* 64.0 */ /**/ {{0x00000000, 0x40600000} }, /* 128.0 */ - }; + }; - static const number -/**/ one = {{0x00000000, 0x3ff00000} }, /* 1 */ -/**/ two = {{0x00000000, 0x40000000} }; /* 2 */ + const number +/**/ __atan_one = {{0x00000000, 0x3ff00000} }, /* 1 */ +/**/ __atan_two = {{0x00000000, 0x40000000} }; /* 2 */ #endif #endif -#define ONE one.d -#define TWO two.d - - static const int - np[33] = { 0, 0, 0, 0, 6, 8,10,11,13,15,17,19,21,23,25,27,28, - 30,32,34,36,38,40,42,43,45,47,49,51,53,55,57,59}; + const int + __atan_np[33] = { 0, 0, 0, 0, 6, 8,10,11,13,15,17,19,21,23,25,27,28, + 30,32,34,36,38,40,42,43,45,47,49,51,53,55,57,59}; #endif +#endif + +#define ONE __atan_one.d +#define TWO __atan_two.d diff --git a/libc/sysdeps/ieee754/dbl-64/mpatan2.c b/libc/sysdeps/ieee754/dbl-64/mpatan2.c index 8977ec904..1deb05641 100644 --- a/libc/sysdeps/ieee754/dbl-64/mpatan2.c +++ b/libc/sysdeps/ieee754/dbl-64/mpatan2.c @@ -1,8 +1,7 @@ - /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -38,18 +37,24 @@ #include "mpa.h" +#ifndef SECTION +# define SECTION +#endif + void __mpsqrt(mp_no *, mp_no *, int); void __mpatan(mp_no *, mp_no *, int); /* Multi-Precision Atan2(y,x) function subroutine, for p >= 4. */ /* y=0 is not permitted if x<=0. No error messages are given. */ -void __mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) { +void +SECTION +__mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) { static const double ZERO = 0.0, ONE = 1.0; mp_no mpone = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; mp_no mpt1,mpt2,mpt3; diff --git a/libc/sysdeps/ieee754/dbl-64/mpexp.c b/libc/sysdeps/ieee754/dbl-64/mpexp.c index e2ab71b2c..b0cffe2fe 100644 --- a/libc/sysdeps/ieee754/dbl-64/mpexp.c +++ b/libc/sysdeps/ieee754/dbl-64/mpexp.c @@ -1,8 +1,7 @@ - /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -34,34 +33,40 @@ #include "mpa.h" #include "mpexp.h" +#ifndef SECTION +# define SECTION +#endif + /* Multi-Precision exponential function subroutine (for p >= 4, */ /* 2**(-55) <= abs(x) <= 1024). */ -void __mpexp(mp_no *x, mp_no *y, int p) { +void +SECTION +__mpexp(mp_no *x, mp_no *y, int p) { int i,j,k,m,m1,m2,n; double a,b; static const int np[33] = {0,0,0,0,3,3,4,4,5,4,4,5,5,5,6,6,6,6,6,6, - 6,6,6,6,7,7,7,7,8,8,8,8,8}; + 6,6,6,6,7,7,7,7,8,8,8,8,8}; static const int m1p[33]= {0,0,0,0,17,23,23,28,27,38,42,39,43,47,43,47,50,54, - 57,60,64,67,71,74,68,71,74,77,70,73,76,78,81}; + 57,60,64,67,71,74,68,71,74,77,70,73,76,78,81}; static const int m1np[7][18] = { - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, - { 0, 0, 0, 0,36,48,60,72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, - { 0, 0, 0, 0,24,32,40,48,56,64,72, 0, 0, 0, 0, 0, 0, 0}, - { 0, 0, 0, 0,17,23,29,35,41,47,53,59,65, 0, 0, 0, 0, 0}, - { 0, 0, 0, 0, 0, 0,23,28,33,38,42,47,52,57,62,66, 0, 0}, - { 0, 0, 0, 0, 0, 0, 0, 0,27, 0, 0,39,43,47,51,55,59,63}, - { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,43,47,50,54}}; + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0,36,48,60,72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0,24,32,40,48,56,64,72, 0, 0, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0,17,23,29,35,41,47,53,59,65, 0, 0, 0, 0, 0}, + { 0, 0, 0, 0, 0, 0,23,28,33,38,42,47,52,57,62,66, 0, 0}, + { 0, 0, 0, 0, 0, 0, 0, 0,27, 0, 0,39,43,47,51,55,59,63}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,43,47,50,54}}; mp_no mpone = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; mp_no mpk = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; mp_no mps,mpak,mpt1,mpt2; /* Choose m,n and compute a=2**(-m) */ - n = np[p]; m1 = m1p[p]; a = twomm1[p].d; + n = np[p]; m1 = m1p[p]; a = __mpexp_twomm1[p].d; for (i=0; i<EX; i++) a *= RADIXI; for ( ; i>EX; i--) a *= RADIX; b = X[1]*RADIXI; m2 = 24*EX; @@ -81,12 +86,12 @@ void __mpexp(mp_no *x, mp_no *y, int p) { /* Evaluate the polynomial. Put result in mpt2 */ mpone.e=1; mpone.d[0]=ONE; mpone.d[1]=ONE; - mpk.e = 1; mpk.d[0] = ONE; mpk.d[1]=nn[n].d; + mpk.e = 1; mpk.d[0] = ONE; mpk.d[1]=__mpexp_nn[n].d; __dvd(&mps,&mpk,&mpt1,p); __add(&mpone,&mpt1,&mpak,p); for (k=n-1; k>1; k--) { __mul(&mps,&mpak,&mpt1,p); - mpk.d[1]=nn[k].d; + mpk.d[1]=__mpexp_nn[k].d; __dvd(&mpt1,&mpk,&mpt2,p); __add(&mpone,&mpt2,&mpak,p); } diff --git a/libc/sysdeps/ieee754/dbl-64/mpexp.h b/libc/sysdeps/ieee754/dbl-64/mpexp.h index a0b08ccb4..7985060a8 100644 --- a/libc/sysdeps/ieee754/dbl-64/mpexp.h +++ b/libc/sysdeps/ieee754/dbl-64/mpexp.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation, Inc. + * Copyright (C) 2001, 2011 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -28,9 +28,20 @@ #ifndef MPEXP_H #define MPEXP_H +extern const number __mpexp_twomm1[33] attribute_hidden; +extern const number __mpexp_nn[9] attribute_hidden; +extern const number __mpexp_radix attribute_hidden; +extern const number __mpexp_radixi attribute_hidden; +extern const number __mpexp_zero attribute_hidden; +extern const number __mpexp_one attribute_hidden; +extern const number __mpexp_two attribute_hidden; +extern const number __mpexp_half attribute_hidden; + + +#ifndef AVOID_MPEXP_H #ifdef BIG_ENDI - static const number - twomm1[33] = { /* 2**-m1 */ + const number + __mpexp_twomm1[33] = { /* 2**-m1 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ @@ -65,8 +76,8 @@ /**/ {{0x3b100000, 0x00000000} }, /* 2**-78 */ /**/ {{0x3ae00000, 0x00000000} }, /* 2**-81 */ }; - static const number - nn[9]={ /* n */ + const number + __mpexp_nn[9]={ /* n */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x3ff00000, 0x00000000} }, /* 1 */ /**/ {{0x40000000, 0x00000000} }, /* 2 */ @@ -78,18 +89,18 @@ /**/ {{0x40200000, 0x00000000} }, /* 8 */ }; - static const number -/**/ radix = {{0x41700000, 0x00000000} }, /* 2**24 */ -/**/ radixi = {{0x3e700000, 0x00000000} }, /* 2**-24 */ -/**/ zero = {{0x00000000, 0x00000000} }, /* 0 */ -/**/ one = {{0x3ff00000, 0x00000000} }, /* 1 */ -/**/ two = {{0x40000000, 0x00000000} }, /* 2 */ -/**/ half = {{0x3fe00000, 0x00000000} }; /* 1/2 */ + const number +/**/ __mpexp_radix = {{0x41700000, 0x00000000} }, /* 2**24 */ +/**/ __mpexp_radixi = {{0x3e700000, 0x00000000} }, /* 2**-24 */ +/**/ __mpexp_zero = {{0x00000000, 0x00000000} }, /* 0 */ +/**/ __mpexp_one = {{0x3ff00000, 0x00000000} }, /* 1 */ +/**/ __mpexp_two = {{0x40000000, 0x00000000} }, /* 2 */ +/**/ __mpexp_half = {{0x3fe00000, 0x00000000} }; /* 1/2 */ #else #ifdef LITTLE_ENDI - static const number - twomm1[33] = { /* 2**-m1 */ + const number + __mpexp_twomm1[33] = { /* 2**-m1 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ @@ -124,8 +135,8 @@ /**/ {{0x00000000, 0x3b100000} }, /* 2**-78 */ /**/ {{0x00000000, 0x3ae00000} }, /* 2**-81 */ }; - static const number - nn[9]={ /* n */ + const number + __mpexp_nn[9]={ /* n */ /**/ {{0x00000000, 0x00000000} }, /* 0 */ /**/ {{0x00000000, 0x3ff00000} }, /* 1 */ /**/ {{0x00000000, 0x40000000} }, /* 2 */ @@ -137,22 +148,23 @@ /**/ {{0x00000000, 0x40200000} }, /* 8 */ }; - static const number -/**/ radix = {{0x00000000, 0x41700000} }, /* 2**24 */ -/**/ radixi = {{0x00000000, 0x3e700000} }, /* 2**-24 */ -/**/ zero = {{0x00000000, 0x00000000} }, /* 0 */ -/**/ one = {{0x00000000, 0x3ff00000} }, /* 1 */ -/**/ two = {{0x00000000, 0x40000000} }, /* 2 */ -/**/ half = {{0x00000000, 0x3fe00000} }; /* 1/2 */ + const number +/**/ __mpexp_radix = {{0x00000000, 0x41700000} }, /* 2**24 */ +/**/ __mpexp_radixi = {{0x00000000, 0x3e700000} }, /* 2**-24 */ +/**/ __mpexp_zero = {{0x00000000, 0x00000000} }, /* 0 */ +/**/ __mpexp_one = {{0x00000000, 0x3ff00000} }, /* 1 */ +/**/ __mpexp_two = {{0x00000000, 0x40000000} }, /* 2 */ +/**/ __mpexp_half = {{0x00000000, 0x3fe00000} }; /* 1/2 */ #endif #endif +#endif -#define RADIX radix.d -#define RADIXI radixi.d -#define ZERO zero.d -#define ONE one.d -#define TWO two.d -#define HALF half.d +#define RADIX __mpexp_radix.d +#define RADIXI __mpexp_radixi.d +#define ZERO __mpexp_zero.d +#define ONE __mpexp_one.d +#define TWO __mpexp_two.d +#define HALF __mpexp_half.d #endif diff --git a/libc/sysdeps/ieee754/dbl-64/mpsqrt.c b/libc/sysdeps/ieee754/dbl-64/mpsqrt.c index 9945de306..d1a80f909 100644 --- a/libc/sysdeps/ieee754/dbl-64/mpsqrt.c +++ b/libc/sysdeps/ieee754/dbl-64/mpsqrt.c @@ -1,8 +1,7 @@ - /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -33,6 +32,12 @@ #include "endian.h" #include "mpa.h" +#ifndef SECTION +# define SECTION +#endif + +#include "mpsqrt.h" + /****************************************************************************/ /* Multi-Precision square root function subroutine for precision p >= 4. */ /* The relative error is bounded by 3.501*r**(1-p), where r=2**24. */ @@ -41,20 +46,20 @@ /* p as integer. Routine computes sqrt(*x) and stores result in *y */ /****************************************************************************/ -double fastiroot(double); - -void __mpsqrt(mp_no *x, mp_no *y, int p) { -#include "mpsqrt.h" +static double fastiroot(double); +void +SECTION +__mpsqrt(mp_no *x, mp_no *y, int p) { int i,m,ex,ey; double dx,dy; mp_no mphalf = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}, mp3halfs = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, - 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, + 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; mp_no mpxn,mpz,mpu,mpt1,mpt2; /* Prepare multi-precision 1/2 and 3/2 */ @@ -65,7 +70,7 @@ void __mpsqrt(mp_no *x, mp_no *y, int p) { __mp_dbl(&mpxn,&dx,p); dy=fastiroot(dx); __dbl_mp(dy,&mpu,p); __mul(&mpxn,&mphalf,&mpz,p); - m=mp[p]; + m=__mpsqrt_mp[p]; for (i=0; i<m; i++) { __mul(&mpu,&mpu,&mpt1,p); __mul(&mpt1,&mpz,&mpt2,p); @@ -82,7 +87,9 @@ void __mpsqrt(mp_no *x, mp_no *y, int p) { /* Compute a double precision approximation for 1/sqrt(x) */ /* with the relative error bounded by 2**-51. */ /***********************************************************/ -double fastiroot(double x) { +static double +SECTION +fastiroot(double x) { union {int i[2]; double d;} p,q; double y,z, t; int n; diff --git a/libc/sysdeps/ieee754/dbl-64/mpsqrt.h b/libc/sysdeps/ieee754/dbl-64/mpsqrt.h index 729d57af2..86fa397b2 100644 --- a/libc/sysdeps/ieee754/dbl-64/mpsqrt.h +++ b/libc/sysdeps/ieee754/dbl-64/mpsqrt.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation, Inc. + * Copyright (C) 2001, 2011 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -28,24 +28,31 @@ #ifndef MPSQRT_H #define MPSQRT_H +extern const number __mpsqrt_one attribute_hidden; +extern const number __mpsqrt_halfrad attribute_hidden; +extern const int __mpsqrt_mp[33] attribute_hidden; + + +#ifndef AVOID_MPSQRT_H #ifdef BIG_ENDI - static const number -/**/ one = {{0x3ff00000, 0x00000000} }, /* 1 */ -/**/ halfrad = {{0x41600000, 0x00000000} }; /* 2**23 */ + const number +/**/ __mpsqrt_one = {{0x3ff00000, 0x00000000} }, /* 1 */ +/**/ __mpsqrt_halfrad = {{0x41600000, 0x00000000} }; /* 2**23 */ #else #ifdef LITTLE_ENDI - static const number -/**/ one = {{0x00000000, 0x3ff00000} }, /* 1 */ -/**/ halfrad = {{0x00000000, 0x41600000} }; /* 2**23 */ + const number +/**/ __mpsqrt_one = {{0x00000000, 0x3ff00000} }, /* 1 */ +/**/ __mpsqrt_halfrad = {{0x00000000, 0x41600000} }; /* 2**23 */ #endif #endif -#define ONE one.d -#define HALFRAD halfrad.d + const int __mpsqrt_mp[33] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4, + 4,4,4,4,4,4,4,4,4}; +#endif - static const int mp[33] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4, - 4,4,4,4,4,4,4,4,4}; +#define ONE __mpsqrt_one.d +#define HALFRAD __mpsqrt_halfrad.d #endif diff --git a/libc/sysdeps/ieee754/dbl-64/mptan.c b/libc/sysdeps/ieee754/dbl-64/mptan.c index 267445a19..e1e5d9b92 100644 --- a/libc/sysdeps/ieee754/dbl-64/mptan.c +++ b/libc/sysdeps/ieee754/dbl-64/mptan.c @@ -1,8 +1,7 @@ - /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -38,10 +37,16 @@ #include "endian.h" #include "mpa.h" +#ifndef SECTION +# define SECTION +#endif + int __mpranred(double, mp_no *, int); void __c32(mp_no *, mp_no *, mp_no *, int); -void __mptan(double x, mp_no *mpy, int p) { +void +SECTION +__mptan(double x, mp_no *mpy, int p) { static const double MONE = -1.0; diff --git a/libc/sysdeps/ieee754/dbl-64/s_atan.c b/libc/sysdeps/ieee754/dbl-64/s_atan.c index 5ea83261a..65369ffb2 100644 --- a/libc/sysdeps/ieee754/dbl-64/s_atan.c +++ b/libc/sysdeps/ieee754/dbl-64/s_atan.c @@ -46,7 +46,13 @@ void __mpatan(mp_no *,mp_no *,int); /* see definition in mpatan.c */ static double atanMp(double,const int[]); -double __signArctan(double,double); + + /* Fix the sign of y and return */ +static double __signArctan(double x,double y){ + return __copysign(y, x); +} + + /* An ultimate atan() routine. Given an IEEE double machine number x, */ /* routine computes the correctly rounded (to nearest) value of atan(x). */ double atan(double x) { @@ -54,7 +60,7 @@ double atan(double x) { double cor,s1,ss1,s2,ss2,t1,t2,t3,t7,t8,t9,t10,u,u2,u3, v,vv,w,ww,y,yy,z,zz; -#ifndef DLA_FMA +#ifndef DLA_FMS double t4,t5,t6; #endif #if 0 @@ -203,14 +209,6 @@ double atan(double x) { } - - /* Fix the sign of y and return */ -double __signArctan(double x,double y){ - - if (x<ZERO) return -y; - else return y; -} - /* Final stages. Compute atan(x) by multiple precision arithmetic */ static double atanMp(double x,const int pr[]){ mp_no mpx,mpy,mpy2,mperr,mpt1,mpy1; diff --git a/libc/sysdeps/ieee754/dbl-64/s_ceil.c b/libc/sysdeps/ieee754/dbl-64/s_ceil.c index 695cae5d5..de50e29bf 100644 --- a/libc/sysdeps/ieee754/dbl-64/s_ceil.c +++ b/libc/sysdeps/ieee754/dbl-64/s_ceil.c @@ -32,18 +32,17 @@ __ceil(double x) EXTRACT_WORDS(i0,i1,x); j0 = ((i0>>20)&0x7ff)-0x3ff; if(j0<20) { - if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0<0) {i0=0x80000000;i1=0;} - else if((i0|i1)!=0) { i0=0x3ff00000;i1=0;} - } + if(j0<0) { /* raise inexact if x != 0 */ + math_force_eval(huge+x); + /* return 0*sign(x) if |x|<1 */ + if(i0<0) {i0=0x80000000;i1=0;} + else if((i0|i1)!=0) { i0=0x3ff00000;i1=0;} } else { i = (0x000fffff)>>j0; if(((i0&i)|i1)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0>0) i0 += (0x00100000)>>j0; - i0 &= (~i); i1=0; - } + math_force_eval(huge+x); /* raise inexact flag */ + if(i0>0) i0 += (0x00100000)>>j0; + i0 &= (~i); i1=0; } } else if (j0>51) { if(j0==0x400) return x+x; /* inf or NaN */ @@ -51,17 +50,16 @@ __ceil(double x) } else { i = ((u_int32_t)(0xffffffff))>>(j0-20); if((i1&i)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0>0) { - if(j0==20) i0+=1; - else { - j = i1 + (1<<(52-j0)); - if(j<i1) i0+=1; /* got a carry */ - i1 = j; - } + math_force_eval(huge+x); /* raise inexact flag */ + if(i0>0) { + if(j0==20) i0+=1; + else { + j = i1 + (1<<(52-j0)); + if(j<i1) i0+=1; /* got a carry */ + i1 = j; } - i1 &= (~i); } + i1 &= (~i); } INSERT_WORDS(x,i0,i1); return x; diff --git a/libc/sysdeps/ieee754/dbl-64/s_expm1.c b/libc/sysdeps/ieee754/dbl-64/s_expm1.c index 324354336..589128c08 100644 --- a/libc/sysdeps/ieee754/dbl-64/s_expm1.c +++ b/libc/sysdeps/ieee754/dbl-64/s_expm1.c @@ -13,10 +13,6 @@ for performance improvement on pipelined processors. */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_expm1.c,v 1.8 1995/05/10 20:47:09 jtc Exp $"; -#endif - /* expm1(x) * Returns exp(x)-1, the exponential of x minus 1. * @@ -40,38 +36,38 @@ static char rcsid[] = "$NetBSD: s_expm1.c,v 1.8 1995/05/10 20:47:09 jtc Exp $"; * = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r)) * = 1 - r^2/60 + r^4/2520 - r^6/100800 + ... * We use a special Reme algorithm on [0,0.347] to generate - * a polynomial of degree 5 in r*r to approximate R1. The + * a polynomial of degree 5 in r*r to approximate R1. The * maximum error of this polynomial approximation is bounded * by 2**-61. In other words, * R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5 - * where Q1 = -1.6666666666666567384E-2, - * Q2 = 3.9682539681370365873E-4, - * Q3 = -9.9206344733435987357E-6, - * Q4 = 2.5051361420808517002E-7, - * Q5 = -6.2843505682382617102E-9; - * (where z=r*r, and the values of Q1 to Q5 are listed below) + * where Q1 = -1.6666666666666567384E-2, + * Q2 = 3.9682539681370365873E-4, + * Q3 = -9.9206344733435987357E-6, + * Q4 = 2.5051361420808517002E-7, + * Q5 = -6.2843505682382617102E-9; + * (where z=r*r, and the values of Q1 to Q5 are listed below) * with error bounded by * | 5 | -61 * | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2 * | | * * expm1(r) = exp(r)-1 is then computed by the following - * specific way which minimize the accumulation rounding error: + * specific way which minimize the accumulation rounding error: * 2 3 * r r [ 3 - (R1 + R1*r/2) ] * expm1(r) = r + --- + --- * [--------------------] - * 2 2 [ 6 - r*(3 - R1*r/2) ] + * 2 2 [ 6 - r*(3 - R1*r/2) ] * * To compensate the error in the argument reduction, we use * expm1(r+c) = expm1(r) + c + expm1(r)*c * ~ expm1(r) + c + r*c * Thus c+r*c will be added in as the correction terms for * expm1(r+c). Now rearrange the term to avoid optimization - * screw up: - * ( 2 2 ) - * ({ ( r [ R1 - (3 - R1*r/2) ] ) } r ) + * screw up: + * ( 2 2 ) + * ({ ( r [ R1 - (3 - R1*r/2) ] ) } r ) * expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- ) - * ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 ) + * ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 ) * ( ) * * = r - E @@ -87,7 +83,7 @@ static char rcsid[] = "$NetBSD: s_expm1.c,v 1.8 1995/05/10 20:47:09 jtc Exp $"; * (ii) if k=0, return r-E * (iii) if k=-1, return 0.5*(r-E)-0.5 * (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E) - * else return 1.0+2.0*(r-E); + * else return 1.0+2.0*(r-E); * (v) if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1) * (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else * (vii) return 2^k(1-((E+2^-k)-r)) @@ -116,11 +112,7 @@ static char rcsid[] = "$NetBSD: s_expm1.c,v 1.8 1995/05/10 20:47:09 jtc Exp $"; #include "math.h" #include "math_private.h" #define one Q[0] -#ifdef __STDC__ static const double -#else -static double -#endif huge = 1.0e+300, tiny = 1.0e-300, o_threshold = 7.09782712893383973096e+02,/* 0x40862E42, 0xFEFA39EF */ @@ -134,12 +126,8 @@ Q[] = {1.0, -3.33333333333331316428e-02, /* BFA11111 111110F4 */ 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */ -2.01099218183624371326e-07}; /* BE8AFDB7 6E09C32D */ -#ifdef __STDC__ - double __expm1(double x) -#else - double __expm1(x) - double x; -#endif +double +__expm1(double x) { double y,hi,lo,c,t,e,hxs,hfx,r1,h2,h4,R1,R2,R3; int32_t k,xsb; @@ -153,20 +141,20 @@ Q[] = {1.0, -3.33333333333331316428e-02, /* BFA11111 111110F4 */ /* filter out huge and non-finite argument */ if(hx >= 0x4043687A) { /* if |x|>=56*ln2 */ if(hx >= 0x40862E42) { /* if |x|>=709.78... */ - if(hx>=0x7ff00000) { + if(hx>=0x7ff00000) { u_int32_t low; GET_LOW_WORD(low,x); if(((hx&0xfffff)|low)!=0) - return x+x; /* NaN */ + return x+x; /* NaN */ else return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */ - } - if(x > o_threshold) { + } + if(x > o_threshold) { __set_errno (ERANGE); return huge*huge; /* overflow */ } } if(xsb!=0) { /* x < -56*ln2, return -1.0 with inexact */ - if(x+tiny<0.0) /* raise inexact */ + math_force_eval(x+tiny); /* raise inexact */ return tiny-one; /* return -1 */ } } @@ -187,7 +175,7 @@ Q[] = {1.0, -3.33333333333331316428e-02, /* BFA11111 111110F4 */ x = hi - lo; c = (hi-x)-lo; } - else if(hx < 0x3c900000) { /* when |x|<2**-54, return x */ + else if(hx < 0x3c900000) { /* when |x|<2**-54, return x */ t = huge+x; /* return x with inexact flags when x!=0 */ return x - (t-(huge+x)); } @@ -212,28 +200,28 @@ Q[] = {1.0, -3.33333333333331316428e-02, /* BFA11111 111110F4 */ e -= hxs; if(k== -1) return 0.5*(x-e)-0.5; if(k==1) { - if(x < -0.25) return -2.0*(e-(x+0.5)); - else return one+2.0*(x-e); + if(x < -0.25) return -2.0*(e-(x+0.5)); + else return one+2.0*(x-e); } if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */ - u_int32_t high; - y = one-(e-x); + u_int32_t high; + y = one-(e-x); GET_HIGH_WORD(high,y); SET_HIGH_WORD(y,high+(k<<20)); /* add k to y's exponent */ - return y-one; + return y-one; } t = one; if(k<20) { - u_int32_t high; - SET_HIGH_WORD(t,0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */ - y = t-(e-x); + u_int32_t high; + SET_HIGH_WORD(t,0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */ + y = t-(e-x); GET_HIGH_WORD(high,y); SET_HIGH_WORD(y,high+(k<<20)); /* add k to y's exponent */ } else { - u_int32_t high; + u_int32_t high; SET_HIGH_WORD(t,((0x3ff-k)<<20)); /* 2^-k */ - y = x-(e+t); - y += one; + y = x-(e+t); + y += one; GET_HIGH_WORD(high,y); SET_HIGH_WORD(y,high+(k<<20)); /* add k to y's exponent */ } diff --git a/libc/sysdeps/ieee754/dbl-64/s_floor.c b/libc/sysdeps/ieee754/dbl-64/s_floor.c index 5b593ca31..8b2038e69 100644 --- a/libc/sysdeps/ieee754/dbl-64/s_floor.c +++ b/libc/sysdeps/ieee754/dbl-64/s_floor.c @@ -41,18 +41,16 @@ static double huge = 1.0e300; j0 = ((i0>>20)&0x7ff)-0x3ff; if(j0<20) { if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0>=0) {i0=i1=0;} - else if(((i0&0x7fffffff)|i1)!=0) - { i0=0xbff00000;i1=0;} - } + math_force_eval(huge+x);/* return 0*sign(x) if |x|<1 */ + if(i0>=0) {i0=i1=0;} + else if(((i0&0x7fffffff)|i1)!=0) + { i0=0xbff00000;i1=0;} } else { i = (0x000fffff)>>j0; if(((i0&i)|i1)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0<0) i0 += (0x00100000)>>j0; - i0 &= (~i); i1=0; - } + math_force_eval(huge+x); /* raise inexact flag */ + if(i0<0) i0 += (0x00100000)>>j0; + i0 &= (~i); i1=0; } } else if (j0>51) { if(j0==0x400) return x+x; /* inf or NaN */ @@ -60,17 +58,16 @@ static double huge = 1.0e300; } else { i = ((u_int32_t)(0xffffffff))>>(j0-20); if((i1&i)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0<0) { - if(j0==20) i0+=1; - else { - j = i1+(1<<(52-j0)); - if(j<i1) i0 +=1 ; /* got a carry */ - i1=j; - } + math_force_eval(huge+x); /* raise inexact flag */ + if(i0<0) { + if(j0==20) i0+=1; + else { + j = i1+(1<<(52-j0)); + if(j<i1) i0 +=1 ; /* got a carry */ + i1=j; } - i1 &= (~i); } + i1 &= (~i); } INSERT_WORDS(x,i0,i1); return x; diff --git a/libc/sysdeps/ieee754/dbl-64/s_log1p.c b/libc/sysdeps/ieee754/dbl-64/s_log1p.c index 0a9801a93..dc79a02bb 100644 --- a/libc/sysdeps/ieee754/dbl-64/s_log1p.c +++ b/libc/sysdeps/ieee754/dbl-64/s_log1p.c @@ -13,10 +13,6 @@ for performance improvement on pipelined processors. */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $"; -#endif - /* double log1p(double x) * * Method : @@ -34,14 +30,14 @@ static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $"; * 2. Approximation of log1p(f). * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) * = 2s + 2/3 s**3 + 2/5 s**5 + ....., - * = 2s + s*R + * = 2s + s*R * We use a special Reme algorithm on [0,0.1716] to generate - * a polynomial of degree 14 to approximate R The maximum error + * a polynomial of degree 14 to approximate R The maximum error * of this polynomial approximation is bounded by 2**-58.45. In * other words, - * 2 4 6 8 10 12 14 + * 2 4 6 8 10 12 14 * R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s +Lp6*s +Lp7*s - * (the values of Lp1 to Lp7 are listed in the program) + * (the values of Lp1 to Lp7 are listed in the program) * and * | 2 14 | -58.45 * | Lp1*s +...+Lp7*s - R(z) | <= 2 @@ -52,7 +48,7 @@ static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $"; * log1p(f) = f - (hfsq - s*(hfsq+R)). * * 3. Finally, log1p(x) = k*ln2 + log1p(f). - * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) + * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) * Here ln2 is split into two floating point number: * ln2_hi + ln2_lo, * where n*ln2_hi is always exact for |n| < 2000. @@ -73,7 +69,7 @@ static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $"; * to produce the hexadecimal values shown. * * Note: Assuming log() return accurate answer, the following - * algorithm can be used to compute log1p(x) to within a few ULP: + * algorithm can be used to compute log1p(x) to within a few ULP: * * u = 1+x; * if(u==1.0) return x ; else @@ -85,11 +81,7 @@ static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $"; #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const double -#else -static double -#endif ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ @@ -101,18 +93,10 @@ Lp[] = {0.0, 6.666666666666735130e-01, /* 3FE55555 55555593 */ 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ 1.479819860511658591e-01}; /* 3FC2F112 DF3E5244 */ -#ifdef __STDC__ static const double zero = 0.0; -#else -static double zero = 0.0; -#endif -#ifdef __STDC__ - double __log1p(double x) -#else - double __log1p(x) - double x; -#endif +double +__log1p(double x) { double hfsq,f,c,s,z,R,u,z2,z4,z6,R1,R2,R3,R4; int32_t k,hx,hu,ax; @@ -127,8 +111,8 @@ static double zero = 0.0; else return (x-x)/(x-x); /* log1p(x<-1)=NaN */ } if(ax<0x3e200000) { /* |x| < 2**-29 */ - if(two54+x>zero /* raise inexact */ - &&ax<0x3c900000) /* |x| < 2**-54 */ + math_force_eval(two54+x); /* raise inexact */ + if (ax<0x3c900000) /* |x| < 2**-54 */ return x; else return x - x*x*0.5; @@ -141,22 +125,22 @@ static double zero = 0.0; if(hx<0x43400000) { u = 1.0+x; GET_HIGH_WORD(hu,u); - k = (hu>>20)-1023; - c = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */ + k = (hu>>20)-1023; + c = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */ c /= u; } else { u = x; GET_HIGH_WORD(hu,u); - k = (hu>>20)-1023; + k = (hu>>20)-1023; c = 0; } hu &= 0x000fffff; if(hu<0x6a09e) { - SET_HIGH_WORD(u,hu|0x3ff00000); /* normalize u */ + SET_HIGH_WORD(u,hu|0x3ff00000); /* normalize u */ } else { - k += 1; + k += 1; SET_HIGH_WORD(u,hu|0x3fe00000); /* normalize u/2 */ - hu = (0x00100000-hu)>>2; + hu = (0x00100000-hu)>>2; } f = u-1.0; } @@ -168,9 +152,9 @@ static double zero = 0.0; } R = hfsq*(1.0-0.66666666666666666*f); if(k==0) return f-R; else - return k*ln2_hi-((R-(k*ln2_lo+c))-f); + return k*ln2_hi-((R-(k*ln2_lo+c))-f); } - s = f/(2.0+f); + s = f/(2.0+f); z = s*s; #ifdef DO_NOT_USE_THIS R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7)))))); diff --git a/libc/sysdeps/ieee754/dbl-64/s_round.c b/libc/sysdeps/ieee754/dbl-64/s_round.c index 94fcde0e4..f8a38163f 100644 --- a/libc/sysdeps/ieee754/dbl-64/s_round.c +++ b/libc/sysdeps/ieee754/dbl-64/s_round.c @@ -1,5 +1,5 @@ /* Round double to integer away from zero. - Copyright (C) 1997 Free Software Foundation, Inc. + Copyright (C) 1997, 2011 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -38,13 +38,12 @@ __round (double x) { if (j0 < 0) { - if (huge + x > 0.0) - { - i0 &= 0x80000000; - if (j0 == -1) - i0 |= 0x3ff00000; - i1 = 0; - } + math_force_eval (huge + x > 0.0); + + i0 &= 0x80000000; + if (j0 == -1) + i0 |= 0x3ff00000; + i1 = 0; } else { @@ -52,13 +51,12 @@ __round (double x) if (((i0 & i) | i1) == 0) /* X is integral. */ return x; - if (huge + x > 0.0) - { - /* Raise inexact if x != 0. */ - i0 += 0x00080000 >> j0; - i0 &= ~i; - i1 = 0; - } + math_force_eval (huge + x > 0.0); + + /* Raise inexact if x != 0. */ + i0 += 0x00080000 >> j0; + i0 &= ~i; + i1 = 0; } } else if (j0 > 51) @@ -76,14 +74,13 @@ __round (double x) /* X is integral. */ return x; - if (huge + x > 0.0) - { - /* Raise inexact if x != 0. */ - u_int32_t j = i1 + (1 << (51 - j0)); - if (j < i1) - i0 += 1; - i1 = j; - } + math_force_eval (huge + x > 0.0); + + /* Raise inexact if x != 0. */ + u_int32_t j = i1 + (1 << (51 - j0)); + if (j < i1) + i0 += 1; + i1 = j; i1 &= ~i; } diff --git a/libc/sysdeps/ieee754/dbl-64/s_sin.c b/libc/sysdeps/ieee754/dbl-64/s_sin.c index b40776f5e..6f19f158f 100644 --- a/libc/sysdeps/ieee754/dbl-64/s_sin.c +++ b/libc/sysdeps/ieee754/dbl-64/s_sin.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001, 2009 Free Software Foundation + * Copyright (C) 2001, 2009, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -53,15 +53,24 @@ #include "mydefs.h" #include "usncs.h" #include "MathLib.h" -#include "sincos.tbl" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + +extern const union +{ + int4 i[880]; + double x[440]; +} __sincostab attribute_hidden; + static const double - sn3 = -1.66666666666664880952546298448555E-01, - sn5 = 8.33333214285722277379541354343671E-03, - cs2 = 4.99999999999999999999950396842453E-01, - cs4 = -4.16666666666664434524222570944589E-02, - cs6 = 1.38888874007937613028114285595617E-03; + sn3 = -1.66666666666664880952546298448555E-01, + sn5 = 8.33333214285722277379541354343671E-03, + cs2 = 4.99999999999999999999950396842453E-01, + cs4 = -4.16666666666664434524222570944589E-02, + cs6 = 1.38888874007937613028114285595617E-03; void __dubsin(double x, double dx, double w[]); void __docos(double x, double dx, double w[]); @@ -87,7 +96,9 @@ static double csloww2(double x, double dx, double orig, int n); /* An ultimate sin routine. Given an IEEE double machine number x */ /* it computes the correctly rounded (to nearest) value of sin(x) */ /*******************************************************************/ -double __sin(double x){ +double +SECTION +__sin(double x){ double xx,res,t,cor,y,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2; #if 0 double w[2]; @@ -120,10 +131,10 @@ double __sin(double x){ s = y + y*xx*(sn3 +xx*sn5); c = xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=(m>0)?sincos.x[k]:-sincos.x[k]; - ssn=(m>0)?sincos.x[k+1]:-sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=(m>0)?__sincostab.x[k]:-__sincostab.x[k]; + ssn=(m>0)?__sincostab.x[k+1]:-__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; cor=(ssn+s*ccs-sn*c)+cs*s; res=sn+cor; cor=(sn-res)+cor; @@ -146,10 +157,10 @@ double __sin(double x){ s = y + y*xx*(sn3 +xx*sn5); c = xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; cor=(ccs-s*ssn-cs*c)-sn*s; res=cs+cor; cor=(cs-res)+cor; @@ -174,7 +185,7 @@ double __sin(double x){ xx = a*a; if (n) {a=-a;da=-da;} if (xx < 0.01588) { - /*Taylor series */ + /*Taylor series */ t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da; res = a+t; cor = (a-res)+t; @@ -192,10 +203,10 @@ double __sin(double x){ s = y + (db+y*xx*(sn3 +xx*sn5)); c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; cor=(ssn+s*ccs-sn*c)+cs*s; res=sn+cor; cor=(sn-res)+cor; @@ -212,10 +223,10 @@ double __sin(double x){ y=a-(u.x-big.x)+da; xx=y*y; k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; s = y + y*xx*(sn3 +xx*sn5); c = xx*(cs2 +xx*(cs4 + xx*cs6)); cor=(ccs-s*ssn-cs*c)-sn*s; @@ -253,7 +264,7 @@ double __sin(double x){ xx = a*a; if (n) {a=-a;da=-da;} if (xx < 0.01588) { - /* Taylor series */ + /* Taylor series */ t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da; res = a+t; cor = (a-res)+t; @@ -269,10 +280,10 @@ double __sin(double x){ s = y + (db+y*xx*(sn3 +xx*sn5)); c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; cor=(ssn+s*ccs-sn*c)+cs*s; res=sn+cor; cor=(sn-res)+cor; @@ -289,10 +300,10 @@ double __sin(double x){ y=a-(u.x-big.x)+da; xx=y*y; k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; s = y + y*xx*(sn3 +xx*sn5); c = xx*(cs2 +xx*(cs4 + xx*cs6)); cor=(ccs-s*ssn-cs*c)-sn*s; @@ -344,7 +355,9 @@ double __sin(double x){ /* it computes the correctly rounded (to nearest) value of cos(x) */ /*******************************************************************/ -double __cos(double x) +double +SECTION +__cos(double x) { double y,xx,res,t,cor,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2; mynumber u,v; @@ -364,10 +377,10 @@ double __cos(double x) s = y + y*xx*(sn3 +xx*sn5); c = xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; cor=(ccs-s*ssn-cs*c)-sn*s; res=cs+cor; cor=(cs-res)+cor; @@ -396,10 +409,10 @@ double __cos(double x) s = y + (db+y*xx*(sn3 +xx*sn5)); c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; cor=(ssn+s*ccs-sn*c)+cs*s; res=sn+cor; cor=(sn-res)+cor; @@ -442,10 +455,10 @@ double __cos(double x) s = y + (db+y*xx*(sn3 +xx*sn5)); c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; cor=(ssn+s*ccs-sn*c)+cs*s; res=sn+cor; cor=(sn-res)+cor; @@ -461,10 +474,10 @@ double __cos(double x) y=a-(u.x-big.x)+da; xx=y*y; k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; s = y + y*xx*(sn3 +xx*sn5); c = xx*(cs2 +xx*(cs4 + xx*cs6)); cor=(ccs-s*ssn-cs*c)-sn*s; @@ -473,7 +486,7 @@ double __cos(double x) cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps; return (res==res+cor)? ((n)?-res:res) : csloww2(a,da,x,n); - break; + break; } @@ -517,10 +530,10 @@ double __cos(double x) s = y + (db+y*xx*(sn3 +xx*sn5)); c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; cor=(ssn+s*ccs-sn*c)+cs*s; res=sn+cor; cor=(sn-res)+cor; @@ -536,10 +549,10 @@ double __cos(double x) y=a-(u.x-big.x)+da; xx=y*y; k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; s = y + y*xx*(sn3 +xx*sn5); c = xx*(cs2 +xx*(cs4 + xx*cs6)); cor=(ccs-s*ssn-cs*c)-sn*s; @@ -591,7 +604,9 @@ double __cos(double x) /* precision and if still doesn't accurate enough by mpsin or dubsin */ /************************************************************************/ -static double slow(double x) { +static double +SECTION +slow(double x) { static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ double y,x1,x2,xx,r,t,res,cor,w[2]; x1=(x+th2_36)-th2_36; @@ -611,11 +626,13 @@ static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ } } /*******************************************************************************/ -/* Routine compute sin(x) for 0.25<|x|< 0.855469 by sincos.tbl and Taylor */ +/* Routine compute sin(x) for 0.25<|x|< 0.855469 by __sincostab.tbl and Taylor */ /* and if result still doesn't accurate enough by mpsin or dubsin */ /*******************************************************************************/ -static double slow1(double x) { +static double +SECTION +slow1(double x) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res; static const double t22 = 6291456.0; @@ -627,10 +644,10 @@ static double slow1(double x) { s = y*xx*(sn3 +xx*sn5); c = xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; /* Data */ - ssn=sincos.x[k+1]; /* from */ - cs=sincos.x[k+2]; /* tables */ - ccs=sincos.x[k+3]; /* sincos.tbl */ + sn=__sincostab.x[k]; /* Data */ + ssn=__sincostab.x[k+1]; /* from */ + cs=__sincostab.x[k+2]; /* tables */ + ccs=__sincostab.x[k+3]; /* __sincostab.tbl */ y1 = (y+t22)-t22; y2 = y - y1; c1 = (cs+t22)-t22; @@ -648,10 +665,12 @@ static double slow1(double x) { } } /**************************************************************************/ -/* Routine compute sin(x) for 0.855469 <|x|<2.426265 by sincos.tbl */ +/* Routine compute sin(x) for 0.855469 <|x|<2.426265 by __sincostab.tbl */ /* and if result still doesn't accurate enough by mpsin or dubsin */ /**************************************************************************/ -static double slow2(double x) { +static double +SECTION +slow2(double x) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res,del; static const double t22 = 6291456.0; @@ -672,10 +691,10 @@ static double slow2(double x) { s = y*xx*(sn3 +xx*sn5); c = y*del+xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; y1 = (y+t22)-t22; y2 = (y - y1)+del; e1 = (sn+t22)-t22; @@ -703,7 +722,9 @@ static double slow2(double x) { /* result.And if result not accurate enough routine calls mpsin1 or dubsin */ /***************************************************************************/ -static double sloww(double x,double dx, double orig) { +static double +SECTION +sloww(double x,double dx, double orig) { static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ double y,x1,x2,xx,r,t,res,cor,w[2],a,da,xn; union {int4 i[2]; double x;} v; @@ -750,7 +771,9 @@ static double sloww(double x,double dx, double orig) { /* accurate enough routine calls mpsin1 or dubsin */ /***************************************************************************/ -static double sloww1(double x, double dx, double orig) { +static double +SECTION +sloww1(double x, double dx, double orig) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res; static const double t22 = 6291456.0; @@ -763,10 +786,10 @@ static double sloww1(double x, double dx, double orig) { s = y*xx*(sn3 +xx*sn5); c = xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; y1 = (y+t22)-t22; y2 = (y - y1)+dx; c1 = (cs+t22)-t22; @@ -792,7 +815,9 @@ static double sloww1(double x, double dx, double orig) { /* accurate enough routine calls mpsin1 or dubsin */ /***************************************************************************/ -static double sloww2(double x, double dx, double orig, int n) { +static double +SECTION +sloww2(double x, double dx, double orig, int n) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res; static const double t22 = 6291456.0; @@ -805,10 +830,10 @@ static double sloww2(double x, double dx, double orig, int n) { s = y*xx*(sn3 +xx*sn5); c = y*dx+xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; y1 = (y+t22)-t22; y2 = (y - y1)+dx; @@ -836,7 +861,9 @@ static double sloww2(double x, double dx, double orig, int n) { /* result.And if result not accurate enough routine calls other routines */ /***************************************************************************/ -static double bsloww(double x,double dx, double orig,int n) { +static double +SECTION +bsloww(double x,double dx, double orig,int n) { static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ double y,x1,x2,xx,r,t,res,cor,w[2]; #if 0 @@ -869,7 +896,9 @@ static double bsloww(double x,double dx, double orig,int n) { /* And if result not accurate enough routine calls other routines */ /***************************************************************************/ -static double bsloww1(double x, double dx, double orig,int n) { +static double +SECTION +bsloww1(double x, double dx, double orig,int n) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res; static const double t22 = 6291456.0; @@ -882,10 +911,10 @@ mynumber u; s = y*xx*(sn3 +xx*sn5); c = xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; y1 = (y+t22)-t22; y2 = (y - y1)+dx; c1 = (cs+t22)-t22; @@ -912,7 +941,9 @@ mynumber u; /* And if result not accurate enough routine calls other routines */ /***************************************************************************/ -static double bsloww2(double x, double dx, double orig, int n) { +static double +SECTION +bsloww2(double x, double dx, double orig, int n) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res; static const double t22 = 6291456.0; @@ -925,10 +956,10 @@ mynumber u; s = y*xx*(sn3 +xx*sn5); c = y*dx+xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; y1 = (y+t22)-t22; y2 = (y - y1)+dx; @@ -954,7 +985,9 @@ mynumber u; /* precision and if still doesn't accurate enough by mpcos or docos */ /************************************************************************/ -static double cslow2(double x) { +static double +SECTION +cslow2(double x) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res; static const double t22 = 6291456.0; @@ -966,10 +999,10 @@ static double cslow2(double x) { s = y*xx*(sn3 +xx*sn5); c = xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; y1 = (y+t22)-t22; y2 = y - y1; e1 = (sn+t22)-t22; @@ -997,7 +1030,9 @@ static double cslow2(double x) { /***************************************************************************/ -static double csloww(double x,double dx, double orig) { +static double +SECTION +csloww(double x,double dx, double orig) { static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ double y,x1,x2,xx,r,t,res,cor,w[2],a,da,xn; union {int4 i[2]; double x;} v; @@ -1046,7 +1081,9 @@ static double csloww(double x,double dx, double orig) { /* accurate enough routine calls other routines */ /***************************************************************************/ -static double csloww1(double x, double dx, double orig) { +static double +SECTION +csloww1(double x, double dx, double orig) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res; static const double t22 = 6291456.0; @@ -1059,10 +1096,10 @@ static double csloww1(double x, double dx, double orig) { s = y*xx*(sn3 +xx*sn5); c = xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; y1 = (y+t22)-t22; y2 = (y - y1)+dx; c1 = (cs+t22)-t22; @@ -1090,7 +1127,9 @@ static double csloww1(double x, double dx, double orig) { /* accurate enough routine calls other routines */ /***************************************************************************/ -static double csloww2(double x, double dx, double orig, int n) { +static double +SECTION +csloww2(double x, double dx, double orig, int n) { mynumber u; double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res; static const double t22 = 6291456.0; @@ -1103,10 +1142,10 @@ static double csloww2(double x, double dx, double orig, int n) { s = y*xx*(sn3 +xx*sn5); c = y*dx+xx*(cs2 +xx*(cs4 + xx*cs6)); k=u.i[LOW_HALF]<<2; - sn=sincos.x[k]; - ssn=sincos.x[k+1]; - cs=sincos.x[k+2]; - ccs=sincos.x[k+3]; + sn=__sincostab.x[k]; + ssn=__sincostab.x[k+1]; + cs=__sincostab.x[k+2]; + ccs=__sincostab.x[k+3]; y1 = (y+t22)-t22; y2 = (y - y1)+dx; @@ -1127,12 +1166,17 @@ static double csloww2(double x, double dx, double orig, int n) { } } +#ifndef __cos weak_alias (__cos, cos) +# ifdef NO_LONG_DOUBLE +strong_alias (__cos, __cosl) +weak_alias (__cos, cosl) +# endif +#endif +#ifndef __sin weak_alias (__sin, sin) - -#ifdef NO_LONG_DOUBLE +# ifdef NO_LONG_DOUBLE strong_alias (__sin, __sinl) weak_alias (__sin, sinl) -strong_alias (__cos, __cosl) -weak_alias (__cos, cosl) +# endif #endif diff --git a/libc/sysdeps/ieee754/dbl-64/s_tan.c b/libc/sysdeps/ieee754/dbl-64/s_tan.c index df8eedd92..f0fcd677a 100644 --- a/libc/sysdeps/ieee754/dbl-64/s_tan.c +++ b/libc/sysdeps/ieee754/dbl-64/s_tan.c @@ -41,17 +41,23 @@ #include "MathLib.h" #include "math.h" +#ifndef SECTION +# define SECTION +#endif + static double tanMp(double); void __mptan(double, mp_no *, int); -double tan(double x) { +double +SECTION +tan(double x) { #include "utan.h" #include "utan.tbl" int ux,i,n; double a,da,a2,b,db,c,dc,c1,cc1,c2,cc2,c3,cc3,fi,ffi,gi,pz,s,sy, t,t1,t2,t3,t4,t7,t8,t9,t10,w,x2,xn,xx2,y,ya,yya,z0,z,zz,z2,zz2; -#ifndef DLA_FMA +#ifndef DLA_FMS double t5,t6; #endif int p; @@ -486,7 +492,9 @@ double tan(double x) { /* multiple precision stage */ /* Convert x to multi precision number,compute tan(x) by mptan() routine */ /* and converts result back to double */ -static double tanMp(double x) +static double +SECTION +tanMp(double x) { int p; double y; diff --git a/libc/sysdeps/ieee754/dbl-64/sincos32.c b/libc/sysdeps/ieee754/dbl-64/sincos32.c index a4f896a46..e39aaeea0 100644 --- a/libc/sysdeps/ieee754/dbl-64/sincos32.c +++ b/libc/sysdeps/ieee754/dbl-64/sincos32.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -45,11 +45,17 @@ #include "sincos32.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + /****************************************************************/ /* Compute Multi-Precision sin() function for given p. Receive */ /* Multi Precision number x and result stored at y */ /****************************************************************/ -static void ss32(mp_no *x, mp_no *y, int p) { +static void +SECTION +ss32(mp_no *x, mp_no *y, int p) { int i; double a; #if 0 @@ -79,7 +85,9 @@ static void ss32(mp_no *x, mp_no *y, int p) { /* Compute Multi-Precision cos() function for given p. Receive Multi */ /* Precision number x and result stored at y */ /**********************************************************************/ -static void cc32(mp_no *x, mp_no *y, int p) { +static void +SECTION +cc32(mp_no *x, mp_no *y, int p) { int i; double a; #if 0 @@ -109,7 +117,9 @@ static void cc32(mp_no *x, mp_no *y, int p) { /***************************************************************************/ /* c32() computes both sin(x), cos(x) as Multi precision numbers */ /***************************************************************************/ -void __c32(mp_no *x, mp_no *y, mp_no *z, int p) { +void +SECTION +__c32(mp_no *x, mp_no *y, mp_no *z, int p) { static const mp_no mpt={1,{1.0,2.0}}, one={1,{1.0,1.0}}; mp_no u,t,t1,t2,c,s; int i; @@ -134,7 +144,9 @@ void __c32(mp_no *x, mp_no *y, mp_no *z, int p) { /*result which is more accurate */ /*Computing sin(x) with multi precision routine c32 */ /************************************************************************/ -double __sin32(double x, double res, double res1) { +double +SECTION +__sin32(double x, double res, double res1) { int p; mp_no a,b,c; p=32; @@ -158,7 +170,9 @@ double __sin32(double x, double res, double res1) { /*result which is more accurate */ /*Computing cos(x) with multi precision routine c32 */ /************************************************************************/ -double __cos32(double x, double res, double res1) { +double +SECTION +__cos32(double x, double res, double res1) { int p; mp_no a,b,c; p=32; @@ -172,12 +186,12 @@ double __cos32(double x, double res, double res1) { } else if (x>0.8) { __sub(&hp,&c,&a,p); - __c32(&a,&c,&b,p); + __c32(&a,&c,&b,p); } else __c32(&c,&b,&a,p); /* b=cos(0.5*(res+res1)) */ __dbl_mp(x,&c,p); /* c = x */ __sub(&b,&c,&a,p); - /* if a>0 return max(res,res1), otherwise return min(res,res1) */ + /* if a>0 return max(res,res1), otherwise return min(res,res1) */ if (a.d[0]>0) return (res>res1)?res:res1; else return (res<res1)?res:res1; } @@ -186,7 +200,9 @@ double __cos32(double x, double res, double res1) { /*Compute sin(x+dx) as Multi Precision number and return result as */ /* double */ /*******************************************************************/ -double __mpsin(double x, double dx) { +double +SECTION +__mpsin(double x, double dx) { int p; double y; mp_no a,b,c; @@ -204,7 +220,9 @@ double __mpsin(double x, double dx) { /* Compute cos()of double-length number (x+dx) as Multi Precision */ /* number and return result as double */ /*******************************************************************/ -double __mpcos(double x, double dx) { +double +SECTION +__mpcos(double x, double dx) { int p; double y; mp_no a,b,c; @@ -227,7 +245,9 @@ double __mpcos(double x, double dx) { /* n=0,+-1,+-2,.... */ /* Return int which indicates in which quarter of circle x is */ /******************************************************************/ -int __mpranred(double x, mp_no *y, int p) +int +SECTION +__mpranred(double x, mp_no *y, int p) { number v; double t,xn; @@ -275,7 +295,9 @@ int __mpranred(double x, mp_no *y, int p) /* Multi-Precision sin() function subroutine, for p=32. It is */ /* based on the routines mpranred() and c32(). */ /*******************************************************************/ -double __mpsin1(double x) +double +SECTION +__mpsin1(double x) { int p; int n; @@ -314,7 +336,9 @@ double __mpsin1(double x) /* based on the routines mpranred() and c32(). */ /*****************************************************************/ -double __mpcos1(double x) +double +SECTION +__mpcos1(double x) { int p; int n; diff --git a/libc/sysdeps/ieee754/dbl-64/sincos.tbl b/libc/sysdeps/ieee754/dbl-64/sincostab.c index 9343f2416..49fccac94 100644 --- a/libc/sysdeps/ieee754/dbl-64/sincos.tbl +++ b/libc/sysdeps/ieee754/dbl-64/sincostab.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001, 2007 Free Software Foundation, Inc. + * Copyright (C) 2001, 2007, 2011 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -18,13 +18,16 @@ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ +#include <mydefs.h> +#include <endian.h> + /****************************************************************/ /* TABLES FOR THE usin() and ucos() FUNCTION */ /****************************************************************/ #ifdef BIG_ENDI -static const union {int4 i[880]; double x[440];}sincos = { .i = { +const union {int4 i[880]; double x[440];}__sincostab = { .i = { /**/ 0x00000000, 0x00000000, /**/ 0x00000000, 0x00000000, /**/ 0x3FF00000, 0x00000000, @@ -467,7 +470,7 @@ static const union {int4 i[880]; double x[440];}sincos = { .i = { /**/ 0x3C747A10, 0x8073C259 } }; #else #ifdef LITTLE_ENDI -static const union {int4 i[880]; double x[440];} sincos = { .i = { +const union {int4 i[880]; double x[440];} __sincostab = { .i = { /**/ 0x00000000, 0x00000000, /**/ 0x00000000, 0x00000000, /**/ 0x00000000, 0x3FF00000, diff --git a/libc/sysdeps/ieee754/dbl-64/slowexp.c b/libc/sysdeps/ieee754/dbl-64/slowexp.c index 78c107f70..6a6bce31b 100644 --- a/libc/sysdeps/ieee754/dbl-64/slowexp.c +++ b/libc/sysdeps/ieee754/dbl-64/slowexp.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -31,10 +31,16 @@ #include "mpa.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + void __mpexp(mp_no *x, mp_no *y, int p); /*Converting from double precision to Multi-precision and calculating e^x */ -double __slowexp(double x) { +double +SECTION +__slowexp(double x) { double w,z,res,eps=3.0e-26; #if 0 double y; @@ -47,7 +53,7 @@ double __slowexp(double x) { p=6; __dbl_mp(x,&mpx,p); /* Convert a double precision number x */ - /* into a multiple precision number mpx with prec. p. */ + /* into a multiple precision number mpx with prec. p. */ __mpexp(&mpx, &mpy, p); /* Multi-Precision exponential function */ __dbl_mp(eps,&mpeps,p); __mul(&mpeps,&mpy,&mpcor,p); diff --git a/libc/sysdeps/ieee754/dbl-64/slowpow.c b/libc/sysdeps/ieee754/dbl-64/slowpow.c index e11a532bf..0c57e6d4f 100644 --- a/libc/sysdeps/ieee754/dbl-64/slowpow.c +++ b/libc/sysdeps/ieee754/dbl-64/slowpow.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation + * Copyright (C) 2001, 2011 Free Software Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -35,12 +35,18 @@ #include "mpa.h" #include "math_private.h" +#ifndef SECTION +# define SECTION +#endif + void __mpexp(mp_no *x, mp_no *y, int p); void __mplog(mp_no *x, mp_no *y, int p); double ulog(double); double __halfulp(double x,double y); -double __slowpow(double x, double y, double z) { +double +SECTION +__slowpow(double x, double y, double z) { double res,res1; mp_no mpx, mpy, mpz,mpw,mpp,mpr,mpr1; static const mp_no eps = {-3,{1.0,4.0}}; @@ -48,7 +54,7 @@ double __slowpow(double x, double y, double z) { res = __halfulp(x,y); /* halfulp() returns -10 or x^y */ if (res >= 0) return res; /* if result was really computed by halfulp */ - /* else, if result was not really computed by halfulp */ + /* else, if result was not really computed by halfulp */ p = 10; /* p=precision */ __dbl_mp(x,&mpx,p); __dbl_mp(y,&mpy,p); diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/e_fmod.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/e_fmod.c new file mode 100644 index 000000000..0e20571a7 --- /dev/null +++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/e_fmod.c @@ -0,0 +1,104 @@ +/* Rewritten for 64-bit machines by Ulrich Drepper <drepper@gmail.com>. */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * __ieee754_fmod(x,y) + * Return x mod y in exact arithmetic + * Method: shift and subtract + */ + +#include "math.h" +#include "math_private.h" + +static const double one = 1.0, Zero[] = {0.0, -0.0,}; + +double +__ieee754_fmod (double x, double y) +{ + int32_t n,i,ix,iy; + int64_t hx,hy,hz,sx; + + EXTRACT_WORDS64(hx,x); + EXTRACT_WORDS64(hy,y); + sx = hx&UINT64_C(0x8000000000000000); /* sign of x */ + hx ^=sx; /* |x| */ + hy &= UINT64_C(0x7fffffffffffffff); /* |y| */ + + /* purge off exception values */ + if(__builtin_expect(hy==0 + || hx >= UINT64_C(0x7ff0000000000000) + || hy > UINT64_C(0x7ff0000000000000), 0)) + /* y=0,or x not finite or y is NaN */ + return (x*y)/(x*y); + if(__builtin_expect(hx<=hy, 0)) { + if(hx<hy) return x; /* |x|<|y| return x */ + return Zero[(uint64_t)sx>>63]; /* |x|=|y| return x*0*/ + } + + /* determine ix = ilogb(x) */ + if(__builtin_expect(hx<UINT64_C(0x0010000000000000), 0)) { + /* subnormal x */ + for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1; + } else ix = (hx>>52)-1023; + + /* determine iy = ilogb(y) */ + if(__builtin_expect(hy<UINT64_C(0x0010000000000000), 0)) { /* subnormal y */ + for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1; + } else iy = (hy>>52)-1023; + + /* set up hx, hy and align y to x */ + if(__builtin_expect(ix >= -1022, 1)) + hx = UINT64_C(0x0010000000000000)|(UINT64_C(0x000fffffffffffff)&hx); + else { /* subnormal x, shift x to normal */ + n = -1022-ix; + hx<<=n; + } + if(__builtin_expect(iy >= -1022, 1)) + hy = UINT64_C(0x0010000000000000)|(UINT64_C(0x000fffffffffffff)&hy); + else { /* subnormal y, shift y to normal */ + n = -1022-iy; + hy<<=n; + } + + /* fix point fmod */ + n = ix - iy; + while(n--) { + hz=hx-hy; + if(hz<0){hx = hx+hx;} + else { + if(hz==0) /* return sign(x)*0 */ + return Zero[(uint64_t)sx>>63]; + hx = hz+hz; + } + } + hz=hx-hy; + if(hz>=0) {hx=hz;} + + /* convert back to floating value and restore the sign */ + if(hx==0) /* return sign(x)*0 */ + return Zero[(uint64_t)sx>>63]; + while(hx<UINT64_C(0x0010000000000000)) { /* normalize x */ + hx = hx+hx; + iy -= 1; + } + if(__builtin_expect(iy>= -1022, 1)) { /* normalize output */ + hx = ((hx-UINT64_C(0x0010000000000000))|((uint64_t)(iy+1023)<<52)); + INSERT_WORDS64(x,hx|sx); + } else { /* subnormal output */ + n = -1022 - iy; + hx>>=n; + INSERT_WORDS64(x,hx|sx); + x *= one; /* create necessary signal */ + } + return x; /* exact output */ +} +strong_alias (__ieee754_fmod, __fmod_finite) diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c index e0e71558f..346dab799 100644 --- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c +++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c @@ -32,18 +32,16 @@ __ceil(double x) EXTRACT_WORDS64(i0,x); j0 = ((i0>>52)&0x7ff)-0x3ff; if(j0<=51) { - if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0<0) {i0=INT64_C(0x8000000000000000);} - else if(i0!=0) { i0=INT64_C(0x3ff0000000000000);} - } + if(j0<0) { /* raise inexact if x != 0 */ + math_force_eval(huge+x);/* return 0*sign(x) if |x|<1 */ + if(i0<0) {i0=INT64_C(0x8000000000000000);} + else if(i0!=0) { i0=INT64_C(0x3ff0000000000000);} } else { i = INT64_C(0x000fffffffffffff)>>j0; if((i0&i)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0>0) i0 += UINT64_C(0x0010000000000000)>>j0; - i0 &= (~i); - } + math_force_eval(huge+x); /* raise inexact flag */ + if(i0>0) i0 += UINT64_C(0x0010000000000000)>>j0; + i0 &= (~i); } } else { if(j0==0x400) return x+x; /* inf or NaN */ diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c index 8b7300bb9..5df4ab52f 100644 --- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c +++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c @@ -54,18 +54,16 @@ __floor (double x) int32_t j0 = ((i0>>52)&0x7ff)-0x3ff; if(__builtin_expect(j0<52, 1)) { if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0>=0) {i0=0;} - else if((i0&0x7fffffffffffffffl)!=0) - { i0=0xbff0000000000000l;} - } + math_force_eval(huge+x);/* return 0*sign(x) if |x|<1 */ + if(i0>=0) {i0=0;} + else if((i0&0x7fffffffffffffffl)!=0) + { i0=0xbff0000000000000l;} } else { uint64_t i = (0x000fffffffffffffl)>>j0; if((i0&i)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0<0) i0 += (0x0010000000000000l)>>j0; - i0 &= (~i); - } + math_force_eval(huge+x); /* raise inexact flag */ + if(i0<0) i0 += (0x0010000000000000l)>>j0; + i0 &= (~i); } INSERT_WORDS64(x,i0); } else if (j0==0x400) diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c new file mode 100644 index 000000000..011401ebc --- /dev/null +++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c @@ -0,0 +1,112 @@ +/* Compute remainder and a congruent to the quotient. + Copyright (C) 1997, 2011 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, write to the Free + Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA + 02111-1307 USA. */ + +#include <math.h> + +#include "math_private.h" + + +static const double zero = 0.0; + + +double +__remquo (double x, double y, int *quo) +{ + int64_t hx, hy; + uint64_t sx, qs; + int cquo; + + EXTRACT_WORDS64 (hx, x); + EXTRACT_WORDS64 (hy, y); + sx = hx & UINT64_C(0x8000000000000000); + qs = sx ^ (hy & UINT64_C(0x8000000000000000)); + hy &= UINT64_C(0x7fffffffffffffff); + hx &= UINT64_C(0x7fffffffffffffff); + + /* Purge off exception values. */ + if (__builtin_expect (hy == 0, 0)) + return (x * y) / (x * y); /* y = 0 */ + if (__builtin_expect (hx >= UINT64_C(0x7ff0000000000000) /* x not finite */ + || hy > UINT64_C(0x7ff0000000000000), 0))/* y is NaN */ + return (x * y) / (x * y); + + if (hy <= UINT64_C(0x7fbfffffffffffff)) + x = __ieee754_fmod (x, 8 * y); /* now x < 8y */ + + if (__builtin_expect (hx == hy, 0)) + { + *quo = qs ? -1 : 1; + return zero * x; + } + + INSERT_WORDS64 (x, hx); + INSERT_WORDS64 (y, hy); + cquo = 0; + + if (x >= 4 * y) + { + x -= 4 * y; + cquo += 4; + } + if (x >= 2 * y) + { + x -= 2 * y; + cquo += 2; + } + + if (hy < UINT64_C(0x0020000000000000)) + { + if (x + x > y) + { + x -= y; + ++cquo; + if (x + x >= y) + { + x -= y; + ++cquo; + } + } + } + else + { + double y_half = 0.5 * y; + if (x > y_half) + { + x -= y; + ++cquo; + if (x >= y_half) + { + x -= y; + ++cquo; + } + } + } + + *quo = qs ? -cquo : cquo; + + if (sx) + x = -x; + return x; +} +weak_alias (__remquo, remquo) +#ifdef NO_LONG_DOUBLE +strong_alias (__remquo, __remquol) +weak_alias (__remquo, remquol) +#endif diff --git a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c index 5bd857910..25ff85928 100644 --- a/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c +++ b/libc/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c @@ -1,5 +1,5 @@ /* Round double to integer away from zero. - Copyright (C) 1997, 2009 Free Software Foundation, Inc. + Copyright (C) 1997, 2009, 2011 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -37,12 +37,11 @@ __round (double x) { if (j0 < 0) { - if (huge + x > 0.0) - { - i0 &= UINT64_C(0x8000000000000000); - if (j0 == -1) - i0 |= UINT64_C(0x3ff0000000000000); - } + math_force_eval (huge + x); + + i0 &= UINT64_C(0x8000000000000000); + if (j0 == -1) + i0 |= UINT64_C(0x3ff0000000000000); } else { @@ -50,12 +49,11 @@ __round (double x) if ((i0 & i) == 0) /* X is integral. */ return x; - if (huge + x > 0.0) - { - /* Raise inexact if x != 0. */ - i0 += UINT64_C(0x0008000000000000) >> j0; - i0 &= ~i; - } + math_force_eval (huge + x); + + /* Raise inexact if x != 0. */ + i0 += UINT64_C(0x0008000000000000) >> j0; + i0 &= ~i; } } else diff --git a/libc/sysdeps/ieee754/flt-32/e_atanhf.c b/libc/sysdeps/ieee754/flt-32/e_atanhf.c index ddd18ab30..d98a11ed6 100644 --- a/libc/sysdeps/ieee754/flt-32/e_atanhf.c +++ b/libc/sysdeps/ieee754/flt-32/e_atanhf.c @@ -49,8 +49,11 @@ __ieee754_atanhf (float x) float t; if (xa < 0.5f) { - if (__builtin_expect (xa < 0x1.0p-28f, 0) && (huge + x) > 0.0f) - return x; + if (__builtin_expect (xa < 0x1.0p-28f, 0)) + { + math_force_eval (huge + x); + return x; + } t = xa + xa; t = 0.5f * __log1pf (t + t * xa / (1.0f - xa)); diff --git a/libc/sysdeps/ieee754/flt-32/e_j0f.c b/libc/sysdeps/ieee754/flt-32/e_j0f.c index d2da43f92..181c2e46d 100644 --- a/libc/sysdeps/ieee754/flt-32/e_j0f.c +++ b/libc/sysdeps/ieee754/flt-32/e_j0f.c @@ -66,10 +66,9 @@ __ieee754_j0f(float x) return z; } if(ix<0x39000000) { /* |x| < 2**-13 */ - if(huge+x>one) { /* raise inexact if x != 0 */ - if(ix<0x32000000) return one; /* |x|<2**-27 */ - else return one - (float)0.25*x*x; - } + math_force_eval(huge+x>one); /* raise inexact if x != 0 */ + if(ix<0x32000000) return one; /* |x|<2**-27 */ + else return one - (float)0.25*x*x; } z = x*x; r = z*(R02+z*(R03+z*(R04+z*R05))); diff --git a/libc/sysdeps/ieee754/flt-32/s_ceilf.c b/libc/sysdeps/ieee754/flt-32/s_ceilf.c index 8a83201c1..af926600b 100644 --- a/libc/sysdeps/ieee754/flt-32/s_ceilf.c +++ b/libc/sysdeps/ieee754/flt-32/s_ceilf.c @@ -29,17 +29,15 @@ __ceilf(float x) j0 = ((i0>>23)&0xff)-0x7f; if(j0<23) { if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0<0) {i0=0x80000000;} - else if(i0!=0) { i0=0x3f800000;} - } + math_force_eval(huge+x);/* return 0*sign(x) if |x|<1 */ + if(i0<0) {i0=0x80000000;} + else if(i0!=0) { i0=0x3f800000;} } else { i = (0x007fffff)>>j0; if((i0&i)==0) return x; /* x is integral */ - if(huge+x>(float)0.0) { /* raise inexact flag */ - if(i0>0) i0 += (0x00800000)>>j0; - i0 &= (~i); - } + math_force_eval(huge+x); /* raise inexact flag */ + if(i0>0) i0 += (0x00800000)>>j0; + i0 &= (~i); } } else { if(__builtin_expect(j0==0x80, 0)) return x+x; /* inf or NaN */ diff --git a/libc/sysdeps/ieee754/flt-32/s_expm1f.c b/libc/sysdeps/ieee754/flt-32/s_expm1f.c index 3f4536b90..1d707120b 100644 --- a/libc/sysdeps/ieee754/flt-32/s_expm1f.c +++ b/libc/sysdeps/ieee754/flt-32/s_expm1f.c @@ -13,22 +13,14 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_expm1f.c,v 1.5 1995/05/10 20:47:11 jtc Exp $"; -#endif - #include <errno.h> #include "math.h" #include "math_private.h" -static const volatile float huge = 1.0e+30; -static const volatile float tiny = 1.0e-30; +static const float huge = 1.0e+30; +static const float tiny = 1.0e-30; -#ifdef __STDC__ static const float -#else -static float -#endif one = 1.0, o_threshold = 8.8721679688e+01,/* 0x42b17180 */ ln2_hi = 6.9313812256e-01,/* 0x3f317180 */ @@ -41,12 +33,8 @@ Q3 = -7.9365076090e-05, /* 0xb8a670cd */ Q4 = 4.0082177293e-06, /* 0x36867e54 */ Q5 = -2.0109921195e-07; /* 0xb457edbb */ -#ifdef __STDC__ - float __expm1f(float x) -#else - float __expm1f(x) - float x; -#endif +float +__expm1f(float x) { float y,hi,lo,c,t,e,hxs,hfx,r1; int32_t k,xsb; @@ -60,17 +48,17 @@ Q5 = -2.0109921195e-07; /* 0xb457edbb */ /* filter out huge and non-finite argument */ if(hx >= 0x4195b844) { /* if |x|>=27*ln2 */ if(hx >= 0x42b17218) { /* if |x|>=88.721... */ - if(hx>0x7f800000) - return x+x; /* NaN */ + if(hx>0x7f800000) + return x+x; /* NaN */ if(hx==0x7f800000) return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */ - if(x > o_threshold) { + if(x > o_threshold) { __set_errno (ERANGE); return huge*huge; /* overflow */ } } if(xsb!=0) { /* x < -27*ln2, return -1.0 with inexact */ - if(x+tiny<(float)0.0) /* raise inexact */ + math_force_eval(x+tiny);/* raise inexact */ return tiny-one; /* return -1 */ } } @@ -91,7 +79,7 @@ Q5 = -2.0109921195e-07; /* 0xb457edbb */ x = hi - lo; c = (hi-x)-lo; } - else if(hx < 0x33000000) { /* when |x|<2**-25, return x */ + else if(hx < 0x33000000) { /* when |x|<2**-25, return x */ t = huge+x; /* return x with inexact flags when x!=0 */ return x - (t-(huge+x)); } @@ -109,28 +97,28 @@ Q5 = -2.0109921195e-07; /* 0xb457edbb */ e -= hxs; if(k== -1) return (float)0.5*(x-e)-(float)0.5; if(k==1) { - if(x < (float)-0.25) return -(float)2.0*(e-(x+(float)0.5)); - else return one+(float)2.0*(x-e); + if(x < (float)-0.25) return -(float)2.0*(e-(x+(float)0.5)); + else return one+(float)2.0*(x-e); } if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */ - int32_t i; - y = one-(e-x); + int32_t i; + y = one-(e-x); GET_FLOAT_WORD(i,y); SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ - return y-one; + return y-one; } t = one; if(k<23) { - int32_t i; - SET_FLOAT_WORD(t,0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */ - y = t-(e-x); + int32_t i; + SET_FLOAT_WORD(t,0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */ + y = t-(e-x); GET_FLOAT_WORD(i,y); SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ } else { - int32_t i; + int32_t i; SET_FLOAT_WORD(t,((0x7f-k)<<23)); /* 2^-k */ - y = x-(e+t); - y += one; + y = x-(e+t); + y += one; GET_FLOAT_WORD(i,y); SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ } diff --git a/libc/sysdeps/ieee754/flt-32/s_floorf.c b/libc/sysdeps/ieee754/flt-32/s_floorf.c index dd19c6bc5..de160d211 100644 --- a/libc/sysdeps/ieee754/flt-32/s_floorf.c +++ b/libc/sysdeps/ieee754/flt-32/s_floorf.c @@ -36,18 +36,16 @@ __floorf(float x) j0 = ((i0>>23)&0xff)-0x7f; if(j0<23) { if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0>=0) {i0=0;} - else if((i0&0x7fffffff)!=0) - { i0=0xbf800000;} - } + math_force_eval(huge+x);/* return 0*sign(x) if |x|<1 */ + if(i0>=0) {i0=0;} + else if((i0&0x7fffffff)!=0) + { i0=0xbf800000;} } else { i = (0x007fffff)>>j0; if((i0&i)==0) return x; /* x is integral */ - if(huge+x>(float)0.0) { /* raise inexact flag */ - if(i0<0) i0 += (0x00800000)>>j0; - i0 &= (~i); - } + math_force_eval(huge+x); /* raise inexact flag */ + if(i0<0) i0 += (0x00800000)>>j0; + i0 &= (~i); } } else { if(__builtin_expect(j0==0x80, 0)) return x+x; /* inf or NaN */ diff --git a/libc/sysdeps/ieee754/flt-32/s_log1pf.c b/libc/sysdeps/ieee754/flt-32/s_log1pf.c index bd3d57635..9e555ed57 100644 --- a/libc/sysdeps/ieee754/flt-32/s_log1pf.c +++ b/libc/sysdeps/ieee754/flt-32/s_log1pf.c @@ -13,18 +13,10 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_log1pf.c,v 1.4 1995/05/10 20:47:48 jtc Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const float -#else -static float -#endif ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ two25 = 3.355443200e+07, /* 0x4c000000 */ @@ -36,18 +28,10 @@ Lp5 = 1.8183572590e-01, /* 3E3A3325 */ Lp6 = 1.5313838422e-01, /* 3E1CD04F */ Lp7 = 1.4798198640e-01; /* 3E178897 */ -#ifdef __STDC__ static const float zero = 0.0; -#else -static float zero = 0.0; -#endif -#ifdef __STDC__ - float __log1pf(float x) -#else - float __log1pf(x) - float x; -#endif +float +__log1pf(float x) { float hfsq,f,c,s,z,R,u; int32_t k,hx,hu,ax; @@ -62,8 +46,8 @@ static float zero = 0.0; else return (x-x)/(x-x); /* log1p(x<-1)=NaN */ } if(ax<0x31000000) { /* |x| < 2**-29 */ - if(two25+x>zero /* raise inexact */ - &&ax<0x24800000) /* |x| < 2**-54 */ + math_force_eval(two25+x); /* raise inexact */ + if (ax<0x24800000) /* |x| < 2**-54 */ return x; else return x - x*x*(float)0.5; @@ -76,37 +60,37 @@ static float zero = 0.0; if(hx<0x5a000000) { u = (float)1.0+x; GET_FLOAT_WORD(hu,u); - k = (hu>>23)-127; + k = (hu>>23)-127; /* correction term */ - c = (k>0)? (float)1.0-(u-x):x-(u-(float)1.0); + c = (k>0)? (float)1.0-(u-x):x-(u-(float)1.0); c /= u; } else { u = x; GET_FLOAT_WORD(hu,u); - k = (hu>>23)-127; + k = (hu>>23)-127; c = 0; } hu &= 0x007fffff; if(hu<0x3504f7) { - SET_FLOAT_WORD(u,hu|0x3f800000);/* normalize u */ + SET_FLOAT_WORD(u,hu|0x3f800000);/* normalize u */ } else { - k += 1; + k += 1; SET_FLOAT_WORD(u,hu|0x3f000000); /* normalize u/2 */ - hu = (0x00800000-hu)>>2; + hu = (0x00800000-hu)>>2; } f = u-(float)1.0; } hfsq=(float)0.5*f*f; if(hu==0) { /* |f| < 2**-20 */ if(f==zero) { - if(k==0) return zero; + if(k==0) return zero; else {c += k*ln2_lo; return k*ln2_hi+c;} } R = hfsq*((float)1.0-(float)0.66666666666666666*f); if(k==0) return f-R; else - return k*ln2_hi-((R-(k*ln2_lo+c))-f); + return k*ln2_hi-((R-(k*ln2_lo+c))-f); } - s = f/((float)2.0+f); + s = f/((float)2.0+f); z = s*s; R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7)))))); if(k==0) return f-(hfsq-s*(hfsq+R)); else diff --git a/libc/sysdeps/ieee754/flt-32/s_roundf.c b/libc/sysdeps/ieee754/flt-32/s_roundf.c index 0b2498779..09b38cdc2 100644 --- a/libc/sysdeps/ieee754/flt-32/s_roundf.c +++ b/libc/sysdeps/ieee754/flt-32/s_roundf.c @@ -1,5 +1,5 @@ /* Round float to integer away from zero. - Copyright (C) 1997 Free Software Foundation, Inc. + Copyright (C) 1997, 2011 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -37,12 +37,11 @@ __roundf (float x) { if (j0 < 0) { - if (huge + x > 0.0F) - { - i0 &= 0x80000000; - if (j0 == -1) - i0 |= 0x3f800000; - } + math_force_eval (huge + x > 0.0F); + + i0 &= 0x80000000; + if (j0 == -1) + i0 |= 0x3f800000; } else { @@ -50,12 +49,11 @@ __roundf (float x) if ((i0 & i) == 0) /* X is integral. */ return x; - if (huge + x > 0.0F) - { - /* Raise inexact if x != 0. */ - i0 += 0x00400000 >> j0; - i0 &= ~i; - } + math_force_eval (huge + x > 0.0F); + + /* Raise inexact if x != 0. */ + i0 += 0x00400000 >> j0; + i0 &= ~i; } } else diff --git a/libc/sysdeps/ieee754/ldbl-96/e_atanhl.c b/libc/sysdeps/ieee754/ldbl-96/e_atanhl.c index 5a2aebef3..0f3c7fb59 100644 --- a/libc/sysdeps/ieee754/ldbl-96/e_atanhl.c +++ b/libc/sysdeps/ieee754/ldbl-96/e_atanhl.c @@ -52,7 +52,10 @@ __ieee754_atanhl(long double x) return (x-x)/(x-x); if(ix==0x3fff) return x/zero; - if(ix<0x3fe3&&(huge+x)>zero) return x; /* x<2**-28 */ + if(ix<0x3fe3) { + math_force_eval(huge+x); + return x; /* x<2**-28 */ + } SET_LDOUBLE_EXP(x,ix); if(ix<0x3ffe) { /* x < 0.5 */ t = x+x; diff --git a/libc/sysdeps/ieee754/ldbl-96/e_j0l.c b/libc/sysdeps/ieee754/ldbl-96/e_j0l.c index ce1f0f756..abf4f109f 100644 --- a/libc/sysdeps/ieee754/ldbl-96/e_j0l.c +++ b/libc/sysdeps/ieee754/ldbl-96/e_j0l.c @@ -144,13 +144,12 @@ __ieee754_j0l (long double x) } if (__builtin_expect (ix < 0x3fef, 0)) /* |x| < 2**-16 */ { - if (huge + x > one) - { /* raise inexact if x != 0 */ - if (ix < 0x3fde) /* |x| < 2^-33 */ - return one; - else - return one - 0.25 * x * x; - } + /* raise inexact if x != 0 */ + math_force_eval (huge + x); + if (ix < 0x3fde) /* |x| < 2^-33 */ + return one; + else + return one - 0.25 * x * x; } z = x * x; r = z * (R[0] + z * (R[1] + z * (R[2] + z * (R[3] + z * R[4])))); diff --git a/libc/sysdeps/ieee754/ldbl-96/s_roundl.c b/libc/sysdeps/ieee754/ldbl-96/s_roundl.c index f1399cc20..833ae0d78 100644 --- a/libc/sysdeps/ieee754/ldbl-96/s_roundl.c +++ b/libc/sysdeps/ieee754/ldbl-96/s_roundl.c @@ -1,5 +1,5 @@ /* Round long double to integer away from zero. - Copyright (C) 1997, 2007 Free Software Foundation, Inc. + Copyright (C) 1997, 2007, 2011 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -38,15 +38,13 @@ __roundl (long double x) { if (j0 < 0) { - if (huge + x > 0.0) + math_force_eval (huge + x); + se &= 0x8000; + i0 = i1 = 0; + if (j0 == -1) { - se &= 0x8000; - i0 = i1 = 0; - if (j0 == -1) - { - se |= 0x3fff; - i0 = 0x80000000; - } + se |= 0x3fff; + i0 = 0x80000000; } } else @@ -55,15 +53,14 @@ __roundl (long double x) if (((i0 & i) | i1) == 0) /* X is integral. */ return x; - if (huge + x > 0.0) - { - /* Raise inexact if x != 0. */ - u_int32_t j = i0 + (0x40000000 >> j0); - if (j < i0) - se += 1; - i0 = (j & ~i) | 0x80000000; - i1 = 0; - } + + /* Raise inexact if x != 0. */ + math_force_eval (huge + x); + u_int32_t j = i0 + (0x40000000 >> j0); + if (j < i0) + se += 1; + i0 = (j & ~i) | 0x80000000; + i1 = 0; } } else if (j0 > 62) @@ -81,22 +78,20 @@ __roundl (long double x) /* X is integral. */ return x; - if (huge + x > 0.0) + math_force_eval (huge + x); + /* Raise inexact if x != 0. */ + u_int32_t j = i1 + (1 << (62 - j0)); + if (j < i1) { - /* Raise inexact if x != 0. */ - u_int32_t j = i1 + (1 << (62 - j0)); - if (j < i1) + u_int32_t k = i0 + 1; + if (k < i0) { - u_int32_t k = i0 + 1; - if (k < i0) - { - se += 1; - k |= 0x80000000; - } - i0 = k; + se += 1; + k |= 0x80000000; } - i1 = j; + i0 = k; } + i1 = j; i1 &= ~i; } diff --git a/libc/sysdeps/x86_64/dla.h b/libc/sysdeps/x86_64/dla.h deleted file mode 100644 index 7aa06e5f6..000000000 --- a/libc/sysdeps/x86_64/dla.h +++ /dev/null @@ -1,6 +0,0 @@ -#if defined __FMA4__ || defined __FMA__ -# define DLA_FMS(x,y,z) \ - __builtin_fma (x, y, -z) -#endif - -#include "sysdeps/ieee754/dbl-64/dla.h" diff --git a/libc/sysdeps/x86_64/fpu/dla.h b/libc/sysdeps/x86_64/fpu/dla.h new file mode 100644 index 000000000..fa2d52bbf --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/dla.h @@ -0,0 +1,17 @@ +#include <features.h> + +#ifdef __FMA4__ +# if __GNUC_PREREQ (4, 6) +# define DLA_FMS(x,y,z) \ + __builtin_fma (x, y, -(z)) +# else +# define DLA_FMS(x,y,z) \ + ({ double __z; \ + asm ("vfmsubsd %3, %2, %1, %0" \ + : "=x" (__z) \ + : "x" ((double) (x)), "xm" ((double) (y)) , "x" ((double) (z))); \ + __z; }) +# endif +#endif + +#include "sysdeps/ieee754/dbl-64/dla.h" diff --git a/libc/sysdeps/x86_64/fpu/math_private.h b/libc/sysdeps/x86_64/fpu/math_private.h index d3d84cfda..7f52d5ee5 100644 --- a/libc/sysdeps/x86_64/fpu/math_private.h +++ b/libc/sysdeps/x86_64/fpu/math_private.h @@ -1,59 +1,67 @@ #ifndef _MATH_PRIVATE_H #define math_opt_barrier(x) \ -({ __typeof(x) __x; \ - if (sizeof (x) <= sizeof (double)) \ - __asm ("" : "=x" (__x) : "0" (x)); \ - else \ - __asm ("" : "=t" (__x) : "0" (x)); \ - __x; }) + ({ __typeof(x) __x; \ + if (sizeof (x) <= sizeof (double)) \ + __asm ("" : "=x" (__x) : "0" (x)); \ + else \ + __asm ("" : "=t" (__x) : "0" (x)); \ + __x; }) #define math_force_eval(x) \ -do \ - { \ - if (sizeof (x) <= sizeof (double)) \ - __asm __volatile ("" : : "x" (x)); \ - else \ - __asm __volatile ("" : : "f" (x)); \ - } \ -while (0) + do { \ + if (sizeof (x) <= sizeof (double)) \ + __asm __volatile ("" : : "x" (x)); \ + else \ + __asm __volatile ("" : : "f" (x)); \ + } while (0) #include <math/math_private.h> /* We can do a few things better on x86-64. */ +#ifdef __AVX__ +# define MOVD "vmovd" +#else +# define MOVD "movd" +#endif + /* Direct movement of float into integer register. */ #undef EXTRACT_WORDS64 -#define EXTRACT_WORDS64(i,d) \ -do { \ - long int i_; \ - asm ("movd %1, %0" : "=rm" (i_) : "x" (d)); \ - (i) = i_; \ -} while (0) +#define EXTRACT_WORDS64(i, d) \ + do { \ + long int i_; \ + asm (MOVD " %1, %0" : "=rm" (i_) : "x" ((double) (d))); \ + (i) = i_; \ + } while (0) /* And the reverse. */ #undef INSERT_WORDS64 -#define INSERT_WORDS64(d,i) \ -do { \ - long int i_ = i; \ - asm ("movd %1, %0" : "=x" (d) : "rm" (i_)); \ -} while (0) +#define INSERT_WORDS64(d, i) \ + do { \ + long int i_ = i; \ + double d__; \ + asm (MOVD " %1, %0" : "=x" (d__) : "rm" (i_)); \ + d = d__; \ + } while (0) /* Direct movement of float into integer register. */ #undef GET_FLOAT_WORD -#define GET_FLOAT_WORD(i,d) \ -do { \ - int i_; \ - asm ("movd %1, %0" : "=rm" (i_) : "x" (d)); \ - (i) = i_; \ -} while (0) +#define GET_FLOAT_WORD(i, d) \ + do { \ + int i_; \ + asm (MOVD " %1, %0" : "=rm" (i_) : "x" ((float) (d))); \ + (i) = i_; \ + } while (0) /* And the reverse. */ #undef SET_FLOAT_WORD -#define SET_FLOAT_WORD(d,i) \ -do { \ - int i_ = i; \ - asm ("movd %1, %0" : "=x" (d) : "rm" (i_)); \ -} while (0) +#define SET_FLOAT_WORD(f, i) \ + do { \ + int i_ = i; \ + float f__; \ + asm (MOVD " %1, %0" : "=x" (f__) : "rm" (i_)); \ + f = f__; \ + } while (0) #endif @@ -78,14 +86,25 @@ do { \ ({ int __di; GET_FLOAT_WORD (__di, (float) d); \ (__di & 0x7fffffff) < 0x7f800000; }) -#define __ieee754_sqrt(d) \ +#ifdef __AVX__ +# define __ieee754_sqrt(d) \ + ({ double __res; \ + asm ("vsqrtsd %1, %0, %0" : "=x" (__res) : "xm" ((double) (d))); \ + __res; }) +# define __ieee754_sqrtf(d) \ + ({ float __res; \ + asm ("vsqrtss %1, %0, %0" : "=x" (__res) : "xm" ((float) (d))); \ + __res; }) +#else +# define __ieee754_sqrt(d) \ ({ double __res; \ asm ("sqrtsd %1, %0" : "=x" (__res) : "xm" ((double) (d))); \ __res; }) -#define __ieee754_sqrtf(d) \ +# define __ieee754_sqrtf(d) \ ({ float __res; \ asm ("sqrtss %1, %0" : "=x" (__res) : "xm" ((float) (d))); \ __res; }) +#endif #define __ieee754_sqrtl(d) \ ({ long double __res; \ asm ("fsqrt" : "=t" (__res) : "0" ((long double) (d))); \ @@ -93,29 +112,57 @@ do { \ #ifdef __SSE4_1__ # ifndef __rint -# define __rint(d) \ +# ifdef __AVX__ +# define __rint(d) \ + ({ double __res; \ + asm ("vroundsd $4, %1, %0, %0" : "=x" (__res) : "xm" ((double) (d))); \ + __res; }) +# else +# define __rint(d) \ ({ double __res; \ asm ("roundsd $4, %1, %0" : "=x" (__res) : "xm" ((double) (d))); \ __res; }) +# endif # endif # ifndef __rintf -# define __rintf(d) \ +# ifdef __AVX__ +# define __rintf(d) \ + ({ float __res; \ + asm ("vroundss $4, %1, %0, %0" : "=x" (__res) : "xm" ((float) (d))); \ + __res; }) +# else +# define __rintf(d) \ ({ float __res; \ asm ("roundss $4, %1, %0" : "=x" (__res) : "xm" ((float) (d))); \ __res; }) +# endif # endif # ifndef __floor -# define __floor(d) \ +# ifdef __AVX__ +# define __floor(d) \ + ({ double __res; \ + asm ("vroundsd $1, %1, %0, %0" : "=x" (__res) : "xm" ((double) (d))); \ + __res; }) +# else +# define __floor(d) \ ({ double __res; \ asm ("roundsd $1, %1, %0" : "=x" (__res) : "xm" ((double) (d))); \ __res; }) +# endif # endif # ifndef __floorf -# define __floorf(d) \ +# ifdef __AVX__ +# define __floorf(d) \ + ({ float __res; \ + asm ("vroundss $1, %1, %0, %0" : "=x" (__res) : "xm" ((float) (d))); \ + __res; }) +# else +# define __floorf(d) \ ({ float __res; \ asm ("roundss $1, %1, %0" : "=x" (__res) : "xm" ((float) (d))); \ __res; }) +# endif # endif #endif @@ -146,7 +193,17 @@ do { \ // #define libc_fesetroundl(r) (void) fesetround (r) #undef libc_feholdexcept -#define libc_feholdexcept(e) \ +#ifdef __AVX__ +# define libc_feholdexcept(e) \ + do { \ + unsigned int mxcsr; \ + asm ("vstmxcsr %0" : "=m" (*&mxcsr)); \ + (e)->__mxcsr = mxcsr; \ + mxcsr = (mxcsr | 0x1f80) & ~0x3f; \ + asm volatile ("vldmxcsr %0" : : "m" (*&mxcsr)); \ + } while (0) +#else +# define libc_feholdexcept(e) \ do { \ unsigned int mxcsr; \ asm ("stmxcsr %0" : "=m" (*&mxcsr)); \ @@ -154,12 +211,23 @@ do { \ mxcsr = (mxcsr | 0x1f80) & ~0x3f; \ asm volatile ("ldmxcsr %0" : : "m" (*&mxcsr)); \ } while (0) +#endif #undef libc_feholdexceptf #define libc_feholdexceptf(e) libc_feholdexcept (e) // #define libc_feholdexceptl(e) (void) feholdexcept (e) #undef libc_feholdexcept_setround -#define libc_feholdexcept_setround(e, r) \ +#ifdef __AVX__ +# define libc_feholdexcept_setround(e, r) \ + do { \ + unsigned int mxcsr; \ + asm ("vstmxcsr %0" : "=m" (*&mxcsr)); \ + (e)->__mxcsr = mxcsr; \ + mxcsr = ((mxcsr | 0x1f80) & ~0x603f) | ((r) << 3); \ + asm volatile ("vldmxcsr %0" : : "m" (*&mxcsr)); \ + } while (0) +#else +# define libc_feholdexcept_setround(e, r) \ do { \ unsigned int mxcsr; \ asm ("stmxcsr %0" : "=m" (*&mxcsr)); \ @@ -167,33 +235,55 @@ do { \ mxcsr = ((mxcsr | 0x1f80) & ~0x603f) | ((r) << 3); \ asm volatile ("ldmxcsr %0" : : "m" (*&mxcsr)); \ } while (0) +#endif #undef libc_feholdexcept_setroundf #define libc_feholdexcept_setroundf(e, r) libc_feholdexcept_setround (e, r) // #define libc_feholdexcept_setroundl(e, r) ... #undef libc_fetestexcept -#define libc_fetestexcept(e) \ - ({ unsigned int mxcsr; asm volatile ("stmxcsr %0" : "=m" (*&mxcsr)); \ +#ifdef __AVX__ +# define libc_fetestexcept(e) \ + ({ unsigned int mxcsr; asm volatile ("vstmxcsr %0" : "=m" (*&mxcsr)); \ mxcsr & (e) & FE_ALL_EXCEPT; }) +#else +# define libc_fetestexcept(e) \ + ({ unsigned int mxcsr; asm volatile ("stmxcsr %0" : "=m" (*&mxcsr)); \ + mxcsr & (e) & FE_ALL_EXCEPT; }) +#endif #undef libc_fetestexceptf #define libc_fetestexceptf(e) libc_fetestexcept (e) // #define libc_fetestexceptl(e) fetestexcept (e) #undef libc_fesetenv -#define libc_fesetenv(e) \ +#ifdef __AVX__ +# define libc_fesetenv(e) \ + asm volatile ("vldmxcsr %0" : : "m" ((e)->__mxcsr)) +#else +# define libc_fesetenv(e) \ asm volatile ("ldmxcsr %0" : : "m" ((e)->__mxcsr)) +#endif #undef libc_fesetenvf #define libc_fesetenvf(e) libc_fesetenv (e) // #define libc_fesetenvl(e) (void) fesetenv (e) #undef libc_feupdateenv -#define libc_feupdateenv(e) \ +#ifdef __AVX__ +# define libc_feupdateenv(e) \ + do { \ + unsigned int mxcsr; \ + asm volatile ("vstmxcsr %0" : "=m" (*&mxcsr)); \ + asm volatile ("vldmxcsr %0" : : "m" ((e)->__mxcsr)); \ + __feraiseexcept (mxcsr & FE_ALL_EXCEPT); \ + } while (0) +#else +# define libc_feupdateenv(e) \ do { \ unsigned int mxcsr; \ asm volatile ("stmxcsr %0" : "=m" (*&mxcsr)); \ asm volatile ("ldmxcsr %0" : : "m" ((e)->__mxcsr)); \ __feraiseexcept (mxcsr & FE_ALL_EXCEPT); \ } while (0) +#endif #undef libc_feupdateenvf #define libc_feupdateenvf(e) libc_feupdateenv (e) // #define libc_feupdateenvl(e) (void) feupdateenv (e) diff --git a/libc/sysdeps/x86_64/fpu/multiarch/Makefile b/libc/sysdeps/x86_64/fpu/multiarch/Makefile index bd07e98e2..70cb740aa 100644 --- a/libc/sysdeps/x86_64/fpu/multiarch/Makefile +++ b/libc/sysdeps/x86_64/fpu/multiarch/Makefile @@ -1,4 +1,36 @@ ifeq ($(subdir),math) libm-sysdep_routines += s_floor-c s_ceil-c s_floorf-c s_ceilf-c \ s_rint-c s_rintf-c s_nearbyint-c s_nearbyintf-c + +ifeq ($(have-mfma4),yes) +libm-sysdep_routines += e_exp-fma4 e_log-fma4 e_pow-fma4 s_atan-fma4 \ + e_asin-fma4 e_atan2-fma4 s_sin-fma4 s_tan-fma4 \ + mplog-fma4 mpa-fma4 slowexp-fma4 slowpow-fma4 \ + sincos32-fma4 doasin-fma4 dosincos-fma4 \ + brandred-fma4 halfulp-fma4 mpexp-fma4 \ + mpatan2-fma4 mpatan-fma4 mpsqrt-fma4 mptan-fma4 + +CFLAGS-brandred-fma4.c = -mfma4 +CFLAGS-doasin-fma4.c = -mfma4 +CFLAGS-dosincos-fma4.c = -mfma4 +CFLAGS-e_asin-fma4.c = -mfma4 +CFLAGS-e_atan2-fma4.c = -mfma4 +CFLAGS-e_exp-fma4.c = -mfma4 +CFLAGS-e_log-fma4.c = -mfma4 +CFLAGS-e_pow-fma4.c = -mfma4 +CFLAGS-halfulp-fma4.c = -mfma4 +CFLAGS-mpa-fma4.c = -mfma4 +CFLAGS-mpatan-fma4.c = -mfma4 +CFLAGS-mpatan2-fma4.c = -mfma4 +CFLAGS-mpexp-fma4.c = -mfma4 +CFLAGS-mplog-fma4.c = -mfma4 +CFLAGS-mpsqrt-fma4.c = -mfma4 +CFLAGS-mptan-fma4.c = -mfma4 +CFLAGS-s_atan-fma4.c = -mfma4 +CFLAGS-sincos32-fma4.c = -mfma4 +CFLAGS-slowexp-fma4.c = -mfma4 +CFLAGS-slowpow-fma4.c = -mfma4 +CLFAGS-s_sin-fma4.c = -mfma4 +CLFAGS-s_tan-fma4.c = -mfma4 +endif endif diff --git a/libc/sysdeps/x86_64/fpu/multiarch/brandred-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/brandred-fma4.c new file mode 100644 index 000000000..f4f68ac96 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/brandred-fma4.c @@ -0,0 +1,4 @@ +#define __branred __branred_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/branred.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/doasin-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/doasin-fma4.c new file mode 100644 index 000000000..53eb41947 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/doasin-fma4.c @@ -0,0 +1,4 @@ +#define __doasin __doasin_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/doasin.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/dosincos-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/dosincos-fma4.c new file mode 100644 index 000000000..1578b2fce --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/dosincos-fma4.c @@ -0,0 +1,6 @@ +#define __docos __docos_fma4 +#define __dubcos __dubcos_fma4 +#define __dubsin __dubsin_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/dosincos.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/e_asin-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/e_asin-fma4.c new file mode 100644 index 000000000..2657c31f4 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/e_asin-fma4.c @@ -0,0 +1,11 @@ +#define __ieee754_acos __ieee754_acos_fma4 +#define __ieee754_asin __ieee754_asin_fma4 +#define __cos32 __cos32_fma4 +#define __doasin __doasin_fma4 +#define __docos __docos_fma4 +#define __dubcos __dubcos_fma4 +#define __dubsin __dubsin_fma4 +#define __sin32 __sin32_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/e_asin.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/e_asin.c b/libc/sysdeps/x86_64/fpu/multiarch/e_asin.c new file mode 100644 index 000000000..8882cead9 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/e_asin.c @@ -0,0 +1,23 @@ +#ifdef HAVE_FMA4_SUPPORT +# include <init-arch.h> +# include <math_private.h> + +extern double __ieee754_acos_sse2 (double); +extern double __ieee754_acos_fma4 (double); +extern double __ieee754_asin_sse2 (double); +extern double __ieee754_asin_fma4 (double); + +libm_ifunc (__ieee754_acos, + HAS_FMA4 ? __ieee754_acos_fma4 : __ieee754_acos_sse2); +strong_alias (__ieee754_acos, __acos_finite) + +libm_ifunc (__ieee754_asin, + HAS_FMA4 ? __ieee754_asin_fma4 : __ieee754_asin_sse2); +strong_alias (__ieee754_asin, __asin_finite) + +# define __ieee754_acos __ieee754_acos_sse2 +# define __ieee754_asin __ieee754_asin_sse2 +#endif + + +#include <sysdeps/ieee754/dbl-64/e_asin.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/e_atan2-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/e_atan2-fma4.c new file mode 100644 index 000000000..f4e986293 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/e_atan2-fma4.c @@ -0,0 +1,10 @@ +#define __ieee754_atan2 __ieee754_atan2_fma4 +#define __add __add_fma4 +#define __dbl_mp __dbl_mp_fma4 +#define __dvd __dvd_fma4 +#define __mpatan2 __mpatan2_fma4 +#define __mul __mul_fma4 +#define __sub __sub_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/e_atan2.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/e_atan2.c b/libc/sysdeps/x86_64/fpu/multiarch/e_atan2.c new file mode 100644 index 000000000..12fc92906 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/e_atan2.c @@ -0,0 +1,16 @@ +#ifdef HAVE_FMA4_SUPPORT +# include <init-arch.h> +# include <math_private.h> + +extern double __ieee754_atan2_sse2 (double, double); +extern double __ieee754_atan2_fma4 (double, double); + +libm_ifunc (__ieee754_atan2, + HAS_FMA4 ? __ieee754_atan2_fma4 : __ieee754_atan2_sse2); +strong_alias (__ieee754_atan2, __atan2_finite) + +# define __ieee754_atan2 __ieee754_atan2_sse2 +#endif + + +#include <sysdeps/ieee754/dbl-64/e_atan2.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/e_exp-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/e_exp-fma4.c new file mode 100644 index 000000000..ae6eb6760 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/e_exp-fma4.c @@ -0,0 +1,6 @@ +#define __ieee754_exp __ieee754_exp_fma4 +#define __exp1 __exp1_fma4 +#define __slowexp __slowexp_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/e_exp.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/e_exp.c b/libc/sysdeps/x86_64/fpu/multiarch/e_exp.c new file mode 100644 index 000000000..fc1096b54 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/e_exp.c @@ -0,0 +1,15 @@ +#ifdef HAVE_FMA4_SUPPORT +# include <init-arch.h> +# include <math_private.h> + +extern double __ieee754_exp_sse2 (double); +extern double __ieee754_exp_fma4 (double); + +libm_ifunc (__ieee754_exp, HAS_FMA4 ? __ieee754_exp_fma4 : __ieee754_exp_sse2); +strong_alias (__ieee754_exp, __exp_finite) + +# define __ieee754_exp __ieee754_exp_sse2 +#endif + + +#include <sysdeps/ieee754/dbl-64/e_exp.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/e_log-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/e_log-fma4.c new file mode 100644 index 000000000..a2346cc61 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/e_log-fma4.c @@ -0,0 +1,8 @@ +#define __ieee754_log __ieee754_log_fma4 +#define __mplog __mplog_fma4 +#define __add __add_fma4 +#define __dbl_mp __dbl_mp_fma4 +#define __sub __sub_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/e_log.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/e_log.c b/libc/sysdeps/x86_64/fpu/multiarch/e_log.c new file mode 100644 index 000000000..c54264609 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/e_log.c @@ -0,0 +1,15 @@ +#ifdef HAVE_FMA4_SUPPORT +# include <init-arch.h> +# include <math_private.h> + +extern double __ieee754_log_sse2 (double); +extern double __ieee754_log_fma4 (double); + +libm_ifunc (__ieee754_log, HAS_FMA4 ? __ieee754_log_fma4 : __ieee754_log_sse2); +strong_alias (__ieee754_log, __log_finite) + +# define __ieee754_log __ieee754_log_sse2 +#endif + + +#include <sysdeps/ieee754/dbl-64/e_log.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/e_pow-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/e_pow-fma4.c new file mode 100644 index 000000000..5b3ea8e10 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/e_pow-fma4.c @@ -0,0 +1,6 @@ +#define __ieee754_pow __ieee754_pow_fma4 +#define __exp1 __exp1_fma4 +#define __slowpow __slowpow_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/e_pow.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/e_pow.c b/libc/sysdeps/x86_64/fpu/multiarch/e_pow.c new file mode 100644 index 000000000..a740b6c44 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/e_pow.c @@ -0,0 +1,15 @@ +#ifdef HAVE_FMA4_SUPPORT +# include <init-arch.h> +# include <math_private.h> + +extern double __ieee754_pow_sse2 (double, double); +extern double __ieee754_pow_fma4 (double, double); + +libm_ifunc (__ieee754_pow, HAS_FMA4 ? __ieee754_pow_fma4 : __ieee754_pow_sse2); +strong_alias (__ieee754_pow, __pow_finite) + +# define __ieee754_pow __ieee754_pow_sse2 +#endif + + +#include <sysdeps/ieee754/dbl-64/e_pow.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/halfulp-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/halfulp-fma4.c new file mode 100644 index 000000000..a00c17c01 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/halfulp-fma4.c @@ -0,0 +1,4 @@ +#define __halfulp __halfulp_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/halfulp.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/mpa-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/mpa-fma4.c new file mode 100644 index 000000000..f8ed8f346 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/mpa-fma4.c @@ -0,0 +1,12 @@ +#define __add __add_fma4 +#define __mul __mul_fma4 +#define __sub __sub_fma4 +#define __dbl_mp __dbl_mp_fma4 +#define __dvd __dvd_fma4 + +#define NO___CPY 1 +#define NO___MP_DBL 1 +#define NO___ACR 1 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/mpa.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/mpatan-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/mpatan-fma4.c new file mode 100644 index 000000000..fbd3bd49a --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/mpatan-fma4.c @@ -0,0 +1,10 @@ +#define __mpatan __mpatan_fma4 +#define __add __add_fma4 +#define __dvd __dvd_fma4 +#define __mpsqrt __mpsqrt_fma4 +#define __mul __mul_fma4 +#define __sub __sub_fma4 +#define AVOID_MPATAN_H 1 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/mpatan.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/mpatan2-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/mpatan2-fma4.c new file mode 100644 index 000000000..e6e44d49b --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/mpatan2-fma4.c @@ -0,0 +1,9 @@ +#define __mpatan2 __mpatan2_fma4 +#define __add __add_fma4 +#define __dvd __dvd_fma4 +#define __mpatan __mpatan_fma4 +#define __mpsqrt __mpsqrt_fma4 +#define __mul __mul_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/mpatan2.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/mpexp-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/mpexp-fma4.c new file mode 100644 index 000000000..07ca6e9ad --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/mpexp-fma4.c @@ -0,0 +1,9 @@ +#define __mpexp __mpexp_fma4 +#define __add __add_fma4 +#define __dbl_mp __dbl_mp_fma4 +#define __dvd __dvd_fma4 +#define __mul __mul_fma4 +#define AVOID_MPEXP_H 1 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/mpexp.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/mplog-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/mplog-fma4.c new file mode 100644 index 000000000..b4733118d --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/mplog-fma4.c @@ -0,0 +1,8 @@ +#define __mplog __mplog_fma4 +#define __add __add_fma4 +#define __mpexp __mpexp_fma4 +#define __mul __mul_fma4 +#define __sub __sub_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/mplog.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/mpsqrt-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/mpsqrt-fma4.c new file mode 100644 index 000000000..f8a1ba2d9 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/mpsqrt-fma4.c @@ -0,0 +1,8 @@ +#define __mpsqrt __mpsqrt_fma4 +#define __dbl_mp __dbl_mp_fma4 +#define __mul __mul_fma4 +#define __sub __sub_fma4 +#define AVOID_MPSQRT_H 1 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/mpsqrt.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/mptan-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/mptan-fma4.c new file mode 100644 index 000000000..fb4a9d48c --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/mptan-fma4.c @@ -0,0 +1,7 @@ +#define __mptan __mptan_fma4 +#define __c32 __c32_fma4 +#define __dvd __dvd_fma4 +#define __mpranred __mpranred_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/mptan.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/s_atan-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/s_atan-fma4.c new file mode 100644 index 000000000..9e83e6cda --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/s_atan-fma4.c @@ -0,0 +1,9 @@ +#define atan __atan_fma4 +#define __add __add_fma4 +#define __dbl_mp __dbl_mp_fma4 +#define __mpatan __mpatan_fma4 +#define __mul __mul_fma4 +#define __sub __sub_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/s_atan.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/s_atan.c b/libc/sysdeps/x86_64/fpu/multiarch/s_atan.c new file mode 100644 index 000000000..ffc4a56fa --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/s_atan.c @@ -0,0 +1,14 @@ +#ifdef HAVE_FMA4_SUPPORT +# include <init-arch.h> +# include <math.h> + +extern double __atan_sse2 (double); +extern double __atan_fma4 (double); + +libm_ifunc (atan, HAS_FMA4 ? __atan_fma4 : __atan_sse2); + +# define atan __atan_sse2 +#endif + + +#include <sysdeps/ieee754/dbl-64/s_atan.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/s_sin-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/s_sin-fma4.c new file mode 100644 index 000000000..2501af981 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/s_sin-fma4.c @@ -0,0 +1,12 @@ +#define __cos __cos_fma4 +#define __sin __sin_fma4 +#define __branred __branred_fma4 +#define __docos __docos_fma4 +#define __dubsin __dubsin_fma4 +#define __mpcos __mpcos_fma4 +#define __mpcos1 __mpcos1_fma4 +#define __mpsin __mpsin_fma4 +#define __mpsin1 __mpsin1_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/s_sin.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/s_sin.c b/libc/sysdeps/x86_64/fpu/multiarch/s_sin.c new file mode 100644 index 000000000..a7c35dc85 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/s_sin.c @@ -0,0 +1,22 @@ +#ifdef HAVE_FMA4_SUPPORT +# include <init-arch.h> +# include <math.h> +# undef NAN + +extern double __cos_sse2 (double); +extern double __cos_fma4 (double); +extern double __sin_sse2 (double); +extern double __sin_fma4 (double); + +libm_ifunc (__cos, HAS_FMA4 ? __cos_fma4 : __cos_sse2); +weak_alias (__cos, cos) + +libm_ifunc (__sin, HAS_FMA4 ? __sin_fma4 : __sin_sse2); +weak_alias (__sin, sin) + +# define __cos __cos_sse2 +# define __sin __sin_sse2 +#endif + + +#include <sysdeps/ieee754/dbl-64/s_sin.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/s_tan-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/s_tan-fma4.c new file mode 100644 index 000000000..d7dab3ca9 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/s_tan-fma4.c @@ -0,0 +1,9 @@ +#define tan __tan_fma4 +#define __branred __branred_fma4 +#define __dbl_mp __dbl_mp_fma4 +#define __mpranred __mpranred_fma4 +#define __mptan __mptan_fma4 +#define __sub __sub_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/s_tan.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/s_tan.c b/libc/sysdeps/x86_64/fpu/multiarch/s_tan.c new file mode 100644 index 000000000..cca02b54d --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/s_tan.c @@ -0,0 +1,14 @@ +#ifdef HAVE_FMA4_SUPPORT +# include <init-arch.h> +# include <math.h> + +extern double __tan_sse2 (double); +extern double __tan_fma4 (double); + +libm_ifunc (tan, HAS_FMA4 ? __tan_fma4 : __tan_sse2); + +# define tan __tan_sse2 +#endif + + +#include <sysdeps/ieee754/dbl-64/s_tan.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/sincos32-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/sincos32-fma4.c new file mode 100644 index 000000000..ebbfa18cc --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/sincos32-fma4.c @@ -0,0 +1,15 @@ +#define __cos32 __cos32_fma4 +#define __sin32 __sin32_fma4 +#define __c32 __c32_fma4 +#define __mpsin __mpsin_fma4 +#define __mpsin1 __mpsin1_fma4 +#define __mpcos __mpcos_fma4 +#define __mpcos1 __mpcos1_fma4 +#define __mpranred __mpranred_fma4 +#define __add __add_fma4 +#define __dbl_mp __dbl_mp_fma4 +#define __mul __mul_fma4 +#define __sub __sub_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/sincos32.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/slowexp-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/slowexp-fma4.c new file mode 100644 index 000000000..3bcde8423 --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/slowexp-fma4.c @@ -0,0 +1,9 @@ +#define __slowexp __slowexp_fma4 +#define __add __add_fma4 +#define __dbl_mp __dbl_mp_fma4 +#define __mpexp __mpexp_fma4 +#define __mul __mul_fma4 +#define __sub __sub_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/slowexp.c> diff --git a/libc/sysdeps/x86_64/fpu/multiarch/slowpow-fma4.c b/libc/sysdeps/x86_64/fpu/multiarch/slowpow-fma4.c new file mode 100644 index 000000000..69d69823b --- /dev/null +++ b/libc/sysdeps/x86_64/fpu/multiarch/slowpow-fma4.c @@ -0,0 +1,11 @@ +#define __slowpow __slowpow_fma4 +#define __add __add_fma4 +#define __dbl_mp __dbl_mp_fma4 +#define __mpexp __mpexp_fma4 +#define __mplog __mplog_fma4 +#define __mul __mul_fma4 +#define __sub __sub_fma4 +#define __halfulp __halfulp_fma4 +#define SECTION __attribute__ ((section (".text.fma4"))) + +#include <sysdeps/ieee754/dbl-64/slowpow.c> diff --git a/libc/wcsmbs/wcscmp.c b/libc/wcsmbs/wcscmp.c index ddbd4aa93..98728359b 100644 --- a/libc/wcsmbs/wcscmp.c +++ b/libc/wcsmbs/wcscmp.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1995, 1996, 1997 Free Software Foundation, Inc. +/* Copyright (C) 1995, 1996, 1997, 2011 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@gnu.ai.mit.edu>, 1995. @@ -31,12 +31,12 @@ WCSCMP (s1, s2) const wchar_t *s1; const wchar_t *s2; { - wint_t c1, c2; + wchar_t c1, c2; do { - c1 = (wint_t) *s1++; - c2 = (wint_t) *s2++; + c1 = *s1++; + c2 = *s2++; if (c2 == L'\0') return c1 - c2; } diff --git a/libc/wcsmbs/wcslen.c b/libc/wcsmbs/wcslen.c index 4d7972b7c..991e15172 100644 --- a/libc/wcsmbs/wcslen.c +++ b/libc/wcsmbs/wcslen.c @@ -20,12 +20,12 @@ #include <wchar.h> /* Return length of string S. */ -#ifndef WCSLEN -# define WCSLEN __wcslen +#ifdef WCSLEN +# define __wcslen WCSLEN #endif size_t -WCSLEN (s) +__wcslen (s) const wchar_t *s; { size_t len = 0; diff --git a/libc/wcsmbs/wmemcmp.c b/libc/wcsmbs/wmemcmp.c index 73c6394f8..9bd556e64 100644 --- a/libc/wcsmbs/wmemcmp.c +++ b/libc/wcsmbs/wmemcmp.c @@ -29,25 +29,25 @@ WMEMCMP (s1, s2, n) const wchar_t *s2; size_t n; { - register wint_t c1; - register wint_t c2; + register wchar_t c1; + register wchar_t c2; while (n >= 4) { - c1 = (wint_t) s1[0]; - c2 = (wint_t) s2[0]; + c1 = s1[0]; + c2 = s2[0]; if (c1 - c2 != 0) return c1 > c2 ? 1 : -1; - c1 = (wint_t) s1[1]; - c2 = (wint_t) s2[1]; + c1 = s1[1]; + c2 = s2[1]; if (c1 - c2 != 0) return c1 > c2 ? 1 : -1; - c1 = (wint_t) s1[2]; - c2 = (wint_t) s2[2]; + c1 = s1[2]; + c2 = s2[2]; if (c1 - c2 != 0) return c1 > c2 ? 1 : -1; - c1 = (wint_t) s1[3]; - c2 = (wint_t) s2[3]; + c1 = s1[3]; + c2 = s2[3]; if (c1 - c2 != 0) return c1 > c2 ? 1 : -1; s1 += 4; @@ -57,8 +57,8 @@ WMEMCMP (s1, s2, n) if (n > 0) { - c1 = (wint_t) s1[0]; - c2 = (wint_t) s2[0]; + c1 = s1[0]; + c2 = s2[0]; if (c1 - c2 != 0) return c1 > c2 ? 1 : -1; ++s1; @@ -67,8 +67,8 @@ WMEMCMP (s1, s2, n) } if (n > 0) { - c1 = (wint_t) s1[0]; - c2 = (wint_t) s2[0]; + c1 = s1[0]; + c2 = s2[0]; if (c1 - c2 != 0) return c1 > c2 ? 1 : -1; ++s1; @@ -77,8 +77,8 @@ WMEMCMP (s1, s2, n) } if (n > 0) { - c1 = (wint_t) s1[0]; - c2 = (wint_t) s2[0]; + c1 = s1[0]; + c2 = s2[0]; if (c1 - c2 != 0) return c1 > c2 ? 1 : -1; } |