X-Git-Url: http://review.tizen.org/git/?a=blobdiff_plain;f=src%2Fsqr.c;h=3cef847e9a4fbcff8e04af039a05d2d806cf6e8a;hb=f34506723f150a03d5e6e389820bd369b89e30b0;hp=f1ce1bac8a40f814ad97c739a90922066ba50fa9;hpb=c9266275a7240e07da227b662cc6afcb0000d7a6;p=platform%2Fupstream%2Fmpc.git diff --git a/src/sqr.c b/src/sqr.c index f1ce1ba..3cef847 100644 --- a/src/sqr.c +++ b/src/sqr.c @@ -1,324 +1,253 @@ /* mpc_sqr -- Square a complex number. -Copyright (C) 2002, 2005, 2008, 2009, 2010, 2011, 2012 INRIA +Copyright (C) 2002, 2005, 2008, 2009 Andreas Enge, Paul Zimmermann, Philippe Th\'eveny -This file is part of GNU MPC. +This file is part of the MPC Library. -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 +The MPC 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. -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 -FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for -more details. +The MPC 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 this program. If not, see http://www.gnu.org/licenses/ . -*/ +along with the MPC Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ #include /* for MPC_ASSERT */ #include "mpc-impl.h" - -static int -mpfr_fsss (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr c, mpfr_rnd_t rnd) -{ - /* Computes z = a^2 - c^2. - Assumes that a and c are finite and non-zero; so a squaring yielding - an infinity is an overflow, and a squaring yielding 0 is an underflow. - Assumes further that z is distinct from a and c. */ - - int inex; - mpfr_t u, v; - - /* u=a^2, v=c^2 exactly */ - mpfr_init2 (u, 2*mpfr_get_prec (a)); - mpfr_init2 (v, 2*mpfr_get_prec (c)); - mpfr_sqr (u, a, GMP_RNDN); - mpfr_sqr (v, c, GMP_RNDN); - - /* tentatively compute z as u-v; here we need z to be distinct - from a and c to not lose the latter */ - inex = mpfr_sub (z, u, v, rnd); - - if (mpfr_inf_p (z)) { - /* replace by "correctly rounded overflow" */ - mpfr_set_si (z, (mpfr_signbit (z) ? -1 : 1), GMP_RNDN); - inex = mpfr_mul_2ui (z, z, mpfr_get_emax (), rnd); - } - else if (mpfr_zero_p (u) && !mpfr_zero_p (v)) { - /* exactly u underflowed, determine inexact flag */ - inex = (mpfr_signbit (u) ? 1 : -1); - } - else if (mpfr_zero_p (v) && !mpfr_zero_p (u)) { - /* exactly v underflowed, determine inexact flag */ - inex = (mpfr_signbit (v) ? -1 : 1); - } - else if (mpfr_nan_p (z) || (mpfr_zero_p (u) && mpfr_zero_p (v))) { - /* In the first case, u and v are +inf. - In the second case, u and v are zeroes; their difference may be 0 - or the least representable number, with a sign to be determined. - Redo the computations with mpz_t exponents */ - mpfr_exp_t ea, ec; - mpz_t eu, ev; - /* cheat to work around the const qualifiers */ - - /* Normalise the input by shifting and keep track of the shifts in - the exponents of u and v */ - ea = mpfr_get_exp (a); - ec = mpfr_get_exp (c); - - mpfr_set_exp ((mpfr_ptr) a, (mpfr_prec_t) 0); - mpfr_set_exp ((mpfr_ptr) c, (mpfr_prec_t) 0); - - mpz_init (eu); - mpz_init (ev); - mpz_set_si (eu, (long int) ea); - mpz_mul_2exp (eu, eu, 1); - mpz_set_si (ev, (long int) ec); - mpz_mul_2exp (ev, ev, 1); - - /* recompute u and v and move exponents to eu and ev */ - mpfr_sqr (u, a, GMP_RNDN); - /* exponent of u is non-positive */ - mpz_sub_ui (eu, eu, (unsigned long int) (-mpfr_get_exp (u))); - mpfr_set_exp (u, (mpfr_prec_t) 0); - mpfr_sqr (v, c, GMP_RNDN); - mpz_sub_ui (ev, ev, (unsigned long int) (-mpfr_get_exp (v))); - mpfr_set_exp (v, (mpfr_prec_t) 0); - if (mpfr_nan_p (z)) { - mpfr_exp_t emax = mpfr_get_emax (); - int overflow; - /* We have a = ma * 2^ea with 1/2 <= |ma| < 1 and ea <= emax. - So eu <= 2*emax, and eu > emax since we have - an overflow. The same holds for ev. Shift u and v by as much as - possible so that one of them has exponent emax and the - remaining exponents in eu and ev are the same. Then carry out - the addition. Shifting u and v prevents an underflow. */ - if (mpz_cmp (eu, ev) >= 0) { - mpfr_set_exp (u, emax); - mpz_sub_ui (eu, eu, (long int) emax); - mpz_sub (ev, ev, eu); - mpfr_set_exp (v, (mpfr_exp_t) mpz_get_ui (ev)); - /* remaining common exponent is now in eu */ - } - else { - mpfr_set_exp (v, emax); - mpz_sub_ui (ev, ev, (long int) emax); - mpz_sub (eu, eu, ev); - mpfr_set_exp (u, (mpfr_exp_t) mpz_get_ui (eu)); - mpz_set (eu, ev); - /* remaining common exponent is now also in eu */ - } - inex = mpfr_sub (z, u, v, rnd); - /* Result is finite since u and v have the same sign. */ - overflow = mpfr_mul_2ui (z, z, mpz_get_ui (eu), rnd); - if (overflow) - inex = overflow; - } - else { - int underflow; - /* Subtraction of two zeroes. We have a = ma * 2^ea - with 1/2 <= |ma| < 1 and ea >= emin and similarly for b. - So 2*emin < 2*emin+1 <= eu < emin < 0, and analogously for v. */ - mpfr_exp_t emin = mpfr_get_emin (); - if (mpz_cmp (eu, ev) <= 0) { - mpfr_set_exp (u, emin); - mpz_add_ui (eu, eu, (unsigned long int) (-emin)); - mpz_sub (ev, ev, eu); - mpfr_set_exp (v, (mpfr_exp_t) mpz_get_si (ev)); - } - else { - mpfr_set_exp (v, emin); - mpz_add_ui (ev, ev, (unsigned long int) (-emin)); - mpz_sub (eu, eu, ev); - mpfr_set_exp (u, (mpfr_exp_t) mpz_get_si (eu)); - mpz_set (eu, ev); - } - inex = mpfr_sub (z, u, v, rnd); - mpz_neg (eu, eu); - underflow = mpfr_div_2ui (z, z, mpz_get_ui (eu), rnd); - if (underflow) - inex = underflow; - } - - mpz_clear (eu); - mpz_clear (ev); - - mpfr_set_exp ((mpfr_ptr) a, ea); - mpfr_set_exp ((mpfr_ptr) c, ec); - /* works also when a == c */ - } - - mpfr_clear (u); - mpfr_clear (v); - - return inex; -} - - int mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { int ok; mpfr_t u, v; mpfr_t x; - /* temporary variable to hold the real part of op, + /* rop temporary variable to hold the real part of op, needed in the case rop==op */ - mpfr_prec_t prec; + mp_prec_t prec; int inex_re, inex_im, inexact; - mpfr_exp_t emin; - int saved_underflow; + mp_exp_t old_emax, old_emin, emin, emax; /* special values: NaN and infinities */ - if (!mpc_fin_p (op)) { - if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op))) { - mpfr_set_nan (mpc_realref (rop)); - mpfr_set_nan (mpc_imagref (rop)); + if (!mpfr_number_p (MPC_RE (op)) || !mpfr_number_p (MPC_IM (op))) { + if (mpfr_nan_p (MPC_RE (op)) || mpfr_nan_p (MPC_IM (op))) { + mpfr_set_nan (MPC_RE (rop)); + mpfr_set_nan (MPC_IM (rop)); } - else if (mpfr_inf_p (mpc_realref (op))) { - if (mpfr_inf_p (mpc_imagref (op))) { - mpfr_set_inf (mpc_imagref (rop), - MPFR_SIGN (mpc_realref (op)) * MPFR_SIGN (mpc_imagref (op))); - mpfr_set_nan (mpc_realref (rop)); + else if (mpfr_inf_p (MPC_RE (op))) { + if (mpfr_inf_p (MPC_IM (op))) { + mpfr_set_nan (MPC_RE (rop)); + mpfr_set_inf (MPC_IM (rop), + MPFR_SIGN (MPC_RE (op)) * MPFR_SIGN (MPC_IM (op))); } else { - if (mpfr_zero_p (mpc_imagref (op))) - mpfr_set_nan (mpc_imagref (rop)); + mpfr_set_inf (MPC_RE (rop), +1); + if (mpfr_zero_p (MPC_IM (op))) + mpfr_set_nan (MPC_IM (rop)); else - mpfr_set_inf (mpc_imagref (rop), - MPFR_SIGN (mpc_realref (op)) * MPFR_SIGN (mpc_imagref (op))); - mpfr_set_inf (mpc_realref (rop), +1); + mpfr_set_inf (MPC_IM (rop), + MPFR_SIGN (MPC_RE (op)) * MPFR_SIGN (MPC_IM (op))); } } else /* IM(op) is infinity, RE(op) is not */ { - if (mpfr_zero_p (mpc_realref (op))) - mpfr_set_nan (mpc_imagref (rop)); + mpfr_set_inf (MPC_RE (rop), -1); + if (mpfr_zero_p (MPC_RE (op))) + mpfr_set_nan (MPC_IM (rop)); else - mpfr_set_inf (mpc_imagref (rop), - MPFR_SIGN (mpc_realref (op)) * MPFR_SIGN (mpc_imagref (op))); - mpfr_set_inf (mpc_realref (rop), -1); + mpfr_set_inf (MPC_IM (rop), + MPFR_SIGN (MPC_RE (op)) * MPFR_SIGN (MPC_IM (op))); } return MPC_INEX (0, 0); /* exact */ } prec = MPC_MAX_PREC(rop); - /* Check for real resp. purely imaginary number */ - if (mpfr_zero_p (mpc_imagref(op))) { - int same_sign = mpfr_signbit (mpc_realref (op)) == mpfr_signbit (mpc_imagref (op)); - inex_re = mpfr_sqr (mpc_realref(rop), mpc_realref(op), MPC_RND_RE(rnd)); - inex_im = mpfr_set_ui (mpc_imagref(rop), 0ul, GMP_RNDN); + /* first check for real resp. purely imaginary number */ + if (mpfr_zero_p (MPC_IM(op))) + { + int same_sign = mpfr_signbit (MPC_RE (op)) == mpfr_signbit (MPC_IM (op)); + inex_re = mpfr_sqr (MPC_RE(rop), MPC_RE(op), MPC_RND_RE(rnd)); + inex_im = mpfr_set_ui (MPC_IM(rop), 0ul, GMP_RNDN); if (!same_sign) mpc_conj (rop, rop, MPC_RNDNN); return MPC_INEX(inex_re, inex_im); } - if (mpfr_zero_p (mpc_realref(op))) { - int same_sign = mpfr_signbit (mpc_realref (op)) == mpfr_signbit (mpc_imagref (op)); - inex_re = -mpfr_sqr (mpc_realref(rop), mpc_imagref(op), INV_RND (MPC_RND_RE(rnd))); - mpfr_neg (mpc_realref(rop), mpc_realref(rop), GMP_RNDN); - inex_im = mpfr_set_ui (mpc_imagref(rop), 0ul, GMP_RNDN); + if (mpfr_zero_p (MPC_RE(op))) + { + int same_sign = mpfr_signbit (MPC_RE (op)) == mpfr_signbit (MPC_IM (op)); + inex_re = -mpfr_sqr (MPC_RE(rop), MPC_IM(op), INV_RND (MPC_RND_RE(rnd))); + mpfr_neg (MPC_RE(rop), MPC_RE(rop), GMP_RNDN); + inex_im = mpfr_set_ui (MPC_IM(rop), 0ul, GMP_RNDN); if (!same_sign) mpc_conj (rop, rop, MPC_RNDNN); return MPC_INEX(inex_re, inex_im); } + /* If the real and imaginary part of the argument have rop very different */ + /* exponent, it is not reasonable to use Karatsuba squaring; compute */ + /* exactly with the standard formulae instead, even if this means an */ + /* additional multiplication. */ + if (SAFE_ABS (mp_exp_t, MPFR_EXP (MPC_RE (op)) - MPFR_EXP (MPC_IM (op))) + > (mp_exp_t) MPC_MAX_PREC (op) / 2) + { + mpfr_init2 (u, 2*MPFR_PREC (MPC_RE (op))); + mpfr_init2 (v, 2*MPFR_PREC (MPC_IM (op))); + + mpfr_sqr (u, MPC_RE (op), GMP_RNDN); + mpfr_sqr (v, MPC_IM (op), GMP_RNDN); /* both are exact */ + inex_im = mpfr_mul (rop->im, op->re, op->im, MPC_RND_IM (rnd)); + mpfr_mul_2exp (rop->im, rop->im, 1, GMP_RNDN); + inex_re = mpfr_sub (rop->re, u, v, MPC_RND_RE (rnd)); + + mpfr_clear (u); + mpfr_clear (v); + return MPC_INEX (inex_re, inex_im); + } + + + mpfr_init (u); + mpfr_init (v); if (rop == op) { - mpfr_init2 (x, MPC_PREC_RE (op)); + mpfr_init2 (x, MPFR_PREC (op->re)); mpfr_set (x, op->re, GMP_RNDN); } else x [0] = op->re [0]; - /* From here on, use x instead of op->re and safely overwrite rop->re. */ - - /* Compute real part of result. */ - if (SAFE_ABS (mpfr_exp_t, - mpfr_get_exp (mpc_realref (op)) - mpfr_get_exp (mpc_imagref (op))) - > (mpfr_exp_t) MPC_MAX_PREC (op) / 2) { - /* If the real and imaginary parts of the argument have very different - exponents, it is not reasonable to use Karatsuba squaring; compute - exactly with the standard formulae instead, even if this means an - additional multiplication. Using the approach copied from mul, over- - and underflows are also handled correctly. */ - - inex_re = mpfr_fsss (rop->re, x, op->im, MPC_RND_RE (rnd)); - } - else { - /* Karatsuba squaring: we compute the real part as (x+y)*(x-y) and the - imaginary part as 2*x*y, with a total of 2M instead of 2S+1M for the - naive algorithm, which computes x^2-y^2 and 2*y*y */ - mpfr_init (u); - mpfr_init (v); - emin = mpfr_get_emin (); + /* store the maximal exponent */ + old_emax = mpfr_get_emax (); + old_emin = mpfr_get_emin (); + mpfr_set_emax (emax = mpfr_get_emax_max ()); + mpfr_set_emin (emin = mpfr_get_emin_min ()); - do + do + { + prec += mpc_ceil_log2 (prec) + 5; + + mpfr_set_prec (u, prec); + mpfr_set_prec (v, prec); + + /* Let op = x + iy. We need u = x+y and v = x-y, rounded away. */ + /* As this is not implemented in mpfr, we round towards zero and */ + /* add one ulp if the result is not exact. The error is bounded */ + /* above by 1 ulp. */ + /* We first let inexact be 1 if the real part is not computed */ + /* exactly and determine the sign later. */ + inexact = 0; + if (mpfr_add (u, x, MPC_IM (op), GMP_RNDZ)) + /* we have to use x here instead of MPC_RE (op), as MPC_RE (rop) + may be modified later in the attempt to assign u to it */ { - prec += mpc_ceil_log2 (prec) + 5; - - mpfr_set_prec (u, prec); - mpfr_set_prec (v, prec); - - /* Let op = x + iy. We need u = x+y and v = x-y, rounded away. */ - /* The error is bounded above by 1 ulp. */ - /* We first let inexact be 1 if the real part is not computed */ - /* exactly and determine the sign later. */ - inexact = ROUND_AWAY (mpfr_add (u, x, mpc_imagref (op), MPFR_RNDA), u) - | ROUND_AWAY (mpfr_sub (v, x, mpc_imagref (op), MPFR_RNDA), v); - - /* compute the real part as u*v, rounded away */ - /* determine also the sign of inex_re */ + mpfr_add_one_ulp (u, GMP_RNDN); + inexact = 1; + } + if (mpfr_sub (v, x, MPC_IM (op), GMP_RNDZ)) + { + mpfr_add_one_ulp (v, GMP_RNDN); + inexact = 1; + } - if (mpfr_sgn (u) == 0 || mpfr_sgn (v) == 0) { - /* as we have rounded away, the result is exact */ - mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); - inex_re = 0; - ok = 1; - } - else { - mpfr_rnd_t rnd_away; - /* FIXME: can be replaced by MPFR_RNDA in mpfr >= 3 */ - rnd_away = (mpfr_sgn (u) * mpfr_sgn (v) > 0 ? GMP_RNDU : GMP_RNDD); - inexact |= ROUND_AWAY (mpfr_mul (u, u, v, MPFR_RNDA), u); /* error 5 */ - if (mpfr_get_exp (u) == emin || mpfr_inf_p (u)) { - /* under- or overflow */ - inex_re = mpfr_fsss (rop->re, x, op->im, MPC_RND_RE (rnd)); - ok = 1; + /* compute the real part as u*v, rounded away */ + /* determine also the sign of inex_re */ + if (mpfr_sgn (u) == 0 || mpfr_sgn (v) == 0) + { + /* as we have rounded away, the result is exact */ + mpfr_set_ui (MPC_RE (rop), 0, GMP_RNDN); + inex_re = 0; + ok = 1; + } + else if (mpfr_sgn (u) * mpfr_sgn (v) > 0) + { + inexact |= mpfr_mul (u, u, v, GMP_RNDU); /* error 5 */ + /* checks that no overflow nor underflow occurs: since u*v > 0 + and we round up, an overflow will give +Inf, but an underflow + will give 0.5*2^emin */ + MPC_ASSERT (mpfr_get_exp (u) != emin); + if (mpfr_inf_p (u)) + { + mp_rnd_t rnd_re = MPC_RND_RE (rnd); + if (rnd_re == GMP_RNDZ || rnd_re == GMP_RNDD) + { + mpfr_set_ui_2exp (MPC_RE (rop), 1, emax, rnd_re); + inex_re = -1; + } + else /* round up or away from zero */ + inex_re = 1; + break; } - else { - ok = (!inexact) | mpfr_can_round (u, prec - 3, - rnd_away, GMP_RNDZ, - MPC_PREC_RE (rop) + (MPC_RND_RE (rnd) == GMP_RNDN)); - if (ok) { - inex_re = mpfr_set (mpc_realref (rop), u, MPC_RND_RE (rnd)); - if (inex_re == 0) - /* remember that u was already rounded */ - inex_re = inexact; - } + ok = (!inexact) | mpfr_can_round (u, prec - 3, GMP_RNDU, GMP_RNDZ, + MPFR_PREC (MPC_RE (rop)) + (MPC_RND_RE (rnd) == GMP_RNDN)); + if (ok) + { + inex_re = mpfr_set (MPC_RE (rop), u, MPC_RND_RE (rnd)); + if (inex_re == 0) + /* remember that u was already rounded */ + inex_re = inexact; } - } - } - while (!ok); - - mpfr_clear (u); - mpfr_clear (v); + } + else + { + inexact |= mpfr_mul (u, u, v, GMP_RNDD); /* error 5 */ + /* checks that no overflow occurs: since u*v < 0 and we round down, + an overflow will give -Inf */ + MPC_ASSERT (mpfr_inf_p (u) == 0); + /* if an underflow happens (i.e., u = -0.5*2^emin since we round + away from zero), the result will be an underflow */ + if (mpfr_get_exp (u) == emin) + { + mp_rnd_t rnd_re = MPC_RND_RE (rnd); + if (rnd_re == GMP_RNDZ || rnd_re == GMP_RNDN || + rnd_re == GMP_RNDU) + { + mpfr_set_ui (MPC_RE (rop), 0, rnd_re); + inex_re = 1; + } + else /* round down or away from zero */ + inex_re = -1; + break; + } + ok = (!inexact) | mpfr_can_round (u, prec - 3, GMP_RNDD, GMP_RNDZ, + MPFR_PREC (MPC_RE (rop)) + (MPC_RND_RE (rnd) == GMP_RNDN)); + if (ok) + { + inex_re = mpfr_set (MPC_RE (rop), u, MPC_RND_RE (rnd)); + if (inex_re == 0) + inex_re = inexact; + } + } } + while (!ok); + + /* compute the imaginary part as 2*x*y, which is always possible */ + if (mpfr_get_exp (x) + mpfr_get_exp(MPC_IM (op)) <= emin + 1) + { + mpfr_mul_2ui (x, x, 1, GMP_RNDN); + inex_im = mpfr_mul (MPC_IM (rop), x, MPC_IM (op), MPC_RND_IM (rnd)); + } + else + { + inex_im = mpfr_mul (MPC_IM (rop), x, MPC_IM (op), MPC_RND_IM (rnd)); + /* checks that no underflow occurs (an overflow can occur here) */ + MPC_ASSERT (mpfr_zero_p (MPC_IM (rop)) == 0); + mpfr_mul_2ui (MPC_IM (rop), MPC_IM (rop), 1, GMP_RNDN); + } - saved_underflow = mpfr_underflow_p (); - mpfr_clear_underflow (); - inex_im = mpfr_mul (rop->im, x, op->im, MPC_RND_IM (rnd)); - if (!mpfr_underflow_p ()) - inex_im |= mpfr_mul_2ui (rop->im, rop->im, 1, MPC_RND_IM (rnd)); - /* We must not multiply by 2 if rop->im has been set to the smallest - representable number. */ - if (saved_underflow) - mpfr_set_underflow (); + mpfr_clear (u); + mpfr_clear (v); if (rop == op) mpfr_clear (x); + /* restore the exponent range */ + mpfr_set_emax (old_emax); + mpfr_set_emin (old_emin); + return MPC_INEX (inex_re, inex_im); }