diff options
| author | (no author) <(no author)@280ebfd0-de03-0410-8827-d642c229c3f4> | 2005-09-19 09:36:55 +0000 |
|---|---|---|
| committer | (no author) <(no author)@280ebfd0-de03-0410-8827-d642c229c3f4> | 2005-09-19 09:36:55 +0000 |
| commit | 18da4e328c2107b75594b979ca0946a96edffc2f (patch) | |
| tree | fad16c5bd12ce7d157751e451220fdb21397eee7 | |
| parent | ac6c37d815d6c7a02f061020553b06103c9a9c87 (diff) | |
| download | mpfr-18da4e328c2107b75594b979ca0946a96edffc2f.tar.gz | |
This commit was manufactured by cvs2svn to create tagmpfr-2-2-0-rel
'mpfr-2-2-0-rel'.
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/tags/mpfr-2-2-0-rel@3866 280ebfd0-de03-0410-8827-d642c229c3f4
| -rw-r--r-- | ChangeLog | 138 | ||||
| -rw-r--r-- | INSTALL | 17 | ||||
| -rw-r--r-- | Makefile.am | 2 | ||||
| -rw-r--r-- | NEWS | 3 | ||||
| -rw-r--r-- | README.dev | 26 | ||||
| -rw-r--r-- | TODO | 1 | ||||
| -rw-r--r-- | acinclude.m4 | 83 | ||||
| -rw-r--r-- | add1sp.c | 10 | ||||
| -rw-r--r-- | algorithms.tex | 60 | ||||
| -rw-r--r-- | configure.in | 7 | ||||
| -rw-r--r-- | const_euler.c | 56 | ||||
| -rw-r--r-- | constant.c | 9 | ||||
| -rw-r--r-- | coth.c | 2 | ||||
| -rw-r--r-- | exp.c | 3 | ||||
| -rwxr-xr-x | fixperm | 8 | ||||
| -rw-r--r-- | gamma.c | 322 | ||||
| -rw-r--r-- | int_ceil_log2.c | 2 | ||||
| -rw-r--r-- | lngamma.c | 458 | ||||
| -rw-r--r-- | log.c | 5 | ||||
| -rw-r--r-- | mpfr-gmp.h | 14 | ||||
| -rw-r--r-- | mpfr-impl.h | 37 | ||||
| -rw-r--r-- | mpfr.h | 1 | ||||
| -rw-r--r-- | mpfr.texi | 20 | ||||
| -rw-r--r-- | mul.c | 239 | ||||
| -rw-r--r-- | set_str_raw.c | 113 | ||||
| -rw-r--r-- | sub1sp.c | 10 | ||||
| -rw-r--r-- | sum.c | 150 | ||||
| -rw-r--r-- | tests/Makefile.am | 4 | ||||
| -rw-r--r-- | tests/tcoth.c | 24 | ||||
| -rw-r--r-- | tests/texp.c | 16 | ||||
| -rw-r--r-- | tests/tgamma.c | 172 | ||||
| -rw-r--r-- | tests/tlngamma.c | 192 | ||||
| -rw-r--r-- | tests/tlog.c | 11 | ||||
| -rw-r--r-- | tests/trint.c | 2 | ||||
| -rw-r--r-- | tests/tset_f.c | 10 | ||||
| -rw-r--r-- | tests/tstckintc.c | 1 | ||||
| -rw-r--r-- | tests/tsub.c | 2 | ||||
| -rw-r--r-- | tests/tversion.c | 4 | ||||
| -rw-r--r-- | tuneup.c | 3 | ||||
| -rw-r--r-- | zeta.c | 2 |
40 files changed, 1667 insertions, 572 deletions
@@ -1,5 +1,143 @@ +2005-09-18 Vincent Lefevre <Vincent.Lefevre@loria.fr> + + * INSTALL: + Mentioned problems with the Tru64 make and other minor changes. + [2.2 branch] + +2005-09-17 Vincent Lefevre <Vincent.Lefevre@loria.fr> + + * configure.in: [2.2 branch] + Fixed the GMP linking test: replaced __gmp_version, which corresponds + to a variable, by __gmpz_init (suggested by the GMP documentation), + which corresponds to a function, as the auto tools assume the symbol + corresponds to a function. It was failing with the AIX linker. + +2005-09-16 Vincent Lefevre <Vincent.Lefevre@loria.fr> + + * zeta.c: Better initial precision. [2.2 branch] + + * const_euler.c: + Fixed overflow problem for large precisions. [2.2 branch] + + * tests/tcoth.c, coth.c: + Fixed stupid bug (coth was computing 1/tan instead of 1/tanh). + [2.2 branch] + + * mpfr.texi: Propagated the changes from the trunk to the 2.2 branch: + - documented the mpfr_get_f and mpfr_pow_z functions; + - fixed typos found by Tomonori Kouya in mpfr_atan2 description; + - updated Contributors section. + +2005-09-13 Vincent Lefevre <Vincent.Lefevre@loria.fr> + + * exp.c, tests/texp.c: + Fixed bug in exp(-eps) for rounding toward zero (test was also wrong). + [2.2 branch] + + * constant.c: Added missing MPFR_THREAD_ATTR. [2.2 branch] + + * NEWS: Minor update. [2.2 branch] + +2005-09-11 Vincent Lefevre <Vincent.Lefevre@loria.fr> + + * mparam_h.in, mul.c, tuneup.c: Fixed tuneup. [2.2 branch] + + * README.dev, mpfr-impl.h, tests/trint.c, tests/tstckintc.c: + [2.2 branch] + Fixed a #ifdef in mpfr-impl.h (when --enable-assert isn't used). + Fixed uninitialized variable in trint.c (with gcc -std=c99). + Added #include <string.h> for memmove in tstckintc.c (see ISO C99). + + * README.dev, int_ceil_log2.c, mpfr-impl.h: + Minor patches from the trunk. [2.2 branch] + +2005-09-09 Vincent Lefevre <Vincent.Lefevre@loria.fr> + + * tests/tset_f.c, algorithms.tex, gamma.c, lngamma.c, mparam_h.in, mul.c, tests/tgamma.c, tests/tlngamma.c: + Ported patches from the trunk to the 2.2 branch. + +2005-09-08 Vincent Lefevre <Vincent.Lefevre@loria.fr> + + * README.dev, fixperm: + Added fixperm script to fix the file permissions broken by CVS. + [2.2 branch] + + * fixperm: + file fixperm was added on branch mpfr-2-2-branch on 2005-09-08 09:05:58 +0000 + + * tests/tlngamma.c: Forgot to add tlngamma.c. [2.2 branch] + + * mpfr.texi: Update (mpfr_lngamma). [2.2 branch] + + * lngamma.c: Forgot to add lngamma.c. [2.2 branch] + + * tests/Makefile.am: Fix: readded dependencies. [2.2 branch] + +2005-09-06 Vincent Lefevre <Vincent.Lefevre@loria.fr> + + * Makefile.am, NEWS, TODO, gamma.c, mpfr-impl.h, mpfr.h, tests/Makefile.am, tests/tgamma.c, tests/tversion.c: + Patches from the trunk: added mpfr_lngamma, fixed mpfr_gamma, + more detailed error message in tversion. [2.2 branch] + +2005-09-06 Vincent Lefevre <Vincent.Lefevre@loria.fr> + + * set_str_raw.c, tests/tsub.c: + Ported patches from the trunk to the 2.2 branch. + + * mpfr-impl.h: + Replaced unsigned long bit-fields (GCC extension) by insigned int + bit-fields. [2.2 branch] + +2005-09-05 Vincent Lefevre <Vincent.Lefevre@loria.fr> + + * mpfr-gmp.h: + MPN_COPY now calls memcpy only if dst != src (otherwise this is an + undefined behavior), and if WANT_ASSERT is defined, it checks that + there is no overlap. [2.2 branch] + +2005-09-02 Vincent Lefevre <Vincent.Lefevre@loria.fr> + + * acinclude.m4, mpfr-impl.h: + Avoid MIPSpro / IRIX64 (incorrect) optimizations for DOUBLE_ISNAN. + [2.2 branch] + + * add1sp.c, sub1sp.c: + Do not put expressions with side effects in MPFR_ASSERTN. [2.2 branch] + +2005-09-01 Vincent Lefevre <Vincent.Lefevre@loria.fr> + + * mpfr.texi: Updated month. [2.2 branch] + + * log.c, tests/tlog.c: Ported fix to the 2.2 branch. + +2005-08-31 Vincent Lefevre <Vincent.Lefevre@loria.fr> + + * sum.c: + Fixed memory leak (missing MPFR_TMP_FREE in a particular case). + [2.2 branch, patch ported from the trunk] + + * mul.c: + Fixed memory leak (missing MPFR_TMP_FREE in a particular case). + [2.2 branch] + +2005-08-31 Vincent Lefevre <Vincent.Lefevre@loria.fr> + + * configure.in: + Fix for the following problem: #error is not sufficient with cc + on IRIX64, and autoconf 2.59 doesn't cope with that. [2.2 branch] + +2005-08-30 Vincent Lefevre <Vincent.Lefevre@loria.fr> + + * mul.c: + Patch from Patrick to solve efficiency problem when one operand is + sparse (e.g. from ui_pow_ui) + other minor changes. [2.2 branch] + 2005-08-24 Vincent Lefevre <Vincent.Lefevre@loria.fr> + * add1sp.c, mpfr-impl.h, mul.c, print_raw.c, sub1sp.c: + When WANT_ASSERT >= 2: the corresponding messages are now output + to stderr instead of stdout. + * set_uj.c: Improved C compliance when mp_limb_t and uintmax_t have the same size (concerning >> BITS_PER_MP_LIMB); some other small improvements. @@ -13,7 +13,7 @@ 3. To check the built library (runs the test files), type: make check -4. To install it (default `/usr/local` | see `--prefix` option), type: +4. To install it (default "/usr/local" | see "--prefix" option), type: make install @@ -45,8 +45,8 @@ Building MPFR with internal GMP header files ./configure options =================== ---prefix=DIR installs MPFR headers and library in DIR/include - and DIR/lib respectively (the default is `/usr/local`). +--prefix=DIR installs MPFR headers and library in DIR/include and + DIR/lib respectively (the default is "/usr/local"). --with-gmp-include=DIR assumes that DIR contains gmp.h [, gmp-impl.h, ...] --with-gmp-lib=DIR assumes that DIR contains libgmp.a @@ -67,8 +67,8 @@ Run "./configure --help" to see the other options (autoconf default options). In case of problem ================== - If the `configure` fails, please check that the C compiler and - its options are the same as those for the GMP build (specially the ABI). + If the "configure" fails, please check that the C compiler and its + options are the same as those for the GMP build (specially the ABI). You can see the latter with the following: grep "^CC\|^CFLAGS" GMPBUILD/Makefile @@ -81,7 +81,8 @@ In case of problem those used by GMP, for example, on powerpc-aix: make AR="ar -X64" - On some architectures, try with `gmake` instead of `make`. + On some architectures, try with "gmake" (GNU make) instead of "make". + Problems have been reported with the Tru64 make. Try to build MPFR with/without GMP internal files. @@ -90,10 +91,10 @@ In case of problem linker. There are two solutions: 1. setenv LD_LIBRARY_PATH /usr/lib:/usr/local/lib [tcsh] or export LD_LIBRARY_PATH=/usr/lib:/usr/local/lib [bash] - 2. Or run MPFR's configure with "--with-gmp=/usr/local/' + 2. Or run MPFR's configure with "--with-gmp=/usr/local/" so that libtool is able to fix this problem. - If you can't solve your problem, you could contact us at mpfr@loria.fr, + If you can't solve your problem, you could contact us at <mpfr@loria.fr>, indicating the machine and operating system used (uname -a), the compiler and version used (gcc -v if you use gcc), the compile options used if any, the version of GMP and MPFR used and a description of the problem diff --git a/Makefile.am b/Makefile.am index f41746f9a..2bbae8961 100644 --- a/Makefile.am +++ b/Makefile.am @@ -7,7 +7,7 @@ include_HEADERS = mpfr.h mpf2mpfr.h lib_LTLIBRARIES = libmpfr.la -libmpfr_la_SOURCES = mpfr.h mpf2mpfr.h mpfr-gmp.h mpfr-impl.h mpfr-longlong.h mpfr-test.h log_b2.h exceptions.c extract.c uceil_exp2.c uceil_log2.c ufloor_log2.c add.c add1.c add_ui.c agm.c clear.c cmp.c cmp_abs.c cmp_si.c cmp_ui.c comparisons.c div_2exp.c div_2si.c div_2ui.c div.c div_ui.c dump.c eq.c exp10.c exp2.c exp3.c exp.c frac.c get_d.c get_exp.c get_str.c init.c inp_str.c isinteger.c isinf.c isnan.c isnum.c const_log2.c log.c mul_2exp.c mul_2si.c mul_2ui.c mul.c mul_ui.c neg.c next.c out_str.c const_pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c round_prec.c set.c setmax.c setmin.c set_d.c set_dfl_prec.c set_exp.c set_rnd.c set_f.c set_prc_raw.c set_prec.c set_q.c set_si.c set_str.c set_str_raw.c set_ui.c set_z.c sqrt.c sqrt_ui.c sub.c sub1.c sub_ui.c rint.c ui_div.c ui_sub.c urandomb.c get_z_exp.c swap.c factorial.c cosh.c sinh.c tanh.c acosh.c asinh.c atanh.c atan.c cmp2.c exp_2.c asin.c const_euler.c cos.c sin.c tan.c fma.c hypot.c log1p.c expm1.c log2.c log10.c ui_pow.c ui_pow_ui.c minmax.c dim.c copysign.c gmp_op.c init2.c acos.c sin_cos.c set_nan.c set_inf.c powerof2.c gamma.c set_ld.c get_ld.c cbrt.c volatile.c fits_s.h fits_sshort.c fits_sint.c fits_slong.c fits_u.h fits_ushort.c fits_uint.c fits_ulong.c fits_uintmax.c fits_intmax.c get_si.c get_ui.c zeta.c cmp_d.c erf.c inits.c inits2.c clears.c sgn.c check.c sub1sp.c version.c mpn_exp.c mpfr-gmp.c mp_clz_tab.c sum.c add1sp.c free_cache.c si_op.c cmp_ld.c set_ui_2exp.c set_si_2exp.c set_uj.c set_sj.c get_sj.c get_uj.c get_z.c iszero.c cache.c sqr.c int_ceil_log2.c isqrt.c strtofr.c pow_z.c logging.c mulders.c get_f.c round_p.c erfc.c atan2.c subnormal.c const_catalan.c root.c gen_inverse.h sec.c csc.c cot.c eint.c sech.c csch.c coth.c round_near_x.c constant.c abort_prec_max.c stack_interface.c +libmpfr_la_SOURCES = mpfr.h mpf2mpfr.h mpfr-gmp.h mpfr-impl.h mpfr-longlong.h mpfr-test.h log_b2.h exceptions.c extract.c uceil_exp2.c uceil_log2.c ufloor_log2.c add.c add1.c add_ui.c agm.c clear.c cmp.c cmp_abs.c cmp_si.c cmp_ui.c comparisons.c div_2exp.c div_2si.c div_2ui.c div.c div_ui.c dump.c eq.c exp10.c exp2.c exp3.c exp.c frac.c get_d.c get_exp.c get_str.c init.c inp_str.c isinteger.c isinf.c isnan.c isnum.c const_log2.c log.c mul_2exp.c mul_2si.c mul_2ui.c mul.c mul_ui.c neg.c next.c out_str.c const_pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c round_prec.c set.c setmax.c setmin.c set_d.c set_dfl_prec.c set_exp.c set_rnd.c set_f.c set_prc_raw.c set_prec.c set_q.c set_si.c set_str.c set_str_raw.c set_ui.c set_z.c sqrt.c sqrt_ui.c sub.c sub1.c sub_ui.c rint.c ui_div.c ui_sub.c urandomb.c get_z_exp.c swap.c factorial.c cosh.c sinh.c tanh.c acosh.c asinh.c atanh.c atan.c cmp2.c exp_2.c asin.c const_euler.c cos.c sin.c tan.c fma.c hypot.c log1p.c expm1.c log2.c log10.c ui_pow.c ui_pow_ui.c minmax.c dim.c copysign.c gmp_op.c init2.c acos.c sin_cos.c set_nan.c set_inf.c powerof2.c gamma.c set_ld.c get_ld.c cbrt.c volatile.c fits_s.h fits_sshort.c fits_sint.c fits_slong.c fits_u.h fits_ushort.c fits_uint.c fits_ulong.c fits_uintmax.c fits_intmax.c get_si.c get_ui.c zeta.c cmp_d.c erf.c inits.c inits2.c clears.c sgn.c check.c sub1sp.c version.c mpn_exp.c mpfr-gmp.c mp_clz_tab.c sum.c add1sp.c free_cache.c si_op.c cmp_ld.c set_ui_2exp.c set_si_2exp.c set_uj.c set_sj.c get_sj.c get_uj.c get_z.c iszero.c cache.c sqr.c int_ceil_log2.c isqrt.c strtofr.c pow_z.c logging.c mulders.c get_f.c round_p.c erfc.c atan2.c subnormal.c const_catalan.c root.c gen_inverse.h sec.c csc.c cot.c eint.c sech.c csch.c coth.c round_near_x.c constant.c abort_prec_max.c stack_interface.c lngamma.c libmpfr_la_LIBADD = @LIBOBJS@ @@ -20,12 +20,13 @@ MA 02110-1301, USA. ############################################################################## -Changes from version 2.1.1 to version 2.2.0: +Changes from versions 2.1.* to version 2.2.0: - Bug fixes. - new functions mpfr_set_overflow, mpfr_set_underflow, mpfr_set_inexflag, mpfr_set_erangeflag, mpfr_set_nanflag, mpfr_erfc, mpfr_atan2, mpfr_pow_z, mpfr_subnormalize, mpfr_const_catalan, mpfr_sec, mpfr_csc, mpfr_cot, mpfr_root, mpfr_eint, mpfr_get_f, mpfr_sech, mpfr_csch, mpfr_coth, + mpfr_lngamma. - new macro: MPFR_VERSION_STRING - Remove the exported MPFR variables from mpfr.h to mpfr-impl.h. (They were undocumented, so programs which respect the API still work). diff --git a/README.dev b/README.dev index 943b26c3a..784a49438 100644 --- a/README.dev +++ b/README.dev @@ -46,17 +46,19 @@ To make a release (for the MPFR team): directly or with emacs), e.g. rcs2log > ChangeLog.2 and edit ChangeLog to insert ChangeLog.2 at the beginning and remove the duplicated lines. - 4) Fix the file permissions that have been broken by CVS. + 4) With fixperm, fix the file permissions that have been broken by CVS. 5) Generate the release version with "make dist". 6) Test the release version on different machines, with and without - the --disable-alloca configure option (or compile gmp with - --enable-alloca=debug), with and without -DXDEBUG in $CFLAGS, - with and without gmp internal files, with and without gmp built - as a shared library, with and without srcdir equals to objdir - (../mpfr/configure --srcdir=/users/spaces/pelissip/mpfr/), - with and without "--std=c99 -O3 -D_XOPEN_SOURCE=500", with and + the --disable-alloca configure option (or compile GMP with + --enable-alloca=debug and MPFR with --with-gmp-build to be able + to get the memory leak errors), with and without -DXDEBUG in + $CFLAGS, with and without gmp internal files, with and without + GMP built as a shared library, with and without srcdir equal to + objdir (../mpfr/configure --srcdir=/users/spaces/pelissip/mpfr/), + with and without "-std=c99 -O3 -D_XOPEN_SOURCE=500", with and without "-ansi" (which allows to turn off features that are - incompatible with ISO C90). + incompatible with ISO C90), and in various locales (LC_ALL=tr_TR + in particular, if installed). 7) If there is no problem, add a tag to the CVS corresponding to the release, e.g. mpfr-2-0-1-rel for mpfr 2.0.1 (note that the periods are replaced by dashes): cvs tag mpfr-2-0-1-rel @@ -199,6 +201,10 @@ use: MPFR_SAVE_EXPO_DECL (expo); mpfr_save_expo_t expo; TMP_DECL (maker); ; + ==================================================================== + +Do not use TMP_DECL / TMP_ALLOC, ... but MPFR_TMP_DECL, MPFR_TMP_ALLOC, ... + =========================================================================== If wou want to use the logging of MPFR, you need to enable it: @@ -260,8 +266,8 @@ Example: the same job as MPFR_LOG_BEGIN and MPFR_LOG_END but it is smatter since it intercepts the return itself to put the end statement. Example - MPFR_LOG_BEGIN (("op[%#R]=%R rnd=%d", op, op"), - ("x[%#R]=%R inexact=%d", x, x, i)); + MPFR_LOG_FUNC (("op[%#R]=%R rnd=%d", op, op"), + ("x[%#R]=%R inexact=%d", x, x, i)); The double brackets "((" and "))" are needed since MPFR must still @@ -110,7 +110,6 @@ New functions to implement: E1(-0) = -Inf DawsonIntegral Psi = LogDerivative - LogGamma GammaD(x) = Gamma(x+1/2) - functions defined in the LIA-2 standard + minimum and maximum (5.2.2): max, min, max_seq, min_seq, mmax_seq diff --git a/acinclude.m4 b/acinclude.m4 index 7592c7457..e8074e443 100644 --- a/acinclude.m4 +++ b/acinclude.m4 @@ -2,19 +2,19 @@ dnl MPFR specific autoconf macros dnl Copyright 2000, 2002, 2003, 2004, 2005 Free Software Foundation. dnl Contributed by the Spaces project, INRIA Lorraine. -dnl +dnl dnl This file is part of the MPFR Library. -dnl +dnl dnl The MPFR Library is free software; you can redistribute it and/or modify dnl it under the terms of the GNU Lesser General Public License as published dnl by the Free Software Foundation; either version 2.1 of the License, or (at dnl your option) any later version. -dnl +dnl dnl The MPFR Library is distributed in the hope that it will be useful, but dnl WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY dnl or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public dnl License for more details. -dnl +dnl dnl You should have received a copy of the GNU Lesser General Public License dnl along with the MPFR Library; see the file COPYING.LIB. If not, write to dnl the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, @@ -28,7 +28,7 @@ dnl to work as expected. AC_PREREQ(2.50) dnl ------------------------------------------------------------ -dnl You must put in MPFR_CONFIGS everything which configure MPFR +dnl You must put in MPFR_CONFIGS everything which configure MPFR dnl except: dnl -everything dealing with CC and CFLAGS in particular the ABI dnl but the IEEE-754 specific flags must be set here. @@ -47,8 +47,8 @@ AC_CHECK_HEADER([float.h],, AC_MSG_ERROR([float.h not found])) dnl Check for stdargs AC_CHECK_HEADER([stdarg.h],[AC_DEFINE([HAVE_STDARG],1,[Define if stdarg])], - [AC_CHECK_HEADER([varargs.h],, - AC_MSG_ERROR([stdarg.h or varargs.h not found]))]) + [AC_CHECK_HEADER([varargs.h],, + AC_MSG_ERROR([stdarg.h or varargs.h not found]))]) dnl sys/fpu.h - MIPS specific AC_CHECK_HEADERS([sys/time.h sys/fpu.h]) @@ -77,8 +77,8 @@ alpha*-*-*) fi esac -AC_CHECK_TYPE( [union fpc_csr], - AC_DEFINE(HAVE_FPC_CSR,1,[Define if union fpc_csr is available]), , +AC_CHECK_TYPE( [union fpc_csr], + AC_DEFINE(HAVE_FPC_CSR,1,[Define if union fpc_csr is available]), , [ #if HAVE_SYS_FPU_H # include <sys/fpu.h> @@ -156,6 +156,29 @@ if test "$mpfr_cv_have_denorms" = "yes"; then AC_DEFINE(HAVE_DENORMS,1,[Define if denormalized floats work.]) fi +dnl Check whether NAN != NAN (as required by the IEEE-754 standard, +dnl but not by the ISO C standard). For instance, this is false with +dnl MIPSpro 7.3.1.3m under IRIX64. By default, assume this is true. +AC_CACHE_CHECK([if NAN == NAN], mpfr_cv_nanisnan, [ +AC_TRY_RUN([ +#include <stdio.h> +#include <math.h> +#ifndef NAN +# define NAN (0.0/0.0) +#endif +int main() { + double d; + d = NAN; + return d != d; +} +], [mpfr_cv_nanisnan="yes"], + [mpfr_cv_nanisnan="no"], + [mpfr_cv_nanisnan="cannot test, assume no"]) +]) +if test "$mpfr_cv_nanisnan" = "yes"; then + AC_DEFINE(MPFR_NANISNAN,1,[Define if NAN == NAN.]) +fi + dnl Check if the chars '0' to '9' are consecutive values AC_MSG_CHECKING([if charset has consecutive values]) AC_RUN_IFELSE(AC_LANG_PROGRAM([[ @@ -192,7 +215,7 @@ int f (double (*func)(double)) { return 0;} a = f (round); return 0; ]])], [ - AC_MSG_RESULT(yes) + AC_MSG_RESULT(yes) AC_DEFINE(HAVE_ROUND, 1,[Have ISO-C99 round function]) ],[AC_MSG_RESULT(no)]) @@ -205,7 +228,7 @@ int f (double (*func)(double)) { return 0;} a = f(trunc); return 0; ]])], [ - AC_MSG_RESULT(yes) + AC_MSG_RESULT(yes) AC_DEFINE(HAVE_TRUNC, 1,[Have ISO-C99 trunc function]) ],[AC_MSG_RESULT(no)]) @@ -218,7 +241,7 @@ int f (double (*func)(double)) { return 0;} a = f(floor); return 0; ]])], [ - AC_MSG_RESULT(yes) + AC_MSG_RESULT(yes) AC_DEFINE(HAVE_FLOOR, 1,[Have ISO-C99 floor function]) ],[AC_MSG_RESULT(no)]) @@ -231,7 +254,7 @@ int f (double (*func)(double)) { return 0;} a = f(ceil); return 0; ]])], [ - AC_MSG_RESULT(yes) + AC_MSG_RESULT(yes) AC_DEFINE(HAVE_CEIL, 1,[Have ISO-C99 ceil function]) ],[AC_MSG_RESULT(no)]) @@ -244,7 +267,7 @@ int f (double (*func)(double)) { return 0;} a = f(nearbyint); return 0; ]])], [ - AC_MSG_RESULT(yes) + AC_MSG_RESULT(yes) AC_DEFINE(HAVE_NEARBYINT, 1,[Have ISO-C99 rint function]) ],[AC_MSG_RESULT(no)]) @@ -426,21 +449,21 @@ BEGIN { } if (got[23] == "300" && \ - got[22] == "031" && \ - got[21] == "326" && \ - got[20] == "363" && \ - got[19] == "105" && \ - got[18] == "100" && \ - got[17] == "000" && \ - got[16] == "000" && \ - got[15] == "000" && \ - got[14] == "000" && \ - got[13] == "000" && \ - got[12] == "000" && \ - got[11] == "000" && \ - got[10] == "000" && \ - got[9] == "000" && \ - got[8] == "000") + got[22] == "031" && \ + got[21] == "326" && \ + got[20] == "363" && \ + got[19] == "105" && \ + got[18] == "100" && \ + got[17] == "000" && \ + got[16] == "000" && \ + got[15] == "000" && \ + got[14] == "000" && \ + got[13] == "000" && \ + got[12] == "000" && \ + got[11] == "000" && \ + got[10] == "000" && \ + got[9] == "000" && \ + got[8] == "000") { print "IEEE quad, big endian" found = 1 @@ -489,7 +512,7 @@ case $mpfr_cv_c_long_double_format in ;; unknown* | "not available") ;; - *) + *) AC_MSG_WARN([oops, unrecognised float format: $mpfr_cv_c_long_double_format]) ;; esac @@ -32,13 +32,17 @@ int mpfr_add1sp2 (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mp_rnd_t); int mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) { mpfr_t tmpa, tmpb, tmpc; - int inexact, inexact2; + int inexb, inexc, inexact, inexact2; mpfr_init2 (tmpa, MPFR_PREC (a)); mpfr_init2 (tmpb, MPFR_PREC (b)); mpfr_init2 (tmpc, MPFR_PREC (c)); - MPFR_ASSERTN (mpfr_set (tmpb, b, GMP_RNDN) == 0); - MPFR_ASSERTN (mpfr_set (tmpc, c, GMP_RNDN) == 0); + + inexb = mpfr_set (tmpb, b, GMP_RNDN); + MPFR_ASSERTN (inexb == 0); + + inexc = mpfr_set (tmpc, c, GMP_RNDN); + MPFR_ASSERTN (inexc == 0); inexact2 = mpfr_add1 (tmpa, tmpb, tmpc, rnd_mode); inexact = mpfr_add1sp2 (a, b, c, rnd_mode); diff --git a/algorithms.tex b/algorithms.tex index 1a5d187ca..f4eb7bc0b 100644 --- a/algorithms.tex +++ b/algorithms.tex @@ -840,6 +840,64 @@ $u'=1 - \cos^2 x$.) Therefore, if $2^{e-1} \leq |s| < 2^e$, the absolute error on $s$ is bounded by $2^{-m} (\frac{5}{2} 2^{1-e}+1) \leq 2^{3-m-e}$. +\subsubsection{An asymptotically fast algorithm for sin and cos.} +We extend here the algorithm proposed by Brent for the exponential function +to the simultaneous computation of sin and cos. The idea is the following. +We first reduce the input $x$ to the range $0 < x < 1/2$. +Then we decompose $x$ as follows: +\[ x = \sum_{i=1}^{k} \frac{r_i}{2^{2^i}}, \] +where $r_i$ is an integer, $0 \leq r_i < 2^{2^{i-1}}$. + +We define $x_j = \sum_{i=j}^{k} \frac{r_i}{2^{2^i}}$; then $x = x_1$, +and we can write $x_j = \frac{r_j}{2^{2^j}} + x_{j+1}$. Thus with +$S_j := \sin \frac{r_j}{2^{2^j}}$ and $C_j := \cos \frac{r_j}{2^{2^j}}$: +\[ \sin x_j = S_j \cos x_{j+1} + C_j \sin x_{j+1}, \quad + \cos x_j = C_j \cos x_{j+1} - S_j \sin x_{j+1}. \] +The $2k$ values $S_j$ and $C_j$ can be computed by a binary splitting +algorithm, each one in $O(M(n) \log n)$. +Then each pair $(\sin x_j, \cos x_j)$ can be computed from +$(\sin x_{j+1}, \sin x_{j+1})$ with four multiplies and two additions or +subtractions. + +\paragraph{Error analysis.} +We use here Higham's method. We assume that the values of $S_j$ +and $C_j$ are approximated up to a multiplicative factor of the form +$(1+u)^3$, where $|u| \leq 2^{-p}$, $p \geq 4$ being the working precision. +We also assume that $\cos x_{j+1}$ and $\sin x_{j+1}$ are +approximated with a factor of the form $(1+u)^{k_j}$. +With rounding to nearest, the values of $S_j \cos x_{j+1}$, +$C_j \sin x_{j+1}$, $C_j \cos x_{j+1}$ and $S_j \sin x_{j+1}$ are thus +approximated with a factor $(1+u)^{k_j+4}$. +The value of $\sin x_j$ is approximated with a factor $(1+u)^{k_j+5}$ since +there all terms are nonnegative. + +We now analyze the effect of the cancellation in +$C_j \cos x_{j+1} - S_j \sin x_{j+1}$. +We have $\frac{r_j}{2^{2^j}} < 2^{-2^{j-1}}$, and for simplicity we define +$l := 2^{j-1}$; +thus $0 \leq S_j \leq 2^{-l}$, and $1-2^{-2l-1} \leq C_j \leq 1$. +Similarly we have $x_{j+1} < 2^{-2l}$, thus +$0 \leq \sin x_{j+1} \leq 2^{-2l}$, and $1-2^{-4l-1} \leq \cos x_{j+1} \leq 1$. +The error is multiplied by a maximal ratio of +\[ \frac{C_j \cos x_{j+1} + S_j \sin x_{j+1}} +{C_j \cos x_{j+1} - S_j \sin x_{j+1}} \leq +\frac{1+2^{-l} \cdot 2^{-2l}}{(1-2^{-2l-1})(1-2^{-4l-1})-2^{-l} \cdot 2^{-2l}}, +\] +which we can bound by +\[ \frac{1+2^{-3l}}{1-2^{-2l}} \leq \frac{1}{(1-2^{-2l})(1-2^{-3l})} +\leq \frac{1}{1-2^{-2l+1}}. \] +The product of all those factors for $j \geq 1$ is bounded by $3$ +(remember $l := 2^{j-1}$). + +In summary, the maximal error is of the form $3 [(1+u)^{5k}-1]$, where +$2^{2^{k-1}} < p \leq 2^{2^k}$. +For $p \geq 4$, $5k \cdot 2^{-p}$ is bounded by $5/16$, and +$(1+2^{-p})^{5k} - 1 \leq e^{5k \cdot 2^{-p}} - 1 \leq \frac{6}{5} +\cdot 5k \cdot 2^{-p} = 6k \cdot 2^{-p}$. +Thus the final relative error bound is $18k \cdot 2^{-p}$. +Since $k \leq 6$ for $p \leq 2^{64}$, this gives a uniform relative error +bound of $2^{-p+7}$. + \subsection{The tangent function} The tangent function is computed from the cosine, @@ -2445,7 +2503,7 @@ $(\bullet\bullet)$ ${\textnormal{error}}(w)$ -$w \leftarrow \o(\frac{u}{v}) $ +$w \leftarrow o(\frac{u}{v}) $ \end{minipage} & \begin{minipage}{7.5cm} diff --git a/configure.in b/configure.in index 42c730e78..173ec8caf 100644 --- a/configure.in +++ b/configure.in @@ -127,6 +127,7 @@ if test -n "$GMP_CFLAGS" && test -z "$user_redefine_cc" ; then AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[ #ifndef __GNUC__ # error "GCC Not Found" +error #endif ]])], [ GCC=yes @@ -141,6 +142,7 @@ AC_MSG_CHECKING(for ICC) AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[ #if !defined(__ICC) # error "ICC Not Found" +error #endif ]], [[]])],[ AC_MSG_RESULT(yes) @@ -212,6 +214,7 @@ case $host in #include "gmp.h" #if !__GMP_LIBGMP_DLL # error "Dead man" +error #endif ]], [[]])],[AC_MSG_RESULT(DLL)],[ AC_MSG_RESULT(static) @@ -221,6 +224,7 @@ case $host in #include "gmp.h" #if __GMP_LIBGMP_DLL # error "Dead man" +error #endif ]], [[]])],[AC_MSG_RESULT(static)],[ AC_MSG_RESULT(DLL) @@ -250,6 +254,7 @@ AC_COMPILE_IFELSE([AC_LANG_SOURCE([[ #include "gmp.h" #if (__GNU_MP_VERSION*100+__GNU_MP_VERSION_MINOR*10 < 410) # error "min GMP version is 4.1.0" +error #endif ]])],[AC_MSG_RESULT(yes)],[ AC_MSG_RESULT(no) @@ -300,7 +305,7 @@ LIBS="-lgmp $LIBS" else dnl Check if we can link with GMP -AC_CHECK_LIB(gmp, __gmp_version, [LIBS="-lgmp $LIBS"], +AC_CHECK_LIB(gmp, __gmpz_init, [LIBS="-lgmp $LIBS"], [AC_MSG_ERROR(libgmp not found)]) dnl Check if we can use mpn_sub_nc diff --git a/const_euler.c b/const_euler.c index d43314819..668cbcd74 100644 --- a/const_euler.c +++ b/const_euler.c @@ -58,7 +58,6 @@ mpfr_const_euler_internal (mpfr_t x, mp_rnd_t rnd) /* since prec >= 1, we have m >= 24 here, which ensures n >= 9 below */ n = 1 + (unsigned long) ((double) m * LOG2 / 2.0); MPFR_ASSERTD (n >= 9); - /* fprintf (stderr, "n=%u\n", n); */ mpfr_const_euler_S2 (y, n); /* error <= 3 ulps */ exp_S = MPFR_EXP(y); mpfr_set_ui (z, n, GMP_RNDN); @@ -100,9 +99,12 @@ mpfr_const_euler_S2_aux (mpz_t P, mpz_t Q, mpz_t T, unsigned long n, { if (a + 1 == b) { - mpz_set_si (P, (a == 1) ? n : - n * (a - 1)); + mpz_set_ui (P, n); + if (a > 1) + mpz_mul_si (P, P, 1 - (long) a); mpz_set (T, P); - mpz_set_ui (Q, a * a); + mpz_set_ui (Q, a); + mpz_mul_ui (Q, Q, a); } else { @@ -142,7 +144,8 @@ mpfr_const_euler_S2_aux (mpz_t P, mpz_t Q, mpz_t T, unsigned long n, } } -/* idem than mpfr_const_euler_S, using binary splitting. +/* computes S(n) = sum(n^k*(-1)^(k-1)/k!/k, k=1..ceil(4.319136566 * n)) + using binary splitting. We have S(n) = sum(f(k), k=1..N) with N=ceil(4.319136566 * n) and f(k) = n^k*(-1)*(k-1)/k!/k, thus f(k)/f(k-1) = -n*(k-1)/k^2 @@ -163,51 +166,6 @@ mpfr_const_euler_S2 (mpfr_t x, unsigned long n) mpz_clear (T); } -#if 0 -/* computes S(n) = sum(n^k*(-1)^(k-1)/k!/k, k=1..ceil(4.319136566 * n)) - with an error of at most ulp(x). - [S(n) >= 2 for n >= 5] - */ -static void -mpfr_const_euler_S (mpfr_t x, unsigned long n) -{ - unsigned long N, k, m; - mpz_t a, s, t; - - N = (long) (ALPHA * (double) n + 1.0); /* ceil(alpha * n) */ - - m = MPFR_PREC(x) + (unsigned long) ((double) n / LOG2) - + MPFR_INT_CEIL_LOG2 (N) + 1; - - mpz_init_set_ui (a, 1); - mpz_mul_2exp (a, a, m); /* a=-2^m where m is the precision of x */ - mpz_init_set_ui (s, 0); - mpz_init (t); - - /* here, a and s are exact */ - for (k = 1; k <= N; k++) - { - mpz_mul_ui (a, a, n); - mpz_div_ui (a, a, k); - mpz_div_ui (t, a, k); - if (k % 2) - mpz_add (s, s, t); - else - mpz_sub (s, s, t); - } - - /* the error on s is at most N (e^n + 1), - thus that the error on x is at most one ulp */ - - mpfr_set_z (x, s, GMP_RNDD); - mpfr_div_2ui (x, x, m, GMP_RNDD); - - mpz_clear (a); - mpz_clear (s); - mpz_clear (t); -} -#endif - /* computes R(n) = exp(-n)/n * sum(k!/(-n)^k, k=0..n-2) with error at most 4*ulp(x). Assumes n>=2. Since x <= exp(-n)/n <= 1/8, then 4*ulp(x) <= ulp(1). diff --git a/constant.c b/constant.c index 7f660b4ef..b5209654c 100644 --- a/constant.c +++ b/constant.c @@ -22,6 +22,9 @@ MA 02110-1301, USA. */ #include "mpfr-impl.h" static const mp_limb_t __gmpfr_limb1[1] = {MPFR_LIMB_HIGHBIT}; -const mpfr_t __gmpfr_one = {{2, MPFR_SIGN_POS, 1, (mp_limb_t*)__gmpfr_limb1}}; -const mpfr_t __gmpfr_two = {{2, MPFR_SIGN_POS, 2, (mp_limb_t*)__gmpfr_limb1}}; -const mpfr_t __gmpfr_four ={{2, MPFR_SIGN_POS, 3, (mp_limb_t*)__gmpfr_limb1}}; +const MPFR_THREAD_ATTR mpfr_t __gmpfr_one = + {{2, MPFR_SIGN_POS, 1, (mp_limb_t*)__gmpfr_limb1}}; +const MPFR_THREAD_ATTR mpfr_t __gmpfr_two = + {{2, MPFR_SIGN_POS, 2, (mp_limb_t*)__gmpfr_limb1}}; +const MPFR_THREAD_ATTR mpfr_t __gmpfr_four = + {{2, MPFR_SIGN_POS, 3, (mp_limb_t*)__gmpfr_limb1}}; @@ -28,7 +28,7 @@ MA 02110-1301, USA. */ */ #define FUNCTION mpfr_coth -#define INVERSE mpfr_tan +#define INVERSE mpfr_tanh #define ACTION_NAN(y) do { MPFR_SET_NAN(y); MPFR_RET_NAN; } while (1) #define ACTION_INF(y) return mpfr_set_si (y, MPFR_IS_POS(x) ? 1 : -1, GMP_RNDN) #define ACTION_ZERO(y,x) do { MPFR_SET_SAME_SIGN(y,x); MPFR_SET_ZERO(y); \ @@ -90,7 +90,8 @@ mpfr_exp (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) int signx = MPFR_SIGN (x); MPFR_SET_POS (y); - if (MPFR_IS_NEG_SIGN (signx) && rnd_mode == GMP_RNDD) + if (MPFR_IS_NEG_SIGN (signx) && (rnd_mode == GMP_RNDD || + rnd_mode == GMP_RNDZ)) { mpfr_setmax (y, 0); /* y = 1 - epsilon */ inexact = -1; diff --git a/fixperm b/fixperm new file mode 100755 index 000000000..ac8b823a0 --- /dev/null +++ b/fixperm @@ -0,0 +1,8 @@ +#!/bin/sh + +# Fix the permissions (this can't be done with CVS). + +for i in mpfr.texi *.c tests/*.c +do + chmod a-x "$i" +done @@ -2,7 +2,7 @@ Copyright 2001, 2002, 2003, 2004, 2005 Free Software Foundation. -This file is part of the MPFR Library, and was contributed by Mathieu Dutour. +This file is part of the MPFR Library. The MPFR 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 @@ -19,16 +19,44 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 51 Franklin Place, Fifth Floor, Boston, MA 02110-1301, USA. */ -/* The error analysis of gamma has been lost. - As a consequence, we can't change the algorithm... - But we may compute the exp(A-k) in the inner loop a lot faster - but less accurate, so it changes the precsion. - All MPFR tests still pass anyway */ -/* #define USE_PRECOMPUTED_EXP */ - #define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" +#define IS_GAMMA +#include "lngamma.c" +#undef IS_GAMMA + +/* return a sufficient precision such that 2-x is exact, assuming x < 0 */ +static mp_prec_t +mpfr_gamma_2_minus_x_exact (mpfr_srcptr x) +{ + /* Since x < 0, 2-x = 2+y with y := -x. + If y < 2, a precision w >= PREC(y) + EXP(2)-EXP(y) = PREC(y) + 2 - EXP(y) + is enough, since no overlap occurs in 2+y, so no carry happens. + If y >= 2, either ULP(y) <= 2, and we need w >= PREC(y)+1 since a + carry can occur, or ULP(y) > 2, and we need w >= EXP(y)-1: + (a) if EXP(y) <= 1, w = PREC(y) + 2 - EXP(y) + (b) if EXP(y) > 1 and EXP(y)-PREC(y) <= 1, w = PREC(y) + 1 + (c) if EXP(y) > 1 and EXP(y)-PREC(y) > 1, w = EXP(y) - 1 */ + return (MPFR_GET_EXP(x) <= 1) ? MPFR_PREC(x) + 2 - MPFR_GET_EXP(x) + : ((MPFR_GET_EXP(x) <= MPFR_PREC(x) + 1) ? MPFR_PREC(x) + 1 + : MPFR_GET_EXP(x) - 1); +} + +/* return a sufficient precision such that 1-x is exact, assuming x < 1 */ +static mp_prec_t +mpfr_gamma_1_minus_x_exact (mpfr_srcptr x) +{ + if (MPFR_IS_POS(x)) + return MPFR_PREC(x) - MPFR_GET_EXP(x); + else if (MPFR_GET_EXP(x) <= 0) + return MPFR_PREC(x) + 1 - MPFR_GET_EXP(x); + else if (MPFR_PREC(x) >= MPFR_GET_EXP(x)) + return MPFR_PREC(x) + 1; + else + return MPFR_GET_EXP(x); +} + /* We use the reflection formula Gamma(1+t) Gamma(1-t) = - Pi t / sin(Pi (1 + t)) in order to treat the case x <= 1, @@ -40,10 +68,7 @@ mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mp_rnd_t rnd_mode) { mpfr_t xp, GammaTrial, tmp, tmp2; mpz_t fact; -#ifdef USE_PRECOMPUTED_EXP - mpfr_t *tab; -#endif - mp_prec_t Prec, realprec, A, k; + mp_prec_t realprec; int compared, inex, is_integer; MPFR_GROUP_DECL (group); MPFR_SAVE_EXPO_DECL (expo); @@ -95,148 +120,172 @@ mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mp_rnd_t rnd_mode) if (compared == 0) return mpfr_set_ui (gamma, 1, rnd_mode); - /* if x is an integer that fits into an unsigned long, use mpfr_fac_ui */ + /* if x is an integer that fits into an unsigned long, use mpfr_fac_ui + if argument is not too large. + If precision is p, fac_ui costs O(u*p), whereas gamma costs O(p*M(p)), + so for u <= M(p), fac_ui should be faster. + We approximate here M(p) by p*log(p)^2, which is not a bad guess. + */ if (is_integer && mpfr_fits_ulong_p (x, GMP_RNDN)) { unsigned long int u; + mp_prec_t p = MPFR_PREC(gamma); u = mpfr_get_ui (x, GMP_RNDN); - return mpfr_fac_ui (gamma, u - 1, rnd_mode); + if (u / (MPFR_INT_CEIL_LOG2 (p) * MPFR_INT_CEIL_LOG2 (p)) <= p) + return mpfr_fac_ui (gamma, u - 1, rnd_mode); } + /* check for overflow: according to (6.1.37) in Abramowitz & Stegun, + gamma(x) >= exp(-x) * x^(x-1/2) * sqrt(2*Pi) + >= 2 * (x/e)^x / x for x >= 1 */ + if (compared > 0) + { + int overflow; + + mpfr_init2 (xp, 53); + mpfr_set_d (xp, EXPM1, GMP_RNDZ); /* 1/e rounded down */ + mpfr_mul (xp, x, xp, GMP_RNDZ); + mpfr_pow (xp, xp, x, GMP_RNDZ); + mpfr_clear_overflow (); + mpfr_mul_2exp (xp, xp, 1, GMP_RNDZ); + overflow = mpfr_overflow_p (); + mpfr_clear (xp); + return (overflow) ? mpfr_overflow (gamma, rnd_mode, 1) + : mpfr_gamma_aux (gamma, x, rnd_mode); + } + + /* now compared < 0 */ + MPFR_SAVE_EXPO_MARK (expo); + /* check for underflow: for x < 1, + gamma(x) = Pi*(x-1)/sin(Pi*(2-x))/gamma(2-x). + Since gamma(2-x) >= 2 * ((2-x)/e)^(2-x) / (2-x), we have + |gamma(x)| <= Pi*(1-x)*(2-x)/2/((2-x)/e)^(2-x) / |sin(Pi*(2-x))| + <= 12 * ((2-x)/e)^x / |sin(Pi*(2-x))|. + To avoid an underflow in ((2-x)/e)^x, we compute the logarithm. + */ + if (MPFR_IS_NEG(x)) + { + int underflow = 0, sgn, ck; + mp_prec_t w; + + mpfr_init2 (xp, 53); + mpfr_init2 (tmp, 53); + mpfr_init2 (tmp2, 53); + /* we want an upper bound for x * [log(2-x)-1]. + since x < 0, we need a lower bound on log(2-x) */ + mpfr_ui_sub (xp, 2, x, GMP_RNDD); + mpfr_log (xp, xp, GMP_RNDD); + mpfr_sub_ui (xp, xp, 1, GMP_RNDD); + mpfr_mul (xp, xp, x, GMP_RNDU); + + /* we need an upper bound on 1/|sin(Pi*(2-x))|, + thus a lower bound on |sin(Pi*(2-x))|. + If 2-x is exact, then the error of Pi*(2-x) is (1+u)^2 with u = 2^(-p) + thus the error on sin(Pi*(2-x)) is less than 1/2ulp + 3Pi(2-x)u, + assuming u <= 1, thus <= u + 3Pi(2-x)u */ + + w = mpfr_gamma_2_minus_x_exact (x); /* 2-x is exact for prec >= w */ + w += 17; /* to get tmp2 small enough */ + mpfr_set_prec (tmp, w); + mpfr_set_prec (tmp2, w); + ck = mpfr_ui_sub (tmp, 2, x, GMP_RNDN); + MPFR_ASSERTD (ck == 0); + mpfr_const_pi (tmp2, GMP_RNDN); + mpfr_mul (tmp2, tmp2, tmp, GMP_RNDN); /* Pi*(2-x) */ + mpfr_sin (tmp, tmp2, GMP_RNDN); /* sin(Pi*(2-x)) */ + mpfr_abs (tmp, tmp, GMP_RNDN); + mpfr_mul_ui (tmp2, tmp2, 3, GMP_RNDU); /* 3Pi(2-x) */ + mpfr_add_ui (tmp2, tmp2, 1, GMP_RNDU); /* 3Pi(2-x)+1 */ + mpfr_div_2exp (tmp2, tmp2, mpfr_get_prec (tmp), GMP_RNDU); + /* if tmp2<|tmp|, we get a lower bound */ + sgn = mpfr_sgn (tmp); + if (mpfr_cmp (tmp2, tmp) < 0) + { + mpfr_sub (tmp, tmp, tmp2, GMP_RNDZ); /* low bnd on |sin(Pi*(2-x))| */ + mpfr_ui_div (tmp, 12, tmp, GMP_RNDU); /* upper bound */ + mpfr_log (tmp, tmp, GMP_RNDU); + mpfr_add (tmp, tmp, xp, GMP_RNDU); + underflow = mpfr_cmp_si (xp, expo.saved_emin - 2) <= 0; + } + + mpfr_clear (xp); + mpfr_clear (tmp); + mpfr_clear (tmp2); + if (underflow) /* the sign is the opposite of that of sin(Pi*(2-x)) */ + { + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_underflow (gamma, (rnd_mode == GMP_RNDN) ? GMP_RNDZ : rnd_mode, -sgn); + } + } + realprec = MPFR_PREC (gamma); - realprec = realprec + MPFR_INT_CEIL_LOG2 (realprec) + 10; + /* we want both 1-x and 2-x to be exact */ + { + mp_prec_t w; + w = mpfr_gamma_1_minus_x_exact (x); + if (realprec < w) + realprec = w; + w = mpfr_gamma_2_minus_x_exact (x); + if (realprec < w) + realprec = w; + } + realprec = realprec + MPFR_INT_CEIL_LOG2 (realprec) + 20; + MPFR_ASSERTD(realprec >= 5); - MPFR_GROUP_INIT_4 (group, realprec + BITS_PER_MP_LIMB, + MPFR_GROUP_INIT_4 (group, realprec + MPFR_INT_CEIL_LOG2 (realprec) + 20, xp, tmp, tmp2, GammaTrial); mpz_init (fact); MPFR_ZIV_INIT (loop, realprec); for (;;) { - /* If compared < 0, we use the reflexion formula */ - /* Precision stuff */ - Prec = realprec + 2 * (compared < 0); - /* A = (prec_nec-0.5)*CST - CST = ln(2)/(ln(2*pi))) = 0.38 - This strange formula is just to avoid any overflow */ - A = (Prec/100)*38 + ((Prec%100)*38+100-38/2)/100 - 1; - /* Estimated_cancel is the amount of bit that will be flushed: - estimated_cancel = A + ecCST * A; - ecCST = {1+sup_{x\in [0,1]} x*ln((1-x)/x)}/ln(2) = 1.84 - This strange formula is just to avoid any overflow */ - Prec += 16 + (A + (A + (A/100)*84 + ((A%100)*84)/100)); - - /* FIXME: for x near 0, we want 1-x to be exact since the error - term does not seem to take into account the possible cancellation - here. Warning: if x < 0, we need one more bit. */ - if (MPFR_EXP(x) < 0 && Prec <= MPFR_PREC(x) - MPFR_EXP(x)) - Prec = MPFR_PREC(x) - MPFR_EXP(x) + 1; - - MPFR_GROUP_REPREC_4 (group, Prec, xp, tmp, tmp2, GammaTrial); - - if (compared < 0) - mpfr_sub (xp, __gmpfr_one, x, GMP_RNDN); - else - mpfr_sub (xp, x, __gmpfr_one, GMP_RNDN); - mpfr_set_ui (GammaTrial, 0, GMP_RNDN); - mpz_set_ui (fact, 1); - - /* It is faster to compute exp from k=1 to A, but - it changes the order of the sum, which change the error - analysis... Don't change it too much until we recompute - the error analysis. - Another trick is to compute factorial k in a MPFR rather than a MPZ - (Once again it changes the error analysis) */ -#ifdef USE_PRECOMPUTED_EXP - tab = (*__gmp_allocate_func) (sizeof (mpfr_t)*(A-1)); - mpfr_init2 (tab[0], Prec); - mpfr_exp (tab[0], __gmpfr_one, GMP_RNDN); - for (k = 1; k < A-1; k++) - { - mpfr_init2 (tab[k], Prec); - mpfr_mul (tab[k], tab[k-1], tab[0], GMP_RNDN); - } -#endif + mp_exp_t err_g; + int ck; + MPFR_GROUP_REPREC_4 (group, realprec, xp, tmp, tmp2, GammaTrial); - for (k = 1; k < A; k++) - { -#ifdef USE_PRECOMPUTED_EXP - mpfr_set (tmp2, tab[A-k-1], GMP_RNDN); -#else - mpfr_set_ui (tmp, A - k, GMP_RNDN); - mpfr_exp (tmp2, tmp, GMP_RNDN); -#endif - /* tmp2 = exp(A-k) */ - mpfr_ui_pow_ui (tmp, A - k, k - 1, GMP_RNDN); - /* tmp = (A-k)^(k-1) */ - mpfr_mul (tmp2, tmp2, tmp, GMP_RNDN); - /* tmp2 = (A-k)^(k-1) * exp(A-k) */ - mpfr_sqrt_ui (tmp, A - k, GMP_RNDN); - mpfr_mul (tmp2, tmp2, tmp, GMP_RNDN); - /* tmp2 = sqrt(A-k) * (A-k)^(k-1) * exp(A-k) */ - if (k >= 3) - { - /* mpfr_fac_ui (tmp, k-1, GMP_RNDN); */ - mpz_mul_ui (fact, fact, k - 1); - /* fact = (k-1)! */ - mpfr_set_z (tmp, fact, GMP_RNDN); - mpfr_div (tmp2, tmp2, tmp, GMP_RNDN); - } - /* tmp2 = sqrt(A-k) * (A-k)^(k-1) * exp(A-k) / (k-1)! */ - mpfr_add_ui (tmp, xp, k, GMP_RNDN); - mpfr_div (tmp2, tmp2, tmp, GMP_RNDN); - /* tmp2 = sqrt(A-k) * (A-k)^(k-1) * exp(A-k) / (k-1)! / (x+k) */ - if ((k & 1) == 0) - MPFR_CHANGE_SIGN (tmp2); - /* (-1)^(k-1)*sqrt(A-k) * (A-k)^(k-1) * exp(A-k) / (k-1)! / (x+k) */ - mpfr_add (GammaTrial, GammaTrial, tmp2, GMP_RNDN); - } -#ifdef DEBUG - printf("GammaTrial ="); - mpfr_out_str (stdout, 10, 0, GammaTrial, GMP_RNDD); - printf ("\n"); -#endif - mpfr_const_pi (tmp, GMP_RNDN); - mpfr_mul_2ui (tmp, tmp, 1, GMP_RNDN); /* 2*Pi */ - mpfr_sqrt (tmp, tmp, GMP_RNDN); /* sqrt(2*Pi) */ - mpfr_add (GammaTrial, GammaTrial, tmp, GMP_RNDN); - /* sqrt(2*Pi) + sum((-1)^(k-1)*sqrt(A-k)*(A-k)^(k-1)*exp(A-k) - /(k-1)!/(xp+k),k=1..A-1) */ - - mpfr_add_ui (tmp2, xp, A, GMP_RNDN); /* xp+A */ - mpfr_set_ui_2exp (tmp, 1, -1, GMP_RNDN); /* tmp= 1/2 */ - mpfr_add (tmp, tmp, xp, GMP_RNDN); /* xp+1/2 */ - mpfr_pow (tmp, tmp2, tmp, GMP_RNDN); /* (xp+A)^(xp+1/2) */ - mpfr_mul (GammaTrial, GammaTrial, tmp, GMP_RNDN); - /* (xp+A)^(xp+1/2) * [sqrt(2*Pi) + - sum((-1)^(k-1)*(A-k)^(k-1/2)*exp(A-k)/(k-1)!/(xp+k),k=1..A-1)] */ + /* reflection formula: gamma(x) = Pi*(x-1)/sin(Pi*(2-x))/gamma(2-x) */ - mpfr_neg (tmp, tmp2, GMP_RNDN); /* -(xp+A) */ - mpfr_exp (tmp, tmp, GMP_RNDN); /* exp(-xp-A) */ + ck = mpfr_ui_sub (xp, 2, x, GMP_RNDN); + MPFR_ASSERTD(ck == 0); /* 2-x, exact */ + mpfr_gamma (tmp, xp, GMP_RNDN); /* gamma(2-x), error (1+u) */ + mpfr_const_pi (tmp2, GMP_RNDN); /* Pi, error (1+u) */ + mpfr_mul (GammaTrial, tmp2, xp, GMP_RNDN); /* Pi*(2-x), error (1+u)^2 */ + err_g = MPFR_GET_EXP(GammaTrial); + mpfr_sin (GammaTrial, GammaTrial, GMP_RNDN); /* sin(Pi*(2-x)) */ + err_g = err_g + 1 - MPFR_GET_EXP(GammaTrial); + /* let g0 the true value of Pi*(2-x), g the computed value. + We have g = g0 + h with |h| <= |(1+u^2)-1|*g. + Thus sin(g) = sin(g0) + h' with |h'| <= |(1+u^2)-1|*g. + The relative error is thus bounded by |(1+u^2)-1|*g/sin(g) + <= |(1+u^2)-1|*2^err_g. <= 2.25*u*2^err_g for |u|<=1/4. + With the rounding error, this gives (0.5 + 2.25*2^err_g)*u. */ + ck = mpfr_sub_ui (xp, x, 1, GMP_RNDN); + MPFR_ASSERTD(ck == 0); /* x-1, exact */ + mpfr_mul (xp, tmp2, xp, GMP_RNDN); /* Pi*(x-1), error (1+u)^2 */ mpfr_mul (GammaTrial, GammaTrial, tmp, GMP_RNDN); + /* [1 + (0.5 + 2.25*2^err_g)*u]*(1+u)^2 = 1 + (2.5 + 2.25*2^err_g)*u + + (0.5 + 2.25*2^err_g)*u*(2u+u^2) + u^2. + For err_g <= realprec-2, we have (0.5 + 2.25*2^err_g)*u <= + 0.5*u + 2.25/4 <= 0.6875 and u^2 <= u/4, thus + (0.5 + 2.25*2^err_g)*u*(2u+u^2) + u^2 <= 0.6875*(2u+u/4) + u/4 + <= 1.8*u, thus the rel. error is bounded by (4.5 + 2.25*2^err_g)*u. */ + mpfr_div (GammaTrial, xp, GammaTrial, GMP_RNDN); + /* the error is of the form (1+u)^3/[1 + (4.5 + 2.25*2^err_g)*u]. + For realprec >= 5 and err_g <= realprec-2, [(4.5 + 2.25*2^err_g)*u]^2 + <= 0.71, and for |y|<=0.71, 1/(1-y) can be written 1+a*y with a<=4. + (1+u)^3 * (1+4*(4.5 + 2.25*2^err_g)*u) + = 1 + (21 + 9*2^err_g)*u + (57+27*2^err_g)*u^2 + (55+27*2^err_g)*u^3 + + (18+9*2^err_g)*u^4 + <= 1 + (21 + 9*2^err_g)*u + (57+27*2^err_g)*u^2 + (56+28*2^err_g)*u^3 + <= 1 + (21 + 9*2^err_g)*u + (59+28*2^err_g)*u^2 + <= 1 + (23 + 10*2^err_g)*u. + The final error is thus bounded by (23 + 10*2^err_g) ulps, + which is <= 2^6 for err_g<=2, and <= 2^(err_g+4) for err_g >= 2. */ + err_g = (err_g <= 2) ? 6 : err_g + 4; - if (compared < 0) - { - mpfr_const_pi (tmp2, GMP_RNDN); - mpfr_sub (tmp, x, __gmpfr_one, GMP_RNDN); - mpfr_mul (tmp, tmp2, tmp, GMP_RNDN); - mpfr_div (GammaTrial, tmp, GammaTrial, GMP_RNDN); - mpfr_sin (tmp, tmp, GMP_RNDN); - mpfr_div (GammaTrial, GammaTrial, tmp, GMP_RNDN); - } -#ifdef DEBUG - printf("GammaTrial ="); - mpfr_out_str (stdout, 10, 0, GammaTrial, GMP_RNDD); - printf ("\n"); -#endif -#ifdef USE_PRECOMPUTED_EXP - for (k = 0 ; k < A-1 ; k++) - mpfr_clear (tab[k]); - (*__gmp_free_func) (tab, sizeof (mpfr_t) * (A-1) ); -#endif - if (mpfr_can_round (GammaTrial, realprec, GMP_RNDD, GMP_RNDZ, - MPFR_PREC(gamma) + (rnd_mode == GMP_RNDN))) + if (MPFR_LIKELY (MPFR_CAN_ROUND (GammaTrial, realprec - err_g, + MPFR_PREC(gamma), rnd_mode))) break; MPFR_ZIV_NEXT (loop, realprec); } @@ -245,6 +294,7 @@ mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mp_rnd_t rnd_mode) inex = mpfr_set (gamma, GammaTrial, rnd_mode); MPFR_GROUP_CLEAR (group); mpz_clear (fact); + MPFR_SAVE_EXPO_FREE (expo); - return mpfr_check_range(gamma, inex, rnd_mode); + return mpfr_check_range (gamma, inex, rnd_mode); } diff --git a/int_ceil_log2.c b/int_ceil_log2.c index 3155e4c03..e4d9dc80b 100644 --- a/int_ceil_log2.c +++ b/int_ceil_log2.c @@ -28,7 +28,7 @@ __gmpfr_int_ceil_log2 (unsigned long n) int b; mp_limb_t limb = n; - MPFR_ASSERTD (limb == n); + MPFR_ASSERTN (limb == n); count_leading_zeros (b, limb); return BITS_PER_MP_LIMB-b; diff --git a/lngamma.c b/lngamma.c new file mode 100644 index 000000000..fcef9bef2 --- /dev/null +++ b/lngamma.c @@ -0,0 +1,458 @@ +/* mpfr_lngamma -- lngamma function + +Copyright 2005 Free Software Foundation. + +This file is part of the MPFR Library. + +The MPFR 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 MPFR 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 MPFR 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. */ + +#define MPFR_NEED_LONGLONG_H +#include "mpfr-impl.h" + +/* assuming b[0]...b[2(n-1)] are computed, computes and stores B[2n]*(2n+1)! + + t/(exp(t)-1) = sum(B[j]*t^j/j!, j=0..infinity) + thus t = (exp(t)-1) * sum(B[j]*t^j/j!, n=0..infinity). + Taking the coefficient of degree n+1 > 1, we get: + 0 = sum(1/(n+1-k)!*B[k]/k!, k=0..n) + which gives: + B[n] = -sum(binomial(n+1,k)*B[k], k=0..n-1)/(n+1). + + Let C[n] = B[n]*(n+1)!. + Then C[n] = -sum(binomial(n+1,k)*C[k]*n!/(k+1)!, k=0..n-1), + which proves that the C[n] are integers. +*/ +static mpz_t* +bernoulli (mpz_t *b, unsigned long n) +{ + if (n == 0) + { + b = (mpz_t *) (*__gmp_allocate_func) (sizeof (mpz_t)); + mpz_init_set_ui (b[0], 1); + } + else + { + mpz_t t; + unsigned long k; + + b = (mpz_t *) (*__gmp_reallocate_func) + (b, n * sizeof (mpz_t), (n + 1) * sizeof (mpz_t)); + mpz_init (b[n]); + /* b[n] = -sum(binomial(2n+1,2k)*C[k]*(2n)!/(2k+1)!, k=0..n-1) */ + mpz_init_set_ui (t, 2 * n + 1); + mpz_mul_ui (t, t, 2 * n - 1); + mpz_mul_ui (t, t, 2 * n); + mpz_mul_ui (t, t, n); + mpz_div_ui (t, t, 3); /* exact: t=binomial(2*n+1,2*k)*(2*n)!/(2*k+1)! + for k=n-1 */ + mpz_mul (b[n], t, b[n-1]); + for (k = n - 1; k-- > 0;) + { + mpz_mul_ui (t, t, 2 * k + 1); + mpz_mul_ui (t, t, 2 * k + 2); + mpz_mul_ui (t, t, 2 * k + 2); + mpz_mul_ui (t, t, 2 * k + 3); + mpz_div_ui (t, t, 2 * (n - k) + 1); + mpz_div_ui (t, t, 2 * (n - k)); + mpz_addmul (b[n], t, b[k]); + } + /* take into account C[1] */ + mpz_mul_ui (t, t, 2 * n + 1); + mpz_div_2exp (t, t, 1); + mpz_sub (b[n], b[n], t); + mpz_neg (b[n], b[n]); + mpz_clear (t); + } + return b; +} + +/* given a precision p, return alpha, such that the argument reduction + will use k = alpha*p*log(2). + + Warning: we should always have alpha >= log(2)/(2Pi) ~ 0.11, + and the smallest value of alpha multiplied by the smallest working + precision should be >= 4. +*/ +static double +mpfr_gamma_alpha (mp_prec_t p) +{ + if (p <= 100) + return 0.6; + else if (p <= 200) + return 0.8; + else if (p <= 500) + return 0.8; + else if (p <= 1000) + return 1.3; + else if (p <= 2000) + return 1.7; + else if (p <= 5000) + return 2.2; + else if (p <= 10000) + return 3.4; + else /* heuristic fit from above */ + return 0.26 * (double) MPFR_INT_CEIL_LOG2 ((unsigned long) p); +} + +/* lngamma(x) = log(gamma(x)). + We use formula [6.1.40] from Abramowitz&Stegun: + lngamma(z) = (z-1/2)*log(z) - z + 1/2*log(2*Pi) + + sum (Bernoulli[2n]/(2m)/(2m-1)/z^(2m-1),m=1..infinity) + According to [6.1.42], if the sum is truncated after m=n, the error + R_n(z) is bounded by |B[2n+2]|*K(z)/(2n+1)/(2n+2)/|z|^(2n+1) + where K(z) = max (z^2/(u^2+z^2)) for u >= 0. + For z real, |K(z)| <= 1 thus R_n(z) is bounded by the first neglected term. + */ +#ifdef IS_GAMMA +static int +#define GAMMA_FUNC mpfr_gamma_aux +#else +int +#define GAMMA_FUNC mpfr_lngamma +#endif +GAMMA_FUNC (mpfr_ptr y, mpfr_srcptr z0, mp_rnd_t rnd) +{ + mp_prec_t precy, w; /* working precision */ + mpfr_t s, t, u, v, z; + unsigned long m, k, maxm; + mpz_t *B; + int inexact, compared; + mp_exp_t err_s, err_t; + unsigned long Bm = 0; /* number of allocated B[] */ + unsigned long oldBm; + double d; + MPFR_SAVE_EXPO_DECL (expo); + +#ifndef IS_GAMMA + /* special cases */ + if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (z0))) + { + if (MPFR_IS_NAN (z0) || MPFR_IS_NEG (z0)) + { + MPFR_SET_NAN (y); + MPFR_RET_NAN; + } + else /* lngamma(+Inf) = lngamma(+0) = +Inf */ + { + MPFR_SET_INF (y); + MPFR_SET_POS (y); + MPFR_RET (0); /* exact */ + } + } + + /* if x < 0 and -2k-1 <= x <= -2k, then lngamma(x) = NaN */ + if (MPFR_IS_NEG(z0) && ((mpfr_get_si (z0, GMP_RNDZ) % 2) == 0 + || mpfr_integer_p (z0))) + { + MPFR_SET_NAN (y); + MPFR_RET_NAN; + } +#endif + + precy = MPFR_PREC(y); + + compared = mpfr_cmp_ui (z0, 1); + +#ifndef IS_GAMMA + if (compared == 0) /* lngamma(1) = +0 */ + return mpfr_set_ui (y, 0, GMP_RNDN); +#endif + + mpfr_init2 (s, MPFR_PREC_MIN); + mpfr_init2 (t, MPFR_PREC_MIN); + mpfr_init2 (u, MPFR_PREC_MIN); + mpfr_init2 (v, MPFR_PREC_MIN); + mpfr_init2 (z, MPFR_PREC_MIN); + + MPFR_SAVE_EXPO_MARK (expo); + + if (compared < 0) + { + mp_exp_t err_u; + + /* use reflection formula: + gamma(x) = Pi*(x-1)/sin(Pi*(2-x))/gamma(2-x) + thus lngamma(x) = log(Pi*(x-1)/sin(Pi*(2-x))) - lngamma(2-x) */ + + w = precy + MPFR_INT_CEIL_LOG2 (precy); + while (1) + { + w += MPFR_INT_CEIL_LOG2 (w) + 13; + MPFR_ASSERTD(w >= 3); + mpfr_set_prec (s, w); + mpfr_set_prec (t, w); + mpfr_set_prec (u, w); + mpfr_set_prec (v, w); + mpfr_ui_sub (s, 2, z0, GMP_RNDD); /* s = (2-z0) * (1+2u) */ + mpfr_const_pi (t, GMP_RNDN); /* t = Pi * (1+u) */ + mpfr_lngamma (u, s, GMP_RNDN); /* lngamma(2-x) */ + /* Let s = (2-z0) + h. By construction, -(2-z0)*(2u) <= h <= 0. + We have lngamma(s) = lngamma(2-z0) + h*Psi(z), z in [2-z0+h,2-z0]. + Since 2-z0+h >= 1 and |Psi(x)| <= max(1,log(x)) for x >= 1, + the error on u is bounded by + ulp(u)/2 + (2-z0)*max(1,log(2-z0))*2^(1-w). */ + d = (double) MPFR_GET_EXP(s) * 0.694; /* upper bound for log(2-z0) */ + err_u = MPFR_GET_EXP(s) + __gmpfr_ceil_log2 (d) + 1 - MPFR_GET_EXP(u); + err_u = (err_u >= 0) ? err_u + 1 : 0; + /* now the error on u is bounded by 2^err_u ulps */ + + mpfr_mul (s, s, t, GMP_RNDN); /* Pi*(2-x), (1+u)^4 */ + err_s = MPFR_GET_EXP(s); /* 2-x <= 2^err_s */ + mpfr_sin (s, s, GMP_RNDN); /* sin(Pi*(2-x)) */ + /* the error on s is bounded by 1/2*ulp(s) + [(1+u)^4-1]*(2-x) + <= 1/2*ulp(s) + 5*2^(-w)*(2-x) for w >= 3 */ + err_s += 3 - MPFR_GET_EXP(s); + err_s = (err_s >= 0) ? err_s + 1 : 0; + /* the error on s is bounded by 2^err_s ulps, thus the relative + error is bounded by 2^(err_s+1) */ + err_s ++; /* relative error */ + + mpfr_sub_ui (v, z0, 1, GMP_RNDN); /* v = (x-1) * (1+u) */ + mpfr_mul (v, v, t, GMP_RNDN); /* v = Pi*(x-1) * (1+u)^3 */ + mpfr_div (v, v, s, GMP_RNDN); /* Pi*(x-1)/sin(Pi*(2-x)) */ + /* (1+u)^(4+2^err_s+1) */ + err_s = (err_s <= 2) ? 3 + (err_s / 2) : err_s + 1; + MPFR_ASSERTD(MPFR_IS_POS(v)); + mpfr_log (v, v, GMP_RNDN); + /* log(v*(1+e)) = log(v)+log(1+e) where |e| <= 2^(err_s-w). + Since |log(1+e)| <= 2*e for |e| <= 1/4, the error on v is + bounded by ulp(v)/2 + 2^(err_s+1-w). */ + if (err_s + 2 > w) + { + w += err_s + 2; + } + else + { + err_s += 1 - MPFR_GET_EXP(v); + err_s = (err_s >= 0) ? err_s + 1 : 0; + /* the error on v is bounded by 2^err_s ulps */ + err_u += MPFR_GET_EXP(u); /* absolute error on u */ + err_s += MPFR_GET_EXP(v); /* absolute error on v */ + mpfr_sub (s, v, u, GMP_RNDN); + /* the total error on s is bounded by ulp(s)/2 + 2^(err_u-w) + + 2^(err_s-w) <= ulp(s)/2 + 2^(max(err_u,err_s)+1-w) */ + err_s = (err_s >= err_u) ? err_s : err_u; + err_s += 1 - MPFR_GET_EXP(s); /* error is 2^err_s ulp(s) */ + err_s = (err_s >= 0) ? err_s + 1 : 0; + if (mpfr_can_round (s, w - err_s, GMP_RNDN, GMP_RNDZ, precy + + (rnd == GMP_RNDN))) + goto end; + } + } + } + + /* now z0 > 1 */ + + MPFR_ASSERTD (compared > 0); + + /* since k is O(w), the value of log(z0*...*(z0+k-1)) is about w*log(w), + so there is a cancellation of ~log(w) in the argument reconstruction */ + w = precy + MPFR_INT_CEIL_LOG2 (precy); + + do + { + w += MPFR_INT_CEIL_LOG2 (w) + 13; + MPFR_ASSERTD (w >= 3); + + mpfr_set_prec (s, 53); + /* we need z >= w*log(2)/(2*Pi) to get an absolute error less than 2^(-w) + but the optimal value is about 0.155665*w */ + mpfr_set_d (s, mpfr_gamma_alpha (precy) * (double) w, GMP_RNDU); + if (mpfr_cmp (z0, s) < 0) + { + mpfr_sub (s, s, z0, GMP_RNDU); + k = mpfr_get_ui (s, GMP_RNDU); + if (k < 3) + k = 3; + } + else + k = 3; + + mpfr_set_prec (s, w); + mpfr_set_prec (t, w); + mpfr_set_prec (u, w); + mpfr_set_prec (v, w); + mpfr_set_prec (z, w); + + mpfr_add_ui (z, z0, k, GMP_RNDN); + /* z = (z0+k)*(1+t1) with |t1| <= 2^(-w) */ + + /* z >= 4 ensures the relative error on log(z) is small, + and also (z-1/2)*log(z)-z >= 0 */ + MPFR_ASSERTD (mpfr_cmp_ui (z, 4) >= 0); + + mpfr_log (s, z, GMP_RNDN); /* log(z) */ + /* we have s = log((z0+k)*(1+t1))*(1+t2) with |t1|, |t2| <= 2^(-w). + Since w >= 2 and z0+k >= 4, we can write log((z0+k)*(1+t1)) + = log(z0+k) * (1+t3) with |t3| <= 2^(-w), thus we have + s = log(z0+k) * (1+t4)^2 with |t4| <= 2^(-w) */ + mpfr_mul_2exp (t, z, 1, GMP_RNDN); /* t = 2z * (1+t5) */ + mpfr_sub_ui (t, t, 1, GMP_RNDN); /* t = 2z-1 * (1+t6)^3 */ + /* since we can write 2z*(1+t5) = (2z-1)*(1+t5') with + t5' = 2z/(2z-1) * t5, thus |t5'| <= 8/7 * t5 */ + mpfr_mul (s, s, t, GMP_RNDN); /* (2z-1)*log(z) * (1+t7)^6 */ + mpfr_div_2exp (s, s, 1, GMP_RNDN); /* (z-1/2)*log(z) * (1+t7)^6 */ + mpfr_sub (s, s, z, GMP_RNDN); /* (z-1/2)*log(z)-z */ + /* s = [(z-1/2)*log(z)-z]*(1+u)^14, s >= 1/2 */ + + mpfr_ui_div (u, 1, z, GMP_RNDN); /* 1/z * (1+u), u <= 1/4 since z >= 4 */ + + /* the first term is B[2]/2/z = 1/12/z: t=1/12/z, C[2]=1 */ + mpfr_div_ui (t, u, 12, GMP_RNDN); /* 1/(12z) * (1+u)^2, t <= 3/128 */ + mpfr_set (v, t, GMP_RNDN); /* (1+u)^2, v < 2^(-5) */ + mpfr_add (s, s, v, GMP_RNDN); /* (1+u)^15 */ + + mpfr_mul (u, u, u, GMP_RNDN); /* 1/z^2 * (1+u)^3 */ + + if (Bm == 0) + { + B = bernoulli ((mpz_t *) 0, 0); + B = bernoulli (B, 1); + Bm = 2; + } + + /* m <= maxm ensures that 2*m*(2*m+1) <= ULONG_MAX */ + maxm = 1UL << (BITS_PER_MP_LIMB / 2 - 1); + + /* s:(1+u)^15, t:(1+u)^2, t <= 3/128 */ + + for (m = 2; MPFR_GET_EXP(v) + (mp_exp_t) w >= MPFR_GET_EXP(s); m++) + { + mpfr_mul (t, t, u, GMP_RNDN); /* (1+u)^(10m-14) */ + if (m <= maxm) + { + mpfr_mul_ui (t, t, 2*(m-1)*(2*m-3), GMP_RNDN); + mpfr_div_ui (t, t, 2*m*(2*m-1), GMP_RNDN); + mpfr_div_ui (t, t, 2*m*(2*m+1), GMP_RNDN); + } + else + { + mpfr_mul_ui (t, t, 2*(m-1), GMP_RNDN); + mpfr_mul_ui (t, t, 2*m-3, GMP_RNDN); + mpfr_div_ui (t, t, 2*m, GMP_RNDN); + mpfr_div_ui (t, t, 2*m-1, GMP_RNDN); + mpfr_div_ui (t, t, 2*m, GMP_RNDN); + mpfr_div_ui (t, t, 2*m+1, GMP_RNDN); + } + /* (1+u)^(10m-8) */ + /* invariant: t=1/(2m)/(2m-1)/z^(2m-1)/(2m+1)! */ + if (Bm <= m) + { + B = bernoulli (B, m); /* B[2m]*(2m+1)!, exact */ + Bm ++; + } + mpfr_mul_z (v, t, B[m], GMP_RNDN); /* (1+u)^(10m-7) */ + MPFR_ASSERTD(MPFR_GET_EXP(v) <= - (2 * m + 3)); + mpfr_add (s, s, v, GMP_RNDN); + } + /* m <= 1/2*Pi*e*z ensures that |v[m]| < 1/2^(2m+3) */ + MPFR_ASSERTD ((double) m <= 4.26 * mpfr_get_d (z, GMP_RNDZ)); + + /* We have sum([(1+u)^(10m-7)-1]*1/2^(2m+3), m=2..infinity) + <= 1.46*u for u <= 2^(-3). + We have 0 < lngamma(z) - [(z - 1/2) ln(z) - z + 1/2 ln(2 Pi)] < 0.021 + for z >= 4, thus since the initial s >= 0.85, the different values of + s differ by at most one binade, and the total rounding error on s + in the for-loop is bounded by 2*(m-1)*ulp(final_s). + The error coming from the v's is bounded by + 1.46*2^(-w) <= 2*ulp(final_s). + Thus the total error so far is bounded by [(1+u)^15-1]*s+2m*ulp(s) + <= (2m+47)*ulp(s). + Taking into account the truncation error (which is bounded by the last + term v[] according to 6.1.42 in A&S), the bound is (2m+48)*ulp(s). + */ + + /* add 1/2*log(2*Pi) and subtract log(z0*(z0+1)*...*(z0+k-1)) */ + mpfr_const_pi (v, GMP_RNDN); /* v = Pi*(1+u) */ + mpfr_mul_2exp (v, v, 1, GMP_RNDN); /* v = 2*Pi * (1+u) */ + if (k) + { + unsigned long l; + mpfr_set (t, z0, GMP_RNDN); /* t = z0*(1+u) */ + for (l = 1; l < k; l++) + { + mpfr_add_ui (u, z0, l, GMP_RNDN); /* u = (z0+l)*(1+u) */ + mpfr_mul (t, t, u, GMP_RNDN); /* (1+u)^(2l+1) */ + } + /* now t: (1+u)^(2k-1) */ + /* instead of computing log(sqrt(2*Pi)/t), we compute + 1/2*log(2*Pi/t^2), which trades a square root for a square */ + mpfr_mul (t, t, t, GMP_RNDN); /* (z0*...*(z0+k-1))^2, (1+u)^(4k-1) */ + mpfr_div (v, v, t, GMP_RNDN); + /* 2*Pi/(z0*...*(z0+k-1))^2 (1+u)^(4k+1) */ + } +#ifdef IS_GAMMA + err_s = MPFR_GET_EXP(s); + mpfr_exp (s, s, GMP_RNDN); + /* before the exponential, we have s = s0 + h where + |h| <= (2m+48)*ulp(s), thus exp(s0) = exp(s) * exp(-h). + For |h| <= 1/4, we have |exp(h)-1| <= 1.2*|h| thus + |exp(s) - exp(s0)| <= 1.2 * exp(s) * (2m+48)* 2^(EXP(s)-w). */ + d = 1.2 * (2.0 * (double) m + 48.0); + /* the error on s is bounded by d*2^err_s * 2^(-w) */ + mpfr_sqrt (t, v, GMP_RNDN); + /* let v0 be the exact value of v. We have v = v0*(1+u)^(4k+1), + thus t = sqrt(v0)*(1+u)^(2k+3/2). */ + mpfr_mul (s, s, t, GMP_RNDN); + /* the error on input s is bounded by (1+u)^(d*2^err_s), + and that on t is (1+u)^(2k+3/2), thus the + total error is (1+u)^(d*2^err_s+2k+5/2) */ + err_s += __gmpfr_ceil_log2 (d); + err_t = __gmpfr_ceil_log2 (2.0 * (double) k + 2.5); + err_s = (err_s >= err_t) ? err_s + 1 : err_t + 1; +#else + mpfr_log (t, v, GMP_RNDN); + /* let v0 be the exact value of v. We have v = v0*(1+u)^(4k+1), + thus log(v) = log(v0) + (4k+1)*log(1+u). Since |log(1+u)/u| <= 1.07 + for |u| <= 2^(-3), the absolute error on log(v) is bounded by + 1.07*(4k+1)*u, and the rounding error by ulp(t). */ + mpfr_div_2exp (t, t, 1, GMP_RNDN); + /* the error on t is now bounded by ulp(t) + 0.54*(4k+1)*2^(-w). + We have sqrt(2*Pi)/(z0*(z0+1)*...*(z0+k-1)) <= sqrt(2*Pi)/k! <= 0.5 + since k>=3, thus t <= -0.5 and ulp(t) >= 2^(-w). + Thus the error on t is bounded by (2.16*k+1.54)*ulp(t). */ + err_t = MPFR_GET_EXP(t) + (mp_exp_t) + __gmpfr_ceil_log2 (2.2 * (double) k + 1.6); + err_s = MPFR_GET_EXP(s) + (mp_exp_t) + __gmpfr_ceil_log2 (2.0 * (double) m + 48.0); + mpfr_add (s, s, t, GMP_RNDN); /* this is a subtraction in fact */ + /* the final error in ulp(s) is + <= 1 + 2^(err_t-EXP(s)) + 2^(err_s-EXP(s)) + <= 2^(1+max(err_t,err_s)-EXP(s)) if err_t <> err_s + <= 2^(2+max(err_t,err_s)-EXP(s)) if err_t = err_s */ + err_s = (err_t == err_s) ? 1 + err_s : ((err_t > err_s) ? err_t : err_s); + err_s += 1 - MPFR_GET_EXP(s); +#endif + } + while (MPFR_UNLIKELY (!MPFR_CAN_ROUND (s, w - err_s, precy, rnd))); + + oldBm = Bm; + while (Bm--) + mpz_clear (B[Bm]); + (*__gmp_free_func) (B, oldBm * sizeof (mpz_t)); + + end: + inexact = mpfr_set (y, s, rnd); + + mpfr_clear (s); + mpfr_clear (t); + mpfr_clear (u); + mpfr_clear (v); + mpfr_clear (z); + + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_check_range (y, inexact, rnd); +} @@ -44,6 +44,7 @@ mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) mp_prec_t p, q; mpfr_t tmp1, tmp2; mp_limb_t *tmp1p, *tmp2p; + MPFR_SAVE_EXPO_DECL (expo); MPFR_ZIV_DECL (loop); MPFR_TMP_DECL(marker); @@ -107,6 +108,7 @@ mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) p += BITS_PER_MP_LIMB - (p%BITS_PER_MP_LIMB); */ MPFR_TMP_MARK(marker); + MPFR_SAVE_EXPO_MARK (expo); MPFR_ZIV_INIT (loop, p); for (;;) @@ -153,5 +155,6 @@ mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) /* We clean */ MPFR_TMP_FREE(marker); - return inexact; /* result is inexact */ + MPFR_SAVE_EXPO_FREE (expo); + return mpfr_check_range (r, inexact, rnd_mode); } diff --git a/mpfr-gmp.h b/mpfr-gmp.h index f68c5d5f7..14a9324c5 100644 --- a/mpfr-gmp.h +++ b/mpfr-gmp.h @@ -94,7 +94,19 @@ extern "C" { #define MPN_ZERO(dst, n) memset((dst), 0, (n)*BYTES_PER_MP_LIMB) #define MPN_COPY_DECR(dst,src,n) memmove((dst),(src),(n)*BYTES_PER_MP_LIMB) #define MPN_COPY_INCR(dst,src,n) memmove((dst),(src),(n)*BYTES_PER_MP_LIMB) -#define MPN_COPY(dst,src,n) memcpy((dst),(src),(n)*BYTES_PER_MP_LIMB) +#define MPN_COPY(dst,src,n) \ + do \ + { \ + if ((dst) != (src)) \ + { \ + MPFR_ASSERTD ((char *) (dst) >= (char *) (src) + \ + (n) * BYTES_PER_MP_LIMB || \ + (char *) (src) >= (char *) (dst) + \ + (n) * BYTES_PER_MP_LIMB); \ + memcpy ((dst), (src), (n) * BYTES_PER_MP_LIMB); \ + } \ + } \ + while (0) /* MPN macros taken from gmp-impl.h */ #define MPN_NORMALIZE(DST, NLIMBS) \ diff --git a/mpfr-impl.h b/mpfr-impl.h index 0f0ffe9f8..11ad48851 100644 --- a/mpfr-impl.h +++ b/mpfr-impl.h @@ -290,6 +290,7 @@ __MPFR_DECLSPEC extern MPFR_THREAD_ATTR const mpfr_t __gmpfr_four; /* Definition of constants */ #define LOG2 0.69314718055994528622 /* log(2) rounded to zero on 53 bits */ #define ALPHA 4.3191365662914471407 /* a+2 = a*log(a), rounded to +infinity */ +#define EXPM1 0.36787944117144227851 /* exp(-1), rounded to zero */ /* MPFR_DOUBLE_SPEC = 1 if the C type 'double' corresponds to IEEE-754 double precision, 0 if it doesn't, and undefined if one doesn't know. @@ -345,7 +346,13 @@ typedef union ieee_double_extract Ieee_double_extract; (((Ieee_double_extract *)&(x))->s.manh != 0))) #else # define DOUBLE_ISINF(x) ((x) > DBL_MAX || (x) < -DBL_MAX) -# define DOUBLE_ISNAN(x) ((x) != (x)) +# if MPFR_NANISNAN +/* Avoid MIPSpro / IRIX64 (incorrect) optimizations. + The + must not be replaced by a ||. */ +# define DOUBLE_ISNAN(x) (!(((x) >= 0.0) + ((x) <= 0.0))) +# else +# define DOUBLE_ISNAN(x) ((x) != (x)) +# endif #endif @@ -388,12 +395,12 @@ typedef union ieee_double_extract Ieee_double_extract; union { \ long double ld; \ struct { \ - unsigned long sign : 1; \ - unsigned long exp : 15; \ - unsigned long man3 : 16; \ - unsigned long man2 : 32; \ - unsigned long man1 : 32; \ - unsigned long man0 : 32; \ + unsigned int sign : 1; \ + unsigned int exp : 15; \ + unsigned int man3 : 16; \ + unsigned int man2 : 32; \ + unsigned int man1 : 32; \ + unsigned int man0 : 32; \ } s; \ } u; \ u.ld = (x); \ @@ -435,11 +442,11 @@ __MPFR_DECLSPEC long double __gmpfr_longdouble_volatile _MPFR_PROTO ((long doubl typedef union { long double ld; struct { - unsigned long manl : 32; - unsigned long manh : 32; - unsigned long expl : 8 ; - unsigned long exph : 7; - unsigned long sign : 1; + unsigned int manl : 32; + unsigned int manh : 32; + unsigned int expl : 8 ; + unsigned int exph : 7; + unsigned int sign : 1; } s; } mpfr_long_double_t; @@ -778,7 +785,7 @@ extern unsigned char *mpfr_stack; #if __MPFR_GNUC(2,95) || __MPFR_ICC(8,1,0) # define MPFR_INT_CEIL_LOG2(x) \ (__extension__ ({int _b; mp_limb_t _limb = (x); \ - MPFR_ASSERTD (_limb == (x)); \ + MPFR_ASSERTN (_limb == (x)); \ count_leading_zeros (_b, _limb); \ (BITS_PER_MP_LIMB - _b); })) #else @@ -1020,7 +1027,7 @@ typedef struct { /* * Round Mantissa (`srcp`, `sprec`) to mpfr_t `dest` using rounding mode `rnd` * assuming dest's sign is `sign`. - * Execute OVERFLOW_HANDLE in case of overflow when rounding (Power 2 case) + * Execute OVERFLOW_HANDLER in case of overflow when rounding (Power 2 case) * Return MPFR_EVEN_INEX in case of EVEN rounding */ #define MPFR_RNDRAW_EVEN(inexact, dest, srcp, sprec, rnd, sign, OVERFLOW_HANDLER)\ @@ -1457,7 +1464,7 @@ __MPFR_DECLSPEC void mpfr_setmin _MPFR_PROTO ((mpfr_ptr, mp_exp_t)); __MPFR_DECLSPEC long mpfr_mpn_exp _MPFR_PROTO ((mp_limb_t *, mp_exp_t *, int, mp_exp_t, size_t)); -#ifdef WANT_ASSERT +#ifdef _MPFR_H_HAVE_FILE __MPFR_DECLSPEC void mpfr_fprint_binary _MPFR_PROTO ((FILE *, mpfr_srcptr)); #endif __MPFR_DECLSPEC void mpfr_print_binary _MPFR_PROTO ((mpfr_srcptr)); @@ -505,6 +505,7 @@ __MPFR_DECLSPEC int mpfr_erfc _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_cbrt _MPFR_PROTO ((mpfr_ptr,mpfr_srcptr,mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_root _MPFR_PROTO ((mpfr_ptr,mpfr_srcptr,unsigned long,mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_gamma _MPFR_PROTO((mpfr_ptr,mpfr_srcptr,mpfr_rnd_t)); +__MPFR_DECLSPEC int mpfr_lngamma _MPFR_PROTO((mpfr_ptr,mpfr_srcptr,mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_zeta _MPFR_PROTO ((mpfr_ptr,mpfr_srcptr,mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_fac_ui _MPFR_PROTO ((mpfr_ptr, unsigned long int, mpfr_rnd_t)); @@ -3,7 +3,7 @@ @setfilename mpfr.info @documentencoding ISO-8859-1 @set VERSION 2.2.0 -@set UPDATED-MONTH August 2005 +@set UPDATED-MONTH September 2005 @settitle MPFR @value{VERSION} @synindex tp fn @iftex @@ -922,6 +922,13 @@ Convert @var{op} to a @code{mpz_t}, after rounding it with respect to @var{rnd}. If @var{op} is NaN or Inf, the result is undefined. @end deftypefun +@deftypefun int mpfr_get_f (mpf_t @var{rop}, mpfr_t @var{op}, mp_rnd_t @var{rnd}) +Convert @var{op} to a @code{mpf_t}, after rounding it with respect to +@var{rnd}. Return zero iff no error occurred, +in particular a non-zero value is returned if +@var{op} is NaN or Inf, which do not exist in @code{mpf}. +@end deftypefun + @deftypefun {char *} mpfr_get_str (char *@var{str}, mp_exp_t *@var{expptr}, int @var{b}, size_t @var{n}, mpfr_t @var{op}, mp_rnd_t @var{rnd}) Convert @var{op} to a string of digits in base @var{b}, with rounding in the direction @var{rnd}. The base may vary from 2 to 36. @@ -1071,6 +1078,7 @@ whatever the parity of @var{k}. @deftypefun int mpfr_pow (mpfr_t @var{rop}, mpfr_t @var{op1}, mpfr_t @var{op2}, mp_rnd_t @var{rnd}) @deftypefunx int mpfr_pow_ui (mpfr_t @var{rop}, mpfr_t @var{op1}, unsigned long int @var{op2}, mp_rnd_t @var{rnd}) @deftypefunx int mpfr_pow_si (mpfr_t @var{rop}, mpfr_t @var{op1}, long int @var{op2}, mp_rnd_t @var{rnd}) +@deftypefunx int mpfr_pow_z (mpfr_t @var{rop}, mpfr_t @var{op1}, mpz_t @var{op2}, mp_rnd_t @var{rnd}) @deftypefunx int mpfr_ui_pow_ui (mpfr_t @var{rop}, unsigned long int @var{op1}, unsigned long int @var{op2}, mp_rnd_t @var{rnd}) @deftypefunx int mpfr_ui_pow (mpfr_t @var{rop}, unsigned long int @var{op1}, mpfr_t @var{op2}, mp_rnd_t @var{rnd}) Set @var{rop} to @m{@var{op1}^{op2}, @var{op1} raised to @var{op2}}, @@ -1293,9 +1301,9 @@ for the @code{atan2} function (note this may change in future versions): @item @code{atan2(+INF, -INF)} returns @m{+3*\pi/4,+3*Pi/4}. @item @code{atan2(-INF, -INF)} returns @m{-3*\pi/4,-3*Pi/4}. @item @code{atan2(+INF, +INF)} returns @m{+\pi/4,+Pi/4}. -@item @code{atan2(-INF, +INF)} returns @m{-\pi/4,-Pi/4} -@item @code{atan2(+INF, x)} returns @m{+\pi/2,+Pi/2} for finite x. -@item @code{atan2(-INF, x)} returns @m{-\pi/2,-Pi/2} for finite x. +@item @code{atan2(-INF, +INF)} returns @m{-\pi/4,-Pi/4}. +@item @code{atan2(+INF, x)} returns @m{+\pi/2,+Pi/2} for finite @math{x}. +@item @code{atan2(-INF, x)} returns @m{-\pi/2,-Pi/2} for finite @math{x}. @item @code{atan2(y, -INF)} returns @m{+\pi,+Pi} for finite @math{y > 0}. @item @code{atan2(y, -INF)} returns @m{-\pi,-Pi} for finite @math{y < 0}. @item @code{atan2(y, +INF)} returns +0 for finite @math{y > 0}. @@ -1355,8 +1363,9 @@ For negative @var{x}, the returned value is NaN. @end deftypefun @deftypefun int mpfr_gamma (mpfr_t @var{rop}, mpfr_t @var{op}, mp_rnd_t @var{rnd}) +@deftypefunx int mpfr_lngamma (mpfr_t @var{rop}, mpfr_t @var{op}, mp_rnd_t @var{rnd}) Set @var{rop} to the value of the Gamma function on @var{op}, -rounded in the direction @var{rnd}. +and its logarithm respectively, rounded in the direction @var{rnd}. When @var{op} is a negative integer, NaN is returned. @end deftypefun @@ -2151,6 +2160,7 @@ a first implementation of the sine and cosine, and improved versions of @code{mpfr_const_log2} and @code{mpfr_const_pi}. Mathieu Dutour contributed the functions @code{mpfr_atan} and @code{mpfr_asin}, +and a previous version of @code{mpfr_gamma}; David Daney contributed the hyperbolic and inverse hyperbolic functions, the base-2 exponential, and the factorial function. Fabrice Rouillier contributed the original version of @file{mul_ui.c}, the @file{gmp_op.c} @@ -107,7 +107,7 @@ mpfr_mul3 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) MPFR_ASSERTD(tn <= k); /* Check for no size_t overflow*/ - MPFR_ASSERTD((size_t) k <= ((size_t) ~0) / BYTES_PER_MP_LIMB); + MPFR_ASSERTD((size_t) k <= ((size_t) -1) / BYTES_PER_MP_LIMB); MPFR_TMP_MARK(marker); tmp = (mp_limb_t *) MPFR_TMP_ALLOC((size_t) k * BYTES_PER_MP_LIMB); @@ -276,8 +276,8 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) if (MPFR_UNLIKELY (ax > __gmpfr_emax + 1)) return mpfr_overflow (a, rnd_mode, sign); if (MPFR_UNLIKELY (ax < __gmpfr_emin - 2)) - return mpfr_underflow (a, rnd_mode == GMP_RNDN ? GMP_RNDZ : rnd_mode, - sign); + return mpfr_underflow (a, rnd_mode == GMP_RNDN ? GMP_RNDZ : rnd_mode, + sign); #endif bq = MPFR_PREC (b); @@ -292,7 +292,7 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) MPFR_ASSERTD (tn <= k); /* tn <= k, thus no int overflow */ /* Check for no size_t overflow*/ - MPFR_ASSERTD ((size_t) k <= ((size_t) ~0) / BYTES_PER_MP_LIMB); + MPFR_ASSERTD ((size_t) k <= ((size_t) -1) / BYTES_PER_MP_LIMB); MPFR_TMP_MARK (marker); tmp = (mp_limb_t *) MPFR_TMP_ALLOC ((size_t) k * BYTES_PER_MP_LIMB); @@ -341,7 +341,6 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) /* Sum those two partial products */ add_ssaaaa (tmp[2], tmp[1], tmp[2], tmp[1], t1, t2); tmp[3] += (tmp[2] < t1); - b1 = tmp[3]; } b1 >>= (BITS_PER_MP_LIMB - 1); @@ -352,100 +351,143 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) else #endif if (MPFR_UNLIKELY (bn > MPFR_MUL_THRESHOLD)) - { - mp_limb_t *bp, *cp; - mp_size_t n; - mp_prec_t p; - - /* Compute estimated precision of mulhigh. - We could use `+ (n < cn) + (n < bn)' instead of `+ 2', - but does it worth it? */ - n = MPFR_LIMB_SIZE (a) + 1; - n = MIN (n, cn); - MPFR_ASSERTD (n >= 1 && 2*n <= k && n <= cn && n <= bn); - p = n*BITS_PER_MP_LIMB - MPFR_INT_CEIL_LOG2 (n + 2); - bp = MPFR_MANT (b) + bn - n; - cp = MPFR_MANT (c) + cn - n; - /* Check if MulHigh can produce a roundable result. - We may lost 1 bit due to RNDN, 1 due to final shift. */ - if (MPFR_UNLIKELY (MPFR_PREC (a) > p - 5)) - { - if (MPFR_UNLIKELY (MPFR_PREC (a) > p - 5 + BITS_PER_MP_LIMB - || bn <= MPFR_MUL_THRESHOLD+1)) - { - /* MulHigh can't produce a roundable result. */ - MPFR_LOG_MSG (("mpfr_mulhigh can't be used (%lu VS %lu)\n", - MPFR_PREC (a), p)); - goto full_multiply; - } - /* Add one extra limb to mantissa of b and c. */ - if (bn > n) - bp --; - else - { - bp = MPFR_TMP_ALLOC ((n+1)*sizeof (mp_limb_t)); - bp[0] = 0; - MPN_COPY (bp+1, MPFR_MANT (b)+bn-n, n); - } - if (cn > n) - cp --; /* FIXME: Could this happen? */ - else - { - cp = MPFR_TMP_ALLOC ((n+1)*sizeof (mp_limb_t)); - cp[0] = 0; - MPN_COPY (cp+1, MPFR_MANT (c)+cn-n, n); - } - /* We will compute with one extra limb */ - n++; - p = n*BITS_PER_MP_LIMB - MPFR_INT_CEIL_LOG2 (n + 2); - /* Due to some nasty reasons we can have only 4 bits */ - MPFR_ASSERTD (MPFR_PREC (a) <= p - 4); - - if (MPFR_LIKELY (k < 2*n)) - { - tmp = MPFR_TMP_ALLOC (2*n*sizeof (mp_limb_t)); - tmp += 2*n-k; /* `tmp' still points to an area of `k' limbs */ - } - } - MPFR_LOG_MSG (("Use mpfr_mulhigh (%lu VS %lu)\n", MPFR_PREC (a), p)); - /* Compute an approximation of the product of b and c */ - mpfr_mulhigh_n (tmp+k-2*n, bp, cp, n); - /* now tmp[0]..tmp[k-1] contains the product of both mantissa, - with tmp[k-1]>=2^(BITS_PER_MP_LIMB-2) */ - b1 = tmp[k-1] >> (BITS_PER_MP_LIMB - 1); /* msb from the product */ - - /* If the mantissas of b and c are uniformly distributed in ]1/2, 1], - then their product is in ]1/4, 1/2] with probability 2*ln(2)-1 ~ 0.386 - and in [1/2, 1] with probability 2-2*ln(2) ~ 0.614 */ - tmp += k - tn; - if (MPFR_UNLIKELY (b1 == 0)) - mpn_lshift (tmp, tmp, tn, 1); - MPFR_ASSERTD (MPFR_LIMB_MSB (tmp[tn-1]) != 0); - - if (MPFR_UNLIKELY (!mpfr_round_p (tmp, tn, p+b1-1, - MPFR_PREC(a)+(rnd_mode==GMP_RNDN)))) - { - tmp -= k-tn; /* tmp may have changed, FIX IT!!!!! */ - goto full_multiply; - } + { + mp_limb_t *bp, *cp; + mp_size_t n; + mp_prec_t p; + + /* Fist check if we can reduce the precision of b or c: + exact values are a nightmare for the short product trick */ + bp = MPFR_MANT (b); + cp = MPFR_MANT (c); + MPFR_ASSERTN (MPFR_MUL_THRESHOLD >= 1); + if (MPFR_UNLIKELY ((bp[0] == 0 && bp[1] == 0) || + (cp[0] == 0 && cp[1] == 0))) + { + mpfr_t b_tmp, c_tmp; + + MPFR_TMP_FREE (marker); + /* Check for b */ + while (*bp == 0) + { + bp++; + bn--; + MPFR_ASSERTD (bn > 0); + } /* This must end since the MSL is != 0 */ + + /* Check for c too */ + while (*cp == 0) + { + cp++; + cn--; + MPFR_ASSERTD (cn > 0); + } /* This must end since the MSL is != 0 */ + + /* It is not the faster way, but it is safer */ + MPFR_SET_SAME_SIGN (b_tmp, b); + MPFR_SET_EXP (b_tmp, MPFR_GET_EXP (b)); + MPFR_PREC (b_tmp) = bn * BITS_PER_MP_LIMB; + MPFR_MANT (b_tmp) = bp; + + MPFR_SET_SAME_SIGN (c_tmp, c); + MPFR_SET_EXP (c_tmp, MPFR_GET_EXP (c)); + MPFR_PREC (c_tmp) = cn * BITS_PER_MP_LIMB; + MPFR_MANT (c_tmp) = cp; + + /* Call again mpfr_mul with the fixed arguments */ + return mpfr_mul (a, b_tmp, c_tmp, rnd_mode); + } + + /* Compute estimated precision of mulhigh. + We could use `+ (n < cn) + (n < bn)' instead of `+ 2', + but does it worth it? */ + n = MPFR_LIMB_SIZE (a) + 1; + n = MIN (n, cn); + MPFR_ASSERTD (n >= 1 && 2*n <= k && n <= cn && n <= bn); + p = n * BITS_PER_MP_LIMB - MPFR_INT_CEIL_LOG2 (n + 2); + bp += bn - n; + cp += cn - n; + + /* Check if MulHigh can produce a roundable result. + We may lost 1 bit due to RNDN, 1 due to final shift. */ + if (MPFR_UNLIKELY (MPFR_PREC (a) > p - 5)) + { + if (MPFR_UNLIKELY (MPFR_PREC (a) > p - 5 + BITS_PER_MP_LIMB + || bn <= MPFR_MUL_THRESHOLD+1)) + { + /* MulHigh can't produce a roundable result. */ + MPFR_LOG_MSG (("mpfr_mulhigh can't be used (%lu VS %lu)\n", + MPFR_PREC (a), p)); + goto full_multiply; + } + /* Add one extra limb to mantissa of b and c. */ + if (bn > n) + bp --; + else + { + bp = MPFR_TMP_ALLOC ((n+1) * sizeof (mp_limb_t)); + bp[0] = 0; + MPN_COPY (bp + 1, MPFR_MANT (b) + bn - n, n); + } + if (cn > n) + cp --; /* FIXME: Could this happen? */ + else + { + cp = MPFR_TMP_ALLOC ((n+1) * sizeof (mp_limb_t)); + cp[0] = 0; + MPN_COPY (cp + 1, MPFR_MANT (c) + cn - n, n); + } + /* We will compute with one extra limb */ + n++; + p = n * BITS_PER_MP_LIMB - MPFR_INT_CEIL_LOG2 (n + 2); + /* Due to some nasty reasons we can have only 4 bits */ + MPFR_ASSERTD (MPFR_PREC (a) <= p - 4); + + if (MPFR_LIKELY (k < 2*n)) + { + tmp = MPFR_TMP_ALLOC (2 * n * sizeof (mp_limb_t)); + tmp += 2*n-k; /* `tmp' still points to an area of `k' limbs */ + } + } + MPFR_LOG_MSG (("Use mpfr_mulhigh (%lu VS %lu)\n", MPFR_PREC (a), p)); + /* Compute an approximation of the product of b and c */ + mpfr_mulhigh_n (tmp+k-2*n, bp, cp, n); + /* now tmp[0]..tmp[k-1] contains the product of both mantissa, + with tmp[k-1]>=2^(BITS_PER_MP_LIMB-2) */ + b1 = tmp[k-1] >> (BITS_PER_MP_LIMB - 1); /* msb from the product */ + + /* If the mantissas of b and c are uniformly distributed in (1/2, 1], + then their product is in (1/4, 1/2] with probability 2*ln(2)-1 + ~ 0.386 and in [1/2, 1] with probability 2-2*ln(2) ~ 0.614 */ + tmp += k - tn; + if (MPFR_UNLIKELY (b1 == 0)) + mpn_lshift (tmp, tmp, tn, 1); + MPFR_ASSERTD (MPFR_LIMB_MSB (tmp[tn-1]) != 0); + + if (MPFR_UNLIKELY (!mpfr_round_p (tmp, tn, p+b1-1, MPFR_PREC(a) + + (rnd_mode == GMP_RNDN)))) + { + tmp -= k - tn; /* tmp may have changed, FIX IT!!!!! */ + goto full_multiply; + } + } + else + { + full_multiply: + MPFR_LOG_MSG (("Use mpn_mul\n", 0)); + b1 = mpn_mul (tmp, MPFR_MANT (b), bn, MPFR_MANT (c), cn); + + /* now tmp[0]..tmp[k-1] contains the product of both mantissa, + with tmp[k-1]>=2^(BITS_PER_MP_LIMB-2) */ + b1 >>= BITS_PER_MP_LIMB - 1; /* msb from the product */ + + /* if the mantissas of b and c are uniformly distributed in (1/2, 1], + then their product is in (1/4, 1/2] with probability 2*ln(2)-1 + ~ 0.386 and in [1/2, 1] with probability 2-2*ln(2) ~ 0.614 */ + tmp += k - tn; + if (MPFR_UNLIKELY (b1 == 0)) + mpn_lshift (tmp, tmp, tn, 1); /* tn <= k, so no stack corruption */ } - else - { - full_multiply: - MPFR_LOG_MSG (("Use mpn_mul\n", 0)); - b1 = mpn_mul (tmp, MPFR_MANT (b), bn, MPFR_MANT (c), cn); - - /* now tmp[0]..tmp[k-1] contains the product of both mantissa, - with tmp[k-1]>=2^(BITS_PER_MP_LIMB-2) */ - b1 >>= BITS_PER_MP_LIMB - 1; /* msb from the product */ - - /* if the mantissas of b and c are uniformly distributed in ]1/2, 1], - then their product is in ]1/4, 1/2] with probability 2*ln(2)-1 ~ 0.386 - and in [1/2, 1] with probability 2-2*ln(2) ~ 0.614 */ - tmp += k - tn; - if (MPFR_UNLIKELY (b1 == 0)) - mpn_lshift (tmp, tmp, tn, 1); /* tn <= k, so no stack corruption */ - } ax2 = ax + (mp_exp_t) (b1 - 1); MPFR_RNDRAW (inexact, a, tmp, bq+cq, rnd_mode, sign, ax2++); @@ -468,4 +510,3 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) } return inexact; } - diff --git a/set_str_raw.c b/set_str_raw.c index 82942b6b7..b9a71655e 100644 --- a/set_str_raw.c +++ b/set_str_raw.c @@ -1,6 +1,6 @@ /* mpfr_set_str_binary -- set a floating-point number from a binary string -Copyright 1999, 2000, 2001, 2002, 2003, 2004 Free Software Foundation, Inc. +Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005 Free Software Foundation, Inc. This file is part of the MPFR Library. @@ -19,11 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 51 Franklin Place, Fifth Floor, Boston, MA 02110-1301, USA. */ -#include <stdlib.h> -#include <string.h> -#include <errno.h> - -#define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" /* Currently the number should be of the form +/- xxxx.xxxxxxEyy, with @@ -33,23 +28,8 @@ MA 02110-1301, USA. */ void mpfr_set_str_binary (mpfr_ptr x, const char *str) { - char *str2, *str0, negative = 0; - unsigned long j, l, k = 0, xsize, cnt, alloc; - mp_limb_t *xp; - long expn = 0, e; - char *endstr2; - - xp = MPFR_MANT(x); - xsize = MPFR_LIMB_SIZE(x); - alloc = strlen(str) + 1; - - if (*str == '-') - { - negative = 1; - str++; - } - else if (*str == '+') - str++; + int has_sign; + int res; if (*str == 'N') { @@ -58,94 +38,17 @@ mpfr_set_str_binary (mpfr_ptr x, const char *str) return; } - if (*str == 'I') + has_sign = *str == '-' || *str == '+'; + if (str[has_sign] == 'I') { MPFR_SET_INF(x); - if (negative) + if (*str == '-') MPFR_SET_NEG(x); else MPFR_SET_POS(x); - /*if (MPFR_IS_STRICTNEG(x) != negative) - MPFR_CHANGE_SIGN(x);*/ return; } - MPFR_CLEAR_FLAGS(x); - - str0 = str2 = (char *) (*__gmp_allocate_func) (alloc); - - while (*str == '0') - str++; - - while (*str == '0' || *str == '1') - { - *(str2++) = *(str++); - k++; - } - - if (*str == '.') - { - str++; - while (*str == '0' || *str == '1') - { - *(str2++) = *(str++); - } - } - - if (*str == 'e' || *str == 'E') - { - char *end; - - e = strtol (++str, &end, 10); /* signed exponent after 'e' or 'E' */ - MPFR_ASSERTN (*end == '\0' && errno != ERANGE); - expn = k + e; - MPFR_ASSERTN (expn >= e); - } - else - expn = k; - - endstr2 = str2; - l = endstr2 - str0; /* length of mantissa */ - if (l == 0) - { /* input is zero */ - MPFR_SET_ZERO(x); - } - else - { - l = (l - 1) / BITS_PER_MP_LIMB + 1; - str2 = str0; - - MPFR_ASSERTN (l <= xsize); - - /* str2[0]..endstr2[-1] contains the mantissa */ - for (k = 1; k <= l; k++) - { - j = 0; - xp[xsize - k] = 0; - while ((str2 < endstr2) && (j < BITS_PER_MP_LIMB)) - { - xp[xsize - k] = (xp[xsize - k] << 1) + (*str2 - '0'); - str2++; j++; - } - xp[xsize - k] <<= (BITS_PER_MP_LIMB - j); - } - - for (; k <= xsize; k++) - xp[xsize - k] = 0; - - count_leading_zeros(cnt, xp[xsize - 1]); - if (cnt) - mpn_lshift(xp, xp, xsize, cnt); - - MPFR_SET_EXP (x, expn - cnt); - if (negative) - MPFR_SET_NEG(x); - else - MPFR_SET_POS(x); - } - - (*__gmp_free_func) (str0, alloc); - - /* May change to take into account : syntax errors, overflow in exponent, - string truncated because of size of x */ + res = mpfr_strtofr (x, str, 0, 2, GMP_RNDZ); + MPFR_ASSERTN (res == 0); } @@ -32,13 +32,17 @@ int mpfr_sub1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode); int mpfr_sub1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) { mpfr_t tmpa, tmpb, tmpc; - int inexact, inexact2; + int inexb, inexc, inexact, inexact2; mpfr_init2 (tmpa, MPFR_PREC (a)); mpfr_init2 (tmpb, MPFR_PREC (b)); mpfr_init2 (tmpc, MPFR_PREC (c)); - MPFR_ASSERTN (mpfr_set (tmpb, b, GMP_RNDN) == 0); - MPFR_ASSERTN (mpfr_set (tmpc, c, GMP_RNDN) == 0); + + inexb = mpfr_set (tmpb, b, GMP_RNDN); + MPFR_ASSERTN (inexb == 0); + + inexc = mpfr_set (tmpc, c, GMP_RNDN); + MPFR_ASSERTN (inexc == 0); inexact2 = mpfr_sub1 (tmpa, tmpb, tmpc, rnd_mode); inexact = mpfr_sub1sp2(a, b, c, rnd_mode); @@ -39,88 +39,89 @@ static void count_sort (mpfr_srcptr *const, unsigned long, mpfr_srcptr *, /* Either sort the tab in perm and returns 0 Or returns 1 for +INF, -1 for -INF and 2 for NAN */ -int mpfr_sum_sort (mpfr_srcptr *const tab, unsigned long n, - mpfr_srcptr *perm) +int +mpfr_sum_sort (mpfr_srcptr *const tab, unsigned long n, mpfr_srcptr *perm) { - mp_exp_t min, max; - mpfr_uexp_t exp_num; - unsigned long i; - int sign_inf; - - sign_inf = 0; - min = MPFR_EMIN_MAX; - max = MPFR_EMAX_MIN; - for (i = 0; i < n; i++) - { - if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tab[i]))) - { - if (MPFR_IS_NAN (tab[i])) - return 2; /* Return NAN code */ - else if (MPFR_IS_INF (tab[i])) - { - if (sign_inf == 0) /* No previous INF */ - sign_inf = MPFR_SIGN (tab[i]); - else if (sign_inf != MPFR_SIGN (tab[i])) - return 2; /* Return NAN */ - } - } - else - { - if (MPFR_GET_EXP (tab[i]) < min) - min = MPFR_GET_EXP(tab[i]); - if (MPFR_GET_EXP (tab[i]) > max) - max = MPFR_GET_EXP(tab[i]); - } - } - if (MPFR_UNLIKELY (sign_inf != 0)) - return sign_inf; - - exp_num = max - min + 1; - /* FIXME : better test */ - if (exp_num > n * MPFR_INT_CEIL_LOG2 (n)) - heap_sort (tab, n, perm); - else - count_sort (tab, n, perm, min, exp_num); - return 0; + mp_exp_t min, max; + mpfr_uexp_t exp_num; + unsigned long i; + int sign_inf; + + sign_inf = 0; + min = MPFR_EMIN_MAX; + max = MPFR_EMAX_MIN; + for (i = 0; i < n; i++) + { + if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tab[i]))) + { + if (MPFR_IS_NAN (tab[i])) + return 2; /* Return NAN code */ + else if (MPFR_IS_INF (tab[i])) + { + if (sign_inf == 0) /* No previous INF */ + sign_inf = MPFR_SIGN (tab[i]); + else if (sign_inf != MPFR_SIGN (tab[i])) + return 2; /* Return NAN */ + } + } + else + { + if (MPFR_GET_EXP (tab[i]) < min) + min = MPFR_GET_EXP(tab[i]); + if (MPFR_GET_EXP (tab[i]) > max) + max = MPFR_GET_EXP(tab[i]); + } + } + if (MPFR_UNLIKELY (sign_inf != 0)) + return sign_inf; + + exp_num = max - min + 1; + /* FIXME : better test */ + if (exp_num > n * MPFR_INT_CEIL_LOG2 (n)) + heap_sort (tab, n, perm); + else + count_sort (tab, n, perm, min, exp_num); + return 0; } #define GET_EXP1(x) (MPFR_IS_ZERO (x) ? min : MPFR_GET_EXP (x)) /* Performs a count sort of the entries */ -static void count_sort (mpfr_srcptr *const tab, unsigned long n, - mpfr_srcptr *perm, mp_exp_t min, mpfr_uexp_t exp_num) +static void +count_sort (mpfr_srcptr *const tab, unsigned long n, + mpfr_srcptr *perm, mp_exp_t min, mpfr_uexp_t exp_num) { - unsigned long *account; - unsigned long target_rank, i; - MPFR_TMP_DECL(marker); - - /* Reserve a place for potential 0 (with EXP min-1) - If there is no zero, we only lose one unused entry */ - min--; - exp_num++; - - /* Performs a counting sort of the entries */ - MPFR_TMP_MARK (marker); - account = (unsigned long *) MPFR_TMP_ALLOC (exp_num * sizeof * account); - for (i = 0; i < exp_num; i++) - account[i] = 0; - for (i = 0; i < n; i++) - account[GET_EXP1 (tab[i]) - min]++; - for (i = exp_num - 1; i >= 1; i--) - account[i - 1] += account[i]; - for (i = 0; i < n; i++) - { - target_rank = --account[GET_EXP1 (tab[i]) - min]; - perm[target_rank] = tab[i]; - } - MPFR_TMP_FREE (marker); + unsigned long *account; + unsigned long target_rank, i; + MPFR_TMP_DECL(marker); + + /* Reserve a place for potential 0 (with EXP min-1) + If there is no zero, we only lose one unused entry */ + min--; + exp_num++; + + /* Performs a counting sort of the entries */ + MPFR_TMP_MARK (marker); + account = (unsigned long *) MPFR_TMP_ALLOC (exp_num * sizeof * account); + for (i = 0; i < exp_num; i++) + account[i] = 0; + for (i = 0; i < n; i++) + account[GET_EXP1 (tab[i]) - min]++; + for (i = exp_num - 1; i >= 1; i--) + account[i - 1] += account[i]; + for (i = 0; i < n; i++) + { + target_rank = --account[GET_EXP1 (tab[i]) - min]; + perm[target_rank] = tab[i]; + } + MPFR_TMP_FREE (marker); } #define GET_EXP2(x) (MPFR_IS_ZERO (x) ? MPFR_EMIN_MIN : MPFR_GET_EXP (x)) /* Performs a heap sort of the entries */ -static void heap_sort (mpfr_srcptr *const tab, unsigned long n, - mpfr_srcptr *perm) +static void +heap_sort (mpfr_srcptr *const tab, unsigned long n, mpfr_srcptr *perm) { unsigned long dernier_traite; unsigned long i, pere; @@ -207,8 +208,8 @@ static void heap_sort (mpfr_srcptr *const tab, unsigned long n, * intermediate size set to F. * Internal use function. */ -static int sum_once (mpfr_ptr ret, mpfr_srcptr *const tab, - unsigned long n, mp_prec_t F) +static int +sum_once (mpfr_ptr ret, mpfr_srcptr *const tab, unsigned long n, mp_prec_t F) { mpfr_t sum; unsigned long i; @@ -230,7 +231,7 @@ static int sum_once (mpfr_ptr ret, mpfr_srcptr *const tab, /* Sum a list of floating-point numbers. * FIXME : add reference to Demmel-Hida's paper. -*/ + */ int mpfr_sum (mpfr_ptr ret, mpfr_ptr *const tab_p, unsigned long n, mp_rnd_t rnd) @@ -257,11 +258,12 @@ mpfr_sum (mpfr_ptr ret, mpfr_ptr *const tab_p, unsigned long n, mp_rnd_t rnd) /* Sort and treat special cases */ MPFR_TMP_MARK (marker); - perm = (mpfr_srcptr *) MPFR_TMP_ALLOC (n * sizeof * perm); + perm = (mpfr_srcptr *) MPFR_TMP_ALLOC (n * sizeof *perm); error_trap = mpfr_sum_sort (tab, n, perm); /* Check if there was a NAN or a INF */ if (MPFR_UNLIKELY (error_trap != 0)) { + MPFR_TMP_FREE (marker); if (error_trap == 2) { MPFR_SET_NAN (ret); diff --git a/tests/Makefile.am b/tests/Makefile.am index 9652a1343..c4306006a 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -1,6 +1,6 @@ -AUTOMAKE_OPTIONS = gnu no-dependencies $(top_builddir)/ansi2knr +AUTOMAKE_OPTIONS = gnu $(top_builddir)/ansi2knr -check_PROGRAMS = tversion tinits tsgn tcheck tisnan texceptions tset_exp tset tabs tset_d tset_f tset_q tset_si tset_str tset_z tset_ld tset_sj tswap tcopysign tcmp tcmp2 tcmpabs tcmp_d tcmp_ld tcomparisons teq tadd tsub tmul tdiv tsub1sp tadd1sp tadd_ui tsub_ui tcmp_ui tdiv_ui tmul_ui tsqrt_ui tui_div tui_sub tgmpop tsi_op tmul_2exp tfma tsum tdim tminmax tnext tfits tget_d tget_d_2exp tget_z tget_str tget_sj tout_str tinp_str toutimpl tcan_round tround_prec tsqrt tconst_log2 tconst_pi tconst_euler trandom ttrunc trint tfrac texp texp2 texpm1 tlog tlog2 tlog10 tlog1p tpow tui_pow tpow3 tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic tasin tacos tcos tatan tsin ttan tsin_cos tagm thypot tfactorial tgamma terf tcbrt tzeta mpf_compat mpfr_compat reuse tsqr tstrtofr tpow_z tget_f tconst_catalan troot tsec tcsc tcot teint tcoth tcsch tsech tstckintc tsubnormal +check_PROGRAMS = tversion tinits tsgn tcheck tisnan texceptions tset_exp tset tabs tset_d tset_f tset_q tset_si tset_str tset_z tset_ld tset_sj tswap tcopysign tcmp tcmp2 tcmpabs tcmp_d tcmp_ld tcomparisons teq tadd tsub tmul tdiv tsub1sp tadd1sp tadd_ui tsub_ui tcmp_ui tdiv_ui tmul_ui tsqrt_ui tui_div tui_sub tgmpop tsi_op tmul_2exp tfma tsum tdim tminmax tnext tfits tget_d tget_d_2exp tget_z tget_str tget_sj tout_str tinp_str toutimpl tcan_round tround_prec tsqrt tconst_log2 tconst_pi tconst_euler trandom ttrunc trint tfrac texp texp2 texpm1 tlog tlog2 tlog10 tlog1p tpow tui_pow tpow3 tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic tasin tacos tcos tatan tsin ttan tsin_cos tagm thypot tfactorial tgamma terf tcbrt tzeta mpf_compat mpfr_compat reuse tsqr tstrtofr tpow_z tget_f tconst_catalan troot tsec tcsc tcot teint tcoth tcsch tsech tstckintc tsubnormal tlngamma EXTRA_DIST = tgeneric.c tgeneric_ui.c mpf_compat.h inp_str.data diff --git a/tests/tcoth.c b/tests/tcoth.c index 74df2d134..80302bfeb 100644 --- a/tests/tcoth.c +++ b/tests/tcoth.c @@ -79,12 +79,36 @@ check_specials (void) mpfr_clear (y); } +static void +check_bugs (void) +{ + mpfr_t x, y; + + mpfr_init (x); + mpfr_init (y); + + /* bug found by Rob (Sisyphus) on 16 Sep 2005 */ + mpfr_set_ui (x, 2, GMP_RNDN); + mpfr_set_prec (y, 2); + mpfr_coth (y, x, GMP_RNDN); + if (mpfr_cmp_ui (y, 1)) + { + printf ("Error for coth(2), expected 1, got "); + mpfr_dump (y); + exit (1); + } + + mpfr_clear (x); + mpfr_clear (y); +} + int main (int argc, char *argv[]) { tests_start_mpfr (); check_specials (); + check_bugs (); test_generic (2, 200, 10); tests_end_mpfr (); diff --git a/tests/texp.c b/tests/texp.c index d398ca2ad..82c0217c3 100644 --- a/tests/texp.c +++ b/tests/texp.c @@ -339,12 +339,16 @@ check_special () if (mpfr_cmp_ui_2exp (y, 3, -2)) { printf ("Error for exp(-1/16), prec=2, RNDD\n"); + printf ("expected 0.11, got "); + mpfr_dump (y); exit (1); } test_exp (y, x, GMP_RNDZ); - if (mpfr_cmp_ui (y, 1)) + if (mpfr_cmp_ui_2exp (y, 3, -2)) { printf ("Error for exp(-1/16), prec=2, RNDZ\n"); + printf ("expected 0.11, got "); + mpfr_dump (y); exit (1); } mpfr_set_str_binary (x, "0.1E-3"); @@ -426,6 +430,16 @@ check_special () mpfr_out_str (stdout, 16, 0, x, GMP_RNDN); putchar ('\n'); } + /* bug found by Guillaume Melquiond, 13 Sep 2005 */ + mpfr_set_prec (x, 53); + mpfr_set_str_binary (x, "-1E-400"); + mpfr_exp (x, x, GMP_RNDZ); + if (mpfr_cmp_ui (x, 1) == 0) + { + printf ("Error for exp(-2^(-400))\n"); + exit (1); + } + mpfr_clear (x); mpfr_clear (y); mpfr_clear (z); diff --git a/tests/tgamma.c b/tests/tgamma.c index f246773d6..176005793 100644 --- a/tests/tgamma.c +++ b/tests/tgamma.c @@ -24,9 +24,8 @@ MA 02110-1301, USA. */ #include "mpfr-test.h" -int mpfr_gamma (mpfr_ptr, mpfr_srcptr, mp_rnd_t); - #define TEST_FUNCTION mpfr_gamma +#define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 1) #include "tgeneric.c" static void @@ -99,25 +98,27 @@ special (void) mpfr_set_prec (y, 53); #define CHECK_X1 "1.0762904832837976166" -#define CHECK_Y1 0.96134843256452096050 +#define CHECK_Y1 "0.96134843256452096050" mpfr_set_str (x, CHECK_X1, 10, GMP_RNDN); mpfr_gamma (y, x, GMP_RNDN); - if (mpfr_get_d (y, GMP_RNDN) != CHECK_Y1 ) + mpfr_set_str (x, CHECK_Y1, 10, GMP_RNDN); + if (mpfr_cmp (y, x)) { printf ("mpfr_gamma("CHECK_X1") is wrong: expected %1.20e, got %1.20e\n", - CHECK_Y1, mpfr_get_d (y, GMP_RNDN)); + mpfr_get_d (x, GMP_RNDN), mpfr_get_d (y, GMP_RNDN)); exit (1); } #define CHECK_X2 "9.23709516716202383435e-01" -#define CHECK_Y2 1.0502315560291053398 +#define CHECK_Y2 "1.0502315560291053398" mpfr_set_str (x, CHECK_X2, 10, GMP_RNDN); mpfr_gamma (y, x, GMP_RNDN); - if (mpfr_get_d (y, GMP_RNDN) != CHECK_Y2) + mpfr_set_str (x, CHECK_Y2, 10, GMP_RNDN); + if (mpfr_cmp (y, x)) { printf ("mpfr_gamma("CHECK_X2") is wrong: expected %1.20e, got %1.20e\n", - CHECK_Y2, mpfr_get_d (y, GMP_RNDN)); + mpfr_get_d (x, GMP_RNDN), mpfr_get_d (y, GMP_RNDN)); exit (1); } @@ -189,6 +190,7 @@ static void special_overflow (void) { mpfr_t x, y; + mp_exp_t emin = mpfr_get_emin (); set_emin (-125); set_emax (128); @@ -199,11 +201,162 @@ special_overflow (void) mpfr_gamma (y, x, GMP_RNDN); if (!mpfr_inf_p (y)) { - printf("Overflow error.\n"); + printf ("Overflow error.\n"); + mpfr_dump (y); + exit (1); + } + + /* problem mentioned by Kenneth Wilder, 18 Aug 2005 */ + mpfr_set_prec (x, 29); + mpfr_set_prec (y, 29); + mpfr_set_str (x, "-200000000.5", 10, GMP_RNDN); /* exact */ + mpfr_gamma (y, x, GMP_RNDN); + if (!(mpfr_zero_p (y) && MPFR_SIGN (y) < 0)) + { + printf ("Error for gamma(-200000000.5)\n"); + printf ("expected -0"); + printf ("got "); + mpfr_dump (y); + exit (1); + } + + mpfr_set_prec (x, 53); + mpfr_set_prec (y, 53); + mpfr_set_str (x, "-200000000.1", 10, GMP_RNDN); + mpfr_gamma (y, x, GMP_RNDN); + if (!(mpfr_zero_p (y) && MPFR_SIGN (y) < 0)) + { + printf ("Error for gamma(-200000000.1), prec=53\n"); + printf ("expected -0"); + printf ("got "); mpfr_dump (y); exit (1); } + /* another problem mentioned by Kenneth Wilder, 29 Aug 2005 */ + mpfr_set_prec (x, 333); + mpfr_set_prec (y, 14); + mpfr_set_str (x, "-2.0000000000000000000000005", 10, GMP_RNDN); + mpfr_gamma (y, x, GMP_RNDN); + mpfr_set_prec (x, 14); + mpfr_set_str_binary (x, "-11010011110001E66"); + if (mpfr_cmp (x, y)) + { + printf ("Error for gamma(-2.0000000000000000000000005)\n"); + printf ("expected "); mpfr_dump (x); + printf ("got "); mpfr_dump (y); + exit (1); + } + + /* another tests from Kenneth Wilder, 31 Aug 2005 */ + set_emax (200); + set_emin (-200); + mpfr_set_prec (x, 38); + mpfr_set_prec (y, 54); + mpfr_set_str_binary (x, "0.11101111011100111101001001010110101001E-166"); + mpfr_gamma (y, x, GMP_RNDN); + mpfr_set_prec (x, 54); + mpfr_set_str_binary (x, "0.100010001101100001110110001010111111010000100101011E167"); + if (mpfr_cmp (x, y)) + { + printf ("Error for gamma (test 1)\n"); + printf ("expected "); mpfr_dump (x); + printf ("got "); mpfr_dump (y); + exit (1); + } + + set_emax (1000); + set_emin (-2000); + mpfr_set_prec (x, 38); + mpfr_set_prec (y, 71); + mpfr_set_str_binary (x, "10101011011100001111111000010111110010E-1034"); + /* 184083777010*2^(-1034) */ + mpfr_gamma (y, x, GMP_RNDN); + mpfr_set_prec (x, 71); + mpfr_set_str_binary (x, "10111111001000011110010001000000000000110011110000000011101011111111100E926"); + /* 1762885132679550982140*2^926 */ + if (mpfr_cmp (x, y)) + { + printf ("Error for gamma (test 2)\n"); + printf ("expected "); mpfr_dump (x); + printf ("got "); mpfr_dump (y); + exit (1); + } + + mpfr_set_prec (x, 38); + mpfr_set_prec (y, 88); + mpfr_set_str_binary (x, "10111100111001010000100001100100100101E-104"); + /* 202824096037*2^(-104) */ + mpfr_gamma (y, x, GMP_RNDN); + mpfr_set_prec (x, 88); + mpfr_set_str_binary (x, "1010110101111000111010111100010110101010100110111000001011000111000011101100001101110010E-21"); + /* 209715199999500283894743922*2^(-21) */ + if (mpfr_cmp (x, y)) + { + printf ("Error for gamma (test 3)\n"); + printf ("expected "); mpfr_dump (x); + printf ("got "); mpfr_dump (y); + exit (1); + } + + mpfr_set_prec (x, 171); + mpfr_set_prec (y, 38); + mpfr_set_str (x, "-2993155353253689176481146537402947624254601559176535", 10, + GMP_RNDN); + mpfr_div_2exp (x, x, 170, GMP_RNDN); + mpfr_gamma (y, x, GMP_RNDN); + mpfr_set_prec (x, 38); + mpfr_set_str (x, "201948391737", 10, GMP_RNDN); + mpfr_mul_2exp (x, x, 92, GMP_RNDN); + if (mpfr_cmp (x, y)) + { + printf ("Error for gamma (test 5)\n"); + printf ("expected "); mpfr_dump (x); + printf ("got "); mpfr_dump (y); + exit (1); + } + + set_emin (-500000); + mpfr_set_prec (x, 337); + mpfr_set_prec (y, 38); + mpfr_set_str (x, "-30000.000000000000000000000000000000000000000000001", 10, + GMP_RNDN); + mpfr_gamma (y, x, GMP_RNDN); + mpfr_set_prec (x, 38); + mpfr_set_str (x, "-3.623795987425E-121243", 10, GMP_RNDN); + if (mpfr_cmp (x, y)) + { + printf ("Error for gamma (test 7)\n"); + printf ("expected "); mpfr_dump (x); + printf ("got "); mpfr_dump (y); + exit (1); + } + + /* was producing infinite loop */ + set_emin (emin); + mpfr_set_prec (x, 71); + mpfr_set_prec (y, 71); + mpfr_set_str (x, "-200000000.1", 10, GMP_RNDN); + mpfr_gamma (y, x, GMP_RNDN); + if (!(mpfr_zero_p (y) && MPFR_SIGN (y) < 0)) + { + printf ("Error for gamma (test 8)\n"); + printf ("expected "); mpfr_dump (x); + printf ("got "); mpfr_dump (y); + exit (1); + } + + set_emax (1073741823); + mpfr_set_prec (x, 29); + mpfr_set_prec (y, 29); + mpfr_set_str (x, "423786866", 10, GMP_RNDN); + mpfr_gamma (y, x, GMP_RNDN); + if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0) + { + printf ("Error for gamma(423786866)\n"); + exit (1); + } + mpfr_clear (y); mpfr_clear (x); set_emin (MPFR_EMIN_MIN); @@ -213,6 +366,7 @@ special_overflow (void) int main (void) { + MPFR_TEST_USE_RANDS (); tests_start_mpfr (); special (); diff --git a/tests/tlngamma.c b/tests/tlngamma.c new file mode 100644 index 000000000..a4e421ba3 --- /dev/null +++ b/tests/tlngamma.c @@ -0,0 +1,192 @@ +/* mpfr_tlngamma -- test file for lngamma function + +Copyright 2005 Free Software Foundation. + +This file is part of the MPFR Library. + +The MPFR 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 MPFR 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 MPFR Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 51 Franklin Place, Fifth Floor, Boston, +MA 02110-1301, USA. */ + +#include <stdio.h> +#include <stdlib.h> + +#include "mpfr-test.h" + +#define TEST_FUNCTION mpfr_lngamma +#include "tgeneric.c" + +static void +special (void) +{ + mpfr_t x, y; + int inex; + + mpfr_init (x); + mpfr_init (y); + + mpfr_set_nan (x); + mpfr_lngamma (y, x, GMP_RNDN); + if (!mpfr_nan_p (y)) + { + printf ("Error for lngamma(NaN)\n"); + exit (1); + } + + mpfr_set_inf (x, -1); + mpfr_lngamma (y, x, GMP_RNDN); + if (!mpfr_nan_p (y)) + { + printf ("Error for lngamma(-Inf)\n"); + exit (1); + } + + mpfr_set_inf (x, 1); + mpfr_lngamma (y, x, GMP_RNDN); + if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0) + { + printf ("Error for lngamma(+Inf)\n"); + exit (1); + } + + mpfr_set_ui (x, 0, GMP_RNDN); + mpfr_lngamma (y, x, GMP_RNDN); + if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0) + { + printf ("Error for lngamma(+0)\n"); + exit (1); + } + + mpfr_set_ui (x, 0, GMP_RNDN); + mpfr_neg (x, x, GMP_RNDN); + mpfr_lngamma (y, x, GMP_RNDN); + if (!mpfr_nan_p (y)) + { + printf ("Error for lngamma(-0)\n"); + exit (1); + } + + mpfr_set_ui (x, 1, GMP_RNDN); + mpfr_lngamma (y, x, GMP_RNDN); + if (mpfr_cmp_ui (y, 0)) + { + printf ("Error for lngamma(1)\n"); + exit (1); + } + + mpfr_set_si (x, -1, GMP_RNDN); + mpfr_lngamma (y, x, GMP_RNDN); + if (!mpfr_nan_p (y)) + { + printf ("Error for lngamma(-1)\n"); + exit (1); + } + + mpfr_set_prec (x, 53); + mpfr_set_prec (y, 53); + +#define CHECK_X1 "1.0762904832837976166" +#define CHECK_Y1 "-0.039418362817587634939" + + mpfr_set_str (x, CHECK_X1, 10, GMP_RNDN); + mpfr_lngamma (y, x, GMP_RNDN); + mpfr_set_str (x, CHECK_Y1, 10, GMP_RNDN); + if (mpfr_cmp (y, x)) + { + printf ("mpfr_lngamma("CHECK_X1") is wrong: expected %1.20e, got %1.20e\n", + mpfr_get_d (x, GMP_RNDN), mpfr_get_d (y, GMP_RNDN)); + exit (1); + } + +#define CHECK_X2 "9.23709516716202383435e-01" +#define CHECK_Y2 "0.049010669407893718563" + mpfr_set_str (x, CHECK_X2, 10, GMP_RNDN); + mpfr_lngamma (y, x, GMP_RNDN); + mpfr_set_str (x, CHECK_Y2, 10, GMP_RNDN); + if (mpfr_cmp (y, x)) + { + printf ("mpfr_lngamma("CHECK_X2") is wrong: expected %1.20e, got %1.20e\n", + mpfr_get_d (x, GMP_RNDN), mpfr_get_d (y, GMP_RNDN)); + exit (1); + } + + mpfr_set_prec (x, 8); + mpfr_set_prec (y, 175); + mpfr_set_ui (x, 33, GMP_RNDN); + mpfr_lngamma (y, x, GMP_RNDU); + mpfr_set_prec (x, 175); + mpfr_set_str_binary (x, "0.1010001100011101101011001101110010100001000001000001110011000001101100001111001001000101011011100100010101011110100111110101010100010011010010000101010111001100011000101111E7"); + if (mpfr_cmp (x, y)) + { + printf ("Error in mpfr_lngamma (1)\n"); + exit (1); + } + + mpfr_set_prec (x, 21); + mpfr_set_prec (y, 8); + mpfr_set_ui (y, 120, GMP_RNDN); + mpfr_lngamma (x, y, GMP_RNDZ); + mpfr_set_prec (y, 21); + mpfr_set_str_binary (y, "0.111000101000001100101E9"); + if (mpfr_cmp (x, y)) + { + printf ("Error in mpfr_lngamma (120)\n"); + printf ("Expected "); mpfr_print_binary (y); puts (""); + printf ("Got "); mpfr_print_binary (x); puts (""); + exit (1); + } + + mpfr_set_prec (x, 3); + mpfr_set_prec (y, 206); + mpfr_set_str_binary (x, "0.110e10"); + inex = mpfr_lngamma (y, x, GMP_RNDN); + mpfr_set_prec (x, 206); + mpfr_set_str_binary (x, "0.10000111011000000011100010101001100110001110000111100011000100100110110010001011011110101001111011110110000001010100111011010000000011100110110101100111000111010011110010000100010111101010001101000110101001E13"); + if (mpfr_cmp (x, y)) + { + printf ("Error in mpfr_lngamma (768)\n"); + exit (1); + } + if (inex >= 0) + { + printf ("Wrong flag for mpfr_lngamma (768)\n"); + exit (1); + } + + mpfr_set_prec (x, 4); + mpfr_set_prec (y, 4); + mpfr_set_str_binary (x, "0.1100E-66"); + mpfr_lngamma (y, x, GMP_RNDN); + mpfr_set_str_binary (x, "0.1100E6"); + if (mpfr_cmp (x, y)) + { + printf ("Error for lngamma(0.1100E-66)\n"); + exit (1); + } + + mpfr_clear (x); + mpfr_clear (y); +} + +int +main (void) +{ + tests_start_mpfr (); + + special (); + test_generic (2, 100, 2); + + tests_end_mpfr (); + return 0; +} diff --git a/tests/tlog.c b/tests/tlog.c index 72d04676d..adf11917e 100644 --- a/tests/tlog.c +++ b/tests/tlog.c @@ -177,6 +177,7 @@ static void special (void) { mpfr_t x, y; + mp_exp_t emax = mpfr_get_emax (); mpfr_init2 (x, 53); mpfr_init2 (y, 53); @@ -213,6 +214,16 @@ special (void) test_log (y, x, GMP_RNDN); MPFR_ASSERTN(mpfr_nan_p (y)); + /* infinite loop when */ + set_emax (128); + mpfr_set_prec (x, 251); + mpfr_set_prec (y, 251); + mpfr_set_str_binary (x, "0.10010111000000000001101E8"); + /* x = 4947981/32768, log(x) ~ 5.017282... */ + test_log (y, x, GMP_RNDN); + + set_emax (emax); + mpfr_clear (x); mpfr_clear (y); } diff --git a/tests/trint.c b/tests/trint.c index d6d2de2ea..8b94e47d0 100644 --- a/tests/trint.c +++ b/tests/trint.c @@ -188,7 +188,7 @@ test_fct (double (*f)(double), int (*g)(), char *s, mp_rnd_t r) static void test_against_libc (void) { - mp_rnd_t r; + mp_rnd_t r = GMP_RNDN; #if HAVE_ROUND TEST_FCT (round); diff --git a/tests/tset_f.c b/tests/tset_f.c index 7c62c71f9..b96f4abbc 100644 --- a/tests/tset_f.c +++ b/tests/tset_f.c @@ -56,19 +56,21 @@ main (void) mpfr_set_str (u, "7.f10872b020c49ba5e353f7ced916872b020c49ba5e353f7ced916872b020c498@2", 16, GMP_RNDN); - mpf_set_str (y, "2033.033", 10); + mpf_set_str (y, "2033033E-3", 10); /* avoid 2033.033 which is + locale-sensitive */ mpfr_set_f (x, y, GMP_RNDN); if (mpfr_cmp (x, u)) { - printf ("mpfr_set_f failed for y=2033.033\n"); + printf ("mpfr_set_f failed for y=2033033E-3\n"); exit (1); } - mpf_set_str (y, "-2033.033", 10); + mpf_set_str (y, "-2033033E-3", 10); /* avoid -2033.033 which is + locale-sensitive */ mpfr_set_f (x, y, GMP_RNDN); mpfr_neg (u, u, GMP_RNDN); if (mpfr_cmp (x, u)) { - printf ("mpfr_set_f failed for y=-2033.033\n"); + printf ("mpfr_set_f failed for y=-2033033E-3\n"); exit (1); } diff --git a/tests/tstckintc.c b/tests/tstckintc.c index 356a0b191..383809cc4 100644 --- a/tests/tstckintc.c +++ b/tests/tstckintc.c @@ -22,6 +22,7 @@ MA 02110-1301, USA. */ #include <stdio.h> #include <stdlib.h> #include <limits.h> +#include <string.h> #include "mpfr-test.h" diff --git a/tests/tsub.c b/tests/tsub.c index abbbf5643..629527784 100644 --- a/tests/tsub.c +++ b/tests/tsub.c @@ -288,7 +288,7 @@ check_diverse (void) mpfr_set_prec (y, 10); mpfr_set_prec (z, 10); mpfr_set_ui (y, 0, GMP_RNDN); - mpfr_set_str_binary (z, "0.100000000000000000000100E15"); + mpfr_set_str_binary (z, "0.10001"); if (test_sub (x, y, z, GMP_RNDN) <= 0) { printf ("Wrong inexact flag in x=mpfr_sub(0,z) for prec(z)>prec(x)\n"); diff --git a/tests/tversion.c b/tests/tversion.c index cf8ed12ea..03a30fd53 100644 --- a/tests/tversion.c +++ b/tests/tversion.c @@ -36,12 +36,12 @@ main (void) version = mpfr_get_version (); if (strcmp (buffer, version) != 0) { - printf ("Incorrect version (%s vs %s)\n", buffer, version); + printf ("Incorrect version [1] (%s vs %s)\n", buffer, version); exit (1); } if (strcmp (MPFR_VERSION_STRING, version) != 0) { - printf ("Incorrect version (%s vs %s)\n", buffer, version); + printf ("Incorrect version [2] (%s vs %s)\n", buffer, version); exit (1); } return 0; @@ -22,6 +22,7 @@ MA 02110-1301, USA. */ #include <stdlib.h> #include <stdio.h> #include <time.h> +#include <limits.h> #define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" @@ -197,7 +198,7 @@ static double domeasure (mp_prec_t *threshold, fprintf (stderr, "Failed to measure function 1!\n"); abort (); } - *threshold = /*MPFR_PREC_MIN*/0; + *threshold = 1; t2 = speed_measure (func, &s); if (t2 == -1.0) { @@ -178,7 +178,7 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) MPFR_FAST_COMPUTE_IF_SMALL_INPUT (z, __gmpfr_one, err, 1, rnd_mode, ); } - d = precz + 11; + d = precz + MPFR_INT_CEIL_LOG2(precz) + 10; /* we want that s1 = s-1 is exact, i.e. we should have PREC(s1) >= EXP(s) */ dint = (mpfr_uexp_t) MPFR_GET_EXP (s); |