diff options
| -rw-r--r-- | BUGS | 9 | ||||
| -rw-r--r-- | NEWS | 1 | ||||
| -rw-r--r-- | TODO | 4 | ||||
| -rw-r--r-- | doc/mpc.texi | 2 | ||||
| -rw-r--r-- | src/mpc-impl.h | 1 | ||||
| -rw-r--r-- | src/norm.c | 149 | ||||
| -rw-r--r-- | tests/norm.dat | 80 |
7 files changed, 176 insertions, 70 deletions
@@ -1,8 +1,7 @@ - Overflow and underflow are not considered in all functions, which might lead to some infinite loops. -- Except for mpc_mul, - no checks are made for intermediate over- or underflows, which may occur - in the middle of the algorithms although the final result may be +- Also, checks are not always made for intermediate over- or underflows, which + may occur in the middle of the algorithms although the final result may be representable. For instance, if in the computation of Im (cos(x+i*y)) = -sin(x)*sinh(y) an overflow occurs in sinh(y), the value sign(Im (cos(x+i*y))) * inf @@ -10,3 +9,7 @@ sin(x)*sinh(y) is representable. If furthermore an underflow occurred in sin(x) (which has not been observed in practice), then the return value would be NaN*(+-inf)=NaN. + +Currently, the following functions have been rewritten to solve these +problems: +mul, norm @@ -1,4 +1,5 @@ Changes in version 1.0: + - First release as a GNU package - License change: LGPLv3+ for code, GFDLv1.3+ for documentation - Bug fixes: - mul now returns correct results even if there are over- or underflows @@ -3,7 +3,7 @@ return the exact ternary value in mpc_div: - for directed rounding modes it is easy (assuming we are able to detect exact cases) - rounding to nearest to n bits is equivalent to directed rounding to n+1 bits -Also for mpc_norm and mpc_abs. +Also for mpc_abs. From Andreas Enge 31 August 2011: implement mul_karatsuba with three multiplications at precision around p, @@ -16,6 +16,8 @@ rewrite fma more efficiently to work with lower (close to target) precision From Andreas Enge 30 August 2011: As soon as dependent on mpfr>=3, remove auxiliary functions from get_version.c and update mpc.h. +Use MPFR_RND? instead of GMP_RND?, and remove workarounds for MPFR_RNDA from +mpc-impl.h. Bench: - from Andreas Enge 9 June 2009: diff --git a/doc/mpc.texi b/doc/mpc.texi index 2388393..f8c6f42 100644 --- a/doc/mpc.texi +++ b/doc/mpc.texi @@ -909,8 +909,6 @@ The returned value is zero iff the result is exact. Set the floating-point number @var{rop} to the norm of @var{op} (i.e., the square of its absolute value), rounded in the direction @var{rnd}. -The returned value is zero iff the result is exact. -Note that the destination is of type @code{mpfr_t}, not @code{mpc_t}. @end deftypefun @deftypefun int mpc_mul_2exp (mpc_t @var{rop}, mpc_t @var{op1}, unsigned long int @var{op2}, mpc_rnd_t @var{rnd}) diff --git a/src/mpc-impl.h b/src/mpc-impl.h index dbd2a9c..baa7a72 100644 --- a/src/mpc-impl.h +++ b/src/mpc-impl.h @@ -32,6 +32,7 @@ along with this program. If not, see http://www.gnu.org/licenses/ . * Miscellaneous useful macros */ +#define MPC_MIN(h,i) ((h) < (i) ? (h) : (i)) #define MPC_MAX(h,i) ((h) > (i) ? (h) : (i)) /* Safe absolute value (to avoid possible integer overflow) */ @@ -1,6 +1,6 @@ /* mpc_norm -- Square of the norm of a complex number. -Copyright (C) 2002, 2005, 2008, 2009, 2010 INRIA +Copyright (C) 2002, 2005, 2008, 2009, 2010, 2011 INRIA This file is part of GNU MPC. @@ -24,67 +24,90 @@ along with this program. If not, see http://www.gnu.org/licenses/ . int mpc_norm (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd) { - mpfr_t u, v; - mpfr_prec_t prec; - int inexact, overflow, underflow; - - prec = mpfr_get_prec (a); - - /* handling of special values; consistent with abs in that - norm = abs^2; so norm (+-inf, nan) = norm (nan, +-inf) = +inf */ - if ( (mpfr_nan_p (MPC_RE (b)) || mpfr_nan_p (MPC_IM (b))) - || (mpfr_inf_p (MPC_RE (b)) || mpfr_inf_p (MPC_IM (b)))) - return mpc_abs (a, b, rnd); - - mpfr_init (u); - mpfr_init (v); - - if (!mpfr_zero_p(MPC_RE(b)) && !mpfr_zero_p(MPC_IM(b)) && - 2 * SAFE_ABS (mpfr_exp_t, - mpfr_get_exp (MPC_RE (b)) - mpfr_get_exp (MPC_IM (b))) - > (mpfr_exp_t)prec) - /* If real and imaginary part have very different magnitudes, then the */ - /* generic code increases the precision too much. Instead, compute the */ - /* squarings _exactly_. */ - { - mpfr_set_prec (u, 2 * MPC_PREC_RE (b)); - mpfr_set_prec (v, 2 * MPC_PREC_IM (b)); - mpfr_sqr (u, MPC_RE (b), GMP_RNDN); - mpfr_sqr (v, MPC_IM (b), GMP_RNDN); - inexact = mpfr_add (a, u, v, rnd); - } - else if (mpfr_zero_p(MPC_RE(b)) && mpfr_zero_p(MPC_IM(b))) - { - inexact = mpfr_set_ui (a, 0, rnd); - } - else - { - do - { - prec += mpc_ceil_log2 (prec) + 3; - - mpfr_set_prec (u, prec); - mpfr_set_prec (v, prec); - - inexact = mpfr_sqr (u, MPC_RE(b), GMP_RNDN); /* err<=1/2ulp */ - inexact |= mpfr_sqr (v, MPC_IM(b), GMP_RNDN); /* err<=1/2ulp*/ - inexact |= mpfr_add (u, u, v, GMP_RNDN); /* err <= 3/2 ulps */ - - overflow = mpfr_inf_p (u); - underflow = mpfr_zero_p (u); + int inexact; + + /* handling of special values; consistent with abs in that + norm = abs^2; so norm (+-inf, nan) = norm (nan, +-inf) = +inf */ + if (!mpc_fin_p (b)) + return mpc_abs (a, b, rnd); + else if (mpfr_zero_p (MPC_RE (b))) { + if (mpfr_zero_p (MPC_IM (b))) + return mpfr_set_ui (a, 0, rnd); /* +0 */ + else + return mpfr_sqr (a, MPC_IM (b), rnd); + } + else if (mpfr_zero_p (MPC_IM (b))) + return mpfr_sqr (a, MPC_RE (b), rnd); + + else /* everything finite and non-zero */ { + mpfr_t u, v, res; + mpfr_prec_t prec, prec_u, prec_v; + int loops; + const int max_loops = 2; + /* switch to exact squarings when loops==max_loops */ + + prec = mpfr_get_prec (a); + + mpfr_init (u); + mpfr_init (v); + mpfr_init (res); + + loops = 0; + mpfr_clear_underflow (); + mpfr_clear_overflow (); + do { + loops++; + prec += mpc_ceil_log2 (prec) + 3; + if (loops >= max_loops) { + prec_u = 2 * MPC_PREC_RE (b); + prec_v = 2 * MPC_PREC_IM (b); + } + else { + prec_u = MPC_MIN (prec, 2 * MPC_PREC_RE (b)); + prec_v = MPC_MIN (prec, 2 * MPC_PREC_IM (b)); + } + + mpfr_set_prec (u, prec_u); + mpfr_set_prec (v, prec_v); + mpfr_set_prec (res, prec); + + inexact = mpfr_sqr (u, MPC_RE(b), GMP_RNDD); /* err <= 1 ulp in prec */ + inexact |= mpfr_sqr (v, MPC_IM(b), GMP_RNDD); /* err <= 1 ulp in prec */ + if (loops < max_loops) + inexact |= mpfr_add (res, u, v, GMP_RNDD); /* err <= 3 ulp in prec */ + + } while (loops < max_loops + && inexact + && !mpfr_can_round (res, prec - 2, GMP_RNDD, GMP_RNDU, + mpfr_get_prec (a) + (rnd == GMP_RNDN))); + + if (mpfr_overflow_p ()) { + /* replace by "correctly rounded overflow" */ + mpfr_set_ui (a, 1ul, GMP_RNDN); + inexact = mpfr_mul_2ui (a, a, mpfr_get_emax (), rnd); } - while (!overflow && !underflow && inexact && - mpfr_can_round (u, prec - 2, GMP_RNDN, rnd, mpfr_get_prec (a)) == 0); - - inexact |= mpfr_set (a, u, rnd); - - if (overflow) - inexact = 1; - if (underflow) - inexact = -1; - } - mpfr_clear (u); - mpfr_clear (v); - - return inexact; + else if (mpfr_underflow_p ()) { + /* squarings were exact except for underflow */ + inexact = mpfr_add (a, u, v, rnd); + if (!inexact) { + if (rnd == GMP_RNDN || rnd == GMP_RNDD || rnd == GMP_RNDZ) + inexact = -1; + else { + mpfr_nextabove (a); + inexact = 1; + } + } + } + else if (loops == max_loops) + /* squarings were exact */ + inexact = mpfr_add (a, u, v, rnd); + else /* no problems, ternary value due to mpfr_can_round trick */ + inexact = mpfr_set (a, res, rnd); + + mpfr_clear (u); + mpfr_clear (v); + mpfr_clear (res); + } + + return inexact; } diff --git a/tests/norm.dat b/tests/norm.dat index 7a99a0f..5a6a2dd 100644 --- a/tests/norm.dat +++ b/tests/norm.dat @@ -1,6 +1,6 @@ # Data file for mpc_norm. # -# Copyright (C) 2008, 2010 INRIA +# Copyright (C) 2008, 2010, 2011 INRIA # # This file is part of GNU MPC. # @@ -77,3 +77,81 @@ # infinite loop reported by E. Thome - 250 +0 250 -0xf.fdda3457c3e69e5841461d505b42987feb42867a4a2d2872179c4efa20054c0@-136943039 250 -0xc.670d1beda685fdd771b6246e32ff49ec4fd70aec25367444e00933e6965d3c0@-136943040 N + +# inexact values: norm (2+i)=5, in the middle between two values at precision 2 +- 2 4 2 2 2 1 D +- 2 4 2 2 2 1 Z ++ 2 6 2 2 2 1 U +- 2 4 2 2 2 1 N + +# over- and underflows ++ 10 inf 10 0 10 0b1@536870912 N ++ 10 inf 10 0 10 0b1@536870912 U +- 10 0b1.111111111@1073741822 10 0 10 0b1@536870912 D +- 10 0b1.111111111@1073741822 10 0 10 0b1@536870912 Z +- 10 0 10 0 10 0b1@-536870913 N +- 10 0 10 0 10 0b1@-536870913 D +- 10 0 10 0 10 0b1@-536870913 Z ++ 10 0b1.000000000e-1073741824 10 0 10 0b1@-536870913 U + ++ 10 inf 10 0b1@536870912 10 0 N ++ 10 inf 10 0b1@536870912 10 0 U +- 10 0b1.111111111@1073741822 10 0b1@536870912 10 0 D +- 10 0b1.111111111@1073741822 10 0b1@536870912 10 0 Z +- 10 0 10 0b1@-536870913 10 0 N +- 10 0 10 0b1@-536870913 10 0 D +- 10 0 10 0b1@-536870913 10 0 Z ++ 10 0b1.000000000e-1073741824 10 0b1@-536870913 10 0 U + ++ 10 inf 10 1 10 0b1@536870912 N ++ 10 inf 10 1 10 0b1@536870912 U +- 10 0b1.111111111@1073741822 10 1 10 0b1@536870912 D +- 10 0b1.111111111@1073741822 10 1 10 0b1@536870912 Z +- 10 1 10 1 10 0b1@-536870913 N +- 10 1 10 1 10 0b1@-536870913 D +- 10 1 10 1 10 0b1@-536870913 Z ++ 10 0b1.000000001 10 1 10 0b1@-536870913 U + ++ 10 inf 10 0b1@536870912 10 1 N ++ 10 inf 10 0b1@536870912 10 1 U +- 10 0b1.111111111@1073741822 10 0b1@536870912 10 1 D +- 10 0b1.111111111@1073741822 10 0b1@536870912 10 1 Z +- 10 1 10 0b1@-536870913 10 1 N +- 10 1 10 0b1@-536870913 10 1 D +- 10 1 10 0b1@-536870913 10 1 Z ++ 10 0b1.000000001 10 0b1@-536870913 10 1 U + ++ 3 inf 10 0b1.1 10 0b1@536870912 N ++ 3 inf 10 0b1.1 10 0b1@536870912 U +- 3 0b1.11@1073741822 10 0b1.1 10 0b1@536870912 D +- 3 0b1.11@1073741822 10 0b1.1 10 0b1@536870912 Z +- 3 2 10 0b1.1 10 0b1@-536870913 N +- 3 2 10 0b1.1 10 0b1@-536870913 D +- 3 2 10 0b1.1 10 0b1@-536870913 Z ++ 3 0b1.01@1 10 0b1.1 10 0b1@-536870913 U + ++ 3 inf 10 0b1@536870912 10 0b1.1 N ++ 3 inf 10 0b1@536870912 10 0b1.1 U +- 3 0b1.11@1073741822 10 0b1@536870912 10 0b1.1 D +- 3 0b1.11@1073741822 10 0b1@536870912 10 0b1.1 Z +- 3 2 10 0b1@-536870913 10 0b1.1 N +- 3 2 10 0b1@-536870913 10 0b1.1 D +- 3 2 10 0b1@-536870913 10 0b1.1 Z ++ 3 0b1.01@1 10 0b1@-536870913 10 0b1.1 U + ++ 10 inf 10 0b1@-536870913 10 0b1@536870912 N ++ 10 inf 10 0b1@-536870913 10 0b1@536870912 U +- 10 0b1.111111111@1073741822 10 0b1@-536870913 10 0b1@536870912 D +- 10 0b1.111111111@1073741822 10 0b1@-536870912 10 0b1@536870912 Z ++ 10 inf 10 0b1@536870912 10 0b1@-536870913 N ++ 10 inf 10 0b1@536870912 10 0b1@-536870913 U +- 10 0b1.111111111@1073741822 10 0b1@536870912 10 0b1@-536870913 D +- 10 0b1.111111111@1073741822 10 0b1@536870912 10 0b1@-536870913 Z ++ 10 inf 10 0b1@536870912 10 0b1@536870912 N ++ 10 inf 10 0b1@536870912 10 0b1@536870912 U +- 10 0b1.111111111@1073741822 10 0b1@536870912 10 0b1@536870912 D +- 10 0b1.111111111@1073741822 10 0b1@536870912 10 0b1@536870912 Z +- 10 0 10 0b1@-536870913 10 0b1@-536870913 N +- 10 0 10 0b1@-536870913 10 0b1@-536870913 D +- 10 0 10 0b1@-536870913 10 0b1@-536870913 Z ++ 10 0b1@-1073741824 10 0b1@-536870913 10 0b1@-536870913 U |