diff options
| author | enge <enge@211d60ee-9f03-0410-a15a-8952a2c7a4e4> | 2011-11-04 18:39:52 +0000 |
|---|---|---|
| committer | enge <enge@211d60ee-9f03-0410-a15a-8952a2c7a4e4> | 2011-11-04 18:39:52 +0000 |
| commit | 5ad0bccecc2c0aafdc23571a491fccc31171c22c (patch) | |
| tree | fbb979252a636e9fc3215e15e12be49ca6d35a95 /src | |
| parent | 18e9d42b59925987d8128096c4a7ea6d49721862 (diff) | |
| download | mpc-5ad0bccecc2c0aafdc23571a491fccc31171c22c.tar.gz | |
replaced MPC_RE by mpc_realref and MPC_IM by mpc_imagref everywhere
git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@1112 211d60ee-9f03-0410-a15a-8952a2c7a4e4
Diffstat (limited to 'src')
| -rw-r--r-- | src/abs.c | 4 | ||||
| -rw-r--r-- | src/acos.c | 100 | ||||
| -rw-r--r-- | src/acosh.c | 26 | ||||
| -rw-r--r-- | src/add.c | 6 | ||||
| -rw-r--r-- | src/add_fr.c | 6 | ||||
| -rw-r--r-- | src/add_si.c | 4 | ||||
| -rw-r--r-- | src/add_ui.c | 6 | ||||
| -rw-r--r-- | src/arg.c | 4 | ||||
| -rw-r--r-- | src/asin.c | 136 | ||||
| -rw-r--r-- | src/asinh.c | 16 | ||||
| -rw-r--r-- | src/atan.c | 96 | ||||
| -rw-r--r-- | src/atanh.c | 16 | ||||
| -rw-r--r-- | src/clear.c | 6 | ||||
| -rw-r--r-- | src/cmp.c | 6 | ||||
| -rw-r--r-- | src/cmp_si_si.c | 6 | ||||
| -rw-r--r-- | src/conj.c | 6 | ||||
| -rw-r--r-- | src/cosh.c | 8 | ||||
| -rw-r--r-- | src/div.c | 162 | ||||
| -rw-r--r-- | src/div_2exp.c | 6 | ||||
| -rw-r--r-- | src/div_fr.c | 10 | ||||
| -rw-r--r-- | src/div_ui.c | 6 | ||||
| -rw-r--r-- | src/exp.c | 68 | ||||
| -rw-r--r-- | src/fma.c | 44 | ||||
| -rw-r--r-- | src/fr_div.c | 10 | ||||
| -rw-r--r-- | src/fr_sub.c | 8 | ||||
| -rw-r--r-- | src/get_x.c | 4 | ||||
| -rw-r--r-- | src/imag.c | 6 | ||||
| -rw-r--r-- | src/init2.c | 6 | ||||
| -rw-r--r-- | src/init3.c | 6 | ||||
| -rw-r--r-- | src/inp_str.c | 6 | ||||
| -rw-r--r-- | src/log.c | 66 | ||||
| -rw-r--r-- | src/mpc-impl.h | 19 | ||||
| -rw-r--r-- | src/mul.c | 160 | ||||
| -rw-r--r-- | src/mul_2exp.c | 6 | ||||
| -rw-r--r-- | src/mul_fr.c | 14 | ||||
| -rw-r--r-- | src/mul_i.c | 32 | ||||
| -rw-r--r-- | src/mul_si.c | 6 | ||||
| -rw-r--r-- | src/mul_ui.c | 6 | ||||
| -rw-r--r-- | src/neg.c | 6 | ||||
| -rw-r--r-- | src/norm.c | 22 | ||||
| -rw-r--r-- | src/out_str.c | 6 | ||||
| -rw-r--r-- | src/pow.c | 144 | ||||
| -rw-r--r-- | src/pow_fr.c | 10 | ||||
| -rw-r--r-- | src/pow_ui.c | 16 | ||||
| -rw-r--r-- | src/proj.c | 4 | ||||
| -rw-r--r-- | src/real.c | 4 | ||||
| -rw-r--r-- | src/set.c | 6 | ||||
| -rw-r--r-- | src/set_prec.c | 6 | ||||
| -rw-r--r-- | src/set_str.c | 6 | ||||
| -rw-r--r-- | src/sin_cos.c | 128 | ||||
| -rw-r--r-- | src/sinh.c | 14 | ||||
| -rw-r--r-- | src/sqr.c | 88 | ||||
| -rw-r--r-- | src/sqrt.c | 102 | ||||
| -rw-r--r-- | src/strtoc.c | 12 | ||||
| -rw-r--r-- | src/sub.c | 6 | ||||
| -rw-r--r-- | src/sub_fr.c | 6 | ||||
| -rw-r--r-- | src/sub_ui.c | 6 | ||||
| -rw-r--r-- | src/swap.c | 6 | ||||
| -rw-r--r-- | src/tan.c | 110 | ||||
| -rw-r--r-- | src/tanh.c | 14 | ||||
| -rw-r--r-- | src/ui_ui_sub.c | 4 | ||||
| -rw-r--r-- | src/urandom.c | 6 |
62 files changed, 916 insertions, 919 deletions
@@ -1,6 +1,6 @@ /* mpc_abs -- Absolute value of a complex number. -Copyright (C) 2008, 2009 INRIA +Copyright (C) 2008, 2009, 2011 INRIA This file is part of GNU MPC. @@ -24,5 +24,5 @@ along with this program. If not, see http://www.gnu.org/licenses/ . int mpc_abs (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd) { - return mpfr_hypot (a, MPC_RE(b), MPC_IM(b), rnd); + return mpfr_hypot (a, mpc_realref(b), mpc_imagref(b), rnd); } @@ -1,6 +1,6 @@ /* mpc_acos -- arccosine of a complex number. -Copyright (C) 2009, 2010 INRIA +Copyright (C) 2009, 2010, 2011 INRIA This file is part of GNU MPC. @@ -36,38 +36,38 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) inex_im = 0; /* special values */ - if (mpfr_nan_p (MPC_RE (op)) || mpfr_nan_p (MPC_IM (op))) + if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op))) { - if (mpfr_inf_p (MPC_RE (op)) || mpfr_inf_p (MPC_IM (op))) + if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op))) { - mpfr_set_inf (MPC_IM (rop), mpfr_signbit (MPC_IM (op)) ? +1 : -1); - mpfr_set_nan (MPC_RE (rop)); + mpfr_set_inf (mpc_imagref (rop), mpfr_signbit (mpc_imagref (op)) ? +1 : -1); + mpfr_set_nan (mpc_realref (rop)); } - else if (mpfr_zero_p (MPC_RE (op))) + else if (mpfr_zero_p (mpc_realref (op))) { - inex_re = set_pi_over_2 (MPC_RE (rop), +1, MPC_RND_RE (rnd)); - mpfr_set_nan (MPC_IM (rop)); + inex_re = set_pi_over_2 (mpc_realref (rop), +1, MPC_RND_RE (rnd)); + mpfr_set_nan (mpc_imagref (rop)); } else { - mpfr_set_nan (MPC_RE (rop)); - mpfr_set_nan (MPC_IM (rop)); + mpfr_set_nan (mpc_realref (rop)); + mpfr_set_nan (mpc_imagref (rop)); } return MPC_INEX (inex_re, 0); } - if (mpfr_inf_p (MPC_RE (op)) || mpfr_inf_p (MPC_IM (op))) + if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op))) { - if (mpfr_inf_p (MPC_RE (op))) + if (mpfr_inf_p (mpc_realref (op))) { - if (mpfr_inf_p (MPC_IM (op))) + if (mpfr_inf_p (mpc_imagref (op))) { - if (mpfr_sgn (MPC_RE (op)) > 0) + if (mpfr_sgn (mpc_realref (op)) > 0) { inex_re = - set_pi_over_2 (MPC_RE (rop), +1, MPC_RND_RE (rnd)); - mpfr_div_2ui (MPC_RE (rop), MPC_RE (rop), 1, GMP_RNDN); + set_pi_over_2 (mpc_realref (rop), +1, MPC_RND_RE (rnd)); + mpfr_div_2ui (mpc_realref (rop), mpc_realref (rop), 1, GMP_RNDN); } else { @@ -81,7 +81,7 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_prec_t prec; int ok; mpfr_init (x); - prec = mpfr_get_prec (MPC_RE (rop)); + prec = mpfr_get_prec (mpc_realref (rop)); p = prec; do @@ -96,61 +96,61 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) } while (ok == 0); inex_re = - mpfr_div_2ui (MPC_RE (rop), x, 2, MPC_RND_RE (rnd)); + mpfr_div_2ui (mpc_realref (rop), x, 2, MPC_RND_RE (rnd)); mpfr_clear (x); } } else { - if (mpfr_sgn (MPC_RE (op)) > 0) - mpfr_set_ui (MPC_RE (rop), 0, GMP_RNDN); + if (mpfr_sgn (mpc_realref (op)) > 0) + mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); else - inex_re = mpfr_const_pi (MPC_RE (rop), MPC_RND_RE (rnd)); + inex_re = mpfr_const_pi (mpc_realref (rop), MPC_RND_RE (rnd)); } } else - inex_re = set_pi_over_2 (MPC_RE (rop), +1, MPC_RND_RE (rnd)); + inex_re = set_pi_over_2 (mpc_realref (rop), +1, MPC_RND_RE (rnd)); - mpfr_set_inf (MPC_IM (rop), mpfr_signbit (MPC_IM (op)) ? +1 : -1); + mpfr_set_inf (mpc_imagref (rop), mpfr_signbit (mpc_imagref (op)) ? +1 : -1); return MPC_INEX (inex_re, 0); } /* pure real argument */ - if (mpfr_zero_p (MPC_IM (op))) + if (mpfr_zero_p (mpc_imagref (op))) { int s_im; - s_im = mpfr_signbit (MPC_IM (op)); + s_im = mpfr_signbit (mpc_imagref (op)); - if (mpfr_cmp_ui (MPC_RE (op), 1) > 0) + if (mpfr_cmp_ui (mpc_realref (op), 1) > 0) { if (s_im) - inex_im = mpfr_acosh (MPC_IM (rop), MPC_RE (op), + inex_im = mpfr_acosh (mpc_imagref (rop), mpc_realref (op), MPC_RND_IM (rnd)); else - inex_im = -mpfr_acosh (MPC_IM (rop), MPC_RE (op), + inex_im = -mpfr_acosh (mpc_imagref (rop), mpc_realref (op), INV_RND (MPC_RND_IM (rnd))); - mpfr_set_ui (MPC_RE (rop), 0, GMP_RNDN); + mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); } - else if (mpfr_cmp_si (MPC_RE (op), -1) < 0) + else if (mpfr_cmp_si (mpc_realref (op), -1) < 0) { mpfr_t minus_op_re; - minus_op_re[0] = MPC_RE (op)[0]; + minus_op_re[0] = mpc_realref (op)[0]; MPFR_CHANGE_SIGN (minus_op_re); if (s_im) - inex_im = mpfr_acosh (MPC_IM (rop), minus_op_re, + inex_im = mpfr_acosh (mpc_imagref (rop), minus_op_re, MPC_RND_IM (rnd)); else - inex_im = -mpfr_acosh (MPC_IM (rop), minus_op_re, + inex_im = -mpfr_acosh (mpc_imagref (rop), minus_op_re, INV_RND (MPC_RND_IM (rnd))); - inex_re = mpfr_const_pi (MPC_RE (rop), MPC_RND_RE (rnd)); + inex_re = mpfr_const_pi (mpc_realref (rop), MPC_RND_RE (rnd)); } else { - inex_re = mpfr_acos (MPC_RE (rop), MPC_RE (op), MPC_RND_RE (rnd)); - mpfr_set_ui (MPC_IM (rop), 0, MPC_RND_IM (rnd)); + inex_re = mpfr_acos (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd)); + mpfr_set_ui (mpc_imagref (rop), 0, MPC_RND_IM (rnd)); } if (!s_im) @@ -160,10 +160,10 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) } /* pure imaginary argument */ - if (mpfr_zero_p (MPC_RE (op))) + if (mpfr_zero_p (mpc_realref (op))) { - inex_re = set_pi_over_2 (MPC_RE (rop), +1, MPC_RND_RE (rnd)); - inex_im = -mpfr_asinh (MPC_IM (rop), MPC_IM (op), + inex_re = set_pi_over_2 (mpc_realref (rop), +1, MPC_RND_RE (rnd)); + inex_im = -mpfr_asinh (mpc_imagref (rop), mpc_imagref (op), INV_RND (MPC_RND_IM (rnd))); mpc_conj (rop,rop, MPC_RNDNN); @@ -171,8 +171,8 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) } /* regular complex argument: acos(z) = Pi/2 - asin(z) */ - p_re = mpfr_get_prec (MPC_RE(rop)); - p_im = mpfr_get_prec (MPC_IM(rop)); + p_re = mpfr_get_prec (mpc_realref(rop)); + p_im = mpfr_get_prec (mpc_imagref(rop)); p = p_re; mpc_init3 (z1, p, p_im); /* we round directly the imaginary part to p_im, with rounding mode opposite to rnd_im */ @@ -181,7 +181,7 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) Im(z) > 0 and rnd_im = RNDZ, we want to round the Im(asin(z)) to -Inf so that -Im(asin(z)) is rounded to zero */ if (rnd_im == GMP_RNDZ) - rnd_im = mpfr_sgn (MPC_IM(op)) > 0 ? GMP_RNDD : GMP_RNDU; + rnd_im = mpfr_sgn (mpc_imagref(op)) > 0 ? GMP_RNDD : GMP_RNDU; else rnd_im = rnd_im == GMP_RNDU ? GMP_RNDD : rnd_im == GMP_RNDD ? GMP_RNDU @@ -192,30 +192,30 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) { p += mpc_ceil_log2 (p) + 3; - mpfr_set_prec (MPC_RE(z1), p); + mpfr_set_prec (mpc_realref(z1), p); mpfr_set_prec (pi_over_2, p); mpfr_const_pi (pi_over_2, GMP_RNDN); mpfr_div_2exp (pi_over_2, pi_over_2, 1, GMP_RNDN); /* Pi/2 */ e1 = 1; /* Exp(pi_over_2) */ inex = mpc_asin (z1, op, rnd1); /* asin(z) */ - MPC_ASSERT (mpfr_sgn (MPC_IM(z1)) * mpfr_sgn (MPC_IM(op)) > 0); + MPC_ASSERT (mpfr_sgn (mpc_imagref(z1)) * mpfr_sgn (mpc_imagref(op)) > 0); inex_im = MPC_INEX_IM(inex); /* inex_im is in {-1, 0, 1} */ - e2 = mpfr_get_exp (MPC_RE(z1)); - mpfr_sub (MPC_RE(z1), pi_over_2, MPC_RE(z1), GMP_RNDN); - if (!mpfr_zero_p (MPC_RE(z1))) + e2 = mpfr_get_exp (mpc_realref(z1)); + mpfr_sub (mpc_realref(z1), pi_over_2, mpc_realref(z1), GMP_RNDN); + if (!mpfr_zero_p (mpc_realref(z1))) { /* the error on x=Re(z1) is bounded by 1/2 ulp(x) + 2^(e1-p-1) + 2^(e2-p-1) */ e1 = e1 >= e2 ? e1 + 1 : e2 + 1; /* the error on x is bounded by 1/2 ulp(x) + 2^(e1-p-1) */ - e1 -= mpfr_get_exp (MPC_RE(z1)); + e1 -= mpfr_get_exp (mpc_realref(z1)); /* the error on x is bounded by 1/2 ulp(x) [1 + 2^e1] */ e1 = e1 <= 0 ? 0 : e1; /* the error on x is bounded by 2^e1 * ulp(x) */ - mpfr_neg (MPC_IM(z1), MPC_IM(z1), GMP_RNDN); /* exact */ + mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), GMP_RNDN); /* exact */ inex_im = -inex_im; - if (mpfr_can_round (MPC_RE(z1), p - e1, GMP_RNDN, GMP_RNDZ, + if (mpfr_can_round (mpc_realref(z1), p - e1, GMP_RNDN, GMP_RNDZ, p_re + (MPC_RND_RE(rnd) == GMP_RNDN))) break; } diff --git a/src/acosh.c b/src/acosh.c index fde4f5e..a4eef4c 100644 --- a/src/acosh.c +++ b/src/acosh.c @@ -1,6 +1,6 @@ /* mpc_acosh -- inverse hyperbolic cosine of a complex number. -Copyright (C) 2009 INRIA +Copyright (C) 2009, 2011 INRIA This file is part of GNU MPC. @@ -33,26 +33,26 @@ mpc_acosh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_t tmp; int inex; - if (mpfr_zero_p (MPC_RE (op)) && mpfr_nan_p (MPC_IM (op))) + if (mpfr_zero_p (mpc_realref (op)) && mpfr_nan_p (mpc_imagref (op))) { - mpfr_set_nan (MPC_RE (rop)); - mpfr_set_nan (MPC_IM (rop)); + mpfr_set_nan (mpc_realref (rop)); + mpfr_set_nan (mpc_imagref (rop)); return 0; } /* Note reversal of precisions due to later multiplication by i or -i */ mpc_init3 (a, MPC_PREC_IM(rop), MPC_PREC_RE(rop)); - if (mpfr_signbit (MPC_IM (op))) + if (mpfr_signbit (mpc_imagref (op))) { inex = mpc_acos (a, op, RNDC (INV_RND (MPC_RND_IM (rnd)), MPC_RND_RE (rnd))); /* change a to -i*a, i.e., -y+i*x to x+i*y */ - tmp[0] = MPC_RE (a)[0]; - MPC_RE (a)[0] = MPC_IM (a)[0]; - MPC_IM (a)[0] = tmp[0]; - MPFR_CHANGE_SIGN (MPC_IM (a)); + tmp[0] = mpc_realref (a)[0]; + mpc_realref (a)[0] = mpc_imagref (a)[0]; + mpc_imagref (a)[0] = tmp[0]; + MPFR_CHANGE_SIGN (mpc_imagref (a)); inex = MPC_INEX (MPC_INEX_IM (inex), -MPC_INEX_RE (inex)); } else @@ -61,10 +61,10 @@ mpc_acosh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) RNDC (MPC_RND_IM (rnd), INV_RND(MPC_RND_RE (rnd)))); /* change a to i*a, i.e., y-i*x to x+i*y */ - tmp[0] = MPC_RE (a)[0]; - MPC_RE (a)[0] = MPC_IM (a)[0]; - MPC_IM (a)[0] = tmp[0]; - MPFR_CHANGE_SIGN (MPC_RE (a)); + tmp[0] = mpc_realref (a)[0]; + mpc_realref (a)[0] = mpc_imagref (a)[0]; + mpc_imagref (a)[0] = tmp[0]; + MPFR_CHANGE_SIGN (mpc_realref (a)); inex = MPC_INEX (-MPC_INEX_IM (inex), MPC_INEX_RE (inex)); } @@ -1,6 +1,6 @@ /* mpc_add -- Add two complex numbers. -Copyright (C) 2002, 2009 INRIA +Copyright (C) 2002, 2009, 2011 INRIA This file is part of GNU MPC. @@ -26,8 +26,8 @@ mpc_add (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_add (MPC_RE(a), MPC_RE(b), MPC_RE(c), MPC_RND_RE(rnd)); - inex_im = mpfr_add (MPC_IM(a), MPC_IM(b), MPC_IM(c), MPC_RND_IM(rnd)); + inex_re = mpfr_add (mpc_realref(a), mpc_realref(b), mpc_realref(c), MPC_RND_RE(rnd)); + inex_im = mpfr_add (mpc_imagref(a), mpc_imagref(b), mpc_imagref(c), MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } diff --git a/src/add_fr.c b/src/add_fr.c index b9f2e9f..ea7b595 100644 --- a/src/add_fr.c +++ b/src/add_fr.c @@ -1,6 +1,6 @@ /* mpc_add_fr -- Add a complex number and a floating-point number. -Copyright (C) 2002, 2009 INRIA +Copyright (C) 2002, 2009, 2011 INRIA This file is part of GNU MPC. @@ -26,8 +26,8 @@ mpc_add_fr (mpc_ptr a, mpc_srcptr b, mpfr_srcptr c, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_add (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); - inex_im = mpfr_set (MPC_IM(a), MPC_IM(b), MPC_RND_IM(rnd)); + inex_re = mpfr_add (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_set (mpc_imagref(a), mpc_imagref(b), MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } diff --git a/src/add_si.c b/src/add_si.c index d48d252..ba14803 100644 --- a/src/add_si.c +++ b/src/add_si.c @@ -25,8 +25,8 @@ mpc_add_si (mpc_ptr rop, mpc_srcptr op1, long int op2, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_add_si (MPC_RE (rop), MPC_RE (op1), op2, MPC_RND_RE (rnd)); - inex_im = mpfr_set (MPC_IM (rop), MPC_IM (op1), MPC_RND_IM (rnd)); + inex_re = mpfr_add_si (mpc_realref (rop), mpc_realref (op1), op2, MPC_RND_RE (rnd)); + inex_im = mpfr_set (mpc_imagref (rop), mpc_imagref (op1), MPC_RND_IM (rnd)); return MPC_INEX (inex_re, inex_im); } diff --git a/src/add_ui.c b/src/add_ui.c index 0ea4aaa..85f4d13 100644 --- a/src/add_ui.c +++ b/src/add_ui.c @@ -1,6 +1,6 @@ /* mpc_add_ui -- Add a complex number and an unsigned long int. -Copyright (C) 2002, 2009 INRIA +Copyright (C) 2002, 2009, 2011 INRIA This file is part of GNU MPC. @@ -26,8 +26,8 @@ mpc_add_ui (mpc_ptr a, mpc_srcptr b, unsigned long int c, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_add_ui (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); - inex_im = mpfr_set (MPC_IM(a), MPC_IM(b), MPC_RND_IM(rnd)); + inex_re = mpfr_add_ui (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_set (mpc_imagref(a), mpc_imagref(b), MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } @@ -1,6 +1,6 @@ /* mpc_arg -- Get the argument of a complex number. -Copyright (C) 2008, 2009 INRIA +Copyright (C) 2008, 2009, 2011 INRIA This file is part of GNU MPC. @@ -23,5 +23,5 @@ along with this program. If not, see http://www.gnu.org/licenses/ . int mpc_arg (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd) { - return mpfr_atan2 (a, MPC_IM (b), MPC_RE (b), rnd); + return mpfr_atan2 (a, mpc_imagref (b), mpc_realref (b), rnd); } @@ -29,117 +29,117 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) int inex; /* special values */ - if (mpfr_nan_p (MPC_RE (op)) || mpfr_nan_p (MPC_IM (op))) + if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op))) { - if (mpfr_inf_p (MPC_RE (op)) || mpfr_inf_p (MPC_IM (op))) + if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op))) { - mpfr_set_nan (MPC_RE (rop)); - mpfr_set_inf (MPC_IM (rop), mpfr_signbit (MPC_IM (op)) ? -1 : +1); + mpfr_set_nan (mpc_realref (rop)); + mpfr_set_inf (mpc_imagref (rop), mpfr_signbit (mpc_imagref (op)) ? -1 : +1); } - else if (mpfr_zero_p (MPC_RE (op))) + else if (mpfr_zero_p (mpc_realref (op))) { - mpfr_set (MPC_RE (rop), MPC_RE (op), GMP_RNDN); - mpfr_set_nan (MPC_IM (rop)); + mpfr_set (mpc_realref (rop), mpc_realref (op), GMP_RNDN); + mpfr_set_nan (mpc_imagref (rop)); } else { - mpfr_set_nan (MPC_RE (rop)); - mpfr_set_nan (MPC_IM (rop)); + mpfr_set_nan (mpc_realref (rop)); + mpfr_set_nan (mpc_imagref (rop)); } return 0; } - if (mpfr_inf_p (MPC_RE (op)) || mpfr_inf_p (MPC_IM (op))) + if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op))) { int inex_re; - if (mpfr_inf_p (MPC_RE (op))) + if (mpfr_inf_p (mpc_realref (op))) { - int inf_im = mpfr_inf_p (MPC_IM (op)); + int inf_im = mpfr_inf_p (mpc_imagref (op)); - inex_re = set_pi_over_2 (MPC_RE (rop), - (mpfr_signbit (MPC_RE (op)) ? -1 : 1), MPC_RND_RE (rnd)); - mpfr_set_inf (MPC_IM (rop), (mpfr_signbit (MPC_IM (op)) ? -1 : 1)); + inex_re = set_pi_over_2 (mpc_realref (rop), + (mpfr_signbit (mpc_realref (op)) ? -1 : 1), MPC_RND_RE (rnd)); + mpfr_set_inf (mpc_imagref (rop), (mpfr_signbit (mpc_imagref (op)) ? -1 : 1)); if (inf_im) - mpfr_div_2ui (MPC_RE (rop), MPC_RE (rop), 1, GMP_RNDN); + mpfr_div_2ui (mpc_realref (rop), mpc_realref (rop), 1, GMP_RNDN); } else { - mpfr_set_zero (MPC_RE (rop), (mpfr_signbit (MPC_RE (op)) ? -1 : 1)); + mpfr_set_zero (mpc_realref (rop), (mpfr_signbit (mpc_realref (op)) ? -1 : 1)); inex_re = 0; - mpfr_set_inf (MPC_IM (rop), (mpfr_signbit (MPC_IM (op)) ? -1 : 1)); + mpfr_set_inf (mpc_imagref (rop), (mpfr_signbit (mpc_imagref (op)) ? -1 : 1)); } return MPC_INEX (inex_re, 0); } /* pure real argument */ - if (mpfr_zero_p (MPC_IM (op))) + if (mpfr_zero_p (mpc_imagref (op))) { int inex_re; int inex_im; int s_im; - s_im = mpfr_signbit (MPC_IM (op)); + s_im = mpfr_signbit (mpc_imagref (op)); - if (mpfr_cmp_ui (MPC_RE (op), 1) > 0) + if (mpfr_cmp_ui (mpc_realref (op), 1) > 0) { if (s_im) - inex_im = -mpfr_acosh (MPC_IM (rop), MPC_RE (op), + inex_im = -mpfr_acosh (mpc_imagref (rop), mpc_realref (op), INV_RND (MPC_RND_IM (rnd))); else - inex_im = mpfr_acosh (MPC_IM (rop), MPC_RE (op), + inex_im = mpfr_acosh (mpc_imagref (rop), mpc_realref (op), MPC_RND_IM (rnd)); - inex_re = set_pi_over_2 (MPC_RE (rop), - (mpfr_signbit (MPC_RE (op)) ? -1 : 1), MPC_RND_RE (rnd)); + inex_re = set_pi_over_2 (mpc_realref (rop), + (mpfr_signbit (mpc_realref (op)) ? -1 : 1), MPC_RND_RE (rnd)); if (s_im) mpc_conj (rop, rop, MPC_RNDNN); } - else if (mpfr_cmp_si (MPC_RE (op), -1) < 0) + else if (mpfr_cmp_si (mpc_realref (op), -1) < 0) { mpfr_t minus_op_re; - minus_op_re[0] = MPC_RE (op)[0]; + minus_op_re[0] = mpc_realref (op)[0]; MPFR_CHANGE_SIGN (minus_op_re); if (s_im) - inex_im = -mpfr_acosh (MPC_IM (rop), minus_op_re, + inex_im = -mpfr_acosh (mpc_imagref (rop), minus_op_re, INV_RND (MPC_RND_IM (rnd))); else - inex_im = mpfr_acosh (MPC_IM (rop), minus_op_re, + inex_im = mpfr_acosh (mpc_imagref (rop), minus_op_re, MPC_RND_IM (rnd)); - inex_re = set_pi_over_2 (MPC_RE (rop), - (mpfr_signbit (MPC_RE (op)) ? -1 : 1), MPC_RND_RE (rnd)); + inex_re = set_pi_over_2 (mpc_realref (rop), + (mpfr_signbit (mpc_realref (op)) ? -1 : 1), MPC_RND_RE (rnd)); if (s_im) mpc_conj (rop, rop, MPC_RNDNN); } else { - inex_im = mpfr_set_ui (MPC_IM (rop), 0, MPC_RND_IM (rnd)); + inex_im = mpfr_set_ui (mpc_imagref (rop), 0, MPC_RND_IM (rnd)); if (s_im) - mpfr_neg (MPC_IM (rop), MPC_IM (rop), GMP_RNDN); - inex_re = mpfr_asin (MPC_RE (rop), MPC_RE (op), MPC_RND_RE (rnd)); + mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), GMP_RNDN); + inex_re = mpfr_asin (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd)); } return MPC_INEX (inex_re, inex_im); } /* pure imaginary argument */ - if (mpfr_zero_p (MPC_RE (op))) + if (mpfr_zero_p (mpc_realref (op))) { int inex_im; int s; - s = mpfr_signbit (MPC_RE (op)); - mpfr_set_ui (MPC_RE (rop), 0, GMP_RNDN); + s = mpfr_signbit (mpc_realref (op)); + mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); if (s) - mpfr_neg (MPC_RE (rop), MPC_RE (rop), GMP_RNDN); - inex_im = mpfr_asinh (MPC_IM (rop), MPC_IM (op), MPC_RND_IM (rnd)); + mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN); + inex_im = mpfr_asinh (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd)); return MPC_INEX (0, inex_im); } /* regular complex: asin(z) = -i*log(i*z+sqrt(1-z^2)) */ - p_re = mpfr_get_prec (MPC_RE(rop)); - p_im = mpfr_get_prec (MPC_IM(rop)); + p_re = mpfr_get_prec (mpc_realref(rop)); + p_im = mpfr_get_prec (mpc_imagref(rop)); rnd_re = MPC_RND_RE(rnd); rnd_im = MPC_RND_IM(rnd); p = p_re >= p_im ? p_re : p_im; @@ -150,71 +150,71 @@ mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) p += mpc_ceil_log2 (p) + 3 + incr_p; /* incr_p is zero initially */ incr_p = p / 2; - mpfr_set_prec (MPC_RE(z1), p); - mpfr_set_prec (MPC_IM(z1), p); + mpfr_set_prec (mpc_realref(z1), p); + mpfr_set_prec (mpc_imagref(z1), p); /* z1 <- z^2 */ mpc_sqr (z1, op, MPC_RNDNN); /* err(x) <= 1/2 ulp(x), err(y) <= 1/2 ulp(y) */ /* z1 <- 1-z1 */ - ex = mpfr_get_exp (MPC_RE(z1)); - mpfr_ui_sub (MPC_RE(z1), 1, MPC_RE(z1), GMP_RNDN); - mpfr_neg (MPC_IM(z1), MPC_IM(z1), GMP_RNDN); - ex = ex - mpfr_get_exp (MPC_RE(z1)); + ex = mpfr_get_exp (mpc_realref(z1)); + mpfr_ui_sub (mpc_realref(z1), 1, mpc_realref(z1), GMP_RNDN); + mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), GMP_RNDN); + ex = ex - mpfr_get_exp (mpc_realref(z1)); ex = (ex <= 0) ? 0 : ex; /* err(x) <= 2^ex * ulp(x) */ - ex = ex + mpfr_get_exp (MPC_RE(z1)) - p; + ex = ex + mpfr_get_exp (mpc_realref(z1)) - p; /* err(x) <= 2^ex */ - ey = mpfr_get_exp (MPC_IM(z1)) - p - 1; + ey = mpfr_get_exp (mpc_imagref(z1)) - p - 1; /* err(y) <= 2^ey */ ex = (ex >= ey) ? ex : ey; /* err(x), err(y) <= 2^ex, i.e., the norm of the error is bounded by |h|<=2^(ex+1/2) */ /* z1 <- sqrt(z1): if z1 = z + h, then sqrt(z1) = sqrt(z) + h/2/sqrt(t) */ - ey = mpfr_get_exp (MPC_RE(z1)) >= mpfr_get_exp (MPC_IM(z1)) - ? mpfr_get_exp (MPC_RE(z1)) : mpfr_get_exp (MPC_IM(z1)); + ey = mpfr_get_exp (mpc_realref(z1)) >= mpfr_get_exp (mpc_imagref(z1)) + ? mpfr_get_exp (mpc_realref(z1)) : mpfr_get_exp (mpc_imagref(z1)); /* we have |z1| >= 2^(ey-1) thus 1/|z1| <= 2^(1-ey) */ mpc_sqrt (z1, z1, MPC_RNDNN); ex = (2 * ex + 1) - 2 - (ey - 1); /* |h^2/4/|t| <= 2^ex */ ex = (ex + 1) / 2; /* ceil(ex/2) */ /* express ex in terms of ulp(z1) */ - ey = mpfr_get_exp (MPC_RE(z1)) <= mpfr_get_exp (MPC_IM(z1)) - ? mpfr_get_exp (MPC_RE(z1)) : mpfr_get_exp (MPC_IM(z1)); + ey = mpfr_get_exp (mpc_realref(z1)) <= mpfr_get_exp (mpc_imagref(z1)) + ? mpfr_get_exp (mpc_realref(z1)) : mpfr_get_exp (mpc_imagref(z1)); ex = ex - ey + p; /* take into account the rounding error in the mpc_sqrt call */ err = (ex <= 0) ? 1 : ex + 1; /* err(x) <= 2^err * ulp(x), err(y) <= 2^err * ulp(y) */ /* z1 <- i*z + z1 */ - ex = mpfr_get_exp (MPC_RE(z1)); - ey = mpfr_get_exp (MPC_IM(z1)); - mpfr_sub (MPC_RE(z1), MPC_RE(z1), MPC_IM(op), GMP_RNDN); - mpfr_add (MPC_IM(z1), MPC_IM(z1), MPC_RE(op), GMP_RNDN); - if (mpfr_cmp_ui (MPC_RE(z1), 0) == 0 || mpfr_cmp_ui (MPC_IM(z1), 0) == 0) + ex = mpfr_get_exp (mpc_realref(z1)); + ey = mpfr_get_exp (mpc_imagref(z1)); + mpfr_sub (mpc_realref(z1), mpc_realref(z1), mpc_imagref(op), GMP_RNDN); + mpfr_add (mpc_imagref(z1), mpc_imagref(z1), mpc_realref(op), GMP_RNDN); + if (mpfr_cmp_ui (mpc_realref(z1), 0) == 0 || mpfr_cmp_ui (mpc_imagref(z1), 0) == 0) continue; - ex -= mpfr_get_exp (MPC_RE(z1)); /* cancellation in x */ - ey -= mpfr_get_exp (MPC_IM(z1)); /* cancellation in y */ + ex -= mpfr_get_exp (mpc_realref(z1)); /* cancellation in x */ + ey -= mpfr_get_exp (mpc_imagref(z1)); /* cancellation in y */ ex = (ex >= ey) ? ex : ey; /* maximum cancellation */ err += ex; err = (err <= 0) ? 1 : err + 1; /* rounding error in sub/add */ /* z1 <- log(z1): if z1 = z + h, then log(z1) = log(z) + h/t with |t| >= min(|z1|,|z|) */ - ex = mpfr_get_exp (MPC_RE(z1)); - ey = mpfr_get_exp (MPC_IM(z1)); + ex = mpfr_get_exp (mpc_realref(z1)); + ey = mpfr_get_exp (mpc_imagref(z1)); ex = (ex >= ey) ? ex : ey; err += ex - p; /* revert to absolute error <= 2^err */ mpc_log (z1, z1, GMP_RNDN); err -= ex - 1; /* 1/|t| <= 1/|z| <= 2^(1-ex) */ /* express err in terms of ulp(z1) */ - ey = mpfr_get_exp (MPC_RE(z1)) <= mpfr_get_exp (MPC_IM(z1)) - ? mpfr_get_exp (MPC_RE(z1)) : mpfr_get_exp (MPC_IM(z1)); + ey = mpfr_get_exp (mpc_realref(z1)) <= mpfr_get_exp (mpc_imagref(z1)) + ? mpfr_get_exp (mpc_realref(z1)) : mpfr_get_exp (mpc_imagref(z1)); err = err - ey + p; /* take into account the rounding error in the mpc_log call */ err = (err <= 0) ? 1 : err + 1; /* z1 <- -i*z1 */ - mpfr_swap (MPC_RE(z1), MPC_IM(z1)); - mpfr_neg (MPC_IM(z1), MPC_IM(z1), GMP_RNDN); - if (mpfr_can_round (MPC_RE(z1), p - err, GMP_RNDN, GMP_RNDZ, + mpfr_swap (mpc_realref(z1), mpc_imagref(z1)); + mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), GMP_RNDN); + if (mpfr_can_round (mpc_realref(z1), p - err, GMP_RNDN, GMP_RNDZ, p_re + (rnd_re == GMP_RNDN)) && - mpfr_can_round (MPC_IM(z1), p - err, GMP_RNDN, GMP_RNDZ, + mpfr_can_round (mpc_imagref(z1), p - err, GMP_RNDN, GMP_RNDZ, p_im + (rnd_im == GMP_RNDN))) break; } diff --git a/src/asinh.c b/src/asinh.c index 9a2cbcb..89477a6 100644 --- a/src/asinh.c +++ b/src/asinh.c @@ -1,6 +1,6 @@ /* mpc_asinh -- inverse hyperbolic sine of a complex number. -Copyright (C) 2009 INRIA +Copyright (C) 2009, 2011 INRIA This file is part of GNU MPC. @@ -29,9 +29,9 @@ mpc_asinh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_t tmp; /* z = i*op */ - MPC_RE (z)[0] = MPC_IM (op)[0]; - MPC_IM (z)[0] = MPC_RE (op)[0]; - MPFR_CHANGE_SIGN (MPC_RE (z)); + mpc_realref (z)[0] = mpc_imagref (op)[0]; + mpc_imagref (z)[0] = mpc_realref (op)[0]; + MPFR_CHANGE_SIGN (mpc_realref (z)); /* Note reversal of precisions due to later multiplication by -i */ mpc_init3 (a, MPC_PREC_IM(rop), MPC_PREC_RE(rop)); @@ -42,10 +42,10 @@ mpc_asinh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* if a = asin(i*op) = x+i*y, and we want y-i*x */ /* change a to -i*a */ - tmp[0] = MPC_RE (a)[0]; - MPC_RE (a)[0] = MPC_IM (a)[0]; - MPC_IM (a)[0] = tmp[0]; - MPFR_CHANGE_SIGN (MPC_IM (a)); + tmp[0] = mpc_realref (a)[0]; + mpc_realref (a)[0] = mpc_imagref (a)[0]; + mpc_imagref (a)[0] = tmp[0]; + MPFR_CHANGE_SIGN (mpc_imagref (a)); mpc_set (rop, a, MPC_RNDNN); /* exact */ @@ -52,45 +52,45 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) inex_re = 0; inex_im = 0; - s_re = mpfr_signbit (MPC_RE (op)); - s_im = mpfr_signbit (MPC_IM (op)); + s_re = mpfr_signbit (mpc_realref (op)); + s_im = mpfr_signbit (mpc_imagref (op)); /* special values */ - if (mpfr_nan_p (MPC_RE (op)) || mpfr_nan_p (MPC_IM (op))) + if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op))) { - if (mpfr_nan_p (MPC_RE (op))) + if (mpfr_nan_p (mpc_realref (op))) { - mpfr_set_nan (MPC_RE (rop)); - if (mpfr_zero_p (MPC_IM (op)) || mpfr_inf_p (MPC_IM (op))) + mpfr_set_nan (mpc_realref (rop)); + if (mpfr_zero_p (mpc_imagref (op)) || mpfr_inf_p (mpc_imagref (op))) { - mpfr_set_ui (MPC_IM (rop), 0, GMP_RNDN); + mpfr_set_ui (mpc_imagref (rop), 0, GMP_RNDN); if (s_im) mpc_conj (rop, rop, MPC_RNDNN); } else - mpfr_set_nan (MPC_IM (rop)); + mpfr_set_nan (mpc_imagref (rop)); } else { - if (mpfr_inf_p (MPC_RE (op))) + if (mpfr_inf_p (mpc_realref (op))) { - inex_re = set_pi_over_2 (MPC_RE (rop), -s_re, MPC_RND_RE (rnd)); - mpfr_set_ui (MPC_IM (rop), 0, GMP_RNDN); + inex_re = set_pi_over_2 (mpc_realref (rop), -s_re, MPC_RND_RE (rnd)); + mpfr_set_ui (mpc_imagref (rop), 0, GMP_RNDN); } else { - mpfr_set_nan (MPC_RE (rop)); - mpfr_set_nan (MPC_IM (rop)); + mpfr_set_nan (mpc_realref (rop)); + mpfr_set_nan (mpc_imagref (rop)); } } return MPC_INEX (inex_re, 0); } - if (mpfr_inf_p (MPC_RE (op)) || mpfr_inf_p (MPC_IM (op))) + if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op))) { - inex_re = set_pi_over_2 (MPC_RE (rop), -s_re, MPC_RND_RE (rnd)); + inex_re = set_pi_over_2 (mpc_realref (rop), -s_re, MPC_RND_RE (rnd)); - mpfr_set_ui (MPC_IM (rop), 0, GMP_RNDN); + mpfr_set_ui (mpc_imagref (rop), 0, GMP_RNDN); if (s_im) mpc_conj (rop, rop, GMP_RNDN); @@ -98,11 +98,11 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) } /* pure real argument */ - if (mpfr_zero_p (MPC_IM (op))) + if (mpfr_zero_p (mpc_imagref (op))) { - inex_re = mpfr_atan (MPC_RE (rop), MPC_RE (op), MPC_RND_RE (rnd)); + inex_re = mpfr_atan (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd)); - mpfr_set_ui (MPC_IM (rop), 0, GMP_RNDN); + mpfr_set_ui (mpc_imagref (rop), 0, GMP_RNDN); if (s_im) mpc_conj (rop, rop, GMP_RNDN); @@ -110,32 +110,32 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) } /* pure imaginary argument */ - if (mpfr_zero_p (MPC_RE (op))) + if (mpfr_zero_p (mpc_realref (op))) { int cmp_1; if (s_im) - cmp_1 = -mpfr_cmp_si (MPC_IM (op), -1); + cmp_1 = -mpfr_cmp_si (mpc_imagref (op), -1); else - cmp_1 = mpfr_cmp_ui (MPC_IM (op), +1); + cmp_1 = mpfr_cmp_ui (mpc_imagref (op), +1); if (cmp_1 < 0) { /* atan(+0+iy) = +0 +i*atanh(y), if |y| < 1 atan(-0+iy) = -0 +i*atanh(y), if |y| < 1 */ - mpfr_set_ui (MPC_RE (rop), 0, GMP_RNDN); + mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); if (s_re) - mpfr_neg (MPC_RE (rop), MPC_RE (rop), GMP_RNDN); + mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN); - inex_im = mpfr_atanh (MPC_IM (rop), MPC_IM (op), MPC_RND_IM (rnd)); + inex_im = mpfr_atanh (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd)); } else if (cmp_1 == 0) { /* atan(+/-0+i) = NaN +i*inf atan(+/-0-i) = NaN -i*inf */ - mpfr_set_nan (MPC_RE (rop)); - mpfr_set_inf (MPC_IM (rop), s_im ? -1 : +1); + mpfr_set_nan (mpc_realref (rop)); + mpfr_set_inf (mpc_imagref (rop), s_im ? -1 : +1); } else { @@ -148,7 +148,7 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) rnd_im = MPC_RND_IM (rnd); mpfr_init (y); - p_im = mpfr_get_prec (MPC_IM (rop)); + p_im = mpfr_get_prec (mpc_imagref (rop)); p = p_im; /* a = o(1/y) with error(a) < 1 ulp(a) @@ -164,7 +164,7 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) p += mpc_ceil_log2 (p) + 2; mpfr_set_prec (y, p); rnd_away = s_im == 0 ? GMP_RNDU : GMP_RNDD; - inex_im = mpfr_ui_div (y, 1, MPC_IM (op), rnd_away); + inex_im = mpfr_ui_div (y, 1, mpc_imagref (op), rnd_away); /* FIXME: should we consider the case with unreasonably huge precision prec(y)>3*exp_min, where atanh(1/Im(op)) could be representable while 1/Im(op) underflows ? @@ -179,8 +179,8 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) p_im + (rnd_im == GMP_RNDN)); } while (ok == 0); - inex_re = set_pi_over_2 (MPC_RE (rop), -s_re, MPC_RND_RE (rnd)); - inex_im = mpfr_set (MPC_IM (rop), y, rnd_im); + inex_re = set_pi_over_2 (mpc_realref (rop), -s_re, MPC_RND_RE (rnd)); + inex_im = mpfr_set (mpc_imagref (rop), y, rnd_im); mpfr_clear (y); } return MPC_INEX (inex_re, inex_im); @@ -199,12 +199,12 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_inits2 (MPFR_PREC_MIN, a, b, x, y, (mpfr_ptr) 0); /* real part: Re(arctan(x+i*y)) = [arctan2(x,1-y) - arctan2(-x,1+y)]/2 */ - minus_op_re[0] = MPC_RE (op)[0]; + minus_op_re[0] = mpc_realref (op)[0]; MPFR_CHANGE_SIGN (minus_op_re); - op_re_exp = mpfr_get_exp (MPC_RE (op)); - op_im_exp = mpfr_get_exp (MPC_IM (op)); + op_re_exp = mpfr_get_exp (mpc_realref (op)); + op_im_exp = mpfr_get_exp (mpc_imagref (op)); - prec = mpfr_get_prec (MPC_RE (rop)); /* result precision */ + prec = mpfr_get_prec (mpc_realref (rop)); /* result precision */ /* a = o(1-y) error(a) < 1 ulp(a) b = o(atan2(x,a)) error(b) < [1+2^{3+Exp(x)-Exp(a)-Exp(b)}] ulp(b) @@ -225,8 +225,8 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* p: working precision */ p = (op_im_exp > 0 || prec > SAFE_ABS (mpfr_prec_t, op_im_exp)) ? prec : (prec - op_im_exp); - rnd1 = mpfr_sgn (MPC_RE (op)) > 0 ? GMP_RNDD : GMP_RNDU; - rnd2 = mpfr_sgn (MPC_RE (op)) < 0 ? GMP_RNDU : GMP_RNDD; + rnd1 = mpfr_sgn (mpc_realref (op)) > 0 ? GMP_RNDD : GMP_RNDU; + rnd2 = mpfr_sgn (mpc_realref (op)) < 0 ? GMP_RNDU : GMP_RNDD; do { @@ -238,27 +238,27 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* x = upper bound for atan (x/(1-y)). Since atan is increasing, we need an upper bound on x/(1-y), i.e., a lower bound on 1-y for x positive, and an upper bound on 1-y for x negative */ - mpfr_ui_sub (a, 1, MPC_IM (op), rnd1); + mpfr_ui_sub (a, 1, mpc_imagref (op), rnd1); if (mpfr_sgn (a) == 0) /* y is near 1, thus 1+y is near 2, and expo will be 1 or 2 below */ { - if (mpfr_cmp_ui (MPC_IM(op), 1) != 0) + if (mpfr_cmp_ui (mpc_imagref(op), 1) != 0) continue; err = 2; /* ensures err will be expo below */ } else err = mpfr_get_exp (a); /* err = Exp(a) with the notations above */ - mpfr_atan2 (x, MPC_RE (op), a, GMP_RNDU); + mpfr_atan2 (x, mpc_realref (op), a, GMP_RNDU); /* b = lower bound for atan (-x/(1+y)): for x negative, we need a lower bound on -x/(1+y), i.e., an upper bound on 1+y */ - mpfr_add_ui (a, MPC_IM(op), 1, rnd2); + mpfr_add_ui (a, mpc_imagref(op), 1, rnd2); /* if a is zero but inexact, try again with a larger precision if a is exactly zero, i.e., Im(op) = -1, then the error on a is 0, and we can simply ignore the terms involving Exp(a) in the error */ if (mpfr_sgn (a) == 0) { - if (mpfr_cmp_si (MPC_IM(op), -1) != 0) + if (mpfr_cmp_si (mpc_imagref(op), -1) != 0) continue; expo = err; /* will leave err unchanged below */ } @@ -283,7 +283,7 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* Imaginary part Im(atan(x+I*y)) = 1/4 * [log(x^2+(1+y)^2) - log (x^2 +(1-y)^2)] */ - prec = mpfr_get_prec (MPC_IM (rop)); /* result precision */ + prec = mpfr_get_prec (mpc_imagref (rop)); /* result precision */ /* a = o(1+y) error(a) < 1 ulp(a) b = o(a^2) error(b) < 5 ulp(b) @@ -304,7 +304,7 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) */ err = 2; p = prec; /* working precision */ - rnd1 = mpfr_cmp_si (MPC_IM (op), -1) > 0 ? GMP_RNDU : GMP_RNDD; + rnd1 = mpfr_cmp_si (mpc_imagref (op), -1) > 0 ? GMP_RNDU : GMP_RNDD; do { @@ -314,16 +314,16 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_set_prec (y, p); /* a = upper bound for log(x^2 + (1+y)^2) */ - mpfr_add_ui (a, MPC_IM (op), 1, rnd1); + mpfr_add_ui (a, mpc_imagref (op), 1, rnd1); mpfr_sqr (a, a, GMP_RNDU); - mpfr_sqr (y, MPC_RE (op), GMP_RNDU); + mpfr_sqr (y, mpc_realref (op), GMP_RNDU); mpfr_add (a, a, y, GMP_RNDU); mpfr_log (a, a, GMP_RNDU); /* b = lower bound for log(x^2 + (1-y)^2) */ - mpfr_ui_sub (b, 1, MPC_IM (op), GMP_RNDZ); + mpfr_ui_sub (b, 1, mpc_imagref (op), GMP_RNDZ); mpfr_sqr (b, b, GMP_RNDU); - /* mpfr_sqr (y, MPC_RE (op), GMP_RNDZ); */ + /* mpfr_sqr (y, mpc_realref (op), GMP_RNDZ); */ mpfr_nextbelow (y); mpfr_add (b, b, y, GMP_RNDZ); mpfr_log (b, b, GMP_RNDZ); diff --git a/src/atanh.c b/src/atanh.c index 0f9d8fd..3e1c448 100644 --- a/src/atanh.c +++ b/src/atanh.c @@ -1,6 +1,6 @@ /* mpc_atanh -- inverse hyperbolic tangent of a complex number. -Copyright (C) 2009 INRIA +Copyright (C) 2009, 2011 INRIA This file is part of GNU MPC. @@ -28,9 +28,9 @@ mpc_atanh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_t tmp; mpc_t z, a; - MPC_RE (z)[0] = MPC_IM (op)[0]; - MPC_IM (z)[0] = MPC_RE (op)[0]; - MPFR_CHANGE_SIGN (MPC_RE (z)); + mpc_realref (z)[0] = mpc_imagref (op)[0]; + mpc_imagref (z)[0] = mpc_realref (op)[0]; + MPFR_CHANGE_SIGN (mpc_realref (z)); /* Note reversal of precisions due to later multiplication by -i */ mpc_init3 (a, MPC_PREC_IM(rop), MPC_PREC_RE(rop)); @@ -39,10 +39,10 @@ mpc_atanh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) RNDC (INV_RND (MPC_RND_IM (rnd)), MPC_RND_RE (rnd))); /* change a to -i*a, i.e., x+i*y to y-i*x */ - tmp[0] = MPC_RE (a)[0]; - MPC_RE (a)[0] = MPC_IM (a)[0]; - MPC_IM (a)[0] = tmp[0]; - MPFR_CHANGE_SIGN (MPC_IM (a)); + tmp[0] = mpc_realref (a)[0]; + mpc_realref (a)[0] = mpc_imagref (a)[0]; + mpc_imagref (a)[0] = tmp[0]; + MPFR_CHANGE_SIGN (mpc_imagref (a)); mpc_set (rop, a, rnd); diff --git a/src/clear.c b/src/clear.c index 69355ca..f76b5db 100644 --- a/src/clear.c +++ b/src/clear.c @@ -1,6 +1,6 @@ /* mpc_clear -- Clear a complex variable. -Copyright (C) 2002, 2009 INRIA +Copyright (C) 2002, 2009, 2011 INRIA This file is part of GNU MPC. @@ -23,6 +23,6 @@ along with this program. If not, see http://www.gnu.org/licenses/ . void mpc_clear (mpc_t x) { - mpfr_clear (MPC_RE(x)); - mpfr_clear (MPC_IM(x)); + mpfr_clear (mpc_realref(x)); + mpfr_clear (mpc_imagref(x)); } @@ -1,6 +1,6 @@ /* mpc_cmp -- Compare two complex numbers. -Copyright (C) 2002, 2009, 2010 INRIA +Copyright (C) 2002, 2009, 2010, 2011 INRIA This file is part of GNU MPC. @@ -26,8 +26,8 @@ mpc_cmp (mpc_srcptr a, mpc_srcptr b) { int cmp_re, cmp_im; - cmp_re = mpfr_cmp (MPC_RE(a), MPC_RE(b)); - cmp_im = mpfr_cmp (MPC_IM(a), MPC_IM(b)); + cmp_re = mpfr_cmp (mpc_realref(a), mpc_realref(b)); + cmp_im = mpfr_cmp (mpc_imagref(a), mpc_imagref(b)); return MPC_INEX(cmp_re, cmp_im); } diff --git a/src/cmp_si_si.c b/src/cmp_si_si.c index eb94e10..a50b758 100644 --- a/src/cmp_si_si.c +++ b/src/cmp_si_si.c @@ -1,7 +1,7 @@ /* mpc_cmp_si_si -- Compare a complex number to a number of the form b+c*i with b and c signed integers. -Copyright (C) 2005, 2009 INRIA +Copyright (C) 2005, 2009, 2011 INRIA This file is part of GNU MPC. @@ -27,8 +27,8 @@ mpc_cmp_si_si (mpc_srcptr a, long int b, long int c) { int cmp_re, cmp_im; - cmp_re = mpfr_cmp_si (MPC_RE(a), b); - cmp_im = mpfr_cmp_si (MPC_IM(a), c); + cmp_re = mpfr_cmp_si (mpc_realref(a), b); + cmp_im = mpfr_cmp_si (mpc_imagref(a), c); return MPC_INEX(cmp_re, cmp_im); } @@ -1,6 +1,6 @@ /* mpc_conj -- Conjugate of a complex number. -Copyright (C) 2002, 2009 INRIA +Copyright (C) 2002, 2009, 2011 INRIA This file is part of GNU MPC. @@ -25,8 +25,8 @@ mpc_conj (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_set (MPC_RE(a), MPC_RE(b), MPC_RND_RE(rnd)); - inex_im = mpfr_neg (MPC_IM(a), MPC_IM(b), MPC_RND_IM(rnd)); + inex_re = mpfr_set (mpc_realref(a), mpc_realref(b), MPC_RND_RE(rnd)); + inex_im = mpfr_neg (mpc_imagref(a), mpc_imagref(b), MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } @@ -1,6 +1,6 @@ /* mpc_cosh -- hyperbolic cosine of a complex number. -Copyright (C) 2008, 2009 INRIA +Copyright (C) 2008, 2009, 2011 INRIA This file is part of GNU MPC. @@ -27,9 +27,9 @@ mpc_cosh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpc_t z; /* z = i*op without copying significand */ - MPC_RE (z)[0] = MPC_IM (op)[0]; - MPC_IM (z)[0] = MPC_RE (op)[0]; - MPFR_CHANGE_SIGN (MPC_RE (z)); + mpc_realref (z)[0] = mpc_imagref (op)[0]; + mpc_imagref (z)[0] = mpc_realref (op)[0]; + MPFR_CHANGE_SIGN (mpc_realref (z)); return mpc_cos (rop, z, rnd); } @@ -24,11 +24,11 @@ static int mpc_div_zero (mpc_ptr a, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd) /* Assumes w==0, implementation according to C99 G.5.1.8 */ { - int sign = MPFR_SIGNBIT (MPC_RE (w)); + int sign = MPFR_SIGNBIT (mpc_realref (w)); mpfr_t infty; mpfr_set_inf (infty, sign); - mpfr_mul (MPC_RE (a), infty, MPC_RE (z), MPC_RND_RE (rnd)); - mpfr_mul (MPC_IM (a), infty, MPC_IM (z), MPC_RND_IM (rnd)); + mpfr_mul (mpc_realref (a), infty, mpc_realref (z), MPC_RND_RE (rnd)); + mpfr_mul (mpc_imagref (a), infty, mpc_imagref (z), MPC_RND_IM (rnd)); return MPC_INEX (0, 0); /* exact */ } @@ -39,14 +39,14 @@ mpc_div_inf_fin (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w) { int a, b, x, y; - a = (mpfr_inf_p (MPC_RE (z)) ? MPFR_SIGNBIT (MPC_RE (z)) : 0); - b = (mpfr_inf_p (MPC_IM (z)) ? MPFR_SIGNBIT (MPC_IM (z)) : 0); + a = (mpfr_inf_p (mpc_realref (z)) ? MPFR_SIGNBIT (mpc_realref (z)) : 0); + b = (mpfr_inf_p (mpc_imagref (z)) ? MPFR_SIGNBIT (mpc_imagref (z)) : 0); - /* x = MPC_MPFR_SIGN (a * MPC_RE (w) + b * MPC_IM (w)) */ - /* y = MPC_MPFR_SIGN (b * MPC_RE (w) - a * MPC_IM (w)) */ + /* x = MPC_MPFR_SIGN (a * mpc_realref (w) + b * mpc_imagref (w)) */ + /* y = MPC_MPFR_SIGN (b * mpc_realref (w) - a * mpc_imagref (w)) */ if (a == 0 || b == 0) { - x = a * MPC_MPFR_SIGN (MPC_RE (w)) + b * MPC_MPFR_SIGN (MPC_IM (w)); - y = b * MPC_MPFR_SIGN (MPC_RE (w)) - a * MPC_MPFR_SIGN (MPC_IM (w)); + x = a * MPC_MPFR_SIGN (mpc_realref (w)) + b * MPC_MPFR_SIGN (mpc_imagref (w)); + y = b * MPC_MPFR_SIGN (mpc_realref (w)) - a * MPC_MPFR_SIGN (mpc_imagref (w)); } else { /* Both parts of z are infinite; x could be determined by sign @@ -58,41 +58,41 @@ mpc_div_inf_fin (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w) if (a == 1) if (b == 1) { - mpfr_add (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN); + mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); x = MPC_MPFR_SIGN (sign); - mpfr_sub (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN); + mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); y = MPC_MPFR_SIGN (sign); } else { /* b == -1 */ - mpfr_sub (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN); + mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); x = MPC_MPFR_SIGN (sign); - mpfr_add (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN); + mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); y = -MPC_MPFR_SIGN (sign); } else /* a == -1 */ if (b == 1) { - mpfr_sub (sign, MPC_IM (w), MPC_RE (w), GMP_RNDN); + mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), GMP_RNDN); x = MPC_MPFR_SIGN (sign); - mpfr_add (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN); + mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); y = MPC_MPFR_SIGN (sign); } else { /* b == -1 */ - mpfr_add (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN); + mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); x = -MPC_MPFR_SIGN (sign); - mpfr_sub (sign, MPC_IM (w), MPC_RE (w), GMP_RNDN); + mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), GMP_RNDN); y = MPC_MPFR_SIGN (sign); } mpfr_clear (sign); } if (x == 0) - mpfr_set_nan (MPC_RE (rop)); + mpfr_set_nan (mpc_realref (rop)); else - mpfr_set_inf (MPC_RE (rop), x); + mpfr_set_inf (mpc_realref (rop), x); if (y == 0) - mpfr_set_nan (MPC_IM (rop)); + mpfr_set_nan (mpc_imagref (rop)); else - mpfr_set_inf (MPC_IM (rop), y); + mpfr_set_inf (mpc_imagref (rop), y); return MPC_INEX (0, 0); /* exact */ } @@ -111,24 +111,24 @@ mpc_div_fin_inf (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w) mpfr_init2 (y, 2); mpfr_init2 (zero, 2); mpfr_set_ui (zero, 0ul, GMP_RNDN); - mpfr_init2 (a, mpfr_get_prec (MPC_RE (z))); - mpfr_init2 (b, mpfr_get_prec (MPC_IM (z))); + mpfr_init2 (a, mpfr_get_prec (mpc_realref (z))); + mpfr_init2 (b, mpfr_get_prec (mpc_imagref (z))); - mpfr_set_ui (c, (mpfr_inf_p (MPC_RE (w)) ? 1 : 0), GMP_RNDN); - MPFR_COPYSIGN (c, c, MPC_RE (w), GMP_RNDN); - mpfr_set_ui (d, (mpfr_inf_p (MPC_IM (w)) ? 1 : 0), GMP_RNDN); - MPFR_COPYSIGN (d, d, MPC_IM (w), GMP_RNDN); + mpfr_set_ui (c, (mpfr_inf_p (mpc_realref (w)) ? 1 : 0), GMP_RNDN); + MPFR_COPYSIGN (c, c, mpc_realref (w), GMP_RNDN); + mpfr_set_ui (d, (mpfr_inf_p (mpc_imagref (w)) ? 1 : 0), GMP_RNDN); + MPFR_COPYSIGN (d, d, mpc_imagref (w), GMP_RNDN); - mpfr_mul (a, MPC_RE (z), c, GMP_RNDN); /* exact */ - mpfr_mul (b, MPC_IM (z), d, GMP_RNDN); + mpfr_mul (a, mpc_realref (z), c, GMP_RNDN); /* exact */ + mpfr_mul (b, mpc_imagref (z), d, GMP_RNDN); mpfr_add (x, a, b, GMP_RNDN); - mpfr_mul (b, MPC_IM (z), c, GMP_RNDN); - mpfr_mul (a, MPC_RE (z), d, GMP_RNDN); + mpfr_mul (b, mpc_imagref (z), c, GMP_RNDN); + mpfr_mul (a, mpc_realref (z), d, GMP_RNDN); mpfr_sub (y, b, a, GMP_RNDN); - MPFR_COPYSIGN (MPC_RE (rop), zero, x, GMP_RNDN); - MPFR_COPYSIGN (MPC_IM (rop), zero, y, GMP_RNDN); + MPFR_COPYSIGN (mpc_realref (rop), zero, x, GMP_RNDN); + MPFR_COPYSIGN (mpc_imagref (rop), zero, y, GMP_RNDN); mpfr_clear (c); mpfr_clear (d); @@ -149,22 +149,22 @@ mpc_div_real (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd) { int inex_re, inex_im; /* save signs of operands in case there are overlaps */ - int zrs = MPFR_SIGNBIT (MPC_RE (z)); - int zis = MPFR_SIGNBIT (MPC_IM (z)); - int wrs = MPFR_SIGNBIT (MPC_RE (w)); - int wis = MPFR_SIGNBIT (MPC_IM (w)); + int zrs = MPFR_SIGNBIT (mpc_realref (z)); + int zis = MPFR_SIGNBIT (mpc_imagref (z)); + int wrs = MPFR_SIGNBIT (mpc_realref (w)); + int wis = MPFR_SIGNBIT (mpc_imagref (w)); /* warning: rop may overlap with z,w so treat the imaginary part first */ - inex_im = mpfr_div (MPC_IM(rop), MPC_IM(z), MPC_RE(w), MPC_RND_IM(rnd)); - inex_re = mpfr_div (MPC_RE(rop), MPC_RE(z), MPC_RE(w), MPC_RND_RE(rnd)); + inex_im = mpfr_div (mpc_imagref(rop), mpc_imagref(z), mpc_realref(w), MPC_RND_IM(rnd)); + inex_re = mpfr_div (mpc_realref(rop), mpc_realref(z), mpc_realref(w), MPC_RND_RE(rnd)); /* correct signs of zeroes if necessary, which does not affect the inexact flags */ - if (mpfr_zero_p (MPC_RE (rop))) - mpfr_setsign (MPC_RE (rop), MPC_RE (rop), (zrs != wrs && zis != wis), + if (mpfr_zero_p (mpc_realref (rop))) + mpfr_setsign (mpc_realref (rop), mpc_realref (rop), (zrs != wrs && zis != wis), GMP_RNDN); /* exact */ - if (mpfr_zero_p (MPC_IM (rop))) - mpfr_setsign (MPC_IM (rop), MPC_IM (rop), (zis != wrs && zrs == wis), + if (mpfr_zero_p (mpc_imagref (rop))) + mpfr_setsign (mpc_imagref (rop), mpc_imagref (rop), (zis != wrs && zrs == wis), GMP_RNDN); return MPC_INEX(inex_re, inex_im); @@ -178,24 +178,24 @@ mpc_div_imag (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd) { int inex_re, inex_im; int overlap = (rop == z) || (rop == w); - int imag_z = mpfr_zero_p (MPC_RE (z)); + int imag_z = mpfr_zero_p (mpc_realref (z)); mpfr_t wloc; mpc_t tmprop; mpc_ptr dest = (overlap) ? tmprop : rop; /* save signs of operands in case there are overlaps */ - int zrs = MPFR_SIGNBIT (MPC_RE (z)); - int zis = MPFR_SIGNBIT (MPC_IM (z)); - int wrs = MPFR_SIGNBIT (MPC_RE (w)); - int wis = MPFR_SIGNBIT (MPC_IM (w)); + int zrs = MPFR_SIGNBIT (mpc_realref (z)); + int zis = MPFR_SIGNBIT (mpc_imagref (z)); + int wrs = MPFR_SIGNBIT (mpc_realref (w)); + int wis = MPFR_SIGNBIT (mpc_imagref (w)); if (overlap) mpc_init3 (tmprop, MPC_PREC_RE (rop), MPC_PREC_IM (rop)); - wloc[0] = MPC_IM(w)[0]; /* copies mpfr struct IM(w) into wloc */ - inex_re = mpfr_div (MPC_RE(dest), MPC_IM(z), wloc, MPC_RND_RE(rnd)); + wloc[0] = mpc_imagref(w)[0]; /* copies mpfr struct IM(w) into wloc */ + inex_re = mpfr_div (mpc_realref(dest), mpc_imagref(z), wloc, MPC_RND_RE(rnd)); mpfr_neg (wloc, wloc, GMP_RNDN); /* changes the sign only in wloc, not in w; no need to correct later */ - inex_im = mpfr_div (MPC_IM(dest), MPC_RE(z), wloc, MPC_RND_IM(rnd)); + inex_im = mpfr_div (mpc_imagref(dest), mpc_realref(z), wloc, MPC_RND_IM(rnd)); if (overlap) { /* Note: we could use mpc_swap here, but this might cause problems @@ -208,11 +208,11 @@ mpc_div_imag (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd) /* correct signs of zeroes if necessary, which does not affect the inexact flags */ - if (mpfr_zero_p (MPC_RE (rop))) - mpfr_setsign (MPC_RE (rop), MPC_RE (rop), (zrs != wrs && zis != wis), + if (mpfr_zero_p (mpc_realref (rop))) + mpfr_setsign (mpc_realref (rop), mpc_realref (rop), (zrs != wrs && zis != wis), GMP_RNDN); /* exact */ if (imag_z) - mpfr_setsign (MPC_IM (rop), MPC_IM (rop), (zis != wrs && zrs == wis), + mpfr_setsign (mpc_imagref (rop), mpc_imagref (rop), (zis != wrs && zrs == wis), GMP_RNDN); return MPC_INEX(inex_re, inex_im); @@ -250,9 +250,9 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) mpc_set_nan (a); return MPC_INEX (0, 0); } - else if (mpfr_zero_p(MPC_IM(c))) + else if (mpfr_zero_p(mpc_imagref(c))) return mpc_div_real (a, b, c, rnd); - else if (mpfr_zero_p(MPC_RE(c))) + else if (mpfr_zero_p(mpc_realref(c))) return mpc_div_imag (a, b, c, rnd); prec = MPC_MAX_PREC(a); @@ -261,9 +261,9 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) mpfr_init (q); /* create the conjugate of c in c_conj without allocating new memory */ - MPC_RE (c_conj)[0] = MPC_RE (c)[0]; - MPC_IM (c_conj)[0] = MPC_IM (c)[0]; - MPFR_CHANGE_SIGN (MPC_IM (c_conj)); + mpc_realref (c_conj)[0] = mpc_realref (c)[0]; + mpc_imagref (c_conj)[0] = mpc_imagref (c)[0]; + MPFR_CHANGE_SIGN (mpc_imagref (c_conj)); do { loops ++; @@ -298,10 +298,10 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) following is not totally rigorous and just a best effort to salvage almost hopeless situations */ if (overflow_prod) { - if (!mpfr_zero_p (MPC_RE (res))) - mpfr_set_inf (MPC_RE (res), mpfr_sgn (MPC_RE (res))); - if (!mpfr_zero_p (MPC_IM (res))) - mpfr_set_inf (MPC_IM (res), mpfr_sgn (MPC_IM (res))); + if (!mpfr_zero_p (mpc_realref (res))) + mpfr_set_inf (mpc_realref (res), mpfr_sgn (mpc_realref (res))); + if (!mpfr_zero_p (mpc_imagref (res))) + mpfr_set_inf (mpc_imagref (res), mpfr_sgn (mpc_imagref (res))); } /* divide the product by the norm */ @@ -311,21 +311,21 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) /* we use the division code of mpfr, which handles the situation. */ mpfr_clear_underflow (); mpfr_clear_overflow (); - inexact_re |= mpfr_div (MPC_RE (res), MPC_RE (res), q, GMP_RNDZ); + inexact_re |= mpfr_div (mpc_realref (res), mpc_realref (res), q, GMP_RNDZ); underflow_re = mpfr_underflow_p (); overflow_re = mpfr_overflow_p (); ok_re = !inexact_re || underflow_re || overflow_re - || mpfr_can_round (MPC_RE (res), prec - 4, GMP_RNDN, + || mpfr_can_round (mpc_realref (res), prec - 4, GMP_RNDN, GMP_RNDZ, MPC_PREC_RE(a) + (rnd_re == GMP_RNDN)); if (ok_re) /* compute imaginary part */ { mpfr_clear_underflow (); mpfr_clear_overflow (); - inexact_im |= mpfr_div (MPC_IM (res), MPC_IM (res), q, GMP_RNDZ); + inexact_im |= mpfr_div (mpc_imagref (res), mpc_imagref (res), q, GMP_RNDZ); underflow_im = mpfr_underflow_p (); overflow_im = mpfr_overflow_p (); ok_im = !inexact_im || underflow_im || overflow_im - || mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDN, + || mpfr_can_round (mpc_imagref (res), prec - 4, GMP_RNDN, GMP_RNDZ, MPC_PREC_IM(a) + (rnd_im == GMP_RNDN)); } } @@ -336,26 +336,26 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) /* if 1/q is inexact, the approximations of the real and imaginary part below will be inexact, unless RE(res) or IM(res) is zero */ - inexact_re |= ~mpfr_zero_p (MPC_RE (res)); - inexact_im |= ~mpfr_zero_p (MPC_IM (res)); + inexact_re |= ~mpfr_zero_p (mpc_realref (res)); + inexact_im |= ~mpfr_zero_p (mpc_imagref (res)); } mpfr_clear_underflow (); mpfr_clear_overflow (); - inexact_re |= mpfr_mul (MPC_RE (res), MPC_RE (res), q, GMP_RNDZ); + inexact_re |= mpfr_mul (mpc_realref (res), mpc_realref (res), q, GMP_RNDZ); underflow_re = mpfr_underflow_p (); overflow_re = mpfr_overflow_p (); ok_re = !inexact_re || underflow_re || overflow_re - || mpfr_can_round (MPC_RE (res), prec - 4, GMP_RNDN, + || mpfr_can_round (mpc_realref (res), prec - 4, GMP_RNDN, GMP_RNDZ, MPC_PREC_RE(a) + (rnd_re == GMP_RNDN)); if (ok_re) /* compute imaginary part */ { mpfr_clear_underflow (); mpfr_clear_overflow (); - inexact_im |= mpfr_mul (MPC_IM (res), MPC_IM (res), q, GMP_RNDZ); + inexact_im |= mpfr_mul (mpc_imagref (res), mpc_imagref (res), q, GMP_RNDZ); underflow_im = mpfr_underflow_p (); overflow_im = mpfr_overflow_p (); ok_im = !inexact_im || underflow_im || overflow_im - || mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDN, + || mpfr_can_round (mpc_imagref (res), prec - 4, GMP_RNDN, GMP_RNDZ, MPC_PREC_IM(a) + (rnd_im == GMP_RNDN)); } } @@ -369,20 +369,20 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) /* fix values and inexact flags in case of overflow/underflow */ /* FIXME: heuristic, certainly does not cover all cases */ if (overflow_re || (underflow_norm && !underflow_prod)) { - mpfr_set_inf (MPC_RE (a), mpfr_sgn (MPC_RE (res))); - inexact_re = mpfr_sgn (MPC_RE (res)); + mpfr_set_inf (mpc_realref (a), mpfr_sgn (mpc_realref (res))); + inexact_re = mpfr_sgn (mpc_realref (res)); } else if (underflow_re || (overflow_norm && !overflow_prod)) { - mpfr_set_zero (MPC_RE (a), mpfr_sgn (MPC_RE (res))); - inexact_re = -mpfr_sgn (MPC_RE (res)); + mpfr_set_zero (mpc_realref (a), mpfr_sgn (mpc_realref (res))); + inexact_re = -mpfr_sgn (mpc_realref (res)); } if (overflow_im || (underflow_norm && !underflow_prod)) { - mpfr_set_inf (MPC_IM (a), mpfr_sgn (MPC_IM (res))); - inexact_im = mpfr_sgn (MPC_IM (res)); + mpfr_set_inf (mpc_imagref (a), mpfr_sgn (mpc_imagref (res))); + inexact_im = mpfr_sgn (mpc_imagref (res)); } else if (underflow_im || (overflow_norm && !overflow_prod)) { - mpfr_set_zero (MPC_IM (a), mpfr_sgn (MPC_IM (res))); - inexact_im = -mpfr_sgn (MPC_IM (res)); + mpfr_set_zero (mpc_imagref (a), mpfr_sgn (mpc_imagref (res))); + inexact_im = -mpfr_sgn (mpc_imagref (res)); } mpc_clear (res); diff --git a/src/div_2exp.c b/src/div_2exp.c index a5ebb12..22e67e4 100644 --- a/src/div_2exp.c +++ b/src/div_2exp.c @@ -1,6 +1,6 @@ /* mpc_div_2exp -- Divide a complex number by 2^e. -Copyright (C) 2002, 2009 INRIA +Copyright (C) 2002, 2009, 2011 INRIA This file is part of GNU MPC. @@ -25,8 +25,8 @@ mpc_div_2exp (mpc_ptr a, mpc_srcptr b, unsigned long int c, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_div_2exp (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); - inex_im = mpfr_div_2exp (MPC_IM(a), MPC_IM(b), c, MPC_RND_IM(rnd)); + inex_re = mpfr_div_2exp (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_div_2exp (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } diff --git a/src/div_fr.c b/src/div_fr.c index 399c3b8..d5ea240 100644 --- a/src/div_fr.c +++ b/src/div_fr.c @@ -1,6 +1,6 @@ /* mpc_div_fr -- Divide a complex number by a floating-point number. -Copyright (C) 2002, 2008, 2009, 2010 INRIA +Copyright (C) 2002, 2008, 2009, 2010, 2011 INRIA This file is part of GNU MPC. @@ -26,12 +26,12 @@ mpc_div_fr (mpc_ptr a, mpc_srcptr b, mpfr_srcptr c, mpc_rnd_t rnd) int inex_re, inex_im; mpfr_t real; - /* We have to use temporary variable in case c=MPC_RE (a). */ + /* We have to use temporary variable in case c=mpc_realref (a). */ mpfr_init2 (real, MPC_PREC_RE (a)); - inex_re = mpfr_div (real, MPC_RE(b), c, MPC_RND_RE(rnd)); - inex_im = mpfr_div (MPC_IM(a), MPC_IM(b), c, MPC_RND_IM(rnd)); - mpfr_set (MPC_RE (a), real, GMP_RNDN); + inex_re = mpfr_div (real, mpc_realref(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_div (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd)); + mpfr_set (mpc_realref (a), real, GMP_RNDN); mpfr_clear (real); diff --git a/src/div_ui.c b/src/div_ui.c index d0b3125..26debf7 100644 --- a/src/div_ui.c +++ b/src/div_ui.c @@ -1,6 +1,6 @@ /* mpc_div_ui -- Divide a complex number by a nonnegative integer. -Copyright (C) 2002, 2009 INRIA +Copyright (C) 2002, 2009, 2011 INRIA This file is part of GNU MPC. @@ -25,8 +25,8 @@ mpc_div_ui (mpc_ptr a, mpc_srcptr b, unsigned long int c, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_div_ui (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); - inex_im = mpfr_div_ui (MPC_IM(a), MPC_IM(b), c, MPC_RND_IM(rnd)); + inex_re = mpfr_div_ui (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_div_ui (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } @@ -1,6 +1,6 @@ /* mpc_exp -- exponential of a complex number. -Copyright (C) 2002, 2009, 2010 INRIA +Copyright (C) 2002, 2009, 2010, 2011 INRIA This file is part of GNU MPC. @@ -29,7 +29,7 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) int inex_re, inex_im; /* special values */ - if (mpfr_nan_p (MPC_RE (op)) || mpfr_nan_p (MPC_IM (op))) + if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op))) /* NaNs exp(nan +i*y) = nan -i*0 if y = -0, nan +i*0 if y = +0, @@ -38,46 +38,46 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) +/-inf +i*nan if x=+inf, nan +i*nan otherwise */ { - if (mpfr_zero_p (MPC_IM (op))) + if (mpfr_zero_p (mpc_imagref (op))) return mpc_set (rop, op, MPC_RNDNN); - if (mpfr_inf_p (MPC_RE (op))) + if (mpfr_inf_p (mpc_realref (op))) { - if (mpfr_signbit (MPC_RE (op))) + if (mpfr_signbit (mpc_realref (op))) return mpc_set_ui_ui (rop, 0, 0, MPC_RNDNN); else { - mpfr_set_inf (MPC_RE (rop), +1); - mpfr_set_nan (MPC_IM (rop)); + mpfr_set_inf (mpc_realref (rop), +1); + mpfr_set_nan (mpc_imagref (rop)); return MPC_INEX(0, 0); /* Inf/NaN are exact */ } } - mpfr_set_nan (MPC_RE (rop)); - mpfr_set_nan (MPC_IM (rop)); + mpfr_set_nan (mpc_realref (rop)); + mpfr_set_nan (mpc_imagref (rop)); return MPC_INEX(0, 0); /* NaN is exact */ } - if (mpfr_zero_p (MPC_IM(op))) + if (mpfr_zero_p (mpc_imagref(op))) /* special case when the input is real exp(x-i*0) = exp(x) -i*0, even if x is NaN exp(x+i*0) = exp(x) +i*0, even if x is NaN */ { - inex_re = mpfr_exp (MPC_RE(rop), MPC_RE(op), MPC_RND_RE(rnd)); - inex_im = mpfr_set (MPC_IM(rop), MPC_IM(op), MPC_RND_IM(rnd)); + inex_re = mpfr_exp (mpc_realref(rop), mpc_realref(op), MPC_RND_RE(rnd)); + inex_im = mpfr_set (mpc_imagref(rop), mpc_imagref(op), MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } - if (mpfr_zero_p (MPC_RE (op))) + if (mpfr_zero_p (mpc_realref (op))) /* special case when the input is imaginary */ { - inex_re = mpfr_cos (MPC_RE (rop), MPC_IM (op), MPC_RND_RE(rnd)); - inex_im = mpfr_sin (MPC_IM (rop), MPC_IM (op), MPC_RND_IM(rnd)); + inex_re = mpfr_cos (mpc_realref (rop), mpc_imagref (op), MPC_RND_RE(rnd)); + inex_im = mpfr_sin (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } - if (mpfr_inf_p (MPC_RE (op))) + if (mpfr_inf_p (mpc_realref (op))) /* real part is an infinity, exp(-inf +i*y) = 0*(cos y +i*sin y) exp(+inf +i*y) = +/-inf +i*nan if y = +/-inf @@ -86,19 +86,19 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_t n; mpfr_init2 (n, 2); - if (mpfr_signbit (MPC_RE (op))) + if (mpfr_signbit (mpc_realref (op))) mpfr_set_ui (n, 0, GMP_RNDN); else mpfr_set_inf (n, +1); - if (mpfr_inf_p (MPC_IM (op))) + if (mpfr_inf_p (mpc_imagref (op))) { - inex_re = mpfr_set (MPC_RE (rop), n, GMP_RNDN); - if (mpfr_signbit (MPC_RE (op))) - inex_im = mpfr_set (MPC_IM (rop), n, GMP_RNDN); + inex_re = mpfr_set (mpc_realref (rop), n, GMP_RNDN); + if (mpfr_signbit (mpc_realref (op))) + inex_im = mpfr_set (mpc_imagref (rop), n, GMP_RNDN); else { - mpfr_set_nan (MPC_IM (rop)); + mpfr_set_nan (mpc_imagref (rop)); inex_im = 0; /* NaN is exact */ } } @@ -108,9 +108,9 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_init2 (c, 2); mpfr_init2 (s, 2); - mpfr_sin_cos (s, c, MPC_IM (op), GMP_RNDN); - inex_re = mpfr_copysign (MPC_RE (rop), n, c, GMP_RNDN); - inex_im = mpfr_copysign (MPC_IM (rop), n, s, GMP_RNDN); + mpfr_sin_cos (s, c, mpc_imagref (op), GMP_RNDN); + inex_re = mpfr_copysign (mpc_realref (rop), n, c, GMP_RNDN); + inex_im = mpfr_copysign (mpc_imagref (rop), n, s, GMP_RNDN); mpfr_clear (s); mpfr_clear (c); @@ -120,11 +120,11 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) return MPC_INEX(inex_re, inex_im); } - if (mpfr_inf_p (MPC_IM (op))) + if (mpfr_inf_p (mpc_imagref (op))) /* real part is finite non-zero number, imaginary part is an infinity */ { - mpfr_set_nan (MPC_RE (rop)); - mpfr_set_nan (MPC_IM (rop)); + mpfr_set_nan (mpc_realref (rop)); + mpfr_set_nan (mpc_imagref (rop)); return MPC_INEX(0, 0); /* NaN is exact */ } @@ -132,8 +132,8 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* from now on, both parts of op are regular numbers */ prec = MPC_MAX_PREC(rop) - + MPC_MAX (MPC_MAX (-mpfr_get_exp (MPC_RE (op)), 0), - -mpfr_get_exp (MPC_IM (op))); + + MPC_MAX (MPC_MAX (-mpfr_get_exp (mpc_realref (op)), 0), + -mpfr_get_exp (mpc_imagref (op))); /* When op is close to 0, then exp is close to 1+Re(op), while cos is close to 1-Im(op); to decide on the ternary value of exp*cos, we need a high enough precision so that none of exp or cos is @@ -154,8 +154,8 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) could be represented in the precision of rop. */ mpfr_clear_overflow (); mpfr_clear_underflow (); - mpfr_exp (x, MPC_RE(op), GMP_RNDN); /* error <= 0.5ulp */ - mpfr_sin_cos (z, y, MPC_IM(op), GMP_RNDN); /* errors <= 0.5ulp */ + mpfr_exp (x, mpc_realref(op), GMP_RNDN); /* error <= 0.5ulp */ + mpfr_sin_cos (z, y, mpc_imagref(op), GMP_RNDN); /* errors <= 0.5ulp */ mpfr_mul (y, y, x, GMP_RNDN); /* error <= 2ulp */ ok = mpfr_overflow_p () || mpfr_zero_p (x) || mpfr_can_round (y, prec - 2, GMP_RNDN, GMP_RNDZ, @@ -170,8 +170,8 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) } while (ok == 0); - inex_re = mpfr_set (MPC_RE(rop), y, MPC_RND_RE(rnd)); - inex_im = mpfr_set (MPC_IM(rop), z, MPC_RND_IM(rnd)); + inex_re = mpfr_set (mpc_realref(rop), y, MPC_RND_RE(rnd)); + inex_im = mpfr_set (mpc_imagref(rop), z, MPC_RND_IM(rnd)); if (mpfr_overflow_p ()) { /* overflow in real exponential, inex is sign of infinite result */ inex_re = mpfr_sgn (y); @@ -46,45 +46,45 @@ mpc_fma (mpc_ptr r, mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) mpfr_prec_t pre12, pre13, pre23, pim12, pim13, pim23; int inex_re, inex_im; - mpfr_init2 (rea_reb, mpfr_get_prec (MPC_RE(a)) + mpfr_get_prec (MPC_RE(b))); - mpfr_init2 (rea_imb, mpfr_get_prec (MPC_RE(a)) + mpfr_get_prec (MPC_IM(b))); - mpfr_init2 (ima_reb, mpfr_get_prec (MPC_IM(a)) + mpfr_get_prec (MPC_RE(b))); - mpfr_init2 (ima_imb, mpfr_get_prec (MPC_IM(a)) + mpfr_get_prec (MPC_IM(b))); + mpfr_init2 (rea_reb, mpfr_get_prec (mpc_realref(a)) + mpfr_get_prec (mpc_realref(b))); + mpfr_init2 (rea_imb, mpfr_get_prec (mpc_realref(a)) + mpfr_get_prec (mpc_imagref(b))); + mpfr_init2 (ima_reb, mpfr_get_prec (mpc_imagref(a)) + mpfr_get_prec (mpc_realref(b))); + mpfr_init2 (ima_imb, mpfr_get_prec (mpc_imagref(a)) + mpfr_get_prec (mpc_imagref(b))); - mpfr_mul (rea_reb, MPC_RE(a), MPC_RE(b), GMP_RNDZ); /* exact */ - mpfr_mul (rea_imb, MPC_RE(a), MPC_IM(b), GMP_RNDZ); /* exact */ - mpfr_mul (ima_reb, MPC_IM(a), MPC_RE(b), GMP_RNDZ); /* exact */ - mpfr_mul (ima_imb, MPC_IM(a), MPC_IM(b), GMP_RNDZ); /* exact */ + mpfr_mul (rea_reb, mpc_realref(a), mpc_realref(b), GMP_RNDZ); /* exact */ + mpfr_mul (rea_imb, mpc_realref(a), mpc_imagref(b), GMP_RNDZ); /* exact */ + mpfr_mul (ima_reb, mpc_imagref(a), mpc_realref(b), GMP_RNDZ); /* exact */ + mpfr_mul (ima_imb, mpc_imagref(a), mpc_imagref(b), GMP_RNDZ); /* exact */ /* Re(r) <- rea_reb - ima_imb + Re(c) */ pre12 = bound_prec_addsub (rea_reb, ima_imb); /* bound on exact precision for rea_reb - ima_imb */ - pre13 = bound_prec_addsub (rea_reb, MPC_RE(c)); + pre13 = bound_prec_addsub (rea_reb, mpc_realref(c)); /* bound for rea_reb + Re(c) */ - pre23 = bound_prec_addsub (ima_imb, MPC_RE(c)); + pre23 = bound_prec_addsub (ima_imb, mpc_realref(c)); /* bound for ima_imb - Re(c) */ if (pre12 <= pre13 && pre12 <= pre23) /* (rea_reb - ima_imb) + Re(c) */ { mpfr_init2 (tmp, pre12); mpfr_sub (tmp, rea_reb, ima_imb, GMP_RNDZ); /* exact */ - inex_re = mpfr_add (MPC_RE(r), tmp, MPC_RE(c), MPC_RND_RE(rnd)); + inex_re = mpfr_add (mpc_realref(r), tmp, mpc_realref(c), MPC_RND_RE(rnd)); /* the only possible bad overlap is between r and c, but since we are only touching the real part of both, it is ok */ } else if (pre13 <= pre23) /* (rea_reb + Re(c)) - ima_imb */ { mpfr_init2 (tmp, pre13); - mpfr_add (tmp, rea_reb, MPC_RE(c), GMP_RNDZ); /* exact */ - inex_re = mpfr_sub (MPC_RE(r), tmp, ima_imb, MPC_RND_RE(rnd)); + mpfr_add (tmp, rea_reb, mpc_realref(c), GMP_RNDZ); /* exact */ + inex_re = mpfr_sub (mpc_realref(r), tmp, ima_imb, MPC_RND_RE(rnd)); /* the only possible bad overlap is between r and c, but since we are only touching the real part of both, it is ok */ } else /* rea_reb + (Re(c) - ima_imb) */ { mpfr_init2 (tmp, pre23); - mpfr_sub (tmp, MPC_RE(c), ima_imb, GMP_RNDZ); /* exact */ - inex_re = mpfr_add (MPC_RE(r), tmp, rea_reb, MPC_RND_RE(rnd)); + mpfr_sub (tmp, mpc_realref(c), ima_imb, GMP_RNDZ); /* exact */ + inex_re = mpfr_add (mpc_realref(r), tmp, rea_reb, MPC_RND_RE(rnd)); /* the only possible bad overlap is between r and c, but since we are only touching the real part of both, it is ok */ } @@ -92,31 +92,31 @@ mpc_fma (mpc_ptr r, mpc_srcptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) /* Im(r) <- rea_imb + ima_reb + Im(c) */ pim12 = bound_prec_addsub (rea_imb, ima_reb); /* bound on exact precision for rea_imb + ima_reb */ - pim13 = bound_prec_addsub (rea_imb, MPC_IM(c)); + pim13 = bound_prec_addsub (rea_imb, mpc_imagref(c)); /* bound for rea_imb + Im(c) */ - pim23 = bound_prec_addsub (ima_reb, MPC_IM(c)); + pim23 = bound_prec_addsub (ima_reb, mpc_imagref(c)); /* bound for ima_reb + Im(c) */ if (pim12 <= pim13 && pim12 <= pim23) /* (rea_imb + ima_reb) + Im(c) */ { mpfr_set_prec (tmp, pim12); mpfr_add (tmp, rea_imb, ima_reb, GMP_RNDZ); /* exact */ - inex_im = mpfr_add (MPC_IM(r), tmp, MPC_IM(c), MPC_RND_IM(rnd)); + inex_im = mpfr_add (mpc_imagref(r), tmp, mpc_imagref(c), MPC_RND_IM(rnd)); /* the only possible bad overlap is between r and c, but since we are only touching the imaginary part of both, it is ok */ } else if (pim13 <= pim23) /* (rea_imb + Im(c)) + ima_reb */ { mpfr_set_prec (tmp, pim13); - mpfr_add (tmp, rea_imb, MPC_IM(c), GMP_RNDZ); /* exact */ - inex_im = mpfr_add (MPC_IM(r), tmp, ima_reb, MPC_RND_IM(rnd)); + mpfr_add (tmp, rea_imb, mpc_imagref(c), GMP_RNDZ); /* exact */ + inex_im = mpfr_add (mpc_imagref(r), tmp, ima_reb, MPC_RND_IM(rnd)); /* the only possible bad overlap is between r and c, but since we are only touching the imaginary part of both, it is ok */ } else /* rea_imb + (Im(c) + ima_reb) */ { mpfr_set_prec (tmp, pre23); - mpfr_add (tmp, MPC_IM(c), ima_reb, GMP_RNDZ); /* exact */ - inex_im = mpfr_add (MPC_IM(r), tmp, rea_imb, MPC_RND_IM(rnd)); + mpfr_add (tmp, mpc_imagref(c), ima_reb, GMP_RNDZ); /* exact */ + inex_im = mpfr_add (mpc_imagref(r), tmp, rea_imb, MPC_RND_IM(rnd)); /* the only possible bad overlap is between r and c, but since we are only touching the imaginary part of both, it is ok */ } diff --git a/src/fr_div.c b/src/fr_div.c index b8201e1..e57eced 100644 --- a/src/fr_div.c +++ b/src/fr_div.c @@ -1,6 +1,6 @@ /* mpc_fr_div -- Divide a floating-point number by a complex number. -Copyright (C) 2008, 2009 INRIA +Copyright (C) 2008, 2009, 2011 INRIA This file is part of GNU MPC. @@ -26,14 +26,14 @@ mpc_fr_div (mpc_ptr a, mpfr_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) mpc_t bc; int inexact; - MPC_RE (bc)[0] = b [0]; - mpfr_init (MPC_IM (bc)); + mpc_realref (bc)[0] = b [0]; + mpfr_init (mpc_imagref (bc)); /* we consider the operand b to have imaginary part +0 */ - mpfr_set_ui (MPC_IM (bc), 0, GMP_RNDN); + mpfr_set_ui (mpc_imagref (bc), 0, GMP_RNDN); inexact = mpc_div (a, bc, c, rnd); - mpfr_clear (MPC_IM (bc)); + mpfr_clear (mpc_imagref (bc)); return inexact; } diff --git a/src/fr_sub.c b/src/fr_sub.c index 9738e4a..b39336c 100644 --- a/src/fr_sub.c +++ b/src/fr_sub.c @@ -1,6 +1,6 @@ /* mpc_fr_sub -- Substract a complex number from a floating-point number. -Copyright (C) 2008, 2009 INRIA +Copyright (C) 2008, 2009, 2011 INRIA This file is part of GNU MPC. @@ -26,9 +26,9 @@ mpc_fr_sub (mpc_ptr a, mpfr_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_sub (MPC_RE(a), b, MPC_RE(c), MPC_RND_RE(rnd)); - inex_im = mpfr_set (MPC_IM (a), MPC_IM (c), INV_RND (MPC_RND_IM (rnd))); - MPFR_CHANGE_SIGN (MPC_IM (a)); + inex_re = mpfr_sub (mpc_realref(a), b, mpc_realref(c), MPC_RND_RE(rnd)); + inex_im = mpfr_set (mpc_imagref (a), mpc_imagref (c), INV_RND (MPC_RND_IM (rnd))); + MPFR_CHANGE_SIGN (mpc_imagref (a)); return MPC_INEX(inex_re, inex_im); } diff --git a/src/get_x.c b/src/get_x.c index 972376d..31610ac 100644 --- a/src/get_x.c +++ b/src/get_x.c @@ -215,8 +215,8 @@ mpc_get_str (int base, size_t n, mpc_srcptr op, mpc_rnd_t rnd) if (base < 2 || base > 36) return NULL; - real_str = get_pretty_str (base, n, MPC_RE (op), MPC_RND_RE (rnd)); - imag_str = get_pretty_str (base, n, MPC_IM (op), MPC_RND_IM (rnd)); + real_str = get_pretty_str (base, n, mpc_realref (op), MPC_RND_RE (rnd)); + imag_str = get_pretty_str (base, n, mpc_imagref (op), MPC_RND_IM (rnd)); needed_size = strlen (real_str) + strlen (imag_str) + 4; @@ -1,6 +1,6 @@ -/* mpc_imag -- Get the real part of a complex number. +/* mpc_imag -- Get the imaginary part of a complex number. -Copyright (C) 2008, 2009 INRIA +Copyright (C) 2008, 2009, 2011 INRIA This file is part of GNU MPC. @@ -23,5 +23,5 @@ along with this program. If not, see http://www.gnu.org/licenses/ . int mpc_imag (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd) { - return mpfr_set (a, MPC_IM (b), rnd); + return mpfr_set (a, mpc_imagref (b), rnd); } diff --git a/src/init2.c b/src/init2.c index dae871a..ce4173e 100644 --- a/src/init2.c +++ b/src/init2.c @@ -1,6 +1,6 @@ /* mpc_init2 -- Initialize a complex variable with a given precision. -Copyright (C) 2002, 2009 INRIA +Copyright (C) 2002, 2009, 2011 INRIA This file is part of GNU MPC. @@ -23,6 +23,6 @@ along with this program. If not, see http://www.gnu.org/licenses/ . void mpc_init2 (mpc_t x, mpfr_prec_t prec) { - mpfr_init2 (MPC_RE(x), prec); - mpfr_init2 (MPC_IM(x), prec); + mpfr_init2 (mpc_realref(x), prec); + mpfr_init2 (mpc_imagref(x), prec); } diff --git a/src/init3.c b/src/init3.c index bd291e8..69f91b2 100644 --- a/src/init3.c +++ b/src/init3.c @@ -1,6 +1,6 @@ /* mpc_init3 -- Initialize a complex variable with given precisions. -Copyright (C) 2002, 2009 INRIA +Copyright (C) 2002, 2009, 2011 INRIA This file is part of GNU MPC. @@ -23,6 +23,6 @@ along with this program. If not, see http://www.gnu.org/licenses/ . void mpc_init3 (mpc_t x, mpfr_prec_t prec_re, mpfr_prec_t prec_im) { - mpfr_init2 (MPC_RE(x), prec_re); - mpfr_init2 (MPC_IM(x), prec_im); + mpfr_init2 (mpc_realref(x), prec_re); + mpfr_init2 (mpc_imagref(x), prec_im); } diff --git a/src/inp_str.c b/src/inp_str.c index 39d8625..695a3ad 100644 --- a/src/inp_str.c +++ b/src/inp_str.c @@ -1,6 +1,6 @@ /* mpc_inp_str -- Input a complex number from a given stream. -Copyright (C) 2009, 2010 INRIA +Copyright (C) 2009, 2010, 2011 INRIA This file is part of GNU MPC. @@ -230,8 +230,8 @@ mpc_rnd_t rnd_mode) error: if (inex == -1) { - mpfr_set_nan (MPC_RE(rop)); - mpfr_set_nan (MPC_IM(rop)); + mpfr_set_nan (mpc_realref(rop)); + mpfr_set_nan (mpc_imagref(rop)); } if (read != NULL) *read = white + nread; @@ -1,6 +1,6 @@ /* mpc_log -- Take the logarithm of a complex number. -Copyright (C) 2008, 2009, 2010 INRIA +Copyright (C) 2008, 2009, 2010, 2011 INRIA This file is part of GNU MPC. @@ -31,58 +31,58 @@ mpc_log (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd){ /* special values: NaN and infinities */ if (!mpc_fin_p (op)) { - if (mpfr_nan_p (MPC_RE (op))) { - if (mpfr_inf_p (MPC_IM (op))) - mpfr_set_inf (MPC_RE (rop), +1); + if (mpfr_nan_p (mpc_realref (op))) { + if (mpfr_inf_p (mpc_imagref (op))) + mpfr_set_inf (mpc_realref (rop), +1); else - mpfr_set_nan (MPC_RE (rop)); - mpfr_set_nan (MPC_IM (rop)); + mpfr_set_nan (mpc_realref (rop)); + mpfr_set_nan (mpc_imagref (rop)); inex_im = 0; /* Inf/NaN is exact */ } - else if (mpfr_nan_p (MPC_IM (op))) { - if (mpfr_inf_p (MPC_RE (op))) - mpfr_set_inf (MPC_RE (rop), +1); + else if (mpfr_nan_p (mpc_imagref (op))) { + if (mpfr_inf_p (mpc_realref (op))) + mpfr_set_inf (mpc_realref (rop), +1); else - mpfr_set_nan (MPC_RE (rop)); - mpfr_set_nan (MPC_IM (rop)); + mpfr_set_nan (mpc_realref (rop)); + mpfr_set_nan (mpc_imagref (rop)); inex_im = 0; /* Inf/NaN is exact */ } else /* We have an infinity in at least one part. */ { - inex_im = mpfr_atan2 (MPC_IM (rop), MPC_IM (op), MPC_RE (op), + inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), mpc_realref (op), MPC_RND_IM (rnd)); - mpfr_set_inf (MPC_RE (rop), +1); + mpfr_set_inf (mpc_realref (rop), +1); } return MPC_INEX(0, inex_im); } /* special cases: real and purely imaginary numbers */ - re_cmp = mpfr_cmp_ui (MPC_RE (op), 0); - im_cmp = mpfr_cmp_ui (MPC_IM (op), 0); + re_cmp = mpfr_cmp_ui (mpc_realref (op), 0); + im_cmp = mpfr_cmp_ui (mpc_imagref (op), 0); if (im_cmp == 0) { if (re_cmp == 0) { - inex_im = mpfr_atan2 (MPC_IM (rop), MPC_IM (op), MPC_RE (op), + inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), mpc_realref (op), MPC_RND_IM (rnd)); - mpfr_set_inf (MPC_RE (rop), -1); + mpfr_set_inf (mpc_realref (rop), -1); inex_re = 0; /* -Inf is exact */ } else if (re_cmp > 0) { - inex_re = mpfr_log (MPC_RE (rop), MPC_RE (op), MPC_RND_RE (rnd)); - inex_im = mpfr_set (MPC_IM (rop), MPC_IM (op), MPC_RND_IM (rnd)); + inex_re = mpfr_log (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd)); + inex_im = mpfr_set (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd)); } else { /* op = x + 0*y; let w = -x = |x| */ int negative_zero; mpfr_rnd_t rnd_im; - negative_zero = mpfr_signbit (MPC_IM (op)); + negative_zero = mpfr_signbit (mpc_imagref (op)); if (negative_zero) rnd_im = INV_RND (MPC_RND_IM (rnd)); else rnd_im = MPC_RND_IM (rnd); - w [0] = *MPC_RE (op); + w [0] = *mpc_realref (op); MPFR_CHANGE_SIGN (w); - inex_re = mpfr_log (MPC_RE (rop), w, MPC_RND_RE (rnd)); - inex_im = mpfr_const_pi (MPC_IM (rop), rnd_im); + inex_re = mpfr_log (mpc_realref (rop), w, MPC_RND_RE (rnd)); + inex_im = mpfr_const_pi (mpc_imagref (rop), rnd_im); if (negative_zero) { mpc_conj (rop, rop, MPC_RNDNN); inex_im = -inex_im; @@ -92,19 +92,19 @@ mpc_log (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd){ } else if (re_cmp == 0) { if (im_cmp > 0) { - inex_re = mpfr_log (MPC_RE (rop), MPC_IM (op), MPC_RND_RE (rnd)); - inex_im = mpfr_const_pi (MPC_IM (rop), MPC_RND_IM (rnd)); + inex_re = mpfr_log (mpc_realref (rop), mpc_imagref (op), MPC_RND_RE (rnd)); + inex_im = mpfr_const_pi (mpc_imagref (rop), MPC_RND_IM (rnd)); /* division by 2 does not change the ternary flag */ - mpfr_div_2ui (MPC_IM (rop), MPC_IM (rop), 1, GMP_RNDN); + mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN); } else { - w [0] = *MPC_IM (op); + w [0] = *mpc_imagref (op); MPFR_CHANGE_SIGN (w); - inex_re = mpfr_log (MPC_RE (rop), w, MPC_RND_RE (rnd)); - inex_im = mpfr_const_pi (MPC_IM (rop), INV_RND (MPC_RND_IM (rnd))); + inex_re = mpfr_log (mpc_realref (rop), w, MPC_RND_RE (rnd)); + inex_im = mpfr_const_pi (mpc_imagref (rop), INV_RND (MPC_RND_IM (rnd))); /* division by 2 does not change the ternary flag */ - mpfr_div_2ui (MPC_IM (rop), MPC_IM (rop), 1, GMP_RNDN); - mpfr_neg (MPC_IM (rop), MPC_IM (rop), GMP_RNDN); + mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN); + mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), GMP_RNDN); inex_im = -inex_im; /* negate the ternary flag */ } return MPC_INEX(inex_re, inex_im); @@ -141,11 +141,11 @@ mpc_log (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd){ } while (ok == 0); /* imaginary part */ - inex_im = mpfr_atan2 (MPC_IM (rop), MPC_IM (op), MPC_RE (op), + inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), mpc_realref (op), MPC_RND_IM (rnd)); /* set the real part; cannot be done before when rop==op */ - inex_re = mpfr_div_2ui (MPC_RE(rop), w, 1ul, MPC_RND_RE (rnd)); + inex_re = mpfr_div_2ui (mpc_realref(rop), w, 1ul, MPC_RND_RE (rnd)); mpfr_clear (w); return MPC_INEX(inex_re, inex_im); } diff --git a/src/mpc-impl.h b/src/mpc-impl.h index 1293273..1b087f2 100644 --- a/src/mpc-impl.h +++ b/src/mpc-impl.h @@ -18,8 +18,8 @@ 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/ . */ -#ifndef __MPC_IMPL_H -#define __MPC_IMPL_H +#ifndef __mpc_imagrefPL_H +#define __mpc_imagrefPL_H /* If mpc-impl.h is included, we are either in the library build or in the tests. Define __MPC_WITHIN_MPC in the first case, since Windows DLL needs @@ -34,9 +34,6 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #endif #include "mpc.h" -#define MPC_RE(x) ((x)->re) -#define MPC_IM(x) ((x)->im) - /* * Miscellaneous useful macros */ @@ -106,22 +103,22 @@ __MPC_DECLSPEC int mpfr_regular_p __MPC_PROTO ((mpfr_srcptr)); * MPC macros */ -#define MPC_PREC_RE(x) (mpfr_get_prec(MPC_RE(x))) -#define MPC_PREC_IM(x) (mpfr_get_prec(MPC_IM(x))) +#define MPC_PREC_RE(x) (mpfr_get_prec(mpc_realref(x))) +#define MPC_PREC_IM(x) (mpfr_get_prec(mpc_imagref(x))) #define MPC_MAX_PREC(x) MPC_MAX(MPC_PREC_RE(x), MPC_PREC_IM(x)) #define INV_RND(r) \ (((r) == GMP_RNDU) ? GMP_RNDD : (((r) == GMP_RNDD) ? GMP_RNDU : (r))) -#define mpc_inf_p(z) (mpfr_inf_p(MPC_RE(z))||mpfr_inf_p(MPC_IM(z))) +#define mpc_inf_p(z) (mpfr_inf_p(mpc_realref(z))||mpfr_inf_p(mpc_imagref(z))) /* Convention in C99 (G.3): z is regarded as an infinity if at least one of its parts is infinite */ -#define mpc_zero_p(z) (mpfr_zero_p(MPC_RE(z))&&mpfr_zero_p(MPC_IM(z))) +#define mpc_zero_p(z) (mpfr_zero_p(mpc_realref(z))&&mpfr_zero_p(mpc_imagref(z))) /* Convention in C99 (G.3): z is regarded as a zero if each of its parts is a zero */ -#define mpc_fin_p(z) (mpfr_number_p(MPC_RE(z))&&mpfr_number_p(MPC_IM(z))) +#define mpc_fin_p(z) (mpfr_number_p(mpc_realref(z))&&mpfr_number_p(mpc_imagref(z))) /* Convention in C99 (G.3): z is regarded as finite if both its parts are */ -#define mpc_nan_p(z) ((mpfr_nan_p(MPC_RE(z)) && !mpfr_inf_p(MPC_IM(z))) || (mpfr_nan_p(MPC_IM(z)) && !mpfr_inf_p(MPC_RE(z)))) +#define mpc_nan_p(z) ((mpfr_nan_p(mpc_realref(z)) && !mpfr_inf_p(mpc_imagref(z))) || (mpfr_nan_p(mpc_imagref(z)) && !mpfr_inf_p(mpc_realref(z)))) /* Consider as NaN all other numbers containing at least one NaN */ @@ -33,10 +33,10 @@ static int mul_infinite (mpc_ptr z, mpc_srcptr x, mpc_srcptr y) { /* Let x=xr+i*xi and y=yr+i*yi; extract the signs of the operands */ - int xrs = mpfr_signbit (MPC_RE (x)) ? -1 : 1; - int xis = mpfr_signbit (MPC_IM (x)) ? -1 : 1; - int yrs = mpfr_signbit (MPC_RE (y)) ? -1 : 1; - int yis = mpfr_signbit (MPC_IM (y)) ? -1 : 1; + int xrs = mpfr_signbit (mpc_realref (x)) ? -1 : 1; + int xis = mpfr_signbit (mpc_imagref (x)) ? -1 : 1; + int yrs = mpfr_signbit (mpc_realref (y)) ? -1 : 1; + int yis = mpfr_signbit (mpc_imagref (y)) ? -1 : 1; int u, v; @@ -44,50 +44,50 @@ mul_infinite (mpc_ptr z, mpc_srcptr x, mpc_srcptr y) u = xrs * yrs * xr * yr - xis * yis * xi * yi v = xrs * yis * xr * yi + xis * yrs * xi * yr +1 if positive, -1 if negatiye, 0 if NaN */ - if ( mpfr_nan_p (MPC_RE (x)) || mpfr_nan_p (MPC_IM (x)) - || mpfr_nan_p (MPC_RE (y)) || mpfr_nan_p (MPC_IM (y))) { + if ( mpfr_nan_p (mpc_realref (x)) || mpfr_nan_p (mpc_imagref (x)) + || mpfr_nan_p (mpc_realref (y)) || mpfr_nan_p (mpc_imagref (y))) { u = 0; v = 0; } - else if (mpfr_inf_p (MPC_RE (x))) { + else if (mpfr_inf_p (mpc_realref (x))) { /* x = (+/-inf) xr + i*xi */ - u = ( mpfr_zero_p (MPC_RE (y)) - || (mpfr_inf_p (MPC_IM (x)) && mpfr_zero_p (MPC_IM (y))) - || (mpfr_zero_p (MPC_IM (x)) && mpfr_inf_p (MPC_IM (y))) - || ( (mpfr_inf_p (MPC_IM (x)) || mpfr_inf_p (MPC_IM (y))) + u = ( mpfr_zero_p (mpc_realref (y)) + || (mpfr_inf_p (mpc_imagref (x)) && mpfr_zero_p (mpc_imagref (y))) + || (mpfr_zero_p (mpc_imagref (x)) && mpfr_inf_p (mpc_imagref (y))) + || ( (mpfr_inf_p (mpc_imagref (x)) || mpfr_inf_p (mpc_imagref (y))) && xrs*yrs == xis*yis) ? 0 : xrs * yrs); - v = ( mpfr_zero_p (MPC_IM (y)) - || (mpfr_inf_p (MPC_IM (x)) && mpfr_zero_p (MPC_RE (y))) - || (mpfr_zero_p (MPC_IM (x)) && mpfr_inf_p (MPC_RE (y))) - || ( (mpfr_inf_p (MPC_IM (x)) || mpfr_inf_p (MPC_IM (x))) + v = ( mpfr_zero_p (mpc_imagref (y)) + || (mpfr_inf_p (mpc_imagref (x)) && mpfr_zero_p (mpc_realref (y))) + || (mpfr_zero_p (mpc_imagref (x)) && mpfr_inf_p (mpc_realref (y))) + || ( (mpfr_inf_p (mpc_imagref (x)) || mpfr_inf_p (mpc_imagref (x))) && xrs*yis != xis*yrs) ? 0 : xrs * yis); } else { /* x = xr + i*(+/-inf) with |xr| != inf */ - u = ( mpfr_zero_p (MPC_IM (y)) - || (mpfr_zero_p (MPC_RE (x)) && mpfr_inf_p (MPC_RE (y))) - || (mpfr_inf_p (MPC_RE (y)) && xrs*yrs == xis*yis) + u = ( mpfr_zero_p (mpc_imagref (y)) + || (mpfr_zero_p (mpc_realref (x)) && mpfr_inf_p (mpc_realref (y))) + || (mpfr_inf_p (mpc_realref (y)) && xrs*yrs == xis*yis) ? 0 : -xis * yis); - v = ( mpfr_zero_p (MPC_RE (y)) - || (mpfr_zero_p (MPC_RE (x)) && mpfr_inf_p (MPC_IM (y))) - || (mpfr_inf_p (MPC_IM (y)) && xrs*yis != xis*yrs) + v = ( mpfr_zero_p (mpc_realref (y)) + || (mpfr_zero_p (mpc_realref (x)) && mpfr_inf_p (mpc_imagref (y))) + || (mpfr_inf_p (mpc_imagref (y)) && xrs*yis != xis*yrs) ? 0 : xis * yrs); } if (u == 0 && v == 0) { /* Naive result is NaN+i*NaN. Obtain an infinity using the algorithm given in Annex G.5.1 of the ISO C99 standard */ - int xr = (mpfr_zero_p (MPC_RE (x)) || mpfr_nan_p (MPC_RE (x)) ? 0 - : (mpfr_inf_p (MPC_RE (x)) ? 1 : 0)); - int xi = (mpfr_zero_p (MPC_IM (x)) || mpfr_nan_p (MPC_IM (x)) ? 0 - : (mpfr_inf_p (MPC_IM (x)) ? 1 : 0)); - int yr = (mpfr_zero_p (MPC_RE (y)) || mpfr_nan_p (MPC_RE (y)) ? 0 : 1); - int yi = (mpfr_zero_p (MPC_IM (y)) || mpfr_nan_p (MPC_IM (y)) ? 0 : 1); + int xr = (mpfr_zero_p (mpc_realref (x)) || mpfr_nan_p (mpc_realref (x)) ? 0 + : (mpfr_inf_p (mpc_realref (x)) ? 1 : 0)); + int xi = (mpfr_zero_p (mpc_imagref (x)) || mpfr_nan_p (mpc_imagref (x)) ? 0 + : (mpfr_inf_p (mpc_imagref (x)) ? 1 : 0)); + int yr = (mpfr_zero_p (mpc_realref (y)) || mpfr_nan_p (mpc_realref (y)) ? 0 : 1); + int yi = (mpfr_zero_p (mpc_imagref (y)) || mpfr_nan_p (mpc_imagref (y)) ? 0 : 1); if (mpc_inf_p (y)) { - yr = mpfr_inf_p (MPC_RE (y)) ? 1 : 0; - yi = mpfr_inf_p (MPC_IM (y)) ? 1 : 0; + yr = mpfr_inf_p (mpc_realref (y)) ? 1 : 0; + yi = mpfr_inf_p (mpc_imagref (y)) ? 1 : 0; } u = xrs * xr * yrs * yr - xis * xi * yis * yi; @@ -95,14 +95,14 @@ mul_infinite (mpc_ptr z, mpc_srcptr x, mpc_srcptr y) } if (u == 0) - mpfr_set_nan (MPC_RE (z)); + mpfr_set_nan (mpc_realref (z)); else - mpfr_set_inf (MPC_RE (z), u); + mpfr_set_inf (mpc_realref (z), u); if (v == 0) - mpfr_set_nan (MPC_IM (z)); + mpfr_set_nan (mpc_imagref (z)); else - mpfr_set_inf (MPC_IM (z), v); + mpfr_set_inf (mpc_imagref (z), v); return MPC_INEX (0, 0); /* exact */ } @@ -116,19 +116,19 @@ mul_real (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) int inex; /* save signs of operands */ - xrs = MPFR_SIGNBIT (MPC_RE (x)); - xis = MPFR_SIGNBIT (MPC_IM (x)); - yrs = MPFR_SIGNBIT (MPC_RE (y)); - yis = MPFR_SIGNBIT (MPC_IM (y)); + xrs = MPFR_SIGNBIT (mpc_realref (x)); + xis = MPFR_SIGNBIT (mpc_imagref (x)); + yrs = MPFR_SIGNBIT (mpc_realref (y)); + yis = MPFR_SIGNBIT (mpc_imagref (y)); - inex = mpc_mul_fr (z, x, MPC_RE (y), rnd); + inex = mpc_mul_fr (z, x, mpc_realref (y), rnd); /* Signs of zeroes may be wrong. Their correction does not change the inexact flag. */ - if (mpfr_zero_p (MPC_RE (z))) - mpfr_setsign (MPC_RE (z), MPC_RE (z), MPC_RND_RE(rnd) == GMP_RNDD + if (mpfr_zero_p (mpc_realref (z))) + mpfr_setsign (mpc_realref (z), mpc_realref (z), MPC_RND_RE(rnd) == GMP_RNDD || (xrs != yrs && xis == yis), GMP_RNDN); - if (mpfr_zero_p (MPC_IM (z))) - mpfr_setsign (MPC_IM (z), MPC_IM (z), MPC_RND_IM (rnd) == GMP_RNDD + if (mpfr_zero_p (mpc_imagref (z))) + mpfr_setsign (mpc_imagref (z), mpc_imagref (z), MPC_RND_IM (rnd) == GMP_RNDD || (xrs != yis && xis != yrs), GMP_RNDN); return inex; @@ -149,19 +149,19 @@ mul_imag (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) else rop [0] = z[0]; - sign = (MPFR_SIGNBIT (MPC_RE (y)) != MPFR_SIGNBIT (MPC_IM (x))) - && (MPFR_SIGNBIT (MPC_IM (y)) != MPFR_SIGNBIT (MPC_RE (x))); + sign = (MPFR_SIGNBIT (mpc_realref (y)) != MPFR_SIGNBIT (mpc_imagref (x))) + && (MPFR_SIGNBIT (mpc_imagref (y)) != MPFR_SIGNBIT (mpc_realref (x))); - inex_re = -mpfr_mul (MPC_RE (rop), MPC_IM (x), MPC_IM (y), + inex_re = -mpfr_mul (mpc_realref (rop), mpc_imagref (x), mpc_imagref (y), INV_RND (MPC_RND_RE (rnd))); - mpfr_neg (MPC_RE (rop), MPC_RE (rop), GMP_RNDN); /* exact */ - inex_im = mpfr_mul (MPC_IM (rop), MPC_RE (x), MPC_IM (y), + mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN); /* exact */ + inex_im = mpfr_mul (mpc_imagref (rop), mpc_realref (x), mpc_imagref (y), MPC_RND_IM (rnd)); mpc_set (z, rop, MPC_RNDNN); /* Sign of zeroes may be wrong (note that Re(z) cannot be zero) */ - if (mpfr_zero_p (MPC_IM (z))) - mpfr_setsign (MPC_IM (z), MPC_IM (z), MPC_RND_IM (rnd) == GMP_RNDD + if (mpfr_zero_p (mpc_imagref (z))) + mpfr_setsign (mpc_imagref (z), mpc_imagref (z), MPC_RND_IM (rnd) == GMP_RNDD || sign, GMP_RNDN); if (overlap) @@ -321,18 +321,18 @@ mpc_mul_naive (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) int overlap, inex; mpc_t rop; - MPC_ASSERT ( mpfr_regular_p (MPC_RE (x)) && mpfr_regular_p (MPC_IM (x)) - && mpfr_regular_p (MPC_RE (y)) && mpfr_regular_p (MPC_IM (y))); + MPC_ASSERT ( mpfr_regular_p (mpc_realref (x)) && mpfr_regular_p (mpc_imagref (x)) + && mpfr_regular_p (mpc_realref (y)) && mpfr_regular_p (mpc_imagref (y))); overlap = (z == x) || (z == y); if (overlap) mpc_init3 (rop, MPC_PREC_RE (z), MPC_PREC_IM (z)); else rop [0] = z [0]; - inex = MPC_INEX (mpfr_fmma (MPC_RE (rop), MPC_RE (x), MPC_RE (y), MPC_IM (x), - MPC_IM (y), -1, MPC_RND_RE (rnd)), - mpfr_fmma (MPC_IM (rop), MPC_RE (x), MPC_IM (y), MPC_IM (x), - MPC_RE (y), +1, MPC_RND_IM (rnd))); + inex = MPC_INEX (mpfr_fmma (mpc_realref (rop), mpc_realref (x), mpc_realref (y), mpc_imagref (x), + mpc_imagref (y), -1, MPC_RND_RE (rnd)), + mpfr_fmma (mpc_imagref (rop), mpc_realref (x), mpc_imagref (y), mpc_imagref (x), + mpc_realref (y), +1, MPC_RND_IM (rnd))); mpc_set (z, rop, MPC_RNDNN); if (overlap) @@ -370,10 +370,10 @@ mpc_mul_karatsuba (mpc_ptr rop, mpc_srcptr op1, mpc_srcptr op2, mpc_rnd_t rnd) else result [0] = rop [0]; - a = MPC_RE(op1); - b = MPC_IM(op1); - c = MPC_RE(op2); - d = MPC_IM(op2); + a = mpc_realref(op1); + b = mpc_imagref(op1); + c = mpc_realref(op2); + d = mpc_imagref(op2); /* (a + i*b) * (c + i*d) = [ac - bd] + i*[ad + bc] */ @@ -528,41 +528,41 @@ mpc_mul_karatsuba (mpc_ptr rop, mpc_srcptr op1, mpc_srcptr op2, mpc_rnd_t rnd) if (mul_i == 0) { - inex_re = mpfr_set (MPC_RE(result), u, MPC_RND_RE(rnd)); + inex_re = mpfr_set (mpc_realref(result), u, MPC_RND_RE(rnd)); if (inex_re == 0) { inex_re = inexact; /* u is rounded away from 0 */ - inex_im = mpfr_add (MPC_IM(result), v, w, MPC_RND_IM(rnd)); + inex_im = mpfr_add (mpc_imagref(result), v, w, MPC_RND_IM(rnd)); } else - inex_im = mpfr_add (MPC_IM(result), v, w, MPC_RND_IM(rnd)); + inex_im = mpfr_add (mpc_imagref(result), v, w, MPC_RND_IM(rnd)); } else if (mul_i == 1) /* (x+i*y)/i = y - i*x */ { - inex_im = mpfr_neg (MPC_IM(result), u, MPC_RND_IM(rnd)); + inex_im = mpfr_neg (mpc_imagref(result), u, MPC_RND_IM(rnd)); if (inex_im == 0) { inex_im = -inexact; /* u is rounded away from 0 */ - inex_re = mpfr_add (MPC_RE(result), v, w, MPC_RND_RE(rnd)); + inex_re = mpfr_add (mpc_realref(result), v, w, MPC_RND_RE(rnd)); } else - inex_re = mpfr_add (MPC_RE(result), v, w, MPC_RND_RE(rnd)); + inex_re = mpfr_add (mpc_realref(result), v, w, MPC_RND_RE(rnd)); } else /* mul_i = 2, z/i^2 = -z */ { - inex_re = mpfr_neg (MPC_RE(result), u, MPC_RND_RE(rnd)); + inex_re = mpfr_neg (mpc_realref(result), u, MPC_RND_RE(rnd)); if (inex_re == 0) { inex_re = -inexact; /* u is rounded away from 0 */ - inex_im = -mpfr_add (MPC_IM(result), v, w, + inex_im = -mpfr_add (mpc_imagref(result), v, w, INV_RND(MPC_RND_IM(rnd))); - mpfr_neg (MPC_IM(result), MPC_IM(result), MPC_RND_IM(rnd)); + mpfr_neg (mpc_imagref(result), mpc_imagref(result), MPC_RND_IM(rnd)); } else { - inex_im = -mpfr_add (MPC_IM(result), v, w, + inex_im = -mpfr_add (mpc_imagref(result), v, w, INV_RND(MPC_RND_IM(rnd))); - mpfr_neg (MPC_IM(result), MPC_IM(result), MPC_RND_IM(rnd)); + mpfr_neg (mpc_imagref(result), mpc_imagref(result), MPC_RND_IM(rnd)); } } @@ -596,23 +596,23 @@ mpc_mul (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) return mul_infinite (a, c, b); /* NaN contamination of both parts in result */ - if (mpfr_nan_p (MPC_RE (b)) || mpfr_nan_p (MPC_IM (b)) - || mpfr_nan_p (MPC_RE (c)) || mpfr_nan_p (MPC_IM (c))) { - mpfr_set_nan (MPC_RE (a)); - mpfr_set_nan (MPC_IM (a)); + if (mpfr_nan_p (mpc_realref (b)) || mpfr_nan_p (mpc_imagref (b)) + || mpfr_nan_p (mpc_realref (c)) || mpfr_nan_p (mpc_imagref (c))) { + mpfr_set_nan (mpc_realref (a)); + mpfr_set_nan (mpc_imagref (a)); return MPC_INEX (0, 0); } /* check for real multiplication */ - if (mpfr_zero_p (MPC_IM (b))) + if (mpfr_zero_p (mpc_imagref (b))) return mul_real (a, c, b, rnd); - if (mpfr_zero_p (MPC_IM (c))) + if (mpfr_zero_p (mpc_imagref (c))) return mul_real (a, b, c, rnd); /* check for purely imaginary multiplication */ - if (mpfr_zero_p (MPC_RE (b))) + if (mpfr_zero_p (mpc_realref (b))) return mul_imag (a, c, b, rnd); - if (mpfr_zero_p (MPC_RE (c))) + if (mpfr_zero_p (mpc_realref (c))) return mul_imag (a, b, c, rnd); /* Check if b==c and call mpc_sqr in this case, to make sure */ @@ -624,10 +624,10 @@ mpc_mul (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) /* If the real and imaginary part of one argument have a very different */ /* exponent, it is not reasonable to use Karatsuba multiplication. */ if ( SAFE_ABS (mpfr_exp_t, - mpfr_get_exp (MPC_RE (b)) - mpfr_get_exp (MPC_IM (b))) + mpfr_get_exp (mpc_realref (b)) - mpfr_get_exp (mpc_imagref (b))) > (mpfr_exp_t) MPC_MAX_PREC (b) / 2 || SAFE_ABS (mpfr_exp_t, - mpfr_get_exp (MPC_RE (c)) - mpfr_get_exp (MPC_IM (c))) + mpfr_get_exp (mpc_realref (c)) - mpfr_get_exp (mpc_imagref (c))) > (mpfr_exp_t) MPC_MAX_PREC (c) / 2) return mpc_mul_naive (a, b, c, rnd); else diff --git a/src/mul_2exp.c b/src/mul_2exp.c index b32001f..ff2efe2 100644 --- a/src/mul_2exp.c +++ b/src/mul_2exp.c @@ -1,6 +1,6 @@ /* mpc_mul_2exp -- Multiply a complex number by 2^e. -Copyright (C) 2002, 2009 INRIA +Copyright (C) 2002, 2009, 2011 INRIA This file is part of GNU MPC. @@ -25,8 +25,8 @@ mpc_mul_2exp (mpc_ptr a, mpc_srcptr b, unsigned long int c, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_mul_2exp (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); - inex_im = mpfr_mul_2exp (MPC_IM(a), MPC_IM(b), c, MPC_RND_IM(rnd)); + inex_re = mpfr_mul_2exp (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_mul_2exp (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } diff --git a/src/mul_fr.c b/src/mul_fr.c index 4b049b6..bd3574d 100644 --- a/src/mul_fr.c +++ b/src/mul_fr.c @@ -1,6 +1,6 @@ /* mpc_mul_fr -- Multiply a complex number by a floating-point number. -Copyright (C) 2002, 2008, 2009, 2010 INRIA +Copyright (C) 2002, 2008, 2009, 2010, 2011 INRIA This file is part of GNU MPC. @@ -26,17 +26,17 @@ mpc_mul_fr (mpc_ptr a, mpc_srcptr b, mpfr_srcptr c, mpc_rnd_t rnd) int inex_re, inex_im; mpfr_t real; - if (c == MPC_RE (a)) + if (c == mpc_realref (a)) /* We have to use a temporary variable. */ mpfr_init2 (real, MPC_PREC_RE (a)); else - real [0] = MPC_RE (a) [0]; + real [0] = mpc_realref (a) [0]; - inex_re = mpfr_mul (real, MPC_RE(b), c, MPC_RND_RE(rnd)); - inex_im = mpfr_mul (MPC_IM(a), MPC_IM(b), c, MPC_RND_IM(rnd)); - mpfr_set (MPC_RE (a), real, GMP_RNDN); /* exact */ + inex_re = mpfr_mul (real, mpc_realref(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_mul (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd)); + mpfr_set (mpc_realref (a), real, GMP_RNDN); /* exact */ - if (c == MPC_RE (a)) + if (c == mpc_realref (a)) mpfr_clear (real); return MPC_INEX(inex_re, inex_im); diff --git a/src/mul_i.c b/src/mul_i.c index 484f7ff..591b0c6 100644 --- a/src/mul_i.c +++ b/src/mul_i.c @@ -1,6 +1,6 @@ /* mpc_mul_i -- Multiply a complex number by plus or minus i. -Copyright (C) 2005, 2009, 2010 INRIA +Copyright (C) 2005, 2009, 2010, 2011 INRIA This file is part of GNU MPC. @@ -32,16 +32,16 @@ mpc_mul_i (mpc_ptr a, mpc_srcptr b, int sign, mpc_rnd_t rnd) && MPC_PREC_IM (b) == MPC_PREC_RE (a)) { if (a == b) - mpfr_swap (MPC_RE (a), MPC_IM (a)); + mpfr_swap (mpc_realref (a), mpc_imagref (a)); else { - mpfr_set (MPC_RE (a), MPC_IM (b), GMP_RNDN); - mpfr_set (MPC_IM (a), MPC_RE (b), GMP_RNDN); + mpfr_set (mpc_realref (a), mpc_imagref (b), GMP_RNDN); + mpfr_set (mpc_imagref (a), mpc_realref (b), GMP_RNDN); } if (sign >= 0) - MPFR_CHANGE_SIGN (MPC_RE (a)); + MPFR_CHANGE_SIGN (mpc_realref (a)); else - MPFR_CHANGE_SIGN (MPC_IM (a)); + MPFR_CHANGE_SIGN (mpc_imagref (a)); inex_re = 0; inex_im = 0; } @@ -52,27 +52,27 @@ mpc_mul_i (mpc_ptr a, mpc_srcptr b, int sign, mpc_rnd_t rnd) mpfr_init2 (tmp, MPC_PREC_RE (a)); if (sign >= 0) { - inex_re = mpfr_neg (tmp, MPC_IM (b), MPC_RND_RE (rnd)); - inex_im = mpfr_set (MPC_IM (a), MPC_RE (b), MPC_RND_IM (rnd)); + inex_re = mpfr_neg (tmp, mpc_imagref (b), MPC_RND_RE (rnd)); + inex_im = mpfr_set (mpc_imagref (a), mpc_realref (b), MPC_RND_IM (rnd)); } else { - inex_re = mpfr_set (tmp, MPC_IM (b), MPC_RND_RE (rnd)); - inex_im = mpfr_neg (MPC_IM (a), MPC_RE (b), MPC_RND_IM (rnd)); + inex_re = mpfr_set (tmp, mpc_imagref (b), MPC_RND_RE (rnd)); + inex_im = mpfr_neg (mpc_imagref (a), mpc_realref (b), MPC_RND_IM (rnd)); } - mpfr_clear (MPC_RE (a)); - MPC_RE (a)[0] = tmp [0]; + mpfr_clear (mpc_realref (a)); + mpc_realref (a)[0] = tmp [0]; } else if (sign >= 0) { - inex_re = mpfr_neg (MPC_RE (a), MPC_IM (b), MPC_RND_RE (rnd)); - inex_im = mpfr_set (MPC_IM (a), MPC_RE (b), MPC_RND_IM (rnd)); + inex_re = mpfr_neg (mpc_realref (a), mpc_imagref (b), MPC_RND_RE (rnd)); + inex_im = mpfr_set (mpc_imagref (a), mpc_realref (b), MPC_RND_IM (rnd)); } else { - inex_re = mpfr_set (MPC_RE (a), MPC_IM (b), MPC_RND_RE (rnd)); - inex_im = mpfr_neg (MPC_IM (a), MPC_RE (b), MPC_RND_IM (rnd)); + inex_re = mpfr_set (mpc_realref (a), mpc_imagref (b), MPC_RND_RE (rnd)); + inex_im = mpfr_neg (mpc_imagref (a), mpc_realref (b), MPC_RND_IM (rnd)); } } diff --git a/src/mul_si.c b/src/mul_si.c index 7cf3bc7..f539d8b 100644 --- a/src/mul_si.c +++ b/src/mul_si.c @@ -1,6 +1,6 @@ /* mpc_mul_si -- Multiply a complex number by a signed integer. -Copyright (C) 2005, 2009 INRIA +Copyright (C) 2005, 2009, 2011 INRIA This file is part of GNU MPC. @@ -25,8 +25,8 @@ mpc_mul_si (mpc_ptr a, mpc_srcptr b, long int c, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_mul_si (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); - inex_im = mpfr_mul_si (MPC_IM(a), MPC_IM(b), c, MPC_RND_IM(rnd)); + inex_re = mpfr_mul_si (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_mul_si (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } diff --git a/src/mul_ui.c b/src/mul_ui.c index 618c36b..922e4b3 100644 --- a/src/mul_ui.c +++ b/src/mul_ui.c @@ -1,6 +1,6 @@ /* mpc_mul_ui -- Multiply a complex number by a nonnegative integer. -Copyright (C) 2002, 2009 INRIA +Copyright (C) 2002, 2009, 2011 INRIA This file is part of GNU MPC. @@ -25,8 +25,8 @@ mpc_mul_ui (mpc_ptr a, mpc_srcptr b, unsigned long int c, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_mul_ui (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); - inex_im = mpfr_mul_ui (MPC_IM(a), MPC_IM(b), c, MPC_RND_IM(rnd)); + inex_re = mpfr_mul_ui (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_mul_ui (mpc_imagref(a), mpc_imagref(b), c, MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } @@ -1,6 +1,6 @@ /* mpc_neg -- Negate a complex number. -Copyright (C) 2002, 2009 INRIA +Copyright (C) 2002, 2009, 2011 INRIA This file is part of GNU MPC. @@ -25,8 +25,8 @@ mpc_neg (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_neg (MPC_RE(a), MPC_RE(b), MPC_RND_RE(rnd)); - inex_im = mpfr_neg (MPC_IM(a), MPC_IM(b), MPC_RND_IM(rnd)); + inex_re = mpfr_neg (mpc_realref(a), mpc_realref(b), MPC_RND_RE(rnd)); + inex_im = mpfr_neg (mpc_imagref(a), mpc_imagref(b), MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } @@ -32,14 +32,14 @@ mpc_norm (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd) norm = abs^2; so norm (+-inf, xxx) = norm (xxx, +-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))) + else if (mpfr_zero_p (mpc_realref (b))) { + if (mpfr_zero_p (mpc_imagref (b))) return mpfr_set_ui (a, 0, rnd); /* +0 */ else - return mpfr_sqr (a, MPC_IM (b), rnd); + return mpfr_sqr (a, mpc_imagref (b), rnd); } - else if (mpfr_zero_p (MPC_IM (b))) - return mpfr_sqr (a, MPC_RE (b), rnd); /* Re(b) <> 0 */ + else if (mpfr_zero_p (mpc_imagref (b))) + return mpfr_sqr (a, mpc_realref (b), rnd); /* Re(b) <> 0 */ else /* everything finite and non-zero */ { mpfr_t u, v, res; @@ -72,8 +72,8 @@ mpc_norm (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd) mpfr_set_prec (u, prec_u); mpfr_set_prec (v, prec_v); - 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 */ + inexact = mpfr_sqr (u, mpc_realref(b), GMP_RNDD); /* err <= 1 ulp in prec */ + inexact |= mpfr_sqr (v, mpc_imagref(b), GMP_RNDD); /* err <= 1 ulp in prec */ /* If loops = max_loops, inexact should be 0 here, except in case of underflow or overflow. @@ -133,20 +133,20 @@ mpc_norm (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd) int inex_underflow; /* scale the input to an average exponent close to 0 */ - exp_re = (unsigned long int) (-mpfr_get_exp (MPC_RE (b))); - exp_im = (unsigned long int) (-mpfr_get_exp (MPC_IM (b))); + exp_re = (unsigned long int) (-mpfr_get_exp (mpc_realref (b))); + exp_im = (unsigned long int) (-mpfr_get_exp (mpc_imagref (b))); scale = exp_re / 2 + exp_im / 2 + (exp_re % 2 + exp_im % 2) / 2; /* (exp_re + exp_im) / 2, computed in a way avoiding integer overflow */ if (mpfr_zero_p (u)) { /* recompute the scaled value exactly */ - mpfr_mul_2ui (u, MPC_RE (b), scale, GMP_RNDN); + mpfr_mul_2ui (u, mpc_realref (b), scale, GMP_RNDN); mpfr_sqr (u, u, GMP_RNDN); } else /* just scale */ mpfr_mul_2ui (u, u, 2*scale, GMP_RNDN); if (mpfr_zero_p (v)) { - mpfr_mul_2ui (v, MPC_IM (b), scale, GMP_RNDN); + mpfr_mul_2ui (v, mpc_imagref (b), scale, GMP_RNDN); mpfr_sqr (v, v, GMP_RNDN); } else diff --git a/src/out_str.c b/src/out_str.c index 1ac9687..87cc823 100644 --- a/src/out_str.c +++ b/src/out_str.c @@ -1,6 +1,6 @@ /* mpc_out_str -- Output a complex number on a given stream. -Copyright (C) 2009 INRIA +Copyright (C) 2009, 2011 INRIA This file is part of GNU MPC. @@ -30,9 +30,9 @@ mpc_out_str (FILE *stream, int base, size_t n, mpc_srcptr op, mpc_rnd_t rnd) { stream = stdout; /* fprintf does not allow NULL as first argument */ fprintf (stream, "("); - size += mpfr_out_str (stream, base, n, MPC_RE(op), MPC_RND_RE(rnd)); + size += mpfr_out_str (stream, base, n, mpc_realref(op), MPC_RND_RE(rnd)); fprintf (stream, " "); - size += mpfr_out_str (stream, base, n, MPC_IM(op), MPC_RND_RE(rnd)); + size += mpfr_out_str (stream, base, n, mpc_imagref(op), MPC_RND_RE(rnd)); fprintf (stream, ")"); return size; @@ -129,7 +129,7 @@ fix_sign (mpc_ptr z, int sign_eps, int sign_a, mpfr_srcptr y) else /* y is not an integer */ goto end; - if (mpfr_zero_p (MPC_RE(z))) + if (mpfr_zero_p (mpc_realref(z))) { /* we assume y is always integer in that case (FIXME: prove it): (eps+I*a)^y = +0 + I*a^y for y = 1 mod 4 and sign_eps = 0 @@ -137,9 +137,9 @@ fix_sign (mpc_ptr z, int sign_eps, int sign_a, mpfr_srcptr y) MPC_ASSERT (ymod4 == 1 || ymod4 == 3); if ((ymod4 == 3 && sign_eps == 0) || (ymod4 == 1 && sign_eps == 1)) - mpfr_neg (MPC_RE(z), MPC_RE(z), GMP_RNDZ); + mpfr_neg (mpc_realref(z), mpc_realref(z), GMP_RNDZ); } - else if (mpfr_zero_p (MPC_IM(z))) + else if (mpfr_zero_p (mpc_imagref(z))) { /* we assume y is always integer in that case (FIXME: prove it): (eps+I*a)^y = a^y - 0*I for y = 0 mod 4 and sign_a = sign_eps @@ -147,7 +147,7 @@ fix_sign (mpc_ptr z, int sign_eps, int sign_a, mpfr_srcptr y) MPC_ASSERT (ymod4 == 0 || ymod4 == 2); if ((ymod4 == 0 && sign_a == sign_eps) || (ymod4 == 2 && sign_a != sign_eps)) - mpfr_neg (MPC_IM(z), MPC_IM(z), GMP_RNDZ); + mpfr_neg (mpc_imagref(z), mpc_imagref(z), GMP_RNDZ); } end: @@ -177,9 +177,9 @@ mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd, mpz_t my, a, b, c, d, u; unsigned long int t; int ret = -2; - int sign_rex = mpfr_signbit (MPC_RE(x)); - int sign_imx = mpfr_signbit (MPC_IM(x)); - int x_imag = mpfr_zero_p (MPC_RE(x)); + int sign_rex = mpfr_signbit (mpc_realref(x)); + int sign_imx = mpfr_signbit (mpc_imagref(x)); + int x_imag = mpfr_zero_p (mpc_realref(x)); int z_is_y = 0; mpfr_t copy_of_y; @@ -210,15 +210,15 @@ mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd, ec = 0; } else - ec = mpfr_get_z_exp (c, MPC_RE(x)); - if (mpfr_zero_p (MPC_IM(x))) + ec = mpfr_get_z_exp (c, mpc_realref(x)); + if (mpfr_zero_p (mpc_imagref(x))) { mpz_set_ui (d, 0); ed = ec; } else { - ed = mpfr_get_z_exp (d, MPC_IM(x)); + ed = mpfr_get_z_exp (d, mpc_imagref(x)); if (x_imag) ec = ed; } @@ -403,10 +403,10 @@ mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd, goto end; } - ret = mpfr_set_z (MPC_RE(z), a, MPC_RND_RE(rnd)); - ret = MPC_INEX(ret, mpfr_set_z (MPC_IM(z), b, MPC_RND_IM(rnd))); - mpfr_mul_2si (MPC_RE(z), MPC_RE(z), ed, MPC_RND_RE(rnd)); - mpfr_mul_2si (MPC_IM(z), MPC_IM(z), ed, MPC_RND_IM(rnd)); + ret = mpfr_set_z (mpc_realref(z), a, MPC_RND_RE(rnd)); + ret = MPC_INEX(ret, mpfr_set_z (mpc_imagref(z), b, MPC_RND_IM(rnd))); + mpfr_mul_2si (mpc_realref(z), mpc_realref(z), ed, MPC_RND_RE(rnd)); + mpfr_mul_2si (mpc_imagref(z), mpc_imagref(z), ed, MPC_RND_IM(rnd)); end: mpz_clear (my); @@ -484,15 +484,15 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) mpc_t t, u; mpfr_prec_t p, pr, pi, maxprec; - x_real = mpfr_zero_p (MPC_IM(x)); - y_real = mpfr_zero_p (MPC_IM(y)); + x_real = mpfr_zero_p (mpc_imagref(x)); + y_real = mpfr_zero_p (mpc_imagref(y)); - if (y_real && mpfr_zero_p (MPC_RE(y))) /* case y zero */ + if (y_real && mpfr_zero_p (mpc_realref(y))) /* case y zero */ { - if (x_real && mpfr_zero_p (MPC_RE(x))) /* 0^0 = NaN +i*NaN */ + if (x_real && mpfr_zero_p (mpc_realref(x))) /* 0^0 = NaN +i*NaN */ { - mpfr_set_nan (MPC_RE(z)); - mpfr_set_nan (MPC_IM(z)); + mpfr_set_nan (mpc_realref(z)); + mpfr_set_nan (mpc_imagref(z)); return 0; } else /* x^0 = 1 +/- i*0 even for x=NaN see algorithms.tex for the @@ -511,10 +511,10 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) if (cx1 == 0 && inex != 0) cx1 = -inex; - sign_zi = (cx1 < 0 && mpfr_signbit (MPC_IM (y)) == 0) + sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0) || (cx1 == 0 - && mpfr_signbit (MPC_IM (x)) != mpfr_signbit (MPC_RE (y))) - || (cx1 > 0 && mpfr_signbit (MPC_IM (y))); + && mpfr_signbit (mpc_imagref (x)) != mpfr_signbit (mpc_realref (y))) + || (cx1 > 0 && mpfr_signbit (mpc_imagref (y))); /* warning: mpc_set_ui_ui does not set Im(z) to -0 if Im(rnd)=RNDD */ ret = mpc_set_ui_ui (z, 1, 0, rnd); @@ -540,7 +540,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) if (x_real) /* case x real */ { - if (mpfr_zero_p (MPC_RE(x))) /* x is zero */ + if (mpfr_zero_p (mpc_realref(x))) /* x is zero */ { /* special values: exp(y*log(x)) */ mpc_init2 (u, 2); @@ -552,11 +552,11 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) } /* Special case 1^y = 1 */ - if (mpfr_cmp_ui (MPC_RE(x), 1) == 0) + if (mpfr_cmp_ui (mpc_realref(x), 1) == 0) { int s1, s2; - s1 = mpfr_signbit (MPC_RE (y)); - s2 = mpfr_signbit (MPC_IM (x)); + s1 = mpfr_signbit (mpc_realref (y)); + s2 = mpfr_signbit (mpc_imagref (x)); ret = mpc_set_ui (z, +1, rnd); /* the sign of the zero imaginary part is known in some cases (see @@ -576,15 +576,15 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) /* x^y is real when: (a) x is real and y is integer (b) x is real non-negative and y is real */ - if (y_real && (mpfr_integer_p (MPC_RE(y)) || - mpfr_cmp_ui (MPC_RE(x), 0) >= 0)) + if (y_real && (mpfr_integer_p (mpc_realref(y)) || + mpfr_cmp_ui (mpc_realref(x), 0) >= 0)) { int s1, s2; - s1 = mpfr_signbit (MPC_RE (y)); - s2 = mpfr_signbit (MPC_IM (x)); + s1 = mpfr_signbit (mpc_realref (y)); + s2 = mpfr_signbit (mpc_imagref (x)); - ret = mpfr_pow (MPC_RE(z), MPC_RE(x), MPC_RE(y), MPC_RND_RE(rnd)); - ret = MPC_INEX(ret, mpfr_set_ui (MPC_IM(z), 0, MPC_RND_IM(rnd))); + ret = mpfr_pow (mpc_realref(z), mpc_realref(x), mpc_realref(y), MPC_RND_RE(rnd)); + ret = MPC_INEX(ret, mpfr_set_ui (mpc_imagref(z), 0, MPC_RND_IM(rnd))); /* the sign of the zero imaginary part is known in some cases (see algorithm.tex). In such cases we have (x +s*0i)^(y+/-0i) @@ -595,20 +595,20 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) because mpfr_set_ui(z_i, 0, rnd) always sets z_i to +0. */ if (MPC_RND_IM(rnd) == GMP_RNDD || s1 != s2) - mpfr_neg (MPC_IM(z), MPC_IM(z), MPC_RND_IM(rnd)); + mpfr_neg (mpc_imagref(z), mpc_imagref(z), MPC_RND_IM(rnd)); goto end; } /* (-1)^(n+I*t) is real for n integer and t real */ - if (mpfr_cmp_si (MPC_RE(x), -1) == 0 && mpfr_integer_p (MPC_RE(y))) + if (mpfr_cmp_si (mpc_realref(x), -1) == 0 && mpfr_integer_p (mpc_realref(y))) z_real = 1; /* for x real, x^y is imaginary when: (a) x is negative and y is half-an-integer (b) x = -1 and Re(y) is half-an-integer */ - if (mpfr_cmp_ui (MPC_RE(x), 0) < 0 && is_odd (MPC_RE(y), 1) && - (y_real || mpfr_cmp_si (MPC_RE(x), -1) == 0)) + if (mpfr_cmp_ui (mpc_realref(x), 0) < 0 && is_odd (mpc_realref(y), 1) && + (y_real || mpfr_cmp_si (mpc_realref(x), -1) == 0)) z_imag = 1; } else /* x non real */ @@ -617,26 +617,26 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) I^(n+t*I) and (-I)^(n+t*I) are imaginary for n odd and t real (s*I)^n is real for n even and imaginary for n odd */ if ((mpc_cmp_si_si (x, 0, 1) == 0 || mpc_cmp_si_si (x, 0, -1) == 0 || - (mpfr_cmp_ui (MPC_RE(x), 0) == 0 && y_real)) && - mpfr_integer_p (MPC_RE(y))) + (mpfr_cmp_ui (mpc_realref(x), 0) == 0 && y_real)) && + mpfr_integer_p (mpc_realref(y))) { /* x is I or -I, and Re(y) is an integer */ - if (is_odd (MPC_RE(y), 0)) + if (is_odd (mpc_realref(y), 0)) z_imag = 1; /* Re(y) odd: z is imaginary */ else z_real = 1; /* Re(y) even: z is real */ } else /* (t+/-t*I)^(2n) is imaginary for n odd and real for n even */ - if (mpfr_cmpabs (MPC_RE(x), MPC_IM(x)) == 0 && y_real && - mpfr_integer_p (MPC_RE(y)) && is_odd (MPC_RE(y), 0) == 0) + if (mpfr_cmpabs (mpc_realref(x), mpc_imagref(x)) == 0 && y_real && + mpfr_integer_p (mpc_realref(y)) && is_odd (mpc_realref(y), 0) == 0) { - if (is_odd (MPC_RE(y), -1)) /* y/2 is odd */ + if (is_odd (mpc_realref(y), -1)) /* y/2 is odd */ z_imag = 1; else z_real = 1; } - pr = mpfr_get_prec (MPC_RE(z)); - pi = mpfr_get_prec (MPC_IM(z)); + pr = mpfr_get_prec (mpc_realref(z)); + pi = mpfr_get_prec (mpc_imagref(z)); p = (pr > pi) ? pr : pi; p += 12; /* experimentally, seems to give less than 10% of failures in Ziv's strategy; probably wrong now since q is not computed */ @@ -647,7 +647,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) pr += MPC_RND_RE(rnd) == GMP_RNDN; pi += MPC_RND_IM(rnd) == GMP_RNDN; maxprec = MPC_MAX_PREC (z); - x_imag = mpfr_zero_p (MPC_RE(x)); + x_imag = mpfr_zero_p (mpc_realref(x)); for (loop = 0;; loop++) { int ret_exp; @@ -660,9 +660,9 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) if (loop == 0) { /* compute q such that |Re (y log x)|, |Im (y log x)| < 2^q */ - q = mpfr_get_exp (MPC_RE(t)) > 0 ? mpfr_get_exp (MPC_RE(t)) : 0; - if (mpfr_get_exp (MPC_IM(t)) > (mpfr_exp_t) q) - q = mpfr_get_exp (MPC_IM(t)); + q = mpfr_get_exp (mpc_realref(t)) > 0 ? mpfr_get_exp (mpc_realref(t)) : 0; + if (mpfr_get_exp (mpc_imagref(t)) > (mpfr_exp_t) q) + q = mpfr_get_exp (mpc_imagref(t)); } mpfr_clear_overflow (); @@ -679,9 +679,9 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) exponent difference between the real and imaginary parts. We assume either the real or the imaginary part of u is not zero. */ - dr = mpfr_zero_p (MPC_RE(u)) ? mpfr_get_exp (MPC_IM(u)) - : mpfr_get_exp (MPC_RE(u)); - di = mpfr_zero_p (MPC_IM(u)) ? dr : mpfr_get_exp (MPC_IM(u)); + dr = mpfr_zero_p (mpc_realref(u)) ? mpfr_get_exp (mpc_imagref(u)) + : mpfr_get_exp (mpc_realref(u)); + di = mpfr_zero_p (mpc_imagref(u)) ? dr : mpfr_get_exp (mpc_imagref(u)); if (dr > di) { di = dr - di; @@ -696,15 +696,15 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) (see algorithms.tex) plus one due to the exponent difference: if z = a + I*b, where the relative error on z is at most 2^(-p), and EXP(a) = EXP(b) + k, the relative error on b is at most 2^(k-p) */ - if ((z_imag || (p > q + 3 + dr && mpfr_can_round (MPC_RE(u), p - q - 3 - dr, GMP_RNDN, GMP_RNDZ, pr))) && - (z_real || (p > q + 3 + di && mpfr_can_round (MPC_IM(u), p - q - 3 - di, GMP_RNDN, GMP_RNDZ, pi)))) + if ((z_imag || (p > q + 3 + dr && mpfr_can_round (mpc_realref(u), p - q - 3 - dr, GMP_RNDN, GMP_RNDZ, pr))) && + (z_real || (p > q + 3 + di && mpfr_can_round (mpc_imagref(u), p - q - 3 - di, GMP_RNDN, GMP_RNDZ, pi)))) break; /* if Re(u) is not known to be zero, assume it is a normal number, i.e., neither zero, Inf or NaN, otherwise we might enter an infinite loop */ - MPC_ASSERT (z_imag || mpfr_number_p (MPC_RE(u))); + MPC_ASSERT (z_imag || mpfr_number_p (mpc_realref(u))); /* idem for Im(u) */ - MPC_ASSERT (z_real || mpfr_number_p (MPC_IM(u))); + MPC_ASSERT (z_real || mpfr_number_p (mpc_imagref(u))); if (ret == -2) /* we did not yet call mpc_pow_exact, or it aborted because intermediate computations had > maxprec bits */ @@ -713,7 +713,7 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) if (y_real) { maxprec *= 2; - ret = mpc_pow_exact (z, x, MPC_RE(y), rnd, maxprec); + ret = mpc_pow_exact (z, x, mpc_realref(y), rnd, maxprec); if (ret != -1 && ret != -2) goto exact; } @@ -741,34 +741,34 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) cx1 > 0 if |x| > 1 */ - sign_rex = mpfr_signbit (MPC_RE (x)); - sign_imx = mpfr_signbit (MPC_IM (x)); + sign_rex = mpfr_signbit (mpc_realref (x)); + sign_imx = mpfr_signbit (mpc_imagref (x)); mpfr_init (n); inex = mpc_norm (n, x, GMP_RNDN); cx1 = mpfr_cmp_ui (n, 1); if (cx1 == 0 && inex != 0) cx1 = -inex; - sign_zi = (cx1 < 0 && mpfr_signbit (MPC_IM (y)) == 0) - || (cx1 == 0 && sign_imx != mpfr_signbit (MPC_RE (y))) - || (cx1 > 0 && mpfr_signbit (MPC_IM (y))); + sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0) + || (cx1 == 0 && sign_imx != mpfr_signbit (mpc_realref (y))) + || (cx1 > 0 && mpfr_signbit (mpc_imagref (y))); /* copy RE(y) to n since if z==y we will destroy Re(y) below */ - mpfr_set_prec (n, mpfr_get_prec (MPC_RE (y))); - mpfr_set (n, MPC_RE (y), GMP_RNDN); - ret = mpfr_set (MPC_RE(z), MPC_RE(u), MPC_RND_RE(rnd)); + mpfr_set_prec (n, mpfr_get_prec (mpc_realref (y))); + mpfr_set (n, mpc_realref (y), GMP_RNDN); + ret = mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd)); if (y_real && (x_real || x_imag)) { /* FIXME: with y_real we assume Im(y) is really 0, which is the case for example when y comes from pow_fr, but in case Im(y) is +0 or -0, we might get different results */ - mpfr_set_ui (MPC_IM (z), 0, MPC_RND_IM (rnd)); + mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd)); fix_sign (z, sign_rex, sign_imx, n); ret = MPC_INEX(ret, 0); /* imaginary part is exact */ } else { - ret = MPC_INEX (ret, mpfr_set_ui (MPC_IM (z), 0, MPC_RND_IM (rnd))); + ret = MPC_INEX (ret, mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd))); /* warning: mpfr_set_ui does not set Im(z) to -0 if Im(rnd) = RNDD */ if (MPC_RND_IM (rnd) == GMP_RNDD || sign_zi) mpc_conj (z, z, MPC_RNDNN); @@ -778,23 +778,23 @@ mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) } else if (z_imag) { - ret = mpfr_set (MPC_IM(z), MPC_IM(u), MPC_RND_IM(rnd)); + ret = mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd)); /* if z is imaginary and y real, then x cannot be real */ if (y_real && x_imag) { - int sign_rex = mpfr_signbit (MPC_RE (x)); + int sign_rex = mpfr_signbit (mpc_realref (x)); /* If z overlaps with y we set Re(z) before checking Re(y) below, but in that case y=0, which was dealt with above. */ - mpfr_set_ui (MPC_RE (z), 0, MPC_RND_RE (rnd)); + mpfr_set_ui (mpc_realref (z), 0, MPC_RND_RE (rnd)); /* Note: fix_sign only does something when y is an integer, then necessarily y = 1 or 3 (mod 4), and in that case the sign of Im(x) is irrelevant. */ - fix_sign (z, sign_rex, 0, MPC_RE (y)); + fix_sign (z, sign_rex, 0, mpc_realref (y)); ret = MPC_INEX(0, ret); } else - ret = MPC_INEX(mpfr_set_ui (MPC_RE(z), 0, MPC_RND_RE(rnd)), ret); + ret = MPC_INEX(mpfr_set_ui (mpc_realref(z), 0, MPC_RND_RE(rnd)), ret); } else ret = mpc_set (z, u, rnd); diff --git a/src/pow_fr.c b/src/pow_fr.c index 96ccbd3..8c5d930 100644 --- a/src/pow_fr.c +++ b/src/pow_fr.c @@ -1,6 +1,6 @@ /* mpc_pow_fr -- Raise a complex number to a floating-point power. -Copyright (C) 2009 INRIA +Copyright (C) 2009, 2011 INRIA This file is part of GNU MPC. @@ -27,11 +27,11 @@ mpc_pow_fr (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd) int inex; /* avoid copying the significand of y by copying only the struct */ - MPC_RE(yy)[0] = y[0]; - mpfr_init2 (MPC_IM(yy), MPFR_PREC_MIN); - mpfr_set_ui (MPC_IM(yy), 0, GMP_RNDN); + mpc_realref(yy)[0] = y[0]; + mpfr_init2 (mpc_imagref(yy), MPFR_PREC_MIN); + mpfr_set_ui (mpc_imagref(yy), 0, GMP_RNDN); inex = mpc_pow (z, x, yy, rnd); - mpfr_clear (MPC_IM(yy)); + mpfr_clear (mpc_imagref(yy)); return inex; } diff --git a/src/pow_ui.c b/src/pow_ui.c index 0ce2e3f..3381028 100644 --- a/src/pow_ui.c +++ b/src/pow_ui.c @@ -1,6 +1,6 @@ /* mpc_pow_ui -- Raise a complex number to an integer power. -Copyright (C) 2009, 2010 INRIA +Copyright (C) 2009, 2010, 2011 INRIA This file is part of GNU MPC. @@ -53,7 +53,7 @@ mpc_pow_usi (mpc_ptr z, mpc_srcptr x, unsigned long y, int sign, int loop, done; /* let mpc_pow deal with special values */ - if (!mpc_fin_p (x) || mpfr_zero_p (MPC_RE (x)) || mpfr_zero_p (MPC_IM(x)) + if (!mpc_fin_p (x) || mpfr_zero_p (mpc_realref (x)) || mpfr_zero_p (mpc_imagref(x)) || y == 0) return mpc_pow_usi_naive (z, x, y, sign, rnd); /* easy special cases */ @@ -67,8 +67,8 @@ mpc_pow_usi (mpc_ptr z, mpc_srcptr x, unsigned long y, int sign, return mpc_sqr (z, x, rnd); /* let mpc_pow treat potential over- and underflows */ else { - mpfr_exp_t exp_r = mpfr_get_exp (MPC_RE (x)), - exp_i = mpfr_get_exp (MPC_IM (x)); + mpfr_exp_t exp_r = mpfr_get_exp (mpc_realref (x)), + exp_i = mpfr_get_exp (mpc_imagref (x)); if ( MPC_MAX (exp_r, exp_i) > mpfr_get_emax () / (mpfr_exp_t) y /* heuristic for overflow */ || MPC_MAX (-exp_r, -exp_i) > (-mpfr_get_emin ()) / (mpfr_exp_t) y @@ -113,7 +113,7 @@ mpc_pow_usi (mpc_ptr z, mpc_srcptr x, unsigned long y, int sign, if (sign < 0) mpc_ui_div (t, 1ul, t, MPC_RNDNN); - if (mpfr_zero_p (MPC_RE(t)) || mpfr_zero_p (MPC_IM(t))) { + if (mpfr_zero_p (mpc_realref(t)) || mpfr_zero_p (mpc_imagref(t))) { inex = mpc_pow_usi_naive (z, x, y, sign, rnd); /* since mpfr_get_exp() is not defined for zero */ done = 1; @@ -124,16 +124,16 @@ mpc_pow_usi (mpc_ptr z, mpc_srcptr x, unsigned long y, int sign, mpfr_exp_t diff; mpfr_prec_t er, ei; - diff = mpfr_get_exp (MPC_RE(t)) - mpfr_get_exp (MPC_IM(t)); + diff = mpfr_get_exp (mpc_realref(t)) - mpfr_get_exp (mpc_imagref(t)); /* the factor on the real part is 2+2^(-diff+2) <= 4 for diff >= 1 and < 2^(-diff+3) for diff <= 0 */ er = (diff >= 1) ? l0 + 3 : l0 + (-diff) + 3; /* the factor on the imaginary part is 2+2^(diff+2) <= 4 for diff <= -1 and < 2^(diff+3) for diff >= 0 */ ei = (diff <= -1) ? l0 + 3 : l0 + diff + 3; - if (mpfr_can_round (MPC_RE(t), p - er, GMP_RNDZ, GMP_RNDZ, + if (mpfr_can_round (mpc_realref(t), p - er, GMP_RNDZ, GMP_RNDZ, MPC_PREC_RE(z) + (MPC_RND_RE(rnd) == GMP_RNDN)) - && mpfr_can_round (MPC_IM(t), p - ei, GMP_RNDZ, GMP_RNDZ, + && mpfr_can_round (mpc_imagref(t), p - ei, GMP_RNDZ, GMP_RNDZ, MPC_PREC_IM(z) + (MPC_RND_IM(rnd) == GMP_RNDN))) { inex = mpc_set (z, t, rnd); done = 1; @@ -25,8 +25,8 @@ mpc_proj (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) { if (mpc_inf_p (b)) { /* infinities project to +Inf +i* copysign(0.0, cimag(z)) */ - mpfr_set_inf (MPC_RE (a), +1); - mpfr_set_zero (MPC_IM (a), (mpfr_signbit (MPC_IM (b)) ? -1 : 1)); + mpfr_set_inf (mpc_realref (a), +1); + mpfr_set_zero (mpc_imagref (a), (mpfr_signbit (mpc_imagref (b)) ? -1 : 1)); return MPC_INEX (0, 0); } else @@ -1,6 +1,6 @@ /* mpc_real -- Get the real part of a complex number. -Copyright (C) 2008, 2009 INRIA +Copyright (C) 2008, 2009, 2011 INRIA This file is part of GNU MPC. @@ -23,5 +23,5 @@ along with this program. If not, see http://www.gnu.org/licenses/ . int mpc_real (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd) { - return mpfr_set (a, MPC_RE (b), rnd); + return mpfr_set (a, mpc_realref (b), rnd); } @@ -1,6 +1,6 @@ /* mpc_set -- Set a complex number from another complex number. -Copyright (C) 2002, 2009 INRIA +Copyright (C) 2002, 2009, 2011 INRIA This file is part of GNU MPC. @@ -25,8 +25,8 @@ mpc_set (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_set (MPC_RE(a), MPC_RE(b), MPC_RND_RE(rnd)); - inex_im = mpfr_set (MPC_IM(a), MPC_IM(b), MPC_RND_IM(rnd)); + inex_re = mpfr_set (mpc_realref(a), mpc_realref(b), MPC_RND_RE(rnd)); + inex_im = mpfr_set (mpc_imagref(a), mpc_imagref(b), MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } diff --git a/src/set_prec.c b/src/set_prec.c index fc3a5cb..c5e6f24 100644 --- a/src/set_prec.c +++ b/src/set_prec.c @@ -1,6 +1,6 @@ /* mpc_set_prec -- reset the precision of a complex variable. -Copyright (C) 2002, 2009 INRIA +Copyright (C) 2002, 2009, 2011 INRIA This file is part of GNU MPC. @@ -23,6 +23,6 @@ along with this program. If not, see http://www.gnu.org/licenses/ . void mpc_set_prec (mpc_t x, mpfr_prec_t prec) { - mpfr_set_prec (MPC_RE(x), prec); - mpfr_set_prec (MPC_IM(x), prec); + mpfr_set_prec (mpc_realref(x), prec); + mpfr_set_prec (mpc_imagref(x), prec); } diff --git a/src/set_str.c b/src/set_str.c index 27a9d89..195b9ac 100644 --- a/src/set_str.c +++ b/src/set_str.c @@ -1,6 +1,6 @@ /* mpc_set_str -- Convert a string into a complex number. -Copyright (C) 2009, 2010 INRIA +Copyright (C) 2009, 2010, 2011 INRIA This file is part of GNU MPC. @@ -36,7 +36,7 @@ mpc_set_str (mpc_t z, const char *str, int base, mpc_rnd_t rnd) return inex; } - mpfr_set_nan (MPC_RE (z)); - mpfr_set_nan (MPC_IM (z)); + mpfr_set_nan (mpc_realref (z)); + mpfr_set_nan (mpc_imagref (z)); return -1; } diff --git a/src/sin_cos.c b/src/sin_cos.c index 99a0772..0cff45a 100644 --- a/src/sin_cos.c +++ b/src/sin_cos.c @@ -37,38 +37,38 @@ mpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, op_loc [0] = op [0]; if (rop_sin != NULL) { - if (mpfr_nan_p (MPC_RE (op_loc)) || mpfr_nan_p (MPC_IM (op_loc))) { + if (mpfr_nan_p (mpc_realref (op_loc)) || mpfr_nan_p (mpc_imagref (op_loc))) { mpc_set (rop_sin, op_loc, rnd_sin); - if (mpfr_nan_p (MPC_IM (op_loc))) { + if (mpfr_nan_p (mpc_imagref (op_loc))) { /* sin(x +i*NaN) = NaN +i*NaN, except for x=0 */ /* sin(-0 +i*NaN) = -0 +i*NaN */ /* sin(+0 +i*NaN) = +0 +i*NaN */ - if (!mpfr_zero_p (MPC_RE (op_loc))) - mpfr_set_nan (MPC_RE (rop_sin)); + if (!mpfr_zero_p (mpc_realref (op_loc))) + mpfr_set_nan (mpc_realref (rop_sin)); } else /* op = NaN + i*y */ - if (!mpfr_inf_p (MPC_IM (op_loc)) && !mpfr_zero_p (MPC_IM (op_loc))) + if (!mpfr_inf_p (mpc_imagref (op_loc)) && !mpfr_zero_p (mpc_imagref (op_loc))) /* sin(NaN -i*Inf) = NaN -i*Inf */ /* sin(NaN -i*0) = NaN -i*0 */ /* sin(NaN +i*0) = NaN +i*0 */ /* sin(NaN +i*Inf) = NaN +i*Inf */ /* sin(NaN +i*y) = NaN +i*NaN, when 0<|y|<Inf */ - mpfr_set_nan (MPC_IM (rop_sin)); + mpfr_set_nan (mpc_imagref (rop_sin)); } - else if (mpfr_inf_p (MPC_RE (op_loc))) { - mpfr_set_nan (MPC_RE (rop_sin)); + else if (mpfr_inf_p (mpc_realref (op_loc))) { + mpfr_set_nan (mpc_realref (rop_sin)); - if (!mpfr_inf_p (MPC_IM (op_loc)) && !mpfr_zero_p (MPC_IM (op_loc))) + if (!mpfr_inf_p (mpc_imagref (op_loc)) && !mpfr_zero_p (mpc_imagref (op_loc))) /* sin(+/-Inf +i*y) = NaN +i*NaN, when 0<|y|<Inf */ - mpfr_set_nan (MPC_IM (rop_sin)); + mpfr_set_nan (mpc_imagref (rop_sin)); else /* sin(+/-Inf -i*Inf) = NaN -i*Inf */ /* sin(+/-Inf +i*Inf) = NaN +i*Inf */ /* sin(+/-Inf -i*0) = NaN -i*0 */ /* sin(+/-Inf +i*0) = NaN +i*0 */ - mpfr_set (MPC_IM (rop_sin), MPC_IM (op_loc), MPC_RND_IM (rnd_sin)); + mpfr_set (mpc_imagref (rop_sin), mpc_imagref (op_loc), MPC_RND_IM (rnd_sin)); } - else if (mpfr_zero_p (MPC_RE (op_loc))) { + else if (mpfr_zero_p (mpc_realref (op_loc))) { /* sin(-0 -i*Inf) = -0 -i*Inf */ /* sin(+0 -i*Inf) = +0 -i*Inf */ /* sin(-0 +i*Inf) = -0 +i*Inf */ @@ -81,46 +81,46 @@ mpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, mpfr_t s, c; mpfr_init2 (s, 2); mpfr_init2 (c, 2); - mpfr_sin_cos (s, c, MPC_RE (op_loc), GMP_RNDZ); - mpfr_set_inf (MPC_RE (rop_sin), MPFR_SIGN (s)); - mpfr_set_inf (MPC_IM (rop_sin), MPFR_SIGN (c)*MPFR_SIGN (MPC_IM (op_loc))); + mpfr_sin_cos (s, c, mpc_realref (op_loc), GMP_RNDZ); + mpfr_set_inf (mpc_realref (rop_sin), MPFR_SIGN (s)); + mpfr_set_inf (mpc_imagref (rop_sin), MPFR_SIGN (c)*MPFR_SIGN (mpc_imagref (op_loc))); mpfr_clear (s); mpfr_clear (c); } } if (rop_cos != NULL) { - if (mpfr_nan_p (MPC_RE (op_loc))) { + if (mpfr_nan_p (mpc_realref (op_loc))) { /* cos(NaN + i * NaN) = NaN + i * NaN */ /* cos(NaN - i * Inf) = +Inf + i * NaN */ /* cos(NaN + i * Inf) = +Inf + i * NaN */ /* cos(NaN - i * 0) = NaN - i * 0 */ /* cos(NaN + i * 0) = NaN + i * 0 */ /* cos(NaN + i * y) = NaN + i * NaN, when y != 0 */ - if (mpfr_inf_p (MPC_IM (op_loc))) - mpfr_set_inf (MPC_RE (rop_cos), +1); + if (mpfr_inf_p (mpc_imagref (op_loc))) + mpfr_set_inf (mpc_realref (rop_cos), +1); else - mpfr_set_nan (MPC_RE (rop_cos)); + mpfr_set_nan (mpc_realref (rop_cos)); - if (mpfr_zero_p (MPC_IM (op_loc))) - mpfr_set (MPC_IM (rop_cos), MPC_IM (op_loc), MPC_RND_IM (rnd_cos)); + if (mpfr_zero_p (mpc_imagref (op_loc))) + mpfr_set (mpc_imagref (rop_cos), mpc_imagref (op_loc), MPC_RND_IM (rnd_cos)); else - mpfr_set_nan (MPC_IM (rop_cos)); + mpfr_set_nan (mpc_imagref (rop_cos)); } - else if (mpfr_nan_p (MPC_IM (op_loc))) { + else if (mpfr_nan_p (mpc_imagref (op_loc))) { /* cos(-Inf + i * NaN) = NaN + i * NaN */ /* cos(+Inf + i * NaN) = NaN + i * NaN */ /* cos(-0 + i * NaN) = NaN - i * 0 */ /* cos(+0 + i * NaN) = NaN + i * 0 */ /* cos(x + i * NaN) = NaN + i * NaN, when x != 0 */ - if (mpfr_zero_p (MPC_RE (op_loc))) - mpfr_set (MPC_IM (rop_cos), MPC_RE (op_loc), MPC_RND_IM (rnd_cos)); + if (mpfr_zero_p (mpc_realref (op_loc))) + mpfr_set (mpc_imagref (rop_cos), mpc_realref (op_loc), MPC_RND_IM (rnd_cos)); else - mpfr_set_nan (MPC_IM (rop_cos)); + mpfr_set_nan (mpc_imagref (rop_cos)); - mpfr_set_nan (MPC_RE (rop_cos)); + mpfr_set_nan (mpc_realref (rop_cos)); } - else if (mpfr_inf_p (MPC_RE (op_loc))) { + else if (mpfr_inf_p (mpc_realref (op_loc))) { /* cos(-Inf -i*Inf) = cos(+Inf +i*Inf) = -Inf +i*NaN */ /* cos(-Inf +i*Inf) = cos(+Inf -i*Inf) = +Inf +i*NaN */ /* cos(-Inf -i*0) = cos(+Inf +i*0) = NaN -i*0 */ @@ -128,28 +128,28 @@ mpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, /* cos(-Inf +i*y) = cos(+Inf +i*y) = NaN +i*NaN, when y != 0 */ const int same_sign = - mpfr_signbit (MPC_RE (op_loc)) == mpfr_signbit (MPC_IM (op_loc)); + mpfr_signbit (mpc_realref (op_loc)) == mpfr_signbit (mpc_imagref (op_loc)); - if (mpfr_inf_p (MPC_IM (op_loc))) - mpfr_set_inf (MPC_RE (rop_cos), (same_sign ? -1 : +1)); + if (mpfr_inf_p (mpc_imagref (op_loc))) + mpfr_set_inf (mpc_realref (rop_cos), (same_sign ? -1 : +1)); else - mpfr_set_nan (MPC_RE (rop_cos)); + mpfr_set_nan (mpc_realref (rop_cos)); - if (mpfr_zero_p (MPC_IM (op_loc))) - mpfr_setsign (MPC_IM (rop_cos), MPC_IM (op_loc), same_sign, + if (mpfr_zero_p (mpc_imagref (op_loc))) + mpfr_setsign (mpc_imagref (rop_cos), mpc_imagref (op_loc), same_sign, MPC_RND_IM(rnd_cos)); else - mpfr_set_nan (MPC_IM (rop_cos)); + mpfr_set_nan (mpc_imagref (rop_cos)); } - else if (mpfr_zero_p (MPC_RE (op_loc))) { + else if (mpfr_zero_p (mpc_realref (op_loc))) { /* cos(-0 -i*Inf) = cos(+0 +i*Inf) = +Inf -i*0 */ /* cos(-0 +i*Inf) = cos(+0 -i*Inf) = +Inf +i*0 */ const int same_sign = - mpfr_signbit (MPC_RE (op_loc)) == mpfr_signbit (MPC_IM (op_loc)); + mpfr_signbit (mpc_realref (op_loc)) == mpfr_signbit (mpc_imagref (op_loc)); - mpfr_setsign (MPC_IM (rop_cos), MPC_RE (op_loc), same_sign, + mpfr_setsign (mpc_imagref (rop_cos), mpc_realref (op_loc), same_sign, MPC_RND_IM (rnd_cos)); - mpfr_set_inf (MPC_RE (rop_cos), +1); + mpfr_set_inf (mpc_realref (rop_cos), +1); } else { /* cos(x -i*Inf) = +Inf*cos(x) +i*Inf*sin(x), when x != 0 */ @@ -157,10 +157,10 @@ mpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, mpfr_t s, c; mpfr_init2 (c, 2); mpfr_init2 (s, 2); - mpfr_sin_cos (s, c, MPC_RE (op_loc), GMP_RNDN); - mpfr_set_inf (MPC_RE (rop_cos), mpfr_sgn (c)); - mpfr_set_inf (MPC_IM (rop_cos), - (mpfr_sgn (MPC_IM (op_loc)) == mpfr_sgn (s) ? -1 : +1)); + mpfr_sin_cos (s, c, mpc_realref (op_loc), GMP_RNDN); + mpfr_set_inf (mpc_realref (rop_cos), mpfr_sgn (c)); + mpfr_set_inf (mpc_imagref (rop_cos), + (mpfr_sgn (mpc_imagref (op_loc)) == mpfr_sgn (s) ? -1 : +1)); mpfr_clear (s); mpfr_clear (c); } @@ -184,7 +184,7 @@ mpc_sin_cos_real (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, by definition, for not computed values. */ mpfr_t s, c; int inex_s, inex_c; - int sign_im = mpfr_signbit (MPC_IM (op)); + int sign_im = mpfr_signbit (mpc_imagref (op)); /* sin(x +-0*i) = sin(x) +-0*i*sign(cos(x)) */ /* cos(x +-i*0) = cos(x) -+i*0*sign(sin(x)) */ @@ -196,24 +196,24 @@ mpc_sin_cos_real (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, mpfr_init2 (c, MPC_PREC_RE (rop_cos)); else mpfr_init2 (c, 2); - inex_s = mpfr_sin (s, MPC_RE (op), MPC_RND_RE (rnd_sin)); - inex_c = mpfr_cos (c, MPC_RE (op), MPC_RND_RE (rnd_cos)); + inex_s = mpfr_sin (s, mpc_realref (op), MPC_RND_RE (rnd_sin)); + inex_c = mpfr_cos (c, mpc_realref (op), MPC_RND_RE (rnd_cos)); /* We cannot use mpfr_sin_cos since we may need two distinct rounding modes and the exact return values. If we need only the sign, an arbitrary rounding mode will work. */ if (rop_sin != NULL) { - mpfr_set (MPC_RE (rop_sin), s, GMP_RNDN); /* exact */ + mpfr_set (mpc_realref (rop_sin), s, GMP_RNDN); /* exact */ inex_sin_re = inex_s; - mpfr_set_zero (MPC_IM (rop_sin), + mpfr_set_zero (mpc_imagref (rop_sin), ( ( sign_im && !mpfr_signbit(c)) || (!sign_im && mpfr_signbit(c)) ? -1 : 1)); } if (rop_cos != NULL) { - mpfr_set (MPC_RE (rop_cos), c, GMP_RNDN); /* exact */ + mpfr_set (mpc_realref (rop_cos), c, GMP_RNDN); /* exact */ inex_cos_re = inex_c; - mpfr_set_zero (MPC_IM (rop_cos), + mpfr_set_zero (mpc_imagref (rop_cos), ( ( sign_im && mpfr_signbit(s)) || (!sign_im && !mpfr_signbit(s)) ? -1 : 1)); } @@ -245,19 +245,19 @@ mpc_sin_cos_imag (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, if (rop_sin != NULL) { /* sin(+-O +i*y) = +-0 +i*sinh(y) */ - mpfr_set (MPC_RE(rop_sin), MPC_RE(op_loc), GMP_RNDN); - inex_sin_im = mpfr_sinh (MPC_IM(rop_sin), MPC_IM(op_loc), MPC_RND_IM(rnd_sin)); + mpfr_set (mpc_realref(rop_sin), mpc_realref(op_loc), GMP_RNDN); + inex_sin_im = mpfr_sinh (mpc_imagref(rop_sin), mpc_imagref(op_loc), MPC_RND_IM(rnd_sin)); } if (rop_cos != NULL) { /* cos(-0 - i * y) = cos(+0 + i * y) = cosh(y) - i * 0, cos(-0 + i * y) = cos(+0 - i * y) = cosh(y) + i * 0, where y > 0 */ - inex_cos_re = mpfr_cosh (MPC_RE (rop_cos), MPC_IM (op_loc), MPC_RND_RE (rnd_cos)); + inex_cos_re = mpfr_cosh (mpc_realref (rop_cos), mpc_imagref (op_loc), MPC_RND_RE (rnd_cos)); - mpfr_set_ui (MPC_IM (rop_cos), 0ul, MPC_RND_IM (rnd_cos)); - if (mpfr_signbit (MPC_RE (op_loc)) == mpfr_signbit (MPC_IM (op_loc))) - MPFR_CHANGE_SIGN (MPC_IM (rop_cos)); + mpfr_set_ui (mpc_imagref (rop_cos), 0ul, MPC_RND_IM (rnd_cos)); + if (mpfr_signbit (mpc_realref (op_loc)) == mpfr_signbit (mpc_imagref (op_loc))) + MPFR_CHANGE_SIGN (mpc_imagref (rop_cos)); } if (overlap) @@ -276,9 +276,9 @@ mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, { if (!mpc_fin_p (op)) return mpc_sin_cos_nonfinite (rop_sin, rop_cos, op, rnd_sin, rnd_cos); - else if (mpfr_zero_p (MPC_IM (op))) + else if (mpfr_zero_p (mpc_imagref (op))) return mpc_sin_cos_real (rop_sin, rop_cos, op, rnd_sin, rnd_cos); - else if (mpfr_zero_p (MPC_RE (op))) + else if (mpfr_zero_p (mpc_realref (op))) return mpc_sin_cos_imag (rop_sin, rop_cos, op, rnd_sin, rnd_cos); else { /* let op = a + i*b, then sin(op) = sin(a)*cosh(b) + i*cos(a)*sinh(b) @@ -325,8 +325,8 @@ mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, mpfr_set_prec (sch, prec); mpfr_set_prec (csh, prec); - mpfr_sin_cos (s, c, MPC_RE(op), GMP_RNDN); - mpfr_sinh_cosh (sh, ch, MPC_IM(op), GMP_RNDN); + mpfr_sin_cos (s, c, mpc_realref(op), GMP_RNDN); + mpfr_sinh_cosh (sh, ch, mpc_imagref(op), GMP_RNDN); if (rop_sin != NULL) { /* real part of sine */ @@ -367,10 +367,10 @@ mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, } while (ok == 0); if (rop_sin != NULL) { - inex_re = mpfr_set (MPC_RE (rop_sin), sch, MPC_RND_RE (rnd_sin)); + inex_re = mpfr_set (mpc_realref (rop_sin), sch, MPC_RND_RE (rnd_sin)); if (mpfr_inf_p (sch)) inex_re = mpfr_sgn (sch); - inex_im = mpfr_set (MPC_IM (rop_sin), csh, MPC_RND_IM (rnd_sin)); + inex_im = mpfr_set (mpc_imagref (rop_sin), csh, MPC_RND_IM (rnd_sin)); if (mpfr_inf_p (csh)) inex_im = mpfr_sgn (csh); inex_sin = MPC_INEX (inex_re, inex_im); @@ -379,10 +379,10 @@ mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, inex_sin = MPC_INEX (0,0); /* return exact if not computed */ if (rop_cos != NULL) { - inex_re = mpfr_set (MPC_RE (rop_cos), c, MPC_RND_RE (rnd_cos)); + inex_re = mpfr_set (mpc_realref (rop_cos), c, MPC_RND_RE (rnd_cos)); if (mpfr_inf_p (c)) inex_re = mpfr_sgn (c); - inex_im = mpfr_set (MPC_IM (rop_cos), s, MPC_RND_IM (rnd_cos)); + inex_im = mpfr_set (mpc_imagref (rop_cos), s, MPC_RND_IM (rnd_cos)); if (mpfr_inf_p (s)) inex_im = mpfr_sgn (s); inex_cos = MPC_INEX (inex_re, inex_im); @@ -1,6 +1,6 @@ /* mpc_sinh -- hyperbolic sine of a complex number. -Copyright (C)2008, 2009 INRIA +Copyright (C)2008, 2009, 2011 INRIA This file is part of GNU MPC. @@ -31,16 +31,16 @@ mpc_sinh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* z := conj(-i * op) and rop = conj(-i * sin(z)), in other words, we have to switch real and imaginary parts. Let us set them without copying significands. */ - MPC_RE (z)[0] = MPC_IM (op)[0]; - MPC_IM (z)[0] = MPC_RE (op)[0]; - MPC_RE (sin_z)[0] = MPC_IM (rop)[0]; - MPC_IM (sin_z)[0] = MPC_RE (rop)[0]; + mpc_realref (z)[0] = mpc_imagref (op)[0]; + mpc_imagref (z)[0] = mpc_realref (op)[0]; + mpc_realref (sin_z)[0] = mpc_imagref (rop)[0]; + mpc_imagref (sin_z)[0] = mpc_realref (rop)[0]; inex = mpc_sin (sin_z, z, RNDC (MPC_RND_IM (rnd), MPC_RND_RE (rnd))); /* sin_z and rop parts share the same significands, copy the rest now. */ - MPC_RE (rop)[0] = MPC_IM (sin_z)[0]; - MPC_IM (rop)[0] = MPC_RE (sin_z)[0]; + mpc_realref (rop)[0] = mpc_imagref (sin_z)[0]; + mpc_imagref (rop)[0] = mpc_realref (sin_z)[0]; /* swap inexact flags for real and imaginary parts */ return MPC_INEX (MPC_INEX_IM (inex), MPC_INEX_RE (inex)); @@ -35,32 +35,32 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* special values: NaN and infinities */ if (!mpc_fin_p (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)); + 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)); } - else if (mpfr_inf_p (MPC_RE (op))) { - if (mpfr_inf_p (MPC_IM (op))) { - mpfr_set_inf (MPC_IM (rop), - MPFR_SIGN (MPC_RE (op)) * MPFR_SIGN (MPC_IM (op))); - mpfr_set_nan (MPC_RE (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_zero_p (MPC_IM (op))) - mpfr_set_nan (MPC_IM (rop)); + if (mpfr_zero_p (mpc_imagref (op))) + mpfr_set_nan (mpc_imagref (rop)); else - mpfr_set_inf (MPC_IM (rop), - MPFR_SIGN (MPC_RE (op)) * MPFR_SIGN (MPC_IM (op))); - mpfr_set_inf (MPC_RE (rop), +1); + mpfr_set_inf (mpc_imagref (rop), + MPFR_SIGN (mpc_realref (op)) * MPFR_SIGN (mpc_imagref (op))); + mpfr_set_inf (mpc_realref (rop), +1); } } else /* IM(op) is infinity, RE(op) is not */ { - if (mpfr_zero_p (MPC_RE (op))) - mpfr_set_nan (MPC_IM (rop)); + if (mpfr_zero_p (mpc_realref (op))) + mpfr_set_nan (mpc_imagref (rop)); else - mpfr_set_inf (MPC_IM (rop), - MPFR_SIGN (MPC_RE (op)) * MPFR_SIGN (MPC_IM (op))); - mpfr_set_inf (MPC_RE (rop), -1); + mpfr_set_inf (mpc_imagref (rop), + MPFR_SIGN (mpc_realref (op)) * MPFR_SIGN (mpc_imagref (op))); + mpfr_set_inf (mpc_realref (rop), -1); } return MPC_INEX (0, 0); /* exact */ } @@ -68,21 +68,21 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) prec = MPC_MAX_PREC(rop); /* first check for real resp. purely imaginary number */ - if (mpfr_zero_p (MPC_IM(op))) + if (mpfr_zero_p (mpc_imagref(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); + 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); if (!same_sign) mpc_conj (rop, rop, MPC_RNDNN); return MPC_INEX(inex_re, inex_im); } - if (mpfr_zero_p (MPC_RE(op))) + if (mpfr_zero_p (mpc_realref(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); + 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 (!same_sign) mpc_conj (rop, rop, MPC_RNDNN); return MPC_INEX(inex_re, inex_im); @@ -92,14 +92,14 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* exactly with the standard formulae instead, even if this means an */ /* additional multiplication. */ if (SAFE_ABS (mpfr_exp_t, - mpfr_get_exp (MPC_RE (op)) - mpfr_get_exp (MPC_IM (op))) + mpfr_get_exp (mpc_realref (op)) - mpfr_get_exp (mpc_imagref (op))) > (mpfr_exp_t) MPC_MAX_PREC (op) / 2) { mpfr_init2 (u, 2*MPC_PREC_RE (op)); mpfr_init2 (v, 2*MPC_PREC_IM (op)); - mpfr_sqr (u, MPC_RE (op), GMP_RNDN); - mpfr_sqr (v, MPC_IM (op), GMP_RNDN); /* both are exact */ + mpfr_sqr (u, mpc_realref (op), GMP_RNDN); + mpfr_sqr (v, mpc_imagref (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)); @@ -135,15 +135,15 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* 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_IM (op), MPFR_RNDA), u) - | ROUND_AWAY (mpfr_sub (v, x, MPC_IM (op), MPFR_RNDA), v); + 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 */ 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); + mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); inex_re = 0; ok = 1; } @@ -156,15 +156,15 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) MPC_ASSERT (mpfr_get_exp (u) != emin); if (mpfr_inf_p (u)) { - /* let MPC_RE(rop) be a "correctly rounded overflow" */ - inex_re = mpfr_set_ui_2exp (MPC_RE (rop), 1, emax, MPC_RND_RE (rnd)); + /* let mpc_realref(rop) be a "correctly rounded overflow" */ + inex_re = mpfr_set_ui_2exp (mpc_realref (rop), 1, emax, MPC_RND_RE (rnd)); break; } ok = (!inexact) | mpfr_can_round (u, prec - 3, GMP_RNDU, GMP_RNDZ, MPC_PREC_RE (rop) + (MPC_RND_RE (rnd) == GMP_RNDN)); if (ok) { - inex_re = mpfr_set (MPC_RE (rop), u, MPC_RND_RE (rnd)); + 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; @@ -184,11 +184,11 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) if (rnd_re == GMP_RNDZ || rnd_re == GMP_RNDN || rnd_re == GMP_RNDU) { - mpfr_set_ui (MPC_RE (rop), 0, rnd_re); + mpfr_set_ui (mpc_realref (rop), 0, rnd_re); inex_re = 1; } else /* round down or away from zero */ { - mpfr_set (MPC_RE (rop), u, rnd_re); + mpfr_set (mpc_realref (rop), u, rnd_re); inex_re = -1; } break; @@ -197,7 +197,7 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) MPC_PREC_RE (rop) + (MPC_RND_RE (rnd) == GMP_RNDN)); if (ok) { - inex_re = mpfr_set (MPC_RE (rop), u, MPC_RND_RE (rnd)); + inex_re = mpfr_set (mpc_realref (rop), u, MPC_RND_RE (rnd)); if (inex_re == 0) inex_re = inexact; } @@ -206,17 +206,17 @@ mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) 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) + if (mpfr_get_exp (x) + mpfr_get_exp(mpc_imagref (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)); + inex_im = mpfr_mul (mpc_imagref (rop), x, mpc_imagref (op), MPC_RND_IM (rnd)); } else { - inex_im = mpfr_mul (MPC_IM (rop), x, MPC_IM (op), MPC_RND_IM (rnd)); + inex_im = mpfr_mul (mpc_imagref (rop), x, mpc_imagref (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); + MPC_ASSERT (mpfr_zero_p (mpc_imagref (rop)) == 0); + mpfr_mul_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN); } mpfr_clear (u); @@ -1,6 +1,6 @@ /* mpc_sqrt -- Take the square root of a complex number. -Copyright (C) 2002, 2008, 2009, 2010 INRIA +Copyright (C) 2002, 2008, 2009, 2010, 2011 INRIA This file is part of GNU MPC. @@ -38,13 +38,13 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) /* be either the real or the imaginary part of a */ mpfr_prec_t prec; int inex_w, inex_t = 1, inex_re, inex_im, loops = 0; - const int re_cmp = mpfr_cmp_ui (MPC_RE (b), 0), - im_cmp = mpfr_cmp_ui (MPC_IM (b), 0); + const int re_cmp = mpfr_cmp_ui (mpc_realref (b), 0), + im_cmp = mpfr_cmp_ui (mpc_imagref (b), 0); /* comparison of the real/imaginary part of b with 0 */ int repr_w, repr_t = 0 /* to avoid gcc warning */ ; /* flag indicating whether the computed value is already representable at the target precision */ - const int im_sgn = mpfr_signbit (MPC_IM (b)) == 0 ? 0 : -1; + const int im_sgn = mpfr_signbit (mpc_imagref (b)) == 0 ? 0 : -1; /* we need to know the sign of Im(b) when it is +/-0 */ const mpfr_rnd_t r = im_sgn ? GMP_RNDD : GMP_RNDU; /* rounding mode used when computing t */ @@ -53,41 +53,41 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) if (!mpc_fin_p (b)) { /* sqrt(x +i*Inf) = +Inf +I*Inf, even if x = NaN */ /* sqrt(x -i*Inf) = +Inf -I*Inf, even if x = NaN */ - if (mpfr_inf_p (MPC_IM (b))) + if (mpfr_inf_p (mpc_imagref (b))) { - mpfr_set_inf (MPC_RE (a), +1); - mpfr_set_inf (MPC_IM (a), im_sgn); + mpfr_set_inf (mpc_realref (a), +1); + mpfr_set_inf (mpc_imagref (a), im_sgn); return MPC_INEX (0, 0); } - if (mpfr_inf_p (MPC_RE (b))) + if (mpfr_inf_p (mpc_realref (b))) { - if (mpfr_signbit (MPC_RE (b))) + if (mpfr_signbit (mpc_realref (b))) { - if (mpfr_number_p (MPC_IM (b))) + if (mpfr_number_p (mpc_imagref (b))) { /* sqrt(-Inf +i*y) = +0 +i*Inf, when y positive */ /* sqrt(-Inf +i*y) = +0 -i*Inf, when y positive */ - mpfr_set_ui (MPC_RE (a), 0, GMP_RNDN); - mpfr_set_inf (MPC_IM (a), im_sgn); + mpfr_set_ui (mpc_realref (a), 0, GMP_RNDN); + mpfr_set_inf (mpc_imagref (a), im_sgn); return MPC_INEX (0, 0); } else { /* sqrt(-Inf +i*NaN) = NaN +/-i*Inf */ - mpfr_set_nan (MPC_RE (a)); - mpfr_set_inf (MPC_IM (a), im_sgn); + mpfr_set_nan (mpc_realref (a)); + mpfr_set_inf (mpc_imagref (a), im_sgn); return MPC_INEX (0, 0); } } else { - if (mpfr_number_p (MPC_IM (b))) + if (mpfr_number_p (mpc_imagref (b))) { /* sqrt(+Inf +i*y) = +Inf +i*0, when y positive */ /* sqrt(+Inf +i*y) = +Inf -i*0, when y positive */ - mpfr_set_inf (MPC_RE (a), +1); - mpfr_set_ui (MPC_IM (a), 0, GMP_RNDN); + mpfr_set_inf (mpc_realref (a), +1); + mpfr_set_ui (mpc_imagref (a), 0, GMP_RNDN); if (im_sgn) mpc_conj (a, a, MPC_RNDNN); return MPC_INEX (0, 0); @@ -104,10 +104,10 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) /* sqrt(x +i*NaN) = NaN +i*NaN, if x is not infinite */ /* sqrt(NaN +i*y) = NaN +i*NaN, if y is not infinite */ - if (mpfr_nan_p (MPC_RE (b)) || mpfr_nan_p (MPC_IM (b))) + if (mpfr_nan_p (mpc_realref (b)) || mpfr_nan_p (mpc_imagref (b))) { - mpfr_set_nan (MPC_RE (a)); - mpfr_set_nan (MPC_IM (a)); + mpfr_set_nan (mpc_realref (a)); + mpfr_set_nan (mpc_imagref (a)); return MPC_INEX (0, 0); } } @@ -124,8 +124,8 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) } else if (re_cmp > 0) { - inex_w = mpfr_sqrt (MPC_RE (a), MPC_RE (b), MPC_RND_RE (rnd)); - mpfr_set_ui (MPC_IM (a), 0, GMP_RNDN); + inex_w = mpfr_sqrt (mpc_realref (a), mpc_realref (b), MPC_RND_RE (rnd)); + mpfr_set_ui (mpc_imagref (a), 0, GMP_RNDN); if (im_sgn) mpc_conj (a, a, MPC_RNDNN); return MPC_INEX (inex_w, 0); @@ -133,16 +133,16 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) else { mpfr_init2 (w, MPC_PREC_RE (b)); - mpfr_neg (w, MPC_RE (b), GMP_RNDN); + mpfr_neg (w, mpc_realref (b), GMP_RNDN); if (im_sgn) { - inex_w = -mpfr_sqrt (MPC_IM (a), w, INV_RND (MPC_RND_IM (rnd))); - mpfr_neg (MPC_IM (a), MPC_IM (a), GMP_RNDN); + inex_w = -mpfr_sqrt (mpc_imagref (a), w, INV_RND (MPC_RND_IM (rnd))); + mpfr_neg (mpc_imagref (a), mpc_imagref (a), GMP_RNDN); } else - inex_w = mpfr_sqrt (MPC_IM (a), w, MPC_RND_IM (rnd)); + inex_w = mpfr_sqrt (mpc_imagref (a), w, MPC_RND_IM (rnd)); - mpfr_set_ui (MPC_RE (a), 0, GMP_RNDN); + mpfr_set_ui (mpc_realref (a), 0, GMP_RNDN); mpfr_clear (w); return MPC_INEX (0, inex_w); } @@ -153,20 +153,20 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) { mpfr_t y; - y[0] = MPC_IM (b)[0]; + y[0] = mpc_imagref (b)[0]; /* If y/2 underflows, so does sqrt(y/2) */ mpfr_div_2ui (y, y, 1, GMP_RNDN); if (im_cmp > 0) { - inex_w = mpfr_sqrt (MPC_RE (a), y, MPC_RND_RE (rnd)); - inex_t = mpfr_sqrt (MPC_IM (a), y, MPC_RND_IM (rnd)); + inex_w = mpfr_sqrt (mpc_realref (a), y, MPC_RND_RE (rnd)); + inex_t = mpfr_sqrt (mpc_imagref (a), y, MPC_RND_IM (rnd)); } else { mpfr_neg (y, y, GMP_RNDN); - inex_w = mpfr_sqrt (MPC_RE (a), y, MPC_RND_RE (rnd)); - inex_t = -mpfr_sqrt (MPC_IM (a), y, INV_RND (MPC_RND_IM (rnd))); - mpfr_neg (MPC_IM (a), MPC_IM (a), GMP_RNDN); + inex_w = mpfr_sqrt (mpc_realref (a), y, MPC_RND_RE (rnd)); + inex_t = -mpfr_sqrt (mpc_imagref (a), y, INV_RND (MPC_RND_IM (rnd))); + mpfr_neg (mpc_imagref (a), mpc_imagref (a), GMP_RNDN); } return MPC_INEX (inex_w, inex_t); } @@ -213,9 +213,9 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) /* total error bounded by 3 ulps */ inex_w = mpc_abs (w, b, GMP_RNDD); if (re_cmp < 0) - inex_w |= mpfr_sub (w, w, MPC_RE (b), GMP_RNDD); + inex_w |= mpfr_sub (w, w, mpc_realref (b), GMP_RNDD); else - inex_w |= mpfr_add (w, w, MPC_RE (b), GMP_RNDD); + inex_w |= mpfr_add (w, w, mpc_realref (b), GMP_RNDD); inex_w |= mpfr_div_2ui (w, w, 1, GMP_RNDD); inex_w |= mpfr_sqrt (w, w, GMP_RNDD); @@ -240,7 +240,7 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) if (!inex_w || ok_w) { /* t = y / 2w, rounded away */ /* total error bounded by 7 ulps */ - inex_t = mpfr_div (t, MPC_IM (b), w, r); + inex_t = mpfr_div (t, mpc_imagref (b), w, r); if (!inex_t && inex_w) /* The division was exact, but w was not. */ inex_t = im_sgn ? -1 : 1; @@ -262,16 +262,16 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) while ((inex_w && !ok_w) || (inex_t && !ok_t)); if (re_cmp > 0) { - inex_re = mpfr_set (MPC_RE (a), w, MPC_RND_RE(rnd)); - inex_im = mpfr_set (MPC_IM (a), t, MPC_RND_IM(rnd)); + inex_re = mpfr_set (mpc_realref (a), w, MPC_RND_RE(rnd)); + inex_im = mpfr_set (mpc_imagref (a), t, MPC_RND_IM(rnd)); } else if (im_cmp > 0) { - inex_re = mpfr_set (MPC_RE(a), t, MPC_RND_RE(rnd)); - inex_im = mpfr_set (MPC_IM(a), w, MPC_RND_IM(rnd)); + inex_re = mpfr_set (mpc_realref(a), t, MPC_RND_RE(rnd)); + inex_im = mpfr_set (mpc_imagref(a), w, MPC_RND_IM(rnd)); } else { - inex_re = mpfr_neg (MPC_RE (a), t, MPC_RND_RE(rnd)); - inex_im = mpfr_neg (MPC_IM (a), w, MPC_RND_IM(rnd)); + inex_re = mpfr_neg (mpc_realref (a), t, MPC_RND_RE(rnd)); + inex_im = mpfr_neg (mpc_imagref (a), w, MPC_RND_IM(rnd)); } if (repr_w && inex_w) { @@ -290,7 +290,7 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) be done now when we are in the target precision */ if (re_cmp > 0) { if (rnd_w == GMP_RNDU) { - MPFR_ADD_ONE_ULP (MPC_RE (a)); + MPFR_ADD_ONE_ULP (mpc_realref (a)); inex_re = +1; } else @@ -298,7 +298,7 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) } else if (im_cmp > 0) { if (rnd_w == GMP_RNDU) { - MPFR_ADD_ONE_ULP (MPC_IM (a)); + MPFR_ADD_ONE_ULP (mpc_imagref (a)); inex_im = +1; } else @@ -306,7 +306,7 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) } else { if (rnd_w == GMP_RNDU) { - MPFR_ADD_ONE_ULP (MPC_IM (a)); + MPFR_ADD_ONE_ULP (mpc_imagref (a)); inex_im = -1; } else @@ -331,9 +331,9 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) inex_im = -inex_t; if ( (im_cmp > 0 && r == GMP_RNDD) || (im_cmp < 0 && r == GMP_RNDU)) - MPFR_ADD_ONE_ULP (MPC_IM (a)); + MPFR_ADD_ONE_ULP (mpc_imagref (a)); else - MPFR_SUB_ONE_ULP (MPC_IM (a)); + MPFR_SUB_ONE_ULP (mpc_imagref (a)); } } else if (im_cmp > 0) { @@ -342,9 +342,9 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) else { inex_re = -inex_t; if (r == GMP_RNDD) - MPFR_ADD_ONE_ULP (MPC_RE (a)); + MPFR_ADD_ONE_ULP (mpc_realref (a)); else - MPFR_SUB_ONE_ULP (MPC_RE (a)); + MPFR_SUB_ONE_ULP (mpc_realref (a)); } } else { @@ -353,9 +353,9 @@ mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd) else { inex_re = inex_t; if (r == GMP_RNDD) - MPFR_SUB_ONE_ULP (MPC_RE (a)); + MPFR_SUB_ONE_ULP (mpc_realref (a)); else - MPFR_ADD_ONE_ULP (MPC_RE (a)); + MPFR_ADD_ONE_ULP (mpc_realref (a)); } } } diff --git a/src/strtoc.c b/src/strtoc.c index 0fb2134..b96ccee 100644 --- a/src/strtoc.c +++ b/src/strtoc.c @@ -1,6 +1,6 @@ /* mpc_strtoc -- Read a complex number from a string. -Copyright (C) 2009, 2010 INRIA +Copyright (C) 2009, 2010, 2011 INRIA This file is part of GNU MPC. @@ -51,20 +51,20 @@ mpc_strtoc (mpc_ptr rop, const char *nptr, char **endptr, int base, mpc_rnd_t rn ++p; } - inex_re = mpfr_strtofr (MPC_RE(rop), p, &end, base, MPC_RND_RE (rnd)); + inex_re = mpfr_strtofr (mpc_realref(rop), p, &end, base, MPC_RND_RE (rnd)); if (end == p) goto error; p = end; if (!bracketed) - inex_im = mpfr_set_ui (MPC_IM (rop), 0ul, GMP_RNDN); + inex_im = mpfr_set_ui (mpc_imagref (rop), 0ul, GMP_RNDN); else { if (!isspace ((unsigned char)*p)) goto error; skip_whitespace (&p); - inex_im = mpfr_strtofr (MPC_IM(rop), p, &end, base, MPC_RND_IM (rnd)); + inex_im = mpfr_strtofr (mpc_imagref(rop), p, &end, base, MPC_RND_IM (rnd)); if (end == p) goto error; p = end; @@ -83,7 +83,7 @@ mpc_strtoc (mpc_ptr rop, const char *nptr, char **endptr, int base, mpc_rnd_t rn error: if (endptr != NULL) *endptr = (char*) nptr; - mpfr_set_nan (MPC_RE (rop)); - mpfr_set_nan (MPC_IM (rop)); + mpfr_set_nan (mpc_realref (rop)); + mpfr_set_nan (mpc_imagref (rop)); return -1; } @@ -1,6 +1,6 @@ /* mpc_sub -- Subtract two complex numbers. -Copyright (C) 2002, 2009 INRIA +Copyright (C) 2002, 2009, 2011 INRIA This file is part of GNU MPC. @@ -25,8 +25,8 @@ mpc_sub (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_sub (MPC_RE(a), MPC_RE(b), MPC_RE(c), MPC_RND_RE(rnd)); - inex_im = mpfr_sub (MPC_IM(a), MPC_IM(b), MPC_IM(c), MPC_RND_IM(rnd)); + inex_re = mpfr_sub (mpc_realref(a), mpc_realref(b), mpc_realref(c), MPC_RND_RE(rnd)); + inex_im = mpfr_sub (mpc_imagref(a), mpc_imagref(b), mpc_imagref(c), MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } diff --git a/src/sub_fr.c b/src/sub_fr.c index 866a0c5..ade33a9 100644 --- a/src/sub_fr.c +++ b/src/sub_fr.c @@ -1,7 +1,7 @@ /* mpc_sub_fr -- Substract a floating-point number to the real part of a complex number. -Copyright (C) 2008, 2009 INRIA +Copyright (C) 2008, 2009, 2011 INRIA This file is part of GNU MPC. @@ -27,8 +27,8 @@ mpc_sub_fr (mpc_ptr a, mpc_srcptr b, mpfr_srcptr c, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_sub (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); - inex_im = mpfr_set (MPC_IM(a), MPC_IM(b), MPC_RND_IM(rnd)); + inex_re = mpfr_sub (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_set (mpc_imagref(a), mpc_imagref(b), MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } diff --git a/src/sub_ui.c b/src/sub_ui.c index c1a9278..561c469 100644 --- a/src/sub_ui.c +++ b/src/sub_ui.c @@ -1,6 +1,6 @@ /* mpc_sub_ui -- Add a complex number and an unsigned long int. -Copyright (C) 2002, 2009 INRIA +Copyright (C) 2002, 2009, 2011 INRIA This file is part of GNU MPC. @@ -26,8 +26,8 @@ mpc_sub_ui (mpc_ptr a, mpc_srcptr b, unsigned long int c, mpc_rnd_t rnd) { int inex_re, inex_im; - inex_re = mpfr_sub_ui (MPC_RE(a), MPC_RE(b), c, MPC_RND_RE(rnd)); - inex_im = mpfr_set (MPC_IM(a), MPC_IM(b), MPC_RND_IM(rnd)); + inex_re = mpfr_sub_ui (mpc_realref(a), mpc_realref(b), c, MPC_RND_RE(rnd)); + inex_im = mpfr_set (mpc_imagref(a), mpc_imagref(b), MPC_RND_IM(rnd)); return MPC_INEX(inex_re, inex_im); } @@ -1,6 +1,6 @@ /* mpc_swap -- Swap two complex numbers. -Copyright (C) 2009 INRIA +Copyright (C) 2009, 2011 INRIA This file is part of GNU MPC. @@ -24,6 +24,6 @@ void mpc_swap (mpc_ptr a, mpc_ptr b) { /* assumes real and imaginary parts do not overlap */ - mpfr_swap (MPC_RE(a), MPC_RE(b)); - mpfr_swap (MPC_IM(a), MPC_IM(b)); + mpfr_swap (mpc_realref(a), mpc_realref(b)); + mpfr_swap (mpc_imagref(a), mpc_imagref(b)); } @@ -34,29 +34,29 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* special values */ if (!mpc_fin_p (op)) { - if (mpfr_nan_p (MPC_RE (op))) + if (mpfr_nan_p (mpc_realref (op))) { - if (mpfr_inf_p (MPC_IM (op))) + if (mpfr_inf_p (mpc_imagref (op))) /* tan(NaN -i*Inf) = +/-0 -i */ /* tan(NaN +i*Inf) = +/-0 +i */ { /* exact unless 1 is not in exponent range */ inex = mpc_set_si_si (rop, 0, - (MPFR_SIGN (MPC_IM (op)) < 0) ? -1 : +1, + (MPFR_SIGN (mpc_imagref (op)) < 0) ? -1 : +1, rnd); } else /* tan(NaN +i*y) = NaN +i*NaN, when y is finite */ /* tan(NaN +i*NaN) = NaN +i*NaN */ { - mpfr_set_nan (MPC_RE (rop)); - mpfr_set_nan (MPC_IM (rop)); + mpfr_set_nan (mpc_realref (rop)); + mpfr_set_nan (mpc_imagref (rop)); inex = MPC_INEX (0, 0); /* always exact */ } } - else if (mpfr_nan_p (MPC_IM (op))) + else if (mpfr_nan_p (mpc_imagref (op))) { - if (mpfr_cmp_ui (MPC_RE (op), 0) == 0) + if (mpfr_cmp_ui (mpc_realref (op), 0) == 0) /* tan(-0 +i*NaN) = -0 +i*NaN */ /* tan(+0 +i*NaN) = +0 +i*NaN */ { @@ -66,28 +66,28 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) else /* tan(x +i*NaN) = NaN +i*NaN, when x != 0 */ { - mpfr_set_nan (MPC_RE (rop)); - mpfr_set_nan (MPC_IM (rop)); + mpfr_set_nan (mpc_realref (rop)); + mpfr_set_nan (mpc_imagref (rop)); inex = MPC_INEX (0, 0); /* always exact */ } } - else if (mpfr_inf_p (MPC_RE (op))) + else if (mpfr_inf_p (mpc_realref (op))) { - if (mpfr_inf_p (MPC_IM (op))) + if (mpfr_inf_p (mpc_imagref (op))) /* tan(-Inf -i*Inf) = -/+0 -i */ /* tan(-Inf +i*Inf) = -/+0 +i */ /* tan(+Inf -i*Inf) = +/-0 -i */ /* tan(+Inf +i*Inf) = +/-0 +i */ { - const int sign_re = mpfr_signbit (MPC_RE (op)); + const int sign_re = mpfr_signbit (mpc_realref (op)); int inex_im; - mpfr_set_ui (MPC_RE (rop), 0, MPC_RND_RE (rnd)); - mpfr_setsign (MPC_RE (rop), MPC_RE (rop), sign_re, GMP_RNDN); + mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd)); + mpfr_setsign (mpc_realref (rop), mpc_realref (rop), sign_re, GMP_RNDN); /* exact, unless 1 is not in exponent range */ - inex_im = mpfr_set_si (MPC_IM (rop), - mpfr_signbit (MPC_IM (op)) ? -1 : +1, + inex_im = mpfr_set_si (mpc_imagref (rop), + mpfr_signbit (mpc_imagref (op)) ? -1 : +1, MPC_RND_IM (rnd)); inex = MPC_INEX (0, inex_im); } @@ -95,8 +95,8 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* tan(-Inf +i*y) = tan(+Inf +i*y) = NaN +i*NaN, when y is finite */ { - mpfr_set_nan (MPC_RE (rop)); - mpfr_set_nan (MPC_IM (rop)); + mpfr_set_nan (mpc_realref (rop)); + mpfr_set_nan (mpc_imagref (rop)); inex = MPC_INEX (0, 0); /* always exact */ } } @@ -111,13 +111,13 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpfr_init (c); mpfr_init (s); - mpfr_sin_cos (s, c, MPC_RE (op), GMP_RNDN); - mpfr_set_ui (MPC_RE (rop), 0, MPC_RND_RE (rnd)); - mpfr_setsign (MPC_RE (rop), MPC_RE (rop), + mpfr_sin_cos (s, c, mpc_realref (op), GMP_RNDN); + mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd)); + mpfr_setsign (mpc_realref (rop), mpc_realref (rop), mpfr_signbit (c) != mpfr_signbit (s), GMP_RNDN); /* exact, unless 1 is not in exponent range */ - inex_im = mpfr_set_si (MPC_IM (rop), - (mpfr_signbit (MPC_IM (op)) ? -1 : +1), + inex_im = mpfr_set_si (mpc_imagref (rop), + (mpfr_signbit (mpc_imagref (op)) ? -1 : +1), MPC_RND_IM (rnd)); inex = MPC_INEX (0, inex_im); @@ -128,26 +128,26 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) return inex; } - if (mpfr_zero_p (MPC_RE (op))) + if (mpfr_zero_p (mpc_realref (op))) /* tan(-0 -i*y) = -0 +i*tanh(y), when y is finite. */ /* tan(+0 +i*y) = +0 +i*tanh(y), when y is finite. */ { int inex_im; - mpfr_set (MPC_RE (rop), MPC_RE (op), MPC_RND_RE (rnd)); - inex_im = mpfr_tanh (MPC_IM (rop), MPC_IM (op), MPC_RND_IM (rnd)); + mpfr_set (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd)); + inex_im = mpfr_tanh (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd)); return MPC_INEX (0, inex_im); } - if (mpfr_zero_p (MPC_IM (op))) + if (mpfr_zero_p (mpc_imagref (op))) /* tan(x -i*0) = tan(x) -i*0, when x is finite. */ /* tan(x +i*0) = tan(x) +i*0, when x is finite. */ { int inex_re; - inex_re = mpfr_tan (MPC_RE (rop), MPC_RE (op), MPC_RND_RE (rnd)); - mpfr_set (MPC_IM (rop), MPC_IM (op), MPC_RND_IM (rnd)); + inex_re = mpfr_tan (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd)); + mpfr_set (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd)); return MPC_INEX (inex_re, 0); } @@ -194,44 +194,44 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) above), sin x and cos y are never exact, so rounding away from 0 is rounding towards 0 and adding one ulp to the absolute value */ mpc_sin_cos (x, y, op, MPC_RNDZZ, MPC_RNDZZ); - MPFR_ADD_ONE_ULP (MPC_RE (x)); - MPFR_ADD_ONE_ULP (MPC_IM (x)); - MPFR_ADD_ONE_ULP (MPC_RE (y)); - MPFR_ADD_ONE_ULP (MPC_IM (y)); - MPC_ASSERT (mpfr_zero_p (MPC_RE (x)) == 0); - - if ( mpfr_inf_p (MPC_RE (x)) || mpfr_inf_p (MPC_IM (x)) - || mpfr_inf_p (MPC_RE (y)) || mpfr_inf_p (MPC_IM (y))) { + MPFR_ADD_ONE_ULP (mpc_realref (x)); + MPFR_ADD_ONE_ULP (mpc_imagref (x)); + MPFR_ADD_ONE_ULP (mpc_realref (y)); + MPFR_ADD_ONE_ULP (mpc_imagref (y)); + MPC_ASSERT (mpfr_zero_p (mpc_realref (x)) == 0); + + if ( mpfr_inf_p (mpc_realref (x)) || mpfr_inf_p (mpc_imagref (x)) + || mpfr_inf_p (mpc_realref (y)) || mpfr_inf_p (mpc_imagref (y))) { /* If the real or imaginary part of x is infinite, it means that Im(op) was large, in which case the result is sign(tan(Re(op)))*0 + sign(Im(op))*I, where sign(tan(Re(op))) = sign(Re(x))*sign(Re(y)). */ int inex_re, inex_im; - mpfr_set_ui (MPC_RE (rop), 0, GMP_RNDN); - if (mpfr_sgn (MPC_RE (x)) * mpfr_sgn (MPC_RE (y)) < 0) + mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); + if (mpfr_sgn (mpc_realref (x)) * mpfr_sgn (mpc_realref (y)) < 0) { - mpfr_neg (MPC_RE (rop), MPC_RE (rop), GMP_RNDN); + mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN); inex_re = 1; } else inex_re = -1; /* +0 is rounded down */ - if (mpfr_sgn (MPC_IM (op)) > 0) + if (mpfr_sgn (mpc_imagref (op)) > 0) { - mpfr_set_ui (MPC_IM (rop), 1, GMP_RNDN); + mpfr_set_ui (mpc_imagref (rop), 1, GMP_RNDN); inex_im = 1; } else { - mpfr_set_si (MPC_IM (rop), -1, GMP_RNDN); + mpfr_set_si (mpc_imagref (rop), -1, GMP_RNDN); inex_im = -1; } inex = MPC_INEX(inex_re, inex_im); goto end; } - exr = mpfr_get_exp (MPC_RE (x)); - eyr = mpfr_get_exp (MPC_RE (y)); - eyi = mpfr_get_exp (MPC_IM (y)); + exr = mpfr_get_exp (mpc_realref (x)); + eyr = mpfr_get_exp (mpc_realref (y)); + eyi = mpfr_get_exp (mpc_imagref (y)); /* some parts of the quotient may be exact */ inex = mpc_div (x, x, y, MPC_RNDZZ); @@ -241,17 +241,17 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) tan(1+14*I) = 1.26e-10 + 1.00*I. For small precision sin(op) and cos(op) differ only by a factor I, thus after mpc_div x = I and its real part is zero. */ - if (mpfr_zero_p (MPC_RE (x)) || mpfr_zero_p (MPC_IM (x))) + if (mpfr_zero_p (mpc_realref (x)) || mpfr_zero_p (mpc_imagref (x))) { err = prec; /* double precision */ continue; } if (MPC_INEX_RE (inex)) - MPFR_ADD_ONE_ULP (MPC_RE (x)); + MPFR_ADD_ONE_ULP (mpc_realref (x)); if (MPC_INEX_IM (inex)) - MPFR_ADD_ONE_ULP (MPC_IM (x)); - MPC_ASSERT (mpfr_zero_p (MPC_RE (x)) == 0); - ezr = mpfr_get_exp (MPC_RE (x)); + MPFR_ADD_ONE_ULP (mpc_imagref (x)); + MPC_ASSERT (mpfr_zero_p (mpc_realref (x)) == 0); + ezr = mpfr_get_exp (mpc_realref (x)); /* FIXME: compute k = Exp(Re(x))+Exp(Re(y))-2min{Exp(Re(y)), Exp(Im(y))}-Exp(Re(x/y)) @@ -260,15 +260,15 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) err = k < 2 ? 7 : (k == 2 ? 8 : (5 + k)); /* Can the real part be rounded? */ - ok = (!mpfr_number_p (MPC_RE (x))) - || mpfr_can_round (MPC_RE(x), prec - err, GMP_RNDN, GMP_RNDZ, + ok = (!mpfr_number_p (mpc_realref (x))) + || mpfr_can_round (mpc_realref(x), prec - err, GMP_RNDN, GMP_RNDZ, MPC_PREC_RE(rop) + (MPC_RND_RE(rnd) == GMP_RNDN)); if (ok) { /* Can the imaginary part be rounded? */ - ok = (!mpfr_number_p (MPC_IM (x))) - || mpfr_can_round (MPC_IM(x), prec - 6, GMP_RNDN, GMP_RNDZ, + ok = (!mpfr_number_p (mpc_imagref (x))) + || mpfr_can_round (mpc_imagref(x), prec - 6, GMP_RNDN, GMP_RNDZ, MPC_PREC_IM(rop) + (MPC_RND_IM(rnd) == GMP_RNDN)); } } @@ -1,6 +1,6 @@ /* mpc_tanh -- hyperbolic tangent of a complex number. -Copyright (C) 2008, 2009 INRIA +Copyright (C) 2008, 2009, 2011 INRIA This file is part of GNU MPC. @@ -31,16 +31,16 @@ mpc_tanh (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) /* z := conj(-i * op) and rop = conj(-i * tan(z)), in other words, we have to switch real and imaginary parts. Let us set them without copying significands. */ - MPC_RE (z)[0] = MPC_IM (op)[0]; - MPC_IM (z)[0] = MPC_RE (op)[0]; - MPC_RE (tan_z)[0] = MPC_IM (rop)[0]; - MPC_IM (tan_z)[0] = MPC_RE (rop)[0]; + mpc_realref (z)[0] = mpc_imagref (op)[0]; + mpc_imagref (z)[0] = mpc_realref (op)[0]; + mpc_realref (tan_z)[0] = mpc_imagref (rop)[0]; + mpc_imagref (tan_z)[0] = mpc_realref (rop)[0]; inex = mpc_tan (tan_z, z, RNDC (MPC_RND_IM (rnd), MPC_RND_RE (rnd))); /* tan_z and rop parts share the same significands, copy the rest now. */ - MPC_RE (rop)[0] = MPC_IM (tan_z)[0]; - MPC_IM (rop)[0] = MPC_RE (tan_z)[0]; + mpc_realref (rop)[0] = mpc_imagref (tan_z)[0]; + mpc_imagref (rop)[0] = mpc_realref (tan_z)[0]; /* swap inexact flags for real and imaginary parts */ return MPC_INEX (MPC_INEX_IM (inex), MPC_INEX_RE (inex)); diff --git a/src/ui_ui_sub.c b/src/ui_ui_sub.c index e1dbdb1..597ee05 100644 --- a/src/ui_ui_sub.c +++ b/src/ui_ui_sub.c @@ -27,8 +27,8 @@ mpc_ui_ui_sub (mpc_ptr rop, unsigned long int re, unsigned long int im, { int inex_re, inex_im; - inex_re = mpfr_ui_sub (MPC_RE (rop), re, MPC_RE (op), MPC_RND_RE (rnd)); - inex_im = mpfr_ui_sub (MPC_IM (rop), im, MPC_IM (op), MPC_RND_IM (rnd)); + inex_re = mpfr_ui_sub (mpc_realref (rop), re, mpc_realref (op), MPC_RND_RE (rnd)); + inex_im = mpfr_ui_sub (mpc_imagref (rop), im, mpc_imagref (op), MPC_RND_IM (rnd)); return MPC_INEX (inex_re, inex_im); } diff --git a/src/urandom.c b/src/urandom.c index 6734ac9..f8e5f90 100644 --- a/src/urandom.c +++ b/src/urandom.c @@ -1,7 +1,7 @@ /* mpc_urandom -- Generate a random complex number uniformly distributed in the interval [0,1[. -Copyright (C) 2008, 2009 INRIA +Copyright (C) 2008, 2009, 2011 INRIA This file is part of GNU MPC. @@ -26,7 +26,7 @@ mpc_urandom (mpc_ptr a, gmp_randstate_t state) { int r, i; - r = mpfr_urandomb (MPC_RE(a), state); - i = mpfr_urandomb (MPC_IM(a), state); + r = mpfr_urandomb (mpc_realref(a), state); + i = mpfr_urandomb (mpc_imagref(a), state); return r && i; } |