From d705269e75ff9846000abea981774b381b076482 Mon Sep 17 00:00:00 2001 From: Ulrich Drepper Date: Thu, 19 Jun 1997 19:11:22 +0000 Subject: Update. 1997-06-19 19:38 Ulrich Drepper * features.h: Define __STDC_IEC_559__ and _STDC_IEC_559_COMPLEX__. * elf/dl-minimal.c (__dcgettext): Remove assertion. * inet/rcmd.c: Correct a few typos. Reported by Eric Troan. * manual/Makefile (distribute): Add dir. * manual/dir: New file. * math/Makefile (libm-support): Rename s_rinttol, s_rinttoll, s_roundtol, and s_roundtoll to s_lrint, s_llrint, s_lround, and s_llround respectively. (libm-calls): Add e_gamma_r. * math/libm-test.c (check_int_exc): New function. (signbit_test): Remove test for sign of NaN. (gamma_test): Clear exception after test of existence. Correct tests and and epsilons. (lgamma_test): Likewise. (ilogb_test): Correct all tests. (scalb_test): Rewrite. (rinttol_test): Rename to lrint_test and correct tests. (rinttoll_test): Likewise. (roundtol_test): Likewise. (roundtoll_test): Likewise. (main): Call lrint/lround functions instead of rinttol/roundtol. * math/math.h: Change prototypes for rinttol/roundtol. * math/mathcalls.h: Rearrange prototypes according to ISO C9X draft. * sysdeps/generic/mathbits.h: Define FP_ILOGB0 and FP_ILOGBNAN. * sysdeps/i386/mathbits.h: Likewise. * sysdeps/libm-i387/e_scalb.S: Handle special cases correctly. * sysdeps/libm-i387/e_scalbf.S: Likewise. * sysdeps/libm-i387/e_scalbl.S: Likewise. * sysdeps/libm-i387/s_asinh.S: Handle -inf correctly. * sysdeps/libm-i387/s_asinhf.S: Likewise. * sysdeps/libm-i387/s_asinhl.S: Likewise. * sysdeps/libm-i387/s_ilogb.S: Optimize. * sysdeps/libm-i387/s_ilogbf.S: Likewise. * sysdeps/libm-i387/s_ilogbl.S: Likewise. * sysdeps/libm-i387/s_rinttol.S: Rename to... * sysdeps/libm-i387/s_lrint.S: ...this. * sysdeps/libm-i387/s_rinttoll.S: Rename to... * sysdeps/libm-i387/s_llrint.S: ...this. * sysdeps/libm-i387/s_remquo.S: Correctly set sign of remainder. * sysdeps/libm-i387/s_remquof.S: Likewise. * sysdeps/libm-i387/s_remquol.S: Likewise. * sysdeps/libm-i387/e_gamma_r.c: New file. Implementation of gamma function according to ISO C. * sysdeps/libm-i387/e_gammaf_r.c: New file. * sysdeps/libm-i387/e_gammal_r.c: New file. * sysdeps/libm-i387/e_lgamma_r.c: Don't let optimize compile the generation of exceptions away. * sysdeps/libm-i387/e_lgammaf_r.c: Likewise. * sysdeps/libm-i387/k_standard.c: Correct return value for infinity points of gamma function when not SVID mode. * sysdeps/libm-i387/s_rinttoll.c: Renamed to... * sysdeps/libm-i387/s_llrint.c: ...this. * sysdeps/libm-i387/s_rinttol.c: Renamed to... * sysdeps/libm-i387/s_lrint.c: ...this. * sysdeps/libm-i387/s_roundtoll.c: Renamed to... * sysdeps/libm-i387/s_llround.c: ...this. * sysdeps/libm-i387/s_roundtol.c: Renamed to... * sysdeps/libm-i387/s_lround.c: ..this. * sysdeps/libm-i387/s_scalbn.c: Change second parameter according to ISO C. * sysdeps/libm-i387/s_scalbnf.c: Likewise. * sysdeps/libm-i387/s_scalbnl.c: Likewise. * sysdeps/libm-i387/w_gamma.c: Call __ieee754_gamma_r if library mode is not _SVID_. * sysdeps/libm-i387/w_gammaf.c: Likewise. * sysdeps/libm-i387/w_gammal.c: Likewise. * sysdeps/m68k/fpu/__math.h: Rename __rinttol to __lrint and rinttol to lrint. * sysdeps/m68k/fpu/s_rinttol.c: Renamed to... * sysdeps/m68k/fpu/s_lrint.c: ...this. * sysdeps/m68k/fpu/s_rinttoll.c: Renamed to... * sysdeps/m68k/fpu/s_llrint.c: ...this. * md5-crypt/Makefile: Link md5test program with md5.o. * stdio-common/temptest.c: Don't use __stdio_gen_tempname which is not exported by the libc.so. * stdio-common/vfscanf.c: Correct scanning of strings after last change. * sysdeps/unix/sysv/linux/i386/sysdep.S: Use .comm to define errno. 1997-06-19 07:37 H.J. Lu * time/tzfile.c (__tzfile_read): Store getc () return in int. 1997-06-13 Miles Bader * argp/argp-parse.c (argp_version_parser): Include `(PROGRAM ERROR)' in the no-version error text to indicate that something's fucked. [!_] (N_): New macro. (argp_default_options, argp_version_options): Wrap doc strings in N_(). 1997-06-12 Miles Bader * argp/argp-parse.c (parser_parse_opt): Detect and report unhandled options here. (parser_parse_arg): Handle ARGP_KEY_ARGS here. Adjust NEXT pointer back if we fail to parse anything. (parser_parse_next): Simplify arg code. Leave state NEXT frobbing to parser_parse_arg. 1997-06-11 Miles Bader * argp/argp.h (ARGP_KEY_ARGS, ARGP_KEY_FINI): New macros. * argp/argp-parse.c (parser_finalize): Do another pass over the parsers with ARGP_KEY_FINI. 1997-06-18 Miles Bader * string/Makefile (routines): Add argz-replace. 1997-06-16 00:16 Miles Bader * manual/string.texi (Argz Functions): Document argz_replace. 1997-06-12 Miles Bader * string/argz.h (__argz_replace, argz_replace): New declarations. * string/argz-replace.c: New file. 1997-06-11 Miles Bader * manual/argp.texi (Argp Special Keys): Document ARGP_KEY_ARGS. 1997-06-16 23:08 Andreas Schwab * Makeconfig (libc-map): Remove definition. * Makerules (libc-map): Define it here, using the full name. (load-map-file): Remove case for empty $(..). ($(common-objpfx)libc.so): Revert last change. 1997-06-17 22:18 Mark Kettenis * login/programs/utmpd.c (main): Improve signal handling. * login/programs/request.c (do_setutent, do_updwtmp): Allow arbitrary length filenames. * login/programs/utmpd.h (setutent_request, updwtmp_request): Get rid of fixed length file field. * login/utmp_daemon.c (do_setutent, do_updwtmp): Allow arbitrary length filenames. * login/programs/request.c (do_pututline): Don't fail if connection->position is -1 on entry. 1997-06-15 16:32 Mark Kettenis * login/utmp_file.c (updwtmp_file): Use the same method for appending an entry as in pututline_file. 1997-06-11 18:59 Mark Kettenis * login/utmp_file.c (getutent_r_file): Use read lock instead of write lock. (getutline_r_file, internal_getut_r): Lock utmp file. (updwtmp_file): Use fcntl to lock file instead of flock. 1997-06-18 00:11 Andreas Schwab * stdio-common/vfscanf.c (inchar, ungetc): Don't count EOF as character read in/put back. * stdio-common/tstscanf.c: Add test case for this. 1997-06-17 22:17 Andreas Schwab * libc.map: Add more symbols. 1997-06-18 12:01 Ulrich Drepper * manual/Makefile (install): Make sure `dir' file exists if we use install-info. 1997-06-17 19:32 Andreas Schwab * manual/Makefile (info): Depend on dir-add.info. 1997-06-17 17:19 Andreas Schwab * elf/ldd.bash.in: Fix spacing in message. 1997-06-17 14:28 Andreas Schwab * elf/do-rel.h (elf_dynamic_do_rel): Always use version information if available. 1997-06-17 11:34 Andreas Schwab * sunrpc/Makefile ($(rpcsvc:%.x=$(objpfx)rpcsvc/%.h)): Make command non-empty to force make to recheck modification time. ($(rpcsvc:%.x=$(objpfx)x%.c)): Likewise. 1997-06-17 00:26 Andreas Schwab * sysdeps/stub/e_acoshl.c: Set errno to ENOSYS. * sysdeps/stub/e_acosl.c: Likewise. * sysdeps/stub/e_asinl.c: Likewise. * sysdeps/stub/e_atan2l.c: Likewise. * sysdeps/stub/e_expl.c: Likewise. * sysdeps/stub/e_fmodl.c: Likewise. * sysdeps/stub/e_j0l.c: Likewise. * sysdeps/stub/e_j1l.c: Likewise. * sysdeps/stub/e_jnl.c: Likewise. * sysdeps/stub/e_lgammal_r.c: Likewise. * sysdeps/stub/e_log10l.c: Likewise. * sysdeps/stub/e_logl.c: Likewise. * sysdeps/stub/e_powl.c: Likewise. * sysdeps/stub/e_rem_pio2l.c: Likewise. * sysdeps/stub/e_sqrtl.c: Likewise. * sysdeps/stub/k_cosl.c: Likewise. * sysdeps/stub/k_rem_pio2l.c: Likewise. * sysdeps/stub/k_sinl.c: Likewise. * sysdeps/stub/k_tanl.c: Likewise. * sysdeps/stub/s_atanl.c: Likewise. * sysdeps/stub/s_erfl.c: Likewise. * sysdeps/stub/s_exp2.c: Likewise. * sysdeps/stub/s_exp2f.c: Likewise. * sysdeps/stub/s_exp2l.c: Likewise. * sysdeps/stub/s_expm1l.c: Likewise. * sysdeps/stub/s_log1pl.c: Likewise. * sysdeps/stub/s_log2l.c: Likewise. 1997-06-18 11:46 Ulrich Drepper * sysdeps/libm-ieee754/w_gamma.c: If _LIB_VERSION is _SVID_ compute result as before last change. * sysdeps/libm-ieee754/w_gammaf.c: Likewise. * sysdeps/libm-ieee754/w_gammal.c: Likewise. 1997-06-16 23:37 Andreas Schwab * sysdeps/libm-ieee754/s_remquo.c: Fix off-by-one when computing quotient. * sysdeps/libm-ieee754/s_remquof.c: Likewise. * sysdeps/libm-ieee754/s_remquol.c: Likewise. * sysdeps/m68k/fpu/s_remquo.c: Remove FIXME and special case for quotient. * sysdeps/libm-ieee754/w_gamma.c: Add missing call to exp function. Don't use global signgam. * sysdeps/libm-ieee754/w_gammaf.c: Likewise. * sysdeps/libm-ieee754/w_gammal.c: Likewise. * math/Makefile (libm-calls): Remove w_gamma_r. * sysdeps/libm-ieee754/w_gamma_r.c: Remove file. * sysdeps/libm-ieee754/w_gammaf_r.c: Remove file. * sysdeps/libm-ieee754/w_gammal_r.c: Remove file. * math/libm-test.c (atanh_test): Declare x only if needed. (signbit_test): Fix typo. (gamma_test): Check whether function is implemented. Add epsilons. (lgamma_test): Likewise. (fmod_test): Add epsilons. (exp2_test): Use right function for existence test. 1997-06-07 09:20 Andreas Schwab * sysdeps/powerpc/Dist: Add fenv_const.c. * sysdeps/unix/sysv/linux/Dist: Add net/if_slip.h. * sysdeps/unix/sysv/linux/powerpc/Dist: Add init-first.h and syscall.h. * sysdeps/unix/sysv/linux/sparc/Dist: Add init-first.h. * string/Makefile (distribute): Add tst-svc.expect. * nis/Makefile (distribute): Add nis_intern.h and Banner. * elf/Makefile (distribute): Add dl-hash.h. * Rules (subdir_echo-distinfo): Add headers from $(distribute). * login/Makefile (others): Add utmpdump. * sysdeps/unix/sysv/linux/Makefile (sysdep_headers): Add net/if_slip.h. * manual/Makefile (dir-add.texi): Also look in indirectly included files. 1997-06-16 23:15 Andreas Schwab * sysdeps/unix/sysv/linux/configure.in: Remove commands to generate stdio_lim.h. * sysdeps/unix/sysv/linux/mk-stdiolim.c: Remove. * sysdeps/unix/sysv/linux/Makefile: Add rules to generate stdio_lim.h here. (common-generated): Add generated files. (inhibit-stdio_lim): Define. * sysdeps/posix/Makefile [$(inhibit-stdio_lim)=yes]: Disable rules to generate stdio_lim.h. * sysdeps/unix/sysv/linux/Makefile: Suppress inclusion of dependecy files if no_deps is set. ($(objpfx)syscall-%.d): Add header file as target to dependency generation. 1997-06-14 19:19 Andreas Schwab * sysdeps/unix/sysv/linux/arm/Dist: New file. 1997-06-14 17:59 Andreas Schwab * extra-lib.mk (others): Depend on versioned shared library, not the unversioned one. * Makerules (build-shlib): Don't make the version link here. ($(common-objpfx)libc.so$(libc.so-version)): New rule for libc version link. 1997-06-16 03:07 Ulrich Drepper * sysdeps/libm-ieee754/k_standard.c: Undo change of Tue Aug 6 01:13:56 1996. * argp/argp-help.c (argp_args_usage): Supply correct argp to filter_doc. * argp/argp-help.c (hol_add_cluster): Initialize CL->depth. * argp/argp-help.c (_help): Supply STATE to argp_args_usage. * argp/argp.h (ARGP_KEY_HELP_ARGS_DOC): New macro. * argp/argp-fmtstream.c: Add casts to prevent warnings. * argp/argp.h (OPTION_NO_USAGE): New macro. * argp/argp-help.c (usage_long_opt, usage_argful_short_opt, * argp/argp-fmtstream.c (__argp_fmtstream_update): Account for case * argp/argp-help.c : New include. * argp/argp.h (argp_state_help, __argp_state_help, argp_usage, * argp/argp.h (argp_program_bug_address): Make const. * argp/argp-parse.c (argp_default_parser): Set STATE->name for OPT_PROGNAME. * argp/argp-help.c (__argp_error, __argp_failure, __argp_state_help): * argp/argp-parse.c (parser_init): Set PARSER->state.flags. --- sysdeps/libm-ieee754/e_gamma_r.c | 47 +++++++ sysdeps/libm-ieee754/e_gammaf_r.c | 46 +++++++ sysdeps/libm-ieee754/e_gammal_r.c | 46 +++++++ sysdeps/libm-ieee754/e_lgamma_r.c | 40 +++--- sysdeps/libm-ieee754/e_lgammaf_r.c | 22 ++-- sysdeps/libm-ieee754/s_llrint.c | 242 +++++++++++++++++++++++++++++++++++++ sysdeps/libm-ieee754/s_llround.c | 179 +++++++++++++++++++++++++++ sysdeps/libm-ieee754/s_lrint.c | 242 +++++++++++++++++++++++++++++++++++++ sysdeps/libm-ieee754/s_lround.c | 177 +++++++++++++++++++++++++++ sysdeps/libm-ieee754/s_remquo.c | 2 + sysdeps/libm-ieee754/s_remquof.c | 2 + sysdeps/libm-ieee754/s_remquol.c | 2 + sysdeps/libm-ieee754/s_rinttol.c | 242 ------------------------------------- sysdeps/libm-ieee754/s_rinttoll.c | 242 ------------------------------------- sysdeps/libm-ieee754/s_roundtol.c | 177 --------------------------- sysdeps/libm-ieee754/s_roundtoll.c | 179 --------------------------- sysdeps/libm-ieee754/s_scalbn.c | 13 +- sysdeps/libm-ieee754/s_scalbnf.c | 20 +-- sysdeps/libm-ieee754/s_scalbnl.c | 34 +++--- sysdeps/libm-ieee754/w_gamma.c | 22 ++-- sysdeps/libm-ieee754/w_gamma_r.c | 51 -------- sysdeps/libm-ieee754/w_gammaf.c | 21 ++-- sysdeps/libm-ieee754/w_gammaf_r.c | 52 -------- sysdeps/libm-ieee754/w_gammal.c | 21 ++-- sysdeps/libm-ieee754/w_gammal_r.c | 52 -------- 25 files changed, 1089 insertions(+), 1084 deletions(-) create mode 100644 sysdeps/libm-ieee754/e_gamma_r.c create mode 100644 sysdeps/libm-ieee754/e_gammaf_r.c create mode 100644 sysdeps/libm-ieee754/e_gammal_r.c create mode 100644 sysdeps/libm-ieee754/s_llrint.c create mode 100644 sysdeps/libm-ieee754/s_llround.c create mode 100644 sysdeps/libm-ieee754/s_lrint.c create mode 100644 sysdeps/libm-ieee754/s_lround.c delete mode 100644 sysdeps/libm-ieee754/s_rinttol.c delete mode 100644 sysdeps/libm-ieee754/s_rinttoll.c delete mode 100644 sysdeps/libm-ieee754/s_roundtol.c delete mode 100644 sysdeps/libm-ieee754/s_roundtoll.c delete mode 100644 sysdeps/libm-ieee754/w_gamma_r.c delete mode 100644 sysdeps/libm-ieee754/w_gammaf_r.c delete mode 100644 sysdeps/libm-ieee754/w_gammal_r.c (limited to 'sysdeps/libm-ieee754') diff --git a/sysdeps/libm-ieee754/e_gamma_r.c b/sysdeps/libm-ieee754/e_gamma_r.c new file mode 100644 index 0000000000..eb80809843 --- /dev/null +++ b/sysdeps/libm-ieee754/e_gamma_r.c @@ -0,0 +1,47 @@ +/* Implementation of gamma function according to ISO C. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper , 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include +#include + + +double +__ieee754_gamma_r (double x, int *signgamp) +{ + /* We don't have a real gamma implementation now. We'll use lgamma + and the exp function. But due to the required boundary + conditions we must check some values separately. */ + int32_t hx; + u_int32_t lx; + + EXTRACT_WORDS (hx, lx, x); + + if ((hx & 0x7fffffff | lx) == 0) + /* Return value for x == 0 is NaN with invalid exception. */ + return x / x; + if (hx < 0 && (u_int32_t) hx < 0xfff00000 && __rint (x) == x) + { + /* Return value for integer x < 0 is NaN with invalid exception. */ + return (x - x) / (x - x); + } + + /* XXX FIXME. */ + return __ieee754_exp (__ieee754_lgamma_r (x, signgamp)); +} diff --git a/sysdeps/libm-ieee754/e_gammaf_r.c b/sysdeps/libm-ieee754/e_gammaf_r.c new file mode 100644 index 0000000000..bf79ffeb57 --- /dev/null +++ b/sysdeps/libm-ieee754/e_gammaf_r.c @@ -0,0 +1,46 @@ +/* Implementation of gamma function according to ISO C. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper , 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include +#include + + +float +__ieee754_gammaf_r (float x, int *signgamp) +{ + /* We don't have a real gamma implementation now. We'll use lgamma + and the exp function. But due to the required boundary + conditions we must check some values separately. */ + int32_t hx; + + GET_FLOAT_WORD (hx, x); + + if ((hx & 0x7fffffff) == 0) + /* Return value for x == 0 is NaN with invalid exception. */ + return x / x; + if (hx < 0 && (u_int32_t) hx < 0xff800000 && __rintf (x) == x) + { + /* Return value for integer x < 0 is NaN with invalid exception. */ + return (x - x) / (x - x); + } + + /* XXX FIXME. */ + return __ieee754_expf (__ieee754_lgammaf_r (x, signgamp)); +} diff --git a/sysdeps/libm-ieee754/e_gammal_r.c b/sysdeps/libm-ieee754/e_gammal_r.c new file mode 100644 index 0000000000..4fd49886d1 --- /dev/null +++ b/sysdeps/libm-ieee754/e_gammal_r.c @@ -0,0 +1,46 @@ +/* Implementation of gamma function according to ISO C. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper , 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include +#include + + +long double +__ieee754_gammal_r (long double x, int *signgamp) +{ + /* We don't have a real gamma implementation now. We'll use lgamma + and the exp function. But due to the required boundary + conditions we must check some values separately. */ + u_int32_t es, hx, lx; + + GET_LDOUBLE_WORDS (es, hx, lx, x); + + if (((es & 0x7fff) | hx | lx) == 0) + /* Return value for x == 0 is NaN with invalid exception. */ + return x / x; + if ((hx & 0x8000) != 0 && (hx & 0x7fff) != 0x7fff && __rintl (x) == x) + { + /* Return value for integer x < 0 is NaN with invalid exception. */ + return (x - x) / (x - x); + } + + /* XXX FIXME. */ + return __ieee754_expl (__ieee754_lgammal_r (x, signgamp)); +} diff --git a/sysdeps/libm-ieee754/e_lgamma_r.c b/sysdeps/libm-ieee754/e_lgamma_r.c index 1be4ddad8f..92e9556568 100644 --- a/sysdeps/libm-ieee754/e_lgamma_r.c +++ b/sysdeps/libm-ieee754/e_lgamma_r.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -15,12 +15,12 @@ static char rcsid[] = "$NetBSD: e_lgamma_r.c,v 1.7 1995/05/10 20:45:42 jtc Exp $ #endif /* __ieee754_lgamma_r(x, signgamp) - * Reentrant version of the logarithm of the Gamma function - * with user provide pointer for the sign of Gamma(x). + * Reentrant version of the logarithm of the Gamma function + * with user provide pointer for the sign of Gamma(x). * * Method: * 1. Argument Reduction for 0 < x <= 8 - * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may + * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may * reduce x to a number in [1.5,2.5] by * lgamma(1+s) = log(s) + lgamma(s) * for example, @@ -58,36 +58,36 @@ static char rcsid[] = "$NetBSD: e_lgamma_r.c,v 1.7 1995/05/10 20:45:42 jtc Exp $ * by * 3 5 11 * w = w0 + w1*z + w2*z + w3*z + ... + w6*z - * where + * where * |w - f(z)| < 2**-58.74 - * + * * 4. For negative x, since (G is gamma function) * -x*G(-x)*G(x) = pi/sin(pi*x), * we have * G(x) = pi/(sin(pi*x)*(-x)*G(-x)) * since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0 - * Hence, for x<0, signgam = sign(sin(pi*x)) and + * Hence, for x<0, signgam = sign(sin(pi*x)) and * lgamma(x) = log(|Gamma(x)|) * = log(pi/(|x*sin(pi*x)|)) - lgamma(-x); - * Note: one should avoid compute pi*(-x) directly in the + * Note: one should avoid compute pi*(-x) directly in the * computation of sin(pi*(-x)). - * + * * 5. Special Cases * lgamma(2+s) ~ s*(1-Euler) for tiny s * lgamma(1)=lgamma(2)=0 * lgamma(x) ~ -log(x) for tiny x * lgamma(0) = lgamma(inf) = inf * lgamma(-integer) = +-inf - * + * */ #include "math.h" #include "math_private.h" #ifdef __STDC__ -static const double +static const double #else -static double +static double #endif two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ @@ -200,9 +200,9 @@ static double zero= 0.00000000000000000000e+00; } switch (n) { case 0: y = __kernel_sin(pi*y,zero,0); break; - case 1: + case 1: case 2: y = __kernel_cos(pi*(0.5-y),zero); break; - case 3: + case 3: case 4: y = __kernel_sin(pi*(one-y),zero,0); break; case 5: case 6: y = -__kernel_cos(pi*(y-1.5),zero); break; @@ -226,9 +226,11 @@ static double zero= 0.00000000000000000000e+00; /* purge off +-inf, NaN, +-0, and negative arguments */ *signgamp = 1; + if ((unsigned int) hx==0xfff00000&&lx==0) + return x-x; ix = hx&0x7fffffff; if(ix>=0x7ff00000) return x*x; - if((ix|lx)==0) return one/zero; + if((ix|lx)==0) return one/fabs(x); if(ix<0x3b900000) { /* |x|<2**-70, return -log(|x|) */ if(hx<0) { *signgamp = -1; @@ -237,9 +239,9 @@ static double zero= 0.00000000000000000000e+00; } if(hx<0) { if(ix>=0x43300000) /* |x|>=2**52, must be -integer */ - return one/zero; + return x/zero; t = sin_pi(x); - if(t==zero) return one/zero; /* -integer */ + if(t==zero) return one/fabsf(t); /* -integer */ nadj = __ieee754_log(pi/fabs(t*x)); if(t=0x7f800000) return x*x; - if(ix==0) return one/zero; + if(ix==0) return one/fabsf(x); if(ix<0x1c800000) { /* |x|<2**-70, return -log(|x|) */ if(hx<0) { *signgamp = -1; @@ -173,9 +175,9 @@ static float zero= 0.0000000000e+00; } if(hx<0) { if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */ - return one/zero; + return x/zero; t = sin_pif(x); - if(t==zero) return one/zero; /* -integer */ + if(t==zero) return one/fabsf(t); /* -integer */ nadj = __ieee754_logf(pi/fabsf(t*x)); if(t, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include + +#include "math_private.h" + +#ifdef NO_LONG_DOUBLE +/* The `long double' is in fact the IEEE `double' type. */ +static long double two52[2] = +{ + 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ + -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ +}; + + +long long int +__llrint (long double x) +{ + int32_t j0,sx; + u_int32_t i0, i1, i; + long double t, w; + long long int result; + + EXTRACT_WORDS (i0, i1, x); + + sx = i0 >> 31; + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + + if (j0 < 20) + { + if (j0 < 0) + { + if (((i0 & 0x7fffffff) | i1) == 0) + /* The number is 0. */ + result = 0; + else + { + i1 |= i0; + i0 &= 0xfffe0000; + i0 |= ((i1 | -i1) >> 12) & 0x80000; + SET_HIGH_WORD (x, i0); + w = two52[sx] + x; + t = w - two52[sx]; + GET_HIGH_WORD (i0, t); + if ((i0 & 0x7fffffff) >= 0x3fff0000) + result = sx ? -1 : 1; + else + result = 0; + } + } + else + { + u_int32_t i = 0x000fffff >> j0; + if (((i0 & i) | i1) == 0) + { + /* X is not integral. */ + i >>= 1; + if (((i0 & i) | i1) != 0) + { + if (j0 == 19) + i1 = 0x40000000; + else + i0 = (i0 & (~i)) | (0x20000 >> j0); + + INSERT_WORDS (x, i0, i1); + w = two52[sx] + x; + x = w - two52[sx]; + EXTRACT_WORDS (i0, i1, x); + + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + } + } + + result = ((i0 >> (20 - j0)) & 0xfffff) | (0x00100000 >> (20 - j0)); + if (sx) + result = -result; + } + } + else if ((unsigned int) j0 < sizeof (long long int) * 8 && j0 < 53) + { + i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); + if ((i1 & i) != 0) + { + /* x is not integral. */ + i >>= 1; + if ((i1 & i) != 0) + i1 = (i1 & (~i)) | (0x40000000 >> (j0 - 20)); + } + + INSERT_WORDS (x, i0, i1); + w = two52[sx] + x; + x = w - two52[sx]; + EXTRACT_WORDS (i0, i1, x); + + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + + result = i0 | 0x00100000; + if (j0 > 20) + { + result <<= j0 - 20; + result |= i1 >> (52 - j0); + } + if (sx) + result = -result; + } + else + /* Too large. The number is either +-inf or NaN or it is too + large to be effected by rounding. The standard leaves it + undefined what to return when the number is too large to fit in + a `long int'. */ + result = (long long int) x; + + return result; +} + +#else +static long double two63[2] = +{ + 9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */ + -9.223372036854775808000000e+18 /* 0xC03E, 0x00000000, 0x00000000 */ +}; + + +long long int +__llrint (long double x) +{ + int32_t se,j0,sx; + u_int32_t i0, i1, i; + long long int result; + long double w, t; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + + sx = (se >> 15) & 1; + j0 = (se & 0x7fff) - 0x3fff; + + if (j0 < 31) + { + if (j0 < 0) + { + if (((se & 0x7fff) | i0 | i1) == 0) + /* The number is 0. */ + result = 0; + else + { + i1 |= i0; + i0 &= 0xe0000000; + i0 |= (i1 | -i1) & 0x80000000; + SET_LDOUBLE_MSW (x, i0); + w = two63[sx] + x; + t = w - two63[sx]; + GET_LDOUBLE_EXP (i0, t); + if ((i0 & 0x7fff) >= 0x3fff) + result = sx ? -1 : 1; + else + result = 0; + } + } + else + { + u_int32_t i = 0x7fffffff >> j0; + if (((i0 & i) | i1) == 0) + { + /* X is not integral. */ + i >>= 1; + if (((i0 & i) | i1) != 0) + { + if (j0 == 31) + i1 = 0x40000000; + else + i0 = (i0 & (~i)) | (0x20000000 >> j0); + + SET_LDOUBLE_WORDS (x, se, i0, i1); + w = two63[sx] + x; + x = w - two63[sx]; + GET_LDOUBLE_WORDS (se, i0, i1, x); + + sx = (se >> 15) & 1; + j0 = (se & 0x7fff) - 0x3fff; + } + } + + + result = i0 >> (31 - j0); + } + } + else if ((unsigned int) j0 < sizeof (long long int) * 8 && j0 < 64) + { + i = ((u_int32_t) (0xffffffff)) >> (j0 - 31); + if ((i1 & i) != 0) + { + /* x is not integral. */ + i >>= 1; + if ((i1 & i) != 0) + i1 = (i1 & (~i)) | (0x40000000 >> (j0 - 31)); + } + + SET_LDOUBLE_WORDS (x, se, i0, i1); + w = two63[sx] + x; + x = w - two63[sx]; + GET_LDOUBLE_WORDS (se, i0, i1, x); + + j0 = (se & 0x7fff) - 0x3fff; + + result = i0; + if (j0 > 31) + { + result <<= j0 - 31; + result |= i1 >> (63 - j0); + } + } + else + /* Too large. The number is either +-inf or NaN or it is too + large to be effected by rounding. The standard leaves it + undefined what to return when the number is too large to fit in + a `long int'. */ + result = (long long int) x; + + return result; +} +#endif + +weak_alias (__llrint, llrint) diff --git a/sysdeps/libm-ieee754/s_llround.c b/sysdeps/libm-ieee754/s_llround.c new file mode 100644 index 0000000000..aee0e31fc0 --- /dev/null +++ b/sysdeps/libm-ieee754/s_llround.c @@ -0,0 +1,179 @@ +/* Round long double value to long long int. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper , 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include + +#include "math_private.h" + + +#ifdef NO_LONG_DOUBLE +/* The `long double' is in fact the IEEE `double' type. */ + +long long int +__llround (long double x) +{ + int32_t j0; + u_int32_t i1, i0; + long long int result; + + EXTRACT_WORDS (i0, i1, x); + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + if (j0 < 20) + { + if (j0 < 0) + result = j0 < -1 ? 0 : ((i0 & 0x80000000) ? -1 : 1); + else + { + u_int32_t i = 0xfffff >> j0; + if (((i0 & i) | i1) == 0) + result = (long long int) ((i0 & 0xfffff) | 0x100000) >> j0; + else + { + /* X is not integral. */ + u_int32_t j = i0 + (0x80000 >> j0); + if (j < i0) + result = (long long int) 0x80000 >> (20 - j0); + else + result = (j | 0x100000) >> (20 - j0); + } + } + } + else if (j0 >= 8 * sizeof (long long int) || j0 > 51) + { + /* The number is too large. It is left implementation defined + what happens. */ + result = (long long int) x; + } + else + { + u_int32_t i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); + if ((i1 & i) != 0) + { + /* x is not integral. */ + u_int32_t j = i1 + (0x80000000 >> (j0 - 20)); + if (j < i1) + { + j = i0 + 1; + if ((j & 0xfffff) == 0) + { + if (sizeof (long long int) <= 4) + /* Overflow. */ + result = (long long int) x; + else + result = 1ll << (j0 + 1); + } + else + result = ((long long int) ((i0 & 0xfffff) | 0x100000) + << (j0 - 31)); + } + else + { + result = ((long long int) ((i0 & 0xfffff) | 0x100000) + << (j0 - 31)); + if (sizeof (long long int) > 4 && j0 > 31) + result |= j >> (63 - j0); + } + } + else + { + result = (long long int) ((i0 & 0xfffff) | 0x100000) << (j0 - 31); + if (sizeof (long long int) > 4 && j0 > 31) + result |= i1 >> (63 - j0); + } + } + + return i0 & 0x80000000 ? -result : result; +} +#else +long long int +__llround (long double x) +{ + int32_t j0; + u_int32_t se, i1, i0; + long long int result; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + j0 = (se & 0x7fff) - 0x3fff; + if (j0 < 31) + { + if (j0 < 0) + result = j0 < -1 ? 0 : 1; + else + { + u_int32_t i = 0x7fffffff >> j0; + if (((i0 & i) | i1) == 0) + result = (long long int) i0 >> j0; + else + { + /* X is not integral. */ + u_int32_t j = i0 + (0x40000000 >> j0); + if (j < i0) + result = 0x80000000l >> (30 - j0); + else + result = j >> (31 - j0); + } + } + } + else if ((unsigned int) j0 >= 8 * sizeof (long long int) || j0 > 62) + { + /* The number is too large. It is left implementation defined + what happens. */ + result = (long long int) x; + } + else + { + u_int32_t i = ((u_int32_t) (0xffffffff)) >> (j0 - 31); + if ((i1 & i) != 0) + { + /* x is not integral. */ + u_int32_t j = i1 + (0x80000000 >> (j0 - 31)); + if (j < i1) + { + j = i0 + 1; + if (j == 0) + { + if (sizeof (long long int) <= 4) + /* Overflow. */ + result = (long long int) x; + else + result = 1ll << (j0 + 1); + } + else + result = (long long int) i0 << (j0 - 31); + } + else + { + result = (long long int) i0 << (j0 - 31); + if (sizeof (long long int) > 4 && j0 > 31) + result |= j >> (63 - j0); + } + } + else + { + result = (long long int) i0 << (j0 - 31); + if (sizeof (long long int) > 4 && j0 > 31) + result |= i1 >> (63 - j0); + } + } + + return se & 0x8000 ? -result : result; +} +#endif +weak_alias (__llround, llround) diff --git a/sysdeps/libm-ieee754/s_lrint.c b/sysdeps/libm-ieee754/s_lrint.c new file mode 100644 index 0000000000..6779f974bd --- /dev/null +++ b/sysdeps/libm-ieee754/s_lrint.c @@ -0,0 +1,242 @@ +/* Round argument to nearest integral value according to current rounding + direction. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper , 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include + +#include "math_private.h" + +#ifdef NO_LONG_DOUBLE +/* The `long double' is in fact the IEEE `double' type. */ +static long double two52[2] = +{ + 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ + -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ +}; + + +long int +__lrint (long double x) +{ + int32_t j0,sx; + u_int32_t i0,i1,i; + long double t, w; + long int result; + + EXTRACT_WORDS (i0, i1, x); + + sx = i0 >> 31; + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + + if (j0 < 20) + { + if (j0 < 0) + { + if (((i0 & 0x7fffffff) | i1) == 0) + /* The number is 0. */ + result = 0; + else + { + i1 |= i0; + i0 &= 0xfffe0000; + i0 |= ((i1 | -i1) >> 12) & 0x80000; + SET_HIGH_WORD (x, i0); + w = two52[sx] + x; + t = w - two52[sx]; + GET_HIGH_WORD (i0, t); + if ((i0 & 0x7fffffff) >= 0x3fff0000) + result = sx ? -1 : 1; + else + result = 0; + } + } + else + { + u_int32_t i = 0x000fffff >> j0; + if (((i0 & i) | i1) == 0) + { + /* X is not integral. */ + i >>= 1; + if (((i0 & i) | i1) != 0) + { + if (j0 == 19) + i1 = 0x40000000; + else + i0 = (i0 & (~i)) | (0x20000 >> j0); + + INSERT_WORDS (x, i0, i1); + w = two52[sx] + x; + x = w - two52[sx]; + EXTRACT_WORDS (i0, i1, x); + + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + } + } + + result = ((i0 >> (20 - j0)) & 0xfffff) | (0x00100000 >> (20 - j0)); + if (sx) + result = -result; + } + } + else if ((unsigned int) j0 < sizeof (long int) * 8 && j0 < 53) + { + i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); + if ((i1 & i) != 0) + { + /* x is not integral. */ + i >>= 1; + if ((i1 & i) != 0) + i1 = (i1 & (~i)) | (0x40000000 >> (j0 - 20)); + } + + INSERT_WORDS (x, i0, i1); + w = two52[sx] + x; + x = w - two52[sx]; + EXTRACT_WORDS (i0, i1, x); + + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + + result = i0 | 0x00100000; + if (j0 > 20) + { + result <<= j0 - 20; + result |= i1 >> (52 - j0); + } + if (sx) + result = -result; + } + else + /* Too large. The number is either +-inf or NaN or it is too + large to be effected by rounding. The standard leaves it + undefined what to return when the number is too large to fit in + a `long int'. */ + result = (long int) x; + + return result; +} + +#else +static long double two63[2] = +{ + 9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */ + -9.223372036854775808000000e+18 /* 0xC03E, 0x00000000, 0x00000000 */ +}; + + +long int +__lrint (long double x) +{ + int32_t se,j0,sx; + u_int32_t i0,i1,i; + long int result; + long double w, t; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + + sx = (se >> 15) & 1; + j0 = (se & 0x7fff) - 0x3fff; + + if (j0 < 31) + { + if (j0 < 0) + { + if (((se & 0x7fff) | i0 | i1) == 0) + /* The number is 0. */ + result = 0; + else + { + i1 |= i0; + i0 &= 0xe0000000; + i0 |= (i1 | -i1) & 0x80000000; + SET_LDOUBLE_MSW (x, i0); + w = two63[sx] + x; + t = w - two63[sx]; + GET_LDOUBLE_EXP (i0, t); + if ((i0 & 0x7fff) >= 0x3fff) + result = sx ? -1 : 1; + else + result = 0; + } + } + else + { + u_int32_t i = 0x7fffffff >> j0; + if (((i0 & i) | i1) == 0) + { + /* X is not integral. */ + i >>= 1; + if (((i0 & i) | i1) != 0) + { + if (j0 == 31) + i1 = 0x40000000; + else + i0 = (i0 & (~i)) | (0x20000000 >> j0); + + SET_LDOUBLE_WORDS (x, se, i0, i1); + w = two63[sx] + x; + x = w - two63[sx]; + GET_LDOUBLE_WORDS (se, i0, i1, x); + + sx = (se >> 15) & 1; + j0 = (se & 0x7fff) - 0x3fff; + } + } + + + result = i0 >> (31 - j0); + } + } + else if ((unsigned int) j0 < sizeof (long int) * 8 && j0 < 64) + { + i = ((u_int32_t) (0xffffffff)) >> (j0 - 31); + if ((i1 & i) != 0) + { + /* x is not integral. */ + i >>= 1; + if ((i1 & i) != 0) + i1 = (i1 & (~i)) | (0x40000000 >> (j0 - 31)); + } + + SET_LDOUBLE_WORDS (x, se, i0, i1); + w = two63[sx] + x; + x = w - two63[sx]; + GET_LDOUBLE_WORDS (se, i0, i1, x); + + j0 = (se & 0x7fff) - 0x3fff; + + result = i0; + if (j0 > 31) + { + result <<= j0 - 31; + result |= i1 >> (63 - j0); + } + } + else + /* Too large. The number is either +-inf or NaN or it is too + large to be effected by rounding. The standard leaves it + undefined what to return when the number is too large to fit in + a `long int'. */ + result = (long int) x; + + return result; +} +#endif + +weak_alias (__lrint, lrint) diff --git a/sysdeps/libm-ieee754/s_lround.c b/sysdeps/libm-ieee754/s_lround.c new file mode 100644 index 0000000000..0f91280a79 --- /dev/null +++ b/sysdeps/libm-ieee754/s_lround.c @@ -0,0 +1,177 @@ +/* Round long double value to long int. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper , 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include + +#include "math_private.h" + + +#ifdef NO_LONG_DOUBLE +/* The `long double' is in fact the IEEE `double' type. */ + +long int +__lround (long double x) +{ + int32_t j0; + u_int32_t i1, i0; + long int result; + + EXTRACT_WORDS (i0, i1, x); + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + if (j0 < 20) + { + if (j0 < 0) + result = j0 < -1 ? 0 : ((i0 & 0x80000000) ? -1 : 1); + else + { + u_int32_t i = 0xfffff >> j0; + if (((i0 & i) | i1) == 0) + result = (long int) ((i0 & 0xfffff) | 0x100000) >> j0; + else + { + /* X is not integral. */ + u_int32_t j = i0 + (0x80000 >> j0); + if (j < i0) + result = (long int) 0x80000 >> (20 - j0); + else + result = (j | 0x100000) >> (20 - j0); + } + } + } + else if (j0 >= 8 * sizeof (long int) || j0 > 51) + { + /* The number is too large. It is left implementation defined + what happens. */ + result = (long int) x; + } + else + { + u_int32_t i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); + if ((i1 & i) != 0) + { + /* x is not integral. */ + u_int32_t j = i1 + (0x80000000 >> (j0 - 20)); + if (j < i1) + { + j = i0 + 1; + if ((j & 0xfffff) == 0) + { + if (sizeof (long int) <= 4) + /* Overflow. */ + result = (long int) x; + else + result = 1l << (j0 + 1); + } + else + result = (long int) ((i0 & 0xfffff) | 0x100000) << (j0 - 31); + } + else + { + result = (long int) ((i0 & 0xfffff) | 0x100000) << (j0 - 31); + if (sizeof (long int) > 4 && j0 > 31) + result |= j >> (63 - j0); + } + } + else + { + result = (long int) ((i0 & 0xfffff) | 0x100000) << (j0 - 31); + if (sizeof (long int) > 4 && j0 > 31) + result |= i1 >> (63 - j0); + } + } + + return i0 & 0x80000000 ? -result : result; +} +#else +long int +__lround (long double x) +{ + int32_t j0; + u_int32_t se, i1, i0; + long int result; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + j0 = (se & 0x7fff) - 0x3fff; + if (j0 < 31) + { + if (j0 < 0) + result = j0 < -1 ? 0 : 1; + else + { + u_int32_t i = 0x7fffffff >> j0; + if (((i0 & i) | i1) == 0) + result = (long int) i0 >> j0; + else + { + /* X is not integral. */ + u_int32_t j = i0 + (0x40000000 >> j0); + if (j < i0) + result = 0x80000000l >> (30 - j0); + else + result = j >> (31 - j0); + } + } + } + else if ((unsigned int) j0 >= 8 * sizeof (long int) || j0 > 62) + { + /* The number is too large. It is left implementation defined + what happens. */ + result = (long int) x; + } + else + { + u_int32_t i = ((u_int32_t) (0xffffffff)) >> (j0 - 31); + if ((i1 & i) != 0) + { + /* x is not integral. */ + u_int32_t j = i1 + (0x80000000 >> (j0 - 31)); + if (j < i1) + { + j = i0 + 1; + if (j == 0) + { + if (sizeof (long int) <= 4) + /* Overflow. */ + result = (long int) x; + else + result = 1l << (j0 + 1); + } + else + result = (long int) i0 << (j0 - 31); + } + else + { + result = (long int) i0 << (j0 - 31); + if (sizeof (long int) > 4 && j0 > 31) + result |= j >> (63 - j0); + } + } + else + { + result = (long int) i0 << (j0 - 31); + if (sizeof (long int) > 4 && j0 > 31) + result |= i1 >> (63 - j0); + } + } + + return se & 0x8000 ? -result : result; +} +#endif +weak_alias (__lround, lround) diff --git a/sysdeps/libm-ieee754/s_remquo.c b/sysdeps/libm-ieee754/s_remquo.c index f0a32130a2..39f6ceab2e 100644 --- a/sysdeps/libm-ieee754/s_remquo.c +++ b/sysdeps/libm-ieee754/s_remquo.c @@ -82,6 +82,7 @@ __remquo (double x, double y, int *quo) if (x + x > y) { x -= y; + ++cquo; if (x + x >= y) { x -= y; @@ -95,6 +96,7 @@ __remquo (double x, double y, int *quo) if (x > y_half) { x -= y; + ++cquo; if (x >= y_half) { x -= y; diff --git a/sysdeps/libm-ieee754/s_remquof.c b/sysdeps/libm-ieee754/s_remquof.c index 6fa02e47b3..b3870f4b02 100644 --- a/sysdeps/libm-ieee754/s_remquof.c +++ b/sysdeps/libm-ieee754/s_remquof.c @@ -81,6 +81,7 @@ __remquof (float x, float y, int *quo) if (x + x > y) { x -= y; + ++cquo; if (x + x >= y) { x -= y; @@ -94,6 +95,7 @@ __remquof (float x, float y, int *quo) if (x > y_half) { x -= y; + ++cquo; if (x >= y_half) { x -= y; diff --git a/sysdeps/libm-ieee754/s_remquol.c b/sysdeps/libm-ieee754/s_remquol.c index 9ef424901b..b7835e6dc7 100644 --- a/sysdeps/libm-ieee754/s_remquol.c +++ b/sysdeps/libm-ieee754/s_remquol.c @@ -82,6 +82,7 @@ __remquol (long double x, long double p, int *quo) if (x + x > p) { x -= p; + ++cquo; if (x + x >= p) { x -= p; @@ -95,6 +96,7 @@ __remquol (long double x, long double p, int *quo) if (x > p_half) { x -= p; + ++cquo; if (x >= p_half) { x -= p; diff --git a/sysdeps/libm-ieee754/s_rinttol.c b/sysdeps/libm-ieee754/s_rinttol.c deleted file mode 100644 index 257bdc744c..0000000000 --- a/sysdeps/libm-ieee754/s_rinttol.c +++ /dev/null @@ -1,242 +0,0 @@ -/* Round argument to nearest integral value according to current rounding - direction. - Copyright (C) 1997 Free Software Foundation, Inc. - This file is part of the GNU C Library. - Contributed by Ulrich Drepper , 1997. - - The GNU C Library is free software; you can redistribute it and/or - modify it under the terms of the GNU Library General Public License as - published by the Free Software Foundation; either version 2 of the - License, or (at your option) any later version. - - 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 - Library General Public License for more details. - - You should have received a copy of the GNU Library General Public - License along with the GNU C Library; see the file COPYING.LIB. If not, - write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, - Boston, MA 02111-1307, USA. */ - -#include - -#include "math_private.h" - -#ifdef NO_LONG_DOUBLE -/* The `long double' is in fact the IEEE `double' type. */ -static long double two52[2] = -{ - 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ - -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ -}; - - -long int -__rinttol (long double x) -{ - int32_t j0,sx; - u_int32_t i0,i1,i; - long double t, w; - long int result; - - EXTRACT_WORDS (i0, i1, x); - - sx = i0 >> 31; - j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; - - if (j0 < 20) - { - if (j0 < 0) - { - if (((i0 & 0x7fffffff) | i1) == 0) - /* The number is 0. */ - result = 0; - else - { - i1 |= i0; - i0 &= 0xfffe0000; - i0 |= ((i1 | -i1) >> 12) & 0x80000; - SET_HIGH_WORD (x, i0); - w = two52[sx] + x; - t = w - two52[sx]; - GET_HIGH_WORD (i0, t); - if ((i0 & 0x7fffffff) >= 0x3fff0000) - result = sx ? -1 : 1; - else - result = 0; - } - } - else - { - u_int32_t i = 0x000fffff >> j0; - if (((i0 & i) | i1) == 0) - { - /* X is not integral. */ - i >>= 1; - if (((i0 & i) | i1) != 0) - { - if (j0 == 19) - i1 = 0x40000000; - else - i0 = (i0 & (~i)) | (0x20000 >> j0); - - INSERT_WORDS (x, i0, i1); - w = two52[sx] + x; - x = w - two52[sx]; - EXTRACT_WORDS (i0, i1, x); - - j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; - } - } - - result = ((i0 >> (20 - j0)) & 0xfffff) | (0x00100000 >> (20 - j0)); - if (sx) - result = -result; - } - } - else if ((unsigned int) j0 < sizeof (long int) * 8 && j0 < 53) - { - i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); - if ((i1 & i) != 0) - { - /* x is not integral. */ - i >>= 1; - if ((i1 & i) != 0) - i1 = (i1 & (~i)) | (0x40000000 >> (j0 - 20)); - } - - INSERT_WORDS (x, i0, i1); - w = two52[sx] + x; - x = w - two52[sx]; - EXTRACT_WORDS (i0, i1, x); - - j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; - - result = i0 | 0x00100000; - if (j0 > 20) - { - result <<= j0 - 20; - result |= i1 >> (52 - j0); - } - if (sx) - result = -result; - } - else - /* Too large. The number is either +-inf or NaN or it is too - large to be effected by rounding. The standard leaves it - undefined what to return when the number is too large to fit in - a `long int'. */ - result = (long int) x; - - return result; -} - -#else -static long double two63[2] = -{ - 9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */ - -9.223372036854775808000000e+18 /* 0xC03E, 0x00000000, 0x00000000 */ -}; - - -long int -__rinttol (long double x) -{ - int32_t se,j0,sx; - u_int32_t i0,i1,i; - long int result; - long double w, t; - - GET_LDOUBLE_WORDS (se, i0, i1, x); - - sx = (se >> 15) & 1; - j0 = (se & 0x7fff) - 0x3fff; - - if (j0 < 31) - { - if (j0 < 0) - { - if (((se & 0x7fff) | i0 | i1) == 0) - /* The number is 0. */ - result = 0; - else - { - i1 |= i0; - i0 &= 0xe0000000; - i0 |= (i1 | -i1) & 0x80000000; - SET_LDOUBLE_MSW (x, i0); - w = two63[sx] + x; - t = w - two63[sx]; - GET_LDOUBLE_EXP (i0, t); - if ((i0 & 0x7fff) >= 0x3fff) - result = sx ? -1 : 1; - else - result = 0; - } - } - else - { - u_int32_t i = 0x7fffffff >> j0; - if (((i0 & i) | i1) == 0) - { - /* X is not integral. */ - i >>= 1; - if (((i0 & i) | i1) != 0) - { - if (j0 == 31) - i1 = 0x40000000; - else - i0 = (i0 & (~i)) | (0x20000000 >> j0); - - SET_LDOUBLE_WORDS (x, se, i0, i1); - w = two63[sx] + x; - x = w - two63[sx]; - GET_LDOUBLE_WORDS (se, i0, i1, x); - - sx = (se >> 15) & 1; - j0 = (se & 0x7fff) - 0x3fff; - } - } - - - result = i0 >> (31 - j0); - } - } - else if ((unsigned int) j0 < sizeof (long int) * 8 && j0 < 64) - { - i = ((u_int32_t) (0xffffffff)) >> (j0 - 31); - if ((i1 & i) != 0) - { - /* x is not integral. */ - i >>= 1; - if ((i1 & i) != 0) - i1 = (i1 & (~i)) | (0x40000000 >> (j0 - 31)); - } - - SET_LDOUBLE_WORDS (x, se, i0, i1); - w = two63[sx] + x; - x = w - two63[sx]; - GET_LDOUBLE_WORDS (se, i0, i1, x); - - j0 = (se & 0x7fff) - 0x3fff; - - result = i0; - if (j0 > 31) - { - result <<= j0 - 31; - result |= i1 >> (63 - j0); - } - } - else - /* Too large. The number is either +-inf or NaN or it is too - large to be effected by rounding. The standard leaves it - undefined what to return when the number is too large to fit in - a `long int'. */ - result = (long int) x; - - return result; -} -#endif - -weak_alias (__rinttol, rinttol) diff --git a/sysdeps/libm-ieee754/s_rinttoll.c b/sysdeps/libm-ieee754/s_rinttoll.c deleted file mode 100644 index b2fccd17b7..0000000000 --- a/sysdeps/libm-ieee754/s_rinttoll.c +++ /dev/null @@ -1,242 +0,0 @@ -/* Round argument to nearest integral value according to current rounding - direction. - Copyright (C) 1997 Free Software Foundation, Inc. - This file is part of the GNU C Library. - Contributed by Ulrich Drepper , 1997. - - The GNU C Library is free software; you can redistribute it and/or - modify it under the terms of the GNU Library General Public License as - published by the Free Software Foundation; either version 2 of the - License, or (at your option) any later version. - - 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 - Library General Public License for more details. - - You should have received a copy of the GNU Library General Public - License along with the GNU C Library; see the file COPYING.LIB. If not, - write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, - Boston, MA 02111-1307, USA. */ - -#include - -#include "math_private.h" - -#ifdef NO_LONG_DOUBLE -/* The `long double' is in fact the IEEE `double' type. */ -static long double two52[2] = -{ - 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ - -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ -}; - - -long long int -__rinttoll (long double x) -{ - int32_t j0,sx; - u_int32_t i0, i1, i; - long double t, w; - long long int result; - - EXTRACT_WORDS (i0, i1, x); - - sx = i0 >> 31; - j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; - - if (j0 < 20) - { - if (j0 < 0) - { - if (((i0 & 0x7fffffff) | i1) == 0) - /* The number is 0. */ - result = 0; - else - { - i1 |= i0; - i0 &= 0xfffe0000; - i0 |= ((i1 | -i1) >> 12) & 0x80000; - SET_HIGH_WORD (x, i0); - w = two52[sx] + x; - t = w - two52[sx]; - GET_HIGH_WORD (i0, t); - if ((i0 & 0x7fffffff) >= 0x3fff0000) - result = sx ? -1 : 1; - else - result = 0; - } - } - else - { - u_int32_t i = 0x000fffff >> j0; - if (((i0 & i) | i1) == 0) - { - /* X is not integral. */ - i >>= 1; - if (((i0 & i) | i1) != 0) - { - if (j0 == 19) - i1 = 0x40000000; - else - i0 = (i0 & (~i)) | (0x20000 >> j0); - - INSERT_WORDS (x, i0, i1); - w = two52[sx] + x; - x = w - two52[sx]; - EXTRACT_WORDS (i0, i1, x); - - j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; - } - } - - result = ((i0 >> (20 - j0)) & 0xfffff) | (0x00100000 >> (20 - j0)); - if (sx) - result = -result; - } - } - else if ((unsigned int) j0 < sizeof (long long int) * 8 && j0 < 53) - { - i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); - if ((i1 & i) != 0) - { - /* x is not integral. */ - i >>= 1; - if ((i1 & i) != 0) - i1 = (i1 & (~i)) | (0x40000000 >> (j0 - 20)); - } - - INSERT_WORDS (x, i0, i1); - w = two52[sx] + x; - x = w - two52[sx]; - EXTRACT_WORDS (i0, i1, x); - - j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; - - result = i0 | 0x00100000; - if (j0 > 20) - { - result <<= j0 - 20; - result |= i1 >> (52 - j0); - } - if (sx) - result = -result; - } - else - /* Too large. The number is either +-inf or NaN or it is too - large to be effected by rounding. The standard leaves it - undefined what to return when the number is too large to fit in - a `long int'. */ - result = (long long int) x; - - return result; -} - -#else -static long double two63[2] = -{ - 9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */ - -9.223372036854775808000000e+18 /* 0xC03E, 0x00000000, 0x00000000 */ -}; - - -long long int -__rinttoll (long double x) -{ - int32_t se,j0,sx; - u_int32_t i0, i1, i; - long long int result; - long double w, t; - - GET_LDOUBLE_WORDS (se, i0, i1, x); - - sx = (se >> 15) & 1; - j0 = (se & 0x7fff) - 0x3fff; - - if (j0 < 31) - { - if (j0 < 0) - { - if (((se & 0x7fff) | i0 | i1) == 0) - /* The number is 0. */ - result = 0; - else - { - i1 |= i0; - i0 &= 0xe0000000; - i0 |= (i1 | -i1) & 0x80000000; - SET_LDOUBLE_MSW (x, i0); - w = two63[sx] + x; - t = w - two63[sx]; - GET_LDOUBLE_EXP (i0, t); - if ((i0 & 0x7fff) >= 0x3fff) - result = sx ? -1 : 1; - else - result = 0; - } - } - else - { - u_int32_t i = 0x7fffffff >> j0; - if (((i0 & i) | i1) == 0) - { - /* X is not integral. */ - i >>= 1; - if (((i0 & i) | i1) != 0) - { - if (j0 == 31) - i1 = 0x40000000; - else - i0 = (i0 & (~i)) | (0x20000000 >> j0); - - SET_LDOUBLE_WORDS (x, se, i0, i1); - w = two63[sx] + x; - x = w - two63[sx]; - GET_LDOUBLE_WORDS (se, i0, i1, x); - - sx = (se >> 15) & 1; - j0 = (se & 0x7fff) - 0x3fff; - } - } - - - result = i0 >> (31 - j0); - } - } - else if ((unsigned int) j0 < sizeof (long long int) * 8 && j0 < 64) - { - i = ((u_int32_t) (0xffffffff)) >> (j0 - 31); - if ((i1 & i) != 0) - { - /* x is not integral. */ - i >>= 1; - if ((i1 & i) != 0) - i1 = (i1 & (~i)) | (0x40000000 >> (j0 - 31)); - } - - SET_LDOUBLE_WORDS (x, se, i0, i1); - w = two63[sx] + x; - x = w - two63[sx]; - GET_LDOUBLE_WORDS (se, i0, i1, x); - - j0 = (se & 0x7fff) - 0x3fff; - - result = i0; - if (j0 > 31) - { - result <<= j0 - 31; - result |= i1 >> (63 - j0); - } - } - else - /* Too large. The number is either +-inf or NaN or it is too - large to be effected by rounding. The standard leaves it - undefined what to return when the number is too large to fit in - a `long int'. */ - result = (long long int) x; - - return result; -} -#endif - -weak_alias (__rinttoll, rinttoll) diff --git a/sysdeps/libm-ieee754/s_roundtol.c b/sysdeps/libm-ieee754/s_roundtol.c deleted file mode 100644 index bc0ceae0f8..0000000000 --- a/sysdeps/libm-ieee754/s_roundtol.c +++ /dev/null @@ -1,177 +0,0 @@ -/* Round long double value to long int. - Copyright (C) 1997 Free Software Foundation, Inc. - This file is part of the GNU C Library. - Contributed by Ulrich Drepper , 1997. - - The GNU C Library is free software; you can redistribute it and/or - modify it under the terms of the GNU Library General Public License as - published by the Free Software Foundation; either version 2 of the - License, or (at your option) any later version. - - 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 - Library General Public License for more details. - - You should have received a copy of the GNU Library General Public - License along with the GNU C Library; see the file COPYING.LIB. If not, - write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, - Boston, MA 02111-1307, USA. */ - -#include - -#include "math_private.h" - - -#ifdef NO_LONG_DOUBLE -/* The `long double' is in fact the IEEE `double' type. */ - -long int -__roundtol (long double x) -{ - int32_t j0; - u_int32_t i1, i0; - long int result; - - EXTRACT_WORDS (i0, i1, x); - j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; - if (j0 < 20) - { - if (j0 < 0) - result = j0 < -1 ? 0 : ((i0 & 0x80000000) ? -1 : 1); - else - { - u_int32_t i = 0xfffff >> j0; - if (((i0 & i) | i1) == 0) - result = (long int) ((i0 & 0xfffff) | 0x100000) >> j0; - else - { - /* X is not integral. */ - u_int32_t j = i0 + (0x80000 >> j0); - if (j < i0) - result = (long int) 0x80000 >> (20 - j0); - else - result = (j | 0x100000) >> (20 - j0); - } - } - } - else if (j0 >= 8 * sizeof (long int) || j0 > 51) - { - /* The number is too large. It is left implementation defined - what happens. */ - result = (long int) x; - } - else - { - u_int32_t i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); - if ((i1 & i) != 0) - { - /* x is not integral. */ - u_int32_t j = i1 + (0x80000000 >> (j0 - 20)); - if (j < i1) - { - j = i0 + 1; - if ((j & 0xfffff) == 0) - { - if (sizeof (long int) <= 4) - /* Overflow. */ - result = (long int) x; - else - result = 1l << (j0 + 1); - } - else - result = (long int) ((i0 & 0xfffff) | 0x100000) << (j0 - 31); - } - else - { - result = (long int) ((i0 & 0xfffff) | 0x100000) << (j0 - 31); - if (sizeof (long int) > 4 && j0 > 31) - result |= j >> (63 - j0); - } - } - else - { - result = (long int) ((i0 & 0xfffff) | 0x100000) << (j0 - 31); - if (sizeof (long int) > 4 && j0 > 31) - result |= i1 >> (63 - j0); - } - } - - return i0 & 0x80000000 ? -result : result; -} -#else -long int -__roundtol (long double x) -{ - int32_t j0; - u_int32_t se, i1, i0; - long int result; - - GET_LDOUBLE_WORDS (se, i0, i1, x); - j0 = (se & 0x7fff) - 0x3fff; - if (j0 < 31) - { - if (j0 < 0) - result = j0 < -1 ? 0 : 1; - else - { - u_int32_t i = 0x7fffffff >> j0; - if (((i0 & i) | i1) == 0) - result = (long int) i0 >> j0; - else - { - /* X is not integral. */ - u_int32_t j = i0 + (0x40000000 >> j0); - if (j < i0) - result = 0x80000000l >> (30 - j0); - else - result = j >> (31 - j0); - } - } - } - else if ((unsigned int) j0 >= 8 * sizeof (long int) || j0 > 62) - { - /* The number is too large. It is left implementation defined - what happens. */ - result = (long int) x; - } - else - { - u_int32_t i = ((u_int32_t) (0xffffffff)) >> (j0 - 31); - if ((i1 & i) != 0) - { - /* x is not integral. */ - u_int32_t j = i1 + (0x80000000 >> (j0 - 31)); - if (j < i1) - { - j = i0 + 1; - if (j == 0) - { - if (sizeof (long int) <= 4) - /* Overflow. */ - result = (long int) x; - else - result = 1l << (j0 + 1); - } - else - result = (long int) i0 << (j0 - 31); - } - else - { - result = (long int) i0 << (j0 - 31); - if (sizeof (long int) > 4 && j0 > 31) - result |= j >> (63 - j0); - } - } - else - { - result = (long int) i0 << (j0 - 31); - if (sizeof (long int) > 4 && j0 > 31) - result |= i1 >> (63 - j0); - } - } - - return se & 0x8000 ? -result : result; -} -#endif -weak_alias (__roundtol, roundtol) diff --git a/sysdeps/libm-ieee754/s_roundtoll.c b/sysdeps/libm-ieee754/s_roundtoll.c deleted file mode 100644 index 49167236a6..0000000000 --- a/sysdeps/libm-ieee754/s_roundtoll.c +++ /dev/null @@ -1,179 +0,0 @@ -/* Round long double value to long long int. - Copyright (C) 1997 Free Software Foundation, Inc. - This file is part of the GNU C Library. - Contributed by Ulrich Drepper , 1997. - - The GNU C Library is free software; you can redistribute it and/or - modify it under the terms of the GNU Library General Public License as - published by the Free Software Foundation; either version 2 of the - License, or (at your option) any later version. - - 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 - Library General Public License for more details. - - You should have received a copy of the GNU Library General Public - License along with the GNU C Library; see the file COPYING.LIB. If not, - write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, - Boston, MA 02111-1307, USA. */ - -#include - -#include "math_private.h" - - -#ifdef NO_LONG_DOUBLE -/* The `long double' is in fact the IEEE `double' type. */ - -long long int -__roundtoll (long double x) -{ - int32_t j0; - u_int32_t i1, i0; - long long int result; - - EXTRACT_WORDS (i0, i1, x); - j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; - if (j0 < 20) - { - if (j0 < 0) - result = j0 < -1 ? 0 : ((i0 & 0x80000000) ? -1 : 1); - else - { - u_int32_t i = 0xfffff >> j0; - if (((i0 & i) | i1) == 0) - result = (long long int) ((i0 & 0xfffff) | 0x100000) >> j0; - else - { - /* X is not integral. */ - u_int32_t j = i0 + (0x80000 >> j0); - if (j < i0) - result = (long long int) 0x80000 >> (20 - j0); - else - result = (j | 0x100000) >> (20 - j0); - } - } - } - else if (j0 >= 8 * sizeof (long long int) || j0 > 51) - { - /* The number is too large. It is left implementation defined - what happens. */ - result = (long long int) x; - } - else - { - u_int32_t i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); - if ((i1 & i) != 0) - { - /* x is not integral. */ - u_int32_t j = i1 + (0x80000000 >> (j0 - 20)); - if (j < i1) - { - j = i0 + 1; - if ((j & 0xfffff) == 0) - { - if (sizeof (long long int) <= 4) - /* Overflow. */ - result = (long long int) x; - else - result = 1ll << (j0 + 1); - } - else - result = ((long long int) ((i0 & 0xfffff) | 0x100000) - << (j0 - 31)); - } - else - { - result = ((long long int) ((i0 & 0xfffff) | 0x100000) - << (j0 - 31)); - if (sizeof (long long int) > 4 && j0 > 31) - result |= j >> (63 - j0); - } - } - else - { - result = (long long int) ((i0 & 0xfffff) | 0x100000) << (j0 - 31); - if (sizeof (long long int) > 4 && j0 > 31) - result |= i1 >> (63 - j0); - } - } - - return i0 & 0x80000000 ? -result : result; -} -#else -long long int -__roundtoll (long double x) -{ - int32_t j0; - u_int32_t se, i1, i0; - long long int result; - - GET_LDOUBLE_WORDS (se, i0, i1, x); - j0 = (se & 0x7fff) - 0x3fff; - if (j0 < 31) - { - if (j0 < 0) - result = j0 < -1 ? 0 : 1; - else - { - u_int32_t i = 0x7fffffff >> j0; - if (((i0 & i) | i1) == 0) - result = (long long int) i0 >> j0; - else - { - /* X is not integral. */ - u_int32_t j = i0 + (0x40000000 >> j0); - if (j < i0) - result = 0x80000000l >> (30 - j0); - else - result = j >> (31 - j0); - } - } - } - else if ((unsigned int) j0 >= 8 * sizeof (long long int) || j0 > 62) - { - /* The number is too large. It is left implementation defined - what happens. */ - result = (long long int) x; - } - else - { - u_int32_t i = ((u_int32_t) (0xffffffff)) >> (j0 - 31); - if ((i1 & i) != 0) - { - /* x is not integral. */ - u_int32_t j = i1 + (0x80000000 >> (j0 - 31)); - if (j < i1) - { - j = i0 + 1; - if (j == 0) - { - if (sizeof (long long int) <= 4) - /* Overflow. */ - result = (long long int) x; - else - result = 1ll << (j0 + 1); - } - else - result = (long long int) i0 << (j0 - 31); - } - else - { - result = (long long int) i0 << (j0 - 31); - if (sizeof (long long int) > 4 && j0 > 31) - result |= j >> (63 - j0); - } - } - else - { - result = (long long int) i0 << (j0 - 31); - if (sizeof (long long int) > 4 && j0 > 31) - result |= i1 >> (63 - j0); - } - } - - return se & 0x8000 ? -result : result; -} -#endif -weak_alias (__roundtoll, roundtoll) diff --git a/sysdeps/libm-ieee754/s_scalbn.c b/sysdeps/libm-ieee754/s_scalbn.c index 439b966e69..0ce0ffdf2f 100644 --- a/sysdeps/libm-ieee754/s_scalbn.c +++ b/sysdeps/libm-ieee754/s_scalbn.c @@ -35,10 +35,10 @@ huge = 1.0e+300, tiny = 1.0e-300; #ifdef __STDC__ - double __scalbn (double x, int n) + double __scalbn (double x, long int n) #else double __scalbn (x,n) - double x; int n; + double x; long int n; #endif { int32_t k,hx,lx; @@ -49,17 +49,16 @@ tiny = 1.0e-300; x *= two54; GET_HIGH_WORD(hx,x); k = ((hx&0x7ff00000)>>20) - 54; - if (n< -50000) return tiny*x; /*underflow*/ } if (k==0x7ff) return x+x; /* NaN or Inf */ k = k+n; - if (k > 0x7fe) return huge*__copysign(huge,x); /* overflow */ + if (n> 50000 || k > 0x7fe) + return huge*__copysign(huge,x); /* overflow */ + if (n< -50000) return tiny*__copysign(tiny,x); /*underflow*/ if (k > 0) /* normal result */ {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;} if (k <= -54) - if (n > 50000) /* in case integer overflow in n+k */ - return huge*__copysign(huge,x); /*overflow*/ - else return tiny*__copysign(tiny,x); /*underflow*/ + return tiny*__copysign(tiny,x); /*underflow*/ k += 54; /* subnormal result */ SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x*twom54; diff --git a/sysdeps/libm-ieee754/s_scalbnf.c b/sysdeps/libm-ieee754/s_scalbnf.c index 3a83e54c56..4799c825c9 100644 --- a/sysdeps/libm-ieee754/s_scalbnf.c +++ b/sysdeps/libm-ieee754/s_scalbnf.c @@ -8,7 +8,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -31,10 +31,10 @@ huge = 1.0e+30, tiny = 1.0e-30; #ifdef __STDC__ - float __scalbnf (float x, int n) + float __scalbnf (float x, long int n) #else float __scalbnf (x,n) - float x; int n; + float x; long int n; #endif { int32_t k,ix; @@ -44,18 +44,18 @@ tiny = 1.0e-30; if ((ix&0x7fffffff)==0) return x; /* +-0 */ x *= two25; GET_FLOAT_WORD(ix,x); - k = ((ix&0x7f800000)>>23) - 25; - if (n< -50000) return tiny*x; /*underflow*/ + k = ((ix&0x7f800000)>>23) - 25; } if (k==0xff) return x+x; /* NaN or Inf */ - k = k+n; - if (k > 0xfe) return huge*copysignf(huge,x); /* overflow */ + k = k+n; + if (n> 50000 || k > 0xfe) + return huge*copysignf(huge,x); /* overflow */ + if (n< -50000) + return tiny*copysignf(tiny,x); /*underflow*/ if (k > 0) /* normal result */ {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;} if (k <= -25) - if (n > 50000) /* in case integer overflow in n+k */ - return huge*copysignf(huge,x); /*overflow*/ - else return tiny*copysignf(tiny,x); /*underflow*/ + return tiny*copysignf(tiny,x); /*underflow*/ k += 25; /* subnormal result */ SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x*twom25; diff --git a/sysdeps/libm-ieee754/s_scalbnl.c b/sysdeps/libm-ieee754/s_scalbnl.c index 6825e9c7b3..3e80d85136 100644 --- a/sysdeps/libm-ieee754/s_scalbnl.c +++ b/sysdeps/libm-ieee754/s_scalbnl.c @@ -33,16 +33,16 @@ static const long double #else static long double #endif -two54 = 1.80143985094819840000e+16, /* 0x4035, 0x00000000, 0x00000000 */ -twom54 = 5.55111512312578270212e-17, /* 0x3FC9, 0x00000000, 0x00000000 */ +two63 = 4.50359962737049600000e+15, +twom63 = 1.08420217248550443400e-19; huge = 1.0e+4900L, tiny = 1.0e-4900L; #ifdef __STDC__ - long double __scalbnl (long double x, int n) + long double __scalbnl (long double x, long int n) #else long double __scalbnl (x,n) - long double x; int n; + long double x; long int n; #endif { int32_t k,es,hx,lx; @@ -50,22 +50,22 @@ tiny = 1.0e-4900L; k = es&0x7fff; /* extract exponent */ if (k==0) { /* 0 or subnormal x */ if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */ - x *= two54; - GET_HIGH_WORD(hx,x); - k = ((hx&0x7ff00000)>>20) - 54; - if (n< -50000) return tiny*x; /*underflow*/ + x *= two63; + GET_LDOUBLE_EXP(es,x); + k = (hx&0x7fff) - 63; } - if (k==0x7ff) return x+x; /* NaN or Inf */ + if (k==0x7fff) return x+x; /* NaN or Inf */ k = k+n; - if (k > 0x7fe) return huge*__copysign(huge,x); /* overflow */ + if (n> 50000 || k > 0x7ffe) + return huge*__copysignl(huge,x); /* overflow */ + if (n< -50000) + return tiny*__copysignl(tiny,x); if (k > 0) /* normal result */ - {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;} - if (k <= -54) - if (n > 50000) /* in case integer overflow in n+k */ - return huge*__copysign(huge,x); /*overflow*/ - else return tiny*__copysign(tiny,x); /*underflow*/ + {SET_LDOUBLE_EXP(x,(es&0x8000)|k); return x;} + if (k <= -63) + return tiny*__copysignl(tiny,x); /*underflow*/ k += 54; /* subnormal result */ - SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); - return x*twom54; + SET_LDOUBLE_EXP(x,(es&0x8000)|k); + return x*twom63; } weak_alias (__scalbnl, scalbnl) diff --git a/sysdeps/libm-ieee754/w_gamma.c b/sysdeps/libm-ieee754/w_gamma.c index 49e4bcf3c3..87a3408c3f 100644 --- a/sysdeps/libm-ieee754/w_gamma.c +++ b/sysdeps/libm-ieee754/w_gamma.c @@ -15,16 +15,13 @@ static char rcsid[] = "$NetBSD: w_gamma.c,v 1.7 1995/11/20 22:06:43 jtc Exp $"; #endif /* double gamma(double x) - * Return the logarithm of the Gamma function of x. - * - * Method: call gamma_r + * Return the logarithm of the Gamma function of x or the Gamma function of x, + * depending on the library mode. */ #include "math.h" #include "math_private.h" -extern int signgam; - #ifdef __STDC__ double __gamma(double x) #else @@ -32,12 +29,19 @@ extern int signgam; double x; #endif { + int signgam; + double y; + if (_LIB_VERSION == _SVID_) + y = __ieee754_lgamma_r(x,&signgam); + else + { + y = __ieee754_gamma_r(x,&signgam); + if (signgam < 0) y = -y; #ifdef _IEEE_LIBM - return __ieee754_lgamma_r(x,&signgam); + return y; #else - double y; - y = __ieee754_lgamma_r(x,&signgam); - if(_LIB_VERSION == _IEEE_) return y; + if(_LIB_VERSION == _IEEE_) return y; + } if(!__finite(y)&&__finite(x)) { if(__floor(x)==x&&x<=0.0) return __kernel_standard(x,x,41); /* gamma pole */ diff --git a/sysdeps/libm-ieee754/w_gamma_r.c b/sysdeps/libm-ieee754/w_gamma_r.c deleted file mode 100644 index f9efc8cf2e..0000000000 --- a/sysdeps/libm-ieee754/w_gamma_r.c +++ /dev/null @@ -1,51 +0,0 @@ -/* @(#)wr_gamma.c 5.1 93/09/24 */ -/* - * ==================================================== - * 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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_gamma_r.c,v 1.7 1995/11/20 22:06:45 jtc Exp $"; -#endif - -/* - * wrapper double gamma_r(double x, int *signgamp) - */ - -#include "math.h" -#include "math_private.h" - - -#ifdef __STDC__ - double __gamma_r(double x, int *signgamp) /* wrapper lgamma_r */ -#else - double __gamma_r(x,signgamp) /* wrapper lgamma_r */ - double x; int *signgamp; -#endif -{ -#ifdef _IEEE_LIBM - return __ieee754_lgamma_r(x,signgamp); -#else - double y; - y = __ieee754_lgamma_r(x,signgamp); - if(_LIB_VERSION == _IEEE_) return y; - if(!__finite(y)&&__finite(x)) { - if(__floor(x)==x&&x<=0.0) - return __kernel_standard(x,x,41); /* gamma pole */ - else - return __kernel_standard(x,x,40); /* gamma overflow */ - } else - return y; -#endif -} -weak_alias (__gamma_r, gamma_r) -#ifdef NO_LONG_DOUBLE -strong_alias (__gamma_r, __gammal_r) -weak_alias (__gamma_r, gammal_r) -#endif diff --git a/sysdeps/libm-ieee754/w_gammaf.c b/sysdeps/libm-ieee754/w_gammaf.c index c2d21daf47..ecde183824 100644 --- a/sysdeps/libm-ieee754/w_gammaf.c +++ b/sysdeps/libm-ieee754/w_gammaf.c @@ -8,7 +8,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -20,8 +20,6 @@ static char rcsid[] = "$NetBSD: w_gammaf.c,v 1.4 1995/11/20 22:06:48 jtc Exp $"; #include "math.h" #include "math_private.h" -extern int signgam; - #ifdef __STDC__ float __gammaf(float x) #else @@ -29,12 +27,19 @@ extern int signgam; float x; #endif { + int signgam; + float y; + if (_LIB_VERSION == _SVID_) + y = __ieee754_lgammaf_r(x,&signgam); + else + { + y = __ieee754_gammaf_r(x,&signgam); + if (signgam < 0) y = -y; #ifdef _IEEE_LIBM - return __ieee754_lgammaf_r(x,&signgam); + return y; #else - float y; - y = __ieee754_lgammaf_r(x,&signgam); - if(_LIB_VERSION == _IEEE_) return y; + if(_LIB_VERSION == _IEEE_) return y; + } if(!__finitef(y)&&__finitef(x)) { if(__floorf(x)==x&&x<=(float)0.0) /* gammaf pole */ @@ -45,5 +50,5 @@ extern int signgam; } else return y; #endif -} +} weak_alias (__gammaf, gammaf) diff --git a/sysdeps/libm-ieee754/w_gammaf_r.c b/sysdeps/libm-ieee754/w_gammaf_r.c deleted file mode 100644 index b15d5e84bb..0000000000 --- a/sysdeps/libm-ieee754/w_gammaf_r.c +++ /dev/null @@ -1,52 +0,0 @@ -/* w_gammaf_r.c -- float version of w_gamma_r.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: w_gammaf_r.c,v 1.4 1995/11/20 22:06:50 jtc Exp $"; -#endif - -/* - * wrapper float gammaf_r(float x, int *signgamp) - */ - -#include "math.h" -#include "math_private.h" - - -#ifdef __STDC__ - float __gammaf_r(float x, int *signgamp) /* wrapper lgammaf_r */ -#else - float __gammaf_r(x,signgamp) /* wrapper lgammaf_r */ - float x; int *signgamp; -#endif -{ -#ifdef _IEEE_LIBM - return __ieee754_lgammaf_r(x,signgamp); -#else - float y; - y = __ieee754_lgammaf_r(x,signgamp); - if(_LIB_VERSION == _IEEE_) return y; - if(!__finitef(y)&&__finitef(x)) { - if(__floorf(x)==x&&x<=(float)0.0) - /* gammaf pole */ - return (float)__kernel_standard((double)x,(double)x,141); - else - /* gamma overflow */ - return (float)__kernel_standard((double)x,(double)x,140); - } else - return y; -#endif -} -weak_alias (__gammaf_r, gammaf_r) diff --git a/sysdeps/libm-ieee754/w_gammal.c b/sysdeps/libm-ieee754/w_gammal.c index f72e2e7e44..605a4167c4 100644 --- a/sysdeps/libm-ieee754/w_gammal.c +++ b/sysdeps/libm-ieee754/w_gammal.c @@ -19,16 +19,12 @@ static char rcsid[] = "$NetBSD: $"; #endif /* long double gammal(double x) - * Return the logarithm of the Gamma function of x. - * - * Method: call gammal_r + * Return the Gamma function of x. */ #include "math.h" #include "math_private.h" -extern int signgam; - #ifdef __STDC__ long double __gammal(long double x) #else @@ -36,12 +32,19 @@ extern int signgam; long double x; #endif { + long double y; + int signgam; + if (_LIB_VERSION == _SVID_) + y = __ieee754_lgammal_r(x,&signgam); + else + { + y = __ieee754_gammal_r(x,&signgam); + if (signgam < 0) y = -y; #ifdef _IEEE_LIBM - return __ieee754_lgammal_r(x,&signgam); + return y; #else - long double y; - y = __ieee754_lgammal_r(x,&signgam); - if(_LIB_VERSION == _IEEE_) return y; + if(_LIB_VERSION == _IEEE_) return y; + } if(!__finitel(y)&&__finitel(x)) { if(__floorl(x)==x&&x<=0.0) return __kernel_standard(x,x,241); /* gamma pole */ diff --git a/sysdeps/libm-ieee754/w_gammal_r.c b/sysdeps/libm-ieee754/w_gammal_r.c deleted file mode 100644 index ad6a43ae9b..0000000000 --- a/sysdeps/libm-ieee754/w_gammal_r.c +++ /dev/null @@ -1,52 +0,0 @@ -/* w_gammal_r.c -- long double version of w_gamma_r.c. - * Conversion to long double by Ulrich Drepper, - * Cygnus Support, drepper@cygnus.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. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: $"; -#endif - -/* - * wrapper long double gammal_r(long double x, int *signgamp) - */ - -#include "math.h" -#include "math_private.h" - - -#ifdef __STDC__ - long double __gammal_r(long double x, int *signgamp) - /* wrapper lgammal_r */ -#else - long double __gammal_r(x,signgamp) /* wrapper lgamma_r */ - long double x; int *signgamp; -#endif -{ -#ifdef _IEEE_LIBM - return __ieee754_lgammal_r(x,signgamp); -#else - long double y; - y = __ieee754_lgammal_r(x,signgamp); - if(_LIB_VERSION == _IEEE_) return y; - if(!__finitel(y)&&__finitel(x)) { - if(__floorl(x)==x&&x<=0.0) - return __kernel_standard(x,x,241); /* gamma pole */ - else - return __kernel_standard(x,x,240); /* gamma overflow */ - } else - return y; -#endif -} -weak_alias (__gammal_r, gammal_r) -- cgit v1.2.3