diff options
author | enge <enge@211d60ee-9f03-0410-a15a-8952a2c7a4e4> | 2011-02-22 15:49:20 +0000 |
---|---|---|
committer | enge <enge@211d60ee-9f03-0410-a15a-8952a2c7a4e4> | 2011-02-22 15:49:20 +0000 |
commit | a568e2d1885f6679342b6f60902a7d2c47b29d8a (patch) | |
tree | 61ff1129f690f33808c58bf203fc6649ab5162c1 /src/mul.c | |
parent | c82b65c57954bd20fff6512484e5df4003026dd8 (diff) | |
download | mpc-a568e2d1885f6679342b6f60902a7d2c47b29d8a.tar.gz |
mul.c: trivial changes in mpc_mul_naive
git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@943 211d60ee-9f03-0410-a15a-8952a2c7a4e4
Diffstat (limited to 'src/mul.c')
-rw-r--r-- | src/mul.c | 57 |
1 files changed, 25 insertions, 32 deletions
@@ -164,46 +164,39 @@ mul_imag (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) } +/* Computes z=x*y by the schoolbook method, where x and y are assumed + to be finite and different. */ int -mpc_mul_naive (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) +mpc_mul_naive (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) { - /* We assume that b and c are different, which is checked in mpc_mul. */ - int overlap, inex_re, inex_im; - mpfr_t v, t; - mpfr_prec_t prec; - - overlap = (a == b) || (a == c); - - prec = MPC_PREC_IM(b) + MPC_MAX_PREC(c); - - mpfr_init2 (v, prec); + int overlap, inex_re, inex_im; + mpfr_t v; + mpc_t rop; - /* Re(a) = Re(b)*Re(c) - Im(b)*Im(c) */ - /* FIXME: this code suffers undue overflows: v can overflow while the result - of the subtraction is representable */ - mpfr_mul (v, MPC_IM(b), MPC_IM(c), GMP_RNDN); /* exact */ + overlap = (z == x) || (z == y); + if (overlap) + mpc_init3 (rop, MPC_PREC_RE (z), MPC_PREC_IM (z)); + else + rop [0] = z [0]; - if (overlap) - { - mpfr_init2 (t, MPC_PREC_RE(a)); - inex_re = mpfr_fms (t, MPC_RE (b), MPC_RE (c), v, MPC_RND_RE (rnd)); - } - else - inex_re = mpfr_fms (MPC_RE (a), MPC_RE (b), MPC_RE (c), v, MPC_RND_RE (rnd)); + mpfr_init2 (v, MPC_PREC_IM (x) + MPC_MAX_PREC (y)); - /* Im(a) = Re(b)*Im(c) + Im(b)*Re(c) */ - mpfr_mul (v, MPC_IM(b), MPC_RE(c), GMP_RNDN); /* exact */ - inex_im = mpfr_fma (MPC_IM(a), MPC_RE (b), MPC_IM (c), v, MPC_RND_IM(rnd)); + /* Re(z) = Re(x)*Re(y) - Im(x)*Im(y) */ + /* FIXME: this code suffers undue overflows: v can overflow while the result + of the subtraction is representable */ + mpfr_mul (v, MPC_IM (x), MPC_IM (y), GMP_RNDN); /* exact */ + inex_re = mpfr_fms (MPC_RE (rop), MPC_RE (x), MPC_RE (y), v, MPC_RND_RE (rnd)); - mpfr_clear (v); + /* Im(z) = Re(x)*Im(y) + Im(x)*Re(y) */ + mpfr_mul (v, MPC_IM (x), MPC_RE (y), GMP_RNDN); /* exact */ + inex_im = mpfr_fma (MPC_IM (rop), MPC_RE (x), MPC_IM (y), v, MPC_RND_IM(rnd)); - if (overlap) - { - mpfr_set (MPC_RE(a), t, GMP_RNDN); /* exact */ - mpfr_clear (t); - } + mpfr_clear (v); + mpc_set (z, rop, MPC_RNDNN); /* exact */ + if (overlap) + mpc_clear (rop); - return MPC_INEX(inex_re, inex_im); + return MPC_INEX (inex_re, inex_im); } |