diff options
author | Philippe Theveny <philippe.theveny@laposte.net> | 2014-01-20 17:05:47 +0000 |
---|---|---|
committer | Philippe Theveny <philippe.theveny@laposte.net> | 2014-01-20 17:05:47 +0000 |
commit | 6bf41d383facbfc69ced85ad17a28f118b764f3a (patch) | |
tree | df89dc30217afed9fb19db5063ae22983934b89c | |
parent | 4167cb23a2b6f74bb4acff38f80e9c777db5571e (diff) | |
download | mpc-git-benchs_tests.tar.gz |
Synchronize with r1413 from trunk.benchs_tests
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/mpc/branches/benchs_tests@1414 211d60ee-9f03-0410-a15a-8952a2c7a4e4
-rw-r--r-- | NEWS | 12 | ||||
-rw-r--r-- | TODO | 5 | ||||
-rw-r--r-- | doc/algorithms.bib | 8 | ||||
-rw-r--r-- | doc/algorithms.tex | 6 | ||||
-rw-r--r-- | m4/mpc.m4 | 3 | ||||
-rw-r--r-- | src/acos.c | 5 | ||||
-rw-r--r-- | src/asin.c | 75 | ||||
-rw-r--r-- | src/sin_cos.c | 24 | ||||
-rw-r--r-- | tests/Makefile.am | 5 | ||||
-rw-r--r-- | tests/asin.dat | 2 | ||||
-rw-r--r-- | tests/tio_str.c | 6 | ||||
-rw-r--r-- | tests/tpow_ld.c | 5 |
12 files changed, 134 insertions, 22 deletions
@@ -1,11 +1,13 @@ Recent changes in the trunk: - Minimally required library version: mpfr 3.0.0 + - Improved speed for corner cases of mpc_asin, mpc_sin, see + http://lists.gforge.inria.fr/pipermail/mpc-discuss/2013-December/001266.html Changes in version 1.0.2: - - Fixed mpc_atan, mpc_atanh for (+-0, +-1), see - http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57994#c7 - - Fixed mpc_log10 for purely imaginary argument, see - http://lists.gforge.inria.fr/pipermail/mpc-discuss/2012-September/001208.html + - Fixed mpc_atan, mpc_atanh for (+-0, +-1), see + http://gcc.gnu.org/bugzilla/show_bug.cgi?id=57994#c7 + - Fixed mpc_log10 for purely imaginary argument, see + http://lists.gforge.inria.fr/pipermail/mpc-discuss/2012-September/001208.html Changes in version 1.0.1: - Switched to automake 1.11.6, see @@ -149,7 +151,7 @@ Changes in version 0.5 ("Aconitum neomontanum"): - mpc_sqrt with directed rounding -Copyright (C) 2008, 2009, 2010, 2011, 2012 INRIA +Copyright (C) 2008, 2009, 2010, 2011, 2012, 2013, 2014 INRIA Copying and distribution of this file, with or without modification, are permitted in any medium without royalty provided the copyright @@ -1,3 +1,7 @@ +From Karim Belabas 9 Jan 2014: +Implement Hurwitz(s,x) -> gives Zeta for x=1. +Cf http://arxiv.org/abs/1309.2877 + From Andreas Enge 27 August 2012: Implement im(atan(x+i*y)) as 1/4 * [log1p (4y / (x^2 +(1-y)^2))] @@ -33,6 +37,7 @@ sage: %timeit asin(MPComplexField()(1,1)) 10000 loops, best of 3: 83.7 us per loop sage: %timeit asin(MPComplexField()(1,1e-1000)) 100 loops, best of 3: 17 ms per loop +-> should be much faster with revision 1402 (check) Same for acos: sage: %timeit acos(MPComplexField()(1,1)) diff --git a/doc/algorithms.bib b/doc/algorithms.bib index 2232a0e..eb4855e 100644 --- a/doc/algorithms.bib +++ b/doc/algorithms.bib @@ -36,3 +36,11 @@ number = 4, pages = {359--367} } + +@Misc{CrTh10, + author = {John E. Cremona and Thotsaphon Thongjunthug}, + title = {The complex {AGM}, periods of elliptic curves over {\C} and complex elliptic logarithms}, + howpublished = {\url{http://arxiv.org/abs/1011.0914}}, + year = 2010, + note = {31 pages}} + diff --git a/doc/algorithms.tex b/doc/algorithms.tex index 669d651..d9acd15 100644 --- a/doc/algorithms.tex +++ b/doc/algorithms.tex @@ -17,6 +17,7 @@ \newcommand {\Ulp}{{\operatorname {ulp}}} \DeclareMathOperator{\Exp}{\operatorname {Exp}} \newcommand {\atantwo}{\operatorname {atan2}} +\newcommand {\asin}{\operatorname {asin}} \newcommand{\error}{\operatorname {error}} \newcommand{\relerror}{\operatorname {relerror}} \newcommand{\Norm}{\operatorname {N}} @@ -1142,6 +1143,11 @@ k_R=\left\{ \right. \end{equation*} +\subsection {\texttt {mpc\_asin}} + +For $0 \leq y \leq 1$, we have $|\Re \asin (1 \pm iy) - \pi/2| +\leq y^{1/2}$, and $|\Im \asin (1 \pm iy) \mp y^{1/2}| \leq (1/12) y^{3/2}$. + \subsection {\texttt {mpc\_pow}} The main issue for the power function is to be able to recognize when the @@ -90,7 +90,7 @@ AC_DEFUN([MPC_C_CHECK_FLAG], [ AC_DEFUN([MPC_C_CHECK_WARNINGCFLAGS], [ AC_REQUIRE([AC_PROG_GREP]) if echo $VERSION | grep -c dev >/dev/null 2>&1 ; then - if test "x$GCC" = "xyes" -a "x$compiler" != "xicc" -a "x$compiler" != "xg++"; then + if test "x$GCC" = "xyes" -a "x$compiler" != "xicc"; then # enable -Werror for myself (Andreas Enge) if test "x$USER" = "xenge"; then MPC_C_CHECK_FLAG(-Werror) @@ -107,6 +107,7 @@ AC_DEFUN([MPC_C_CHECK_WARNINGCFLAGS], [ MPC_C_CHECK_FLAG(-Wstrict-prototypes) MPC_C_CHECK_FLAG(-Wmissing-prototypes) MPC_C_CHECK_FLAG(-Wno-unused-value) + MPC_C_CHECK_FLAG(-Wlogical-op) fi fi ]) @@ -24,7 +24,7 @@ along with this program. If not, see http://www.gnu.org/licenses/ . int mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { - int inex_re, inex_im, inex; + int inex_re, inex_im, inex, loop = 0; mpfr_prec_t p_re, p_im, p; mpc_t z1; mpfr_t pi_over_2; @@ -190,7 +190,8 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_init2 (pi_over_2, p); for (;;) { - p += mpc_ceil_log2 (p) + 3; + loop ++; + p += (loop <= 2) ? mpc_ceil_log2 (p) + 3 : p / 2; mpfr_set_prec (mpc_realref(z1), p); mpfr_set_prec (pi_over_2, p); @@ -1,6 +1,6 @@ /* mpc_asin -- arcsine of a complex number. -Copyright (C) 2009, 2010, 2011, 2012 INRIA +Copyright (C) 2009, 2010, 2011, 2012, 2013 INRIA This file is part of GNU MPC. @@ -20,13 +20,74 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #include "mpc-impl.h" +/* Special case op = 1 + i*y for tiny y (see algorithms.tex). + Return 0 if special formula fails, otherwise put in z1 the approximate + value which needs to be converted to rop. + z1 is a temporary variable with enough precision. + */ +static int +mpc_asin_special (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd, mpc_ptr z1) +{ + mpfr_exp_t ey = mpfr_get_exp (mpc_imagref (op)); + mpfr_t abs_y; + mpfr_prec_t p; + int inex; + + /* |Re(asin(1+i*y)) - pi/2| <= y^(1/2) */ + if (ey >= 0 || ((-ey) / 2 < mpfr_get_prec (mpc_realref (z1)))) + return 0; + if (mpfr_get_prec (mpc_imagref (z1)) > (-ey) + 2) /* p > -ey + 2 */ + return 0; + + mpfr_const_pi (mpc_realref (z1), MPFR_RNDN); + mpfr_div_2exp (mpc_realref (z1), mpc_realref (z1), 1, MPFR_RNDN); /* exact */ + p = mpfr_get_prec (mpc_realref (z1)); + /* if z1 has precision p, the error on z1 is 1/2*ulp(z1) = 2^(-p) so far, + and since ey <= -2p, then y^(1/2) <= 1/2*ulp(z1) too, thus the total + error is bounded by ulp(z1) */ + if (!mpfr_can_round (mpc_realref(z1), p, MPFR_RNDN, MPFR_RNDZ, + mpfr_get_prec (mpc_realref(rop)) + + (MPC_RND_RE(rnd) == MPFR_RNDN))) + return 0; + + /* |Im(asin(1+i*y)) - y^(1/2)| <= (1/12) * y^(3/2) for y >= 0 (err >= 0) + |Im(asin(1-i*y)) + y^(1/2)| <= (1/12) * y^(3/2) for y >= 0 (err <= 0) */ + abs_y[0] = mpc_imagref (op)[0]; + if (mpfr_signbit (mpc_imagref (op))) + MPFR_CHANGE_SIGN (abs_y); + inex = mpfr_sqrt (mpc_imagref (z1), abs_y, MPFR_RNDN); /* error <= 1/2 ulp */ + if (mpfr_signbit (mpc_imagref (op))) + MPFR_CHANGE_SIGN (mpc_imagref (z1)); + /* If z1 has precision p, the error on z1 is 1/2*ulp(z1) = 2^(-p) so far, + and (1/12) * y^(3/2) <= (1/8) * y * y^(1/2) <= 2^(ey-3)*2^p*ulp(y^(1/2)) + thus for p+ey-3 <= -1 we have (1/12) * y^(3/2) <= (1/2) * ulp(y^(1/2)), + and the total error is bounded by ulp(z1). + Note: if y^(1/2) is exactly representable, and ends with many zeroes, + then mpfr_can_round below will fail; however in that case we know that + Im(asin(1+i*y)) is away from +/-y^(1/2) wrt zero. */ + if (inex == 0) /* enlarge im(z1) so that the inexact flag is correct */ + { + if (mpfr_signbit (mpc_imagref (op))) + mpfr_nextbelow (mpc_imagref (z1)); + else + mpfr_nextabove (mpc_imagref (z1)); + return 1; + } + p = mpfr_get_prec (mpc_imagref (z1)); + if (!mpfr_can_round (mpc_imagref(z1), p - 1, MPFR_RNDA, MPFR_RNDZ, + mpfr_get_prec (mpc_imagref(rop)) + + (MPC_RND_IM(rnd) == MPFR_RNDN))) + return 0; + return 1; +} + int mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { - mpfr_prec_t p, p_re, p_im, incr_p = 0; + mpfr_prec_t p, p_re, p_im; mpfr_rnd_t rnd_re, rnd_im; mpc_t z1; - int inex; + int inex, loop = 0; /* special values */ if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op))) @@ -148,11 +209,15 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { mpfr_exp_t ex, ey, err; - p += mpc_ceil_log2 (p) + 3 + incr_p; /* incr_p is zero initially */ - incr_p = p / 2; + loop ++; + p += (loop <= 2) ? mpc_ceil_log2 (p) + 3 : p / 2; mpfr_set_prec (mpc_realref(z1), p); mpfr_set_prec (mpc_imagref(z1), p); + /* try special code for 1+i*y with tiny y */ + if (loop == 1 && mpc_asin_special (rop, op, rnd, z1)) + break; + /* z1 <- z^2 */ mpc_sqr (z1, op, MPC_RNDNN); /* err(x) <= 1/2 ulp(x), err(y) <= 1/2 ulp(y) */ diff --git a/src/sin_cos.c b/src/sin_cos.c index 0bb00cf..7051708 100644 --- a/src/sin_cos.c +++ b/src/sin_cos.c @@ -299,11 +299,28 @@ mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, mpfr_t s, c, sh, ch, sch, csh; mpfr_prec_t prec; int ok; - int inex_re, inex_im, inex_sin, inex_cos; + int inex_re, inex_im, inex_sin, inex_cos, loop = 0; prec = 2; if (rop_sin != NULL) - prec = MPC_MAX (prec, MPC_MAX_PREC (rop_sin)); + { + mp_prec_t er, ei; + prec = MPC_MAX (prec, MPC_MAX_PREC (rop_sin)); + /* since the Taylor expansion of sin(x) at x=0 is x - x^3/6 + O(x^5), + if x <= 2^(-p), then the second term/x is about 2^(-2p)/6, thus we + need at least 2p+3 bits of precision. This is true only when x is + exactly representable in the target precision. */ + if (MPC_MAX_PREC (op) <= prec) + { + er = mpfr_get_exp (mpc_realref (op)); + ei = mpfr_get_exp (mpc_imagref (op)); + /* consider the maximal exponent only */ + er = (er < ei) ? ei : er; + if (er < 0) + if (prec < 2 * (mp_prec_t) (-er) + 3) + prec = 2 * (mp_prec_t) (-er) + 3; + } + } if (rop_cos != NULL) prec = MPC_MAX (prec, MPC_MAX_PREC (rop_cos)); @@ -315,8 +332,9 @@ mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, mpfr_init2 (csh, 2); do { + loop ++; ok = 1; - prec += mpc_ceil_log2 (prec) + 5; + prec += (loop <= 2) ? mpc_ceil_log2 (prec) + 5 : prec / 2; mpfr_set_prec (s, prec); mpfr_set_prec (c, prec); diff --git a/tests/Makefile.am b/tests/Makefile.am index 77e5866..61e5583 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -25,8 +25,9 @@ AM_LDFLAGS = -no-install # LOADLIBES (documented in the "GNU make" manual and equivalent to LDLIBS) # enables to compile a program foo.c in the test directory by simply doing # "make foo". -LOADLIBES=$(DEFS) -I$(top_srcdir)/src -I$(top_builddir) $(CPPFLAGS) \ - $(CFLAGS) -L$(top_builddir)/tests/.libs -lmpc-tests -lmpc -lm $(LIBS) +LOADLIBES=$(DEFS) $(AM_CPPFLAGS) $(CPPFLAGS) $(CFLAGS) \ + $(top_builddir)/tests/.libs/libmpc-tests.a \ + $(top_builddir)/src/.libs/libmpc.a $(LIBS) check_PROGRAMS = tabs tacos tacosh tadd tadd_fr tadd_si tadd_ui targ \ tasin tasinh tatan tatanh tconj tcos tcosh tdiv tdiv_2si tdiv_2ui \ diff --git a/tests/asin.dat b/tests/asin.dat index 9793fcc..e93e8c5 100644 --- a/tests/asin.dat +++ b/tests/asin.dat @@ -124,3 +124,5 @@ + - 53 0x189BF9EC7FCD5Bp-54 53 0x1206ECFA94614Bp-50 53 17 53 42 N N - + 2 1.5 2 6 2 96 2 0x1p-8 N N - - 8 0xC9p-7 8 0x15p-2 2 96 2 0x1p-8 N N +- - 53 0x3243f6a8885a3p-49 53 0x1p-500 53 1 53 0x1p-1000 N N +- + 53 0x3243f6a8885a3p-49 53 -0x1p-500 53 1 53 -0x1p-1000 N N diff --git a/tests/tio_str.c b/tests/tio_str.c index b7daf9d..b4e0c34 100644 --- a/tests/tio_str.c +++ b/tests/tio_str.c @@ -1,6 +1,6 @@ /* tio_str-- Test file for mpc_inp_str and mpc_out_str. -Copyright (C) 2009, 2011, 2013 INRIA +Copyright (C) 2009, 2011, 2013, 2014 INRIA This file is part of GNU MPC. @@ -21,9 +21,9 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #include "mpc-tests.h" #ifdef HAVE_UNISTD_H -#if !defined _POSIX_C_SOURCE +#ifndef _POSIX_C_SOURCE #define _POSIX_C_SOURCE 1 /* apparently needed on Darwin */ -#endif /* _POSIX_C_SOURCE */ +#endif #include <unistd.h> /* for dup, dup2, STDIN_FILENO and STDOUT_FILENO */ #else #define STDIN_FILENO 0 diff --git a/tests/tpow_ld.c b/tests/tpow_ld.c index a02ea6e..a4e3ab9 100644 --- a/tests/tpow_ld.c +++ b/tests/tpow_ld.c @@ -1,6 +1,6 @@ /* tpow_ld -- test file for mpc_pow_ld. -Copyright (C) 2009 INRIA +Copyright (C) 2009, 2014 INRIA This file is part of GNU MPC. @@ -33,6 +33,9 @@ main (void) if (mpc_cmp_si_si (z, -9, 46) != 0) { printf ("Error for mpc_pow_ld (1)\n"); + printf ("expected (-9 46)\ngot "); + mpc_out_str (stdout, 10, 0, z, MPC_RNDNN); + printf ("\n"); exit (1); } mpc_clear (z); |