diff options
Diffstat (limited to 'src/sqr.c')
| -rw-r--r-- | src/sqr.c | 146 |
1 files changed, 0 insertions, 146 deletions
@@ -21,152 +21,6 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #include <stdio.h> /* for MPC_ASSERT */ #include "mpc-impl.h" -#define mpz_add_si(z,x,y) do { \ - if (y >= 0) \ - mpz_add_ui (z, x, (long int) y); \ - else \ - mpz_sub_ui (z, x, (long int) (-y)); \ - } while (0); - -static int -mpfr_fmma (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, - mpfr_srcptr d, int sign, mpfr_rnd_t rnd) -{ - /* Computes z = ab+cd if sign >= 0, or z = ab-cd if sign < 0. - Assumes that a, b, c, d are finite and non-zero; so any multiplication - of two of them yielding an infinity is an overflow, and a - multiplication yielding 0 is an underflow. - Assumes further that z is distinct from a, b, c, d. */ - - int inex; - mpfr_t u, v; - - /* u=a*b, v=sign*c*d exactly */ - mpfr_init2 (u, mpfr_get_prec (a) + mpfr_get_prec (b)); - mpfr_init2 (v, mpfr_get_prec (c) + mpfr_get_prec (d)); - mpfr_mul (u, a, b, GMP_RNDN); - mpfr_mul (v, c, d, GMP_RNDN); - if (sign < 0) - mpfr_neg (v, v, GMP_RNDN); - - /* tentatively compute z as u+v; here we need z to be distinct - from a, b, c, d to not lose the latter */ - inex = mpfr_add (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_inf_p (u) || mpfr_inf_p (v) - || mpfr_zero_p (u) || mpfr_zero_p (v)) { - /* There is at least one over- or underflow. - Redo the computations with mpz_t exponents */ - mpfr_exp_t ea, eb, ec, ed; - 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); - eb = mpfr_get_exp (b); - ec = mpfr_get_exp (c); - ed = mpfr_get_exp (d); - - mpfr_set_exp ((mpfr_ptr) a, (mpfr_prec_t) 0); - mpfr_set_exp ((mpfr_ptr) b, (mpfr_prec_t) 0); - mpfr_set_exp ((mpfr_ptr) c, (mpfr_prec_t) 0); - mpfr_set_exp ((mpfr_ptr) d, (mpfr_prec_t) 0); - - mpz_init (eu); - mpz_init (ev); - mpz_set_si (eu, (long int) ea); - mpz_add_si (eu, eu, (long int) eb); - mpz_set_si (ev, (long int) ec); - mpz_add_si (ev, ev, (long int) ed); - - /* recompute u and v and move exponents to eu and ev */ - mpfr_mul (u, a, b, 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_mul (v, c, d, GMP_RNDN); - if (sign < 0) - mpfr_neg (v, v, 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, and - analogously for b. 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_add (z, u, v, rnd); - /* Result is finite since u and v have different signs. */ - overflow = mpfr_mul_2ui (z, z, mpz_get_ui (eu), rnd); - if (overflow) - inex = overflow; - } - else { - int underflow; - /* Addition of two zeroes with same sign. 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_add (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) b, eb); - mpfr_set_exp ((mpfr_ptr) c, ec); - mpfr_set_exp ((mpfr_ptr) d, ed); - /* works also when some of a, b, c, d are not all distinct */ - } - - mpfr_clear (u); - mpfr_clear (v); - - return inex; -} int mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) |