diff options
| -rw-r--r-- | TODO | 4 | ||||
| -rw-r--r-- | src/log.c | 7 | ||||
| -rw-r--r-- | tests/div.dat | 2 | ||||
| -rw-r--r-- | tests/log.dat | 11 |
4 files changed, 16 insertions, 8 deletions
@@ -1,3 +1,7 @@ +From Andreas Enge 27 June 2012: +implement over- and underflow resistant log(x^2+y^2) in log.c; +uncomment corresponding tests in log.dat + From Andreas Enge 31 August 2011: implement mul_karatsuba with three multiplications at precision around p, instead of two at precision 2*p and one at precision p @@ -28,7 +28,7 @@ mpc_log (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd){ mpfr_prec_t prec; int loops = 0; int re_cmp, im_cmp; - int inex_re, inex_im; + int inex_re, inex_im, inex; /* special values: NaN and infinities */ if (!mpc_fin_p (op)) { @@ -121,10 +121,13 @@ mpc_log (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd){ mpfr_set_prec (w, prec); /* w is rounded down */ - mpc_norm (w, op, GMP_RNDD); + inex = mpc_norm (w, op, GMP_RNDD); /* error 1 ulp */ MPC_ASSERT (!mpfr_inf_p (w)); /* FIXME: intermediate overflow; the logarithm may be representable */ + MPC_ASSERT (!(inex && mpfr_cmp_ui (w, 1) == 0)); + /* FIXME: this is like an underflow; the following call to log will + compute 0 instead of a positive result */ mpfr_log (w, w, GMP_RNDD); /* generic error of log: (2^(2 - exp(w)) + 1) ulp */ diff --git a/tests/div.dat b/tests/div.dat index 2cffc7e..dbdc2b1 100644 --- a/tests/div.dat +++ b/tests/div.dat @@ -1,6 +1,6 @@ # Data file for mpc_div. # -# Copyright (C) 2008, 2009, 2010, 2011 INRIA +# Copyright (C) 2008, 2009, 2010, 2011, 2012 INRIA # # This file is part of GNU MPC. # diff --git a/tests/log.dat b/tests/log.dat index 57e3b35..0f590d0 100644 --- a/tests/log.dat +++ b/tests/log.dat @@ -1,13 +1,13 @@ # Data test file for mpc_log. # -# Copyright (C) 2008, 2009, 2010 INRIA +# Copyright (C) 2008, 2009, 2010, 2012 INRIA # # This file is part of GNU MPC. # # GNU MPC 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 3 of the License, or (at your -#o ption) any later version. +# option) any later version. # # GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS @@ -170,6 +170,7 @@ # instead of the correct result. Since this may happen in other parts of the # library as well, we do not consider it a bug for the moment. # + + 53 0x58B90BFD4BCBFp-22 2 0x1p0 2 0x1p536870912 2 0x1p536870912 U U - -# log(-1 + i*eps) : infinite loop ? -# 0 0 2 0 2 0b11 2 -1 2 0x1p-1073741813 N N +# Due to intermediate "underflow" (see FIXME in log.c), the following example +# loops. +# log(-1 + i*eps) +# - - 2 0 2 3 2 -1 2 0x1p-1073741813 N N |