|
Packit Service |
2e9770 |
/* mpc_asin -- arcsine of a complex number.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
Copyright (C) 2009, 2010, 2011, 2012, 2013, 2014 INRIA
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
This file is part of GNU MPC.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
GNU MPC is free software; you can redistribute it and/or modify it under
|
|
Packit Service |
2e9770 |
the terms of the GNU Lesser General Public License as published by the
|
|
Packit Service |
2e9770 |
Free Software Foundation; either version 3 of the License, or (at your
|
|
Packit Service |
2e9770 |
option) any later version.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
Packit Service |
2e9770 |
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
Packit Service |
2e9770 |
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
|
|
Packit Service |
2e9770 |
more details.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
You should have received a copy of the GNU Lesser General Public License
|
|
Packit Service |
2e9770 |
along with this program. If not, see http://www.gnu.org/licenses/ .
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
#include <stdio.h>
|
|
Packit Service |
2e9770 |
#include "mpc-impl.h"
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* Special case op = 1 + i*y for tiny y (see algorithms.tex).
|
|
Packit Service |
2e9770 |
Return 0 if special formula fails, otherwise put in z1 the approximate
|
|
Packit Service |
2e9770 |
value which needs to be converted to rop.
|
|
Packit Service |
2e9770 |
z1 is a temporary variable with enough precision.
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
static int
|
|
Packit Service |
2e9770 |
mpc_asin_special (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd, mpc_ptr z1)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_exp_t ey = mpfr_get_exp (mpc_imagref (op));
|
|
Packit Service |
2e9770 |
mpfr_t abs_y;
|
|
Packit Service |
2e9770 |
mpfr_prec_t p;
|
|
Packit Service |
2e9770 |
int inex;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* |Re(asin(1+i*y)) - pi/2| <= y^(1/2) */
|
|
Packit Service |
2e9770 |
if (ey >= 0 || ((-ey) / 2 < mpfr_get_prec (mpc_realref (z1))))
|
|
Packit Service |
2e9770 |
return 0;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_const_pi (mpc_realref (z1), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
mpfr_div_2exp (mpc_realref (z1), mpc_realref (z1), 1, MPFR_RNDN); /* exact */
|
|
Packit Service |
2e9770 |
p = mpfr_get_prec (mpc_realref (z1));
|
|
Packit Service |
2e9770 |
/* if z1 has precision p, the error on z1 is 1/2*ulp(z1) = 2^(-p) so far,
|
|
Packit Service |
2e9770 |
and since ey <= -2p, then y^(1/2) <= 1/2*ulp(z1) too, thus the total
|
|
Packit Service |
2e9770 |
error is bounded by ulp(z1) */
|
|
Packit Service |
2e9770 |
if (!mpfr_can_round (mpc_realref(z1), p, MPFR_RNDN, MPFR_RNDZ,
|
|
Packit Service |
2e9770 |
mpfr_get_prec (mpc_realref(rop)) +
|
|
Packit Service |
2e9770 |
(MPC_RND_RE(rnd) == MPFR_RNDN)))
|
|
Packit Service |
2e9770 |
return 0;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* |Im(asin(1+i*y)) - y^(1/2)| <= (1/12) * y^(3/2) for y >= 0 (err >= 0)
|
|
Packit Service |
2e9770 |
|Im(asin(1-i*y)) + y^(1/2)| <= (1/12) * y^(3/2) for y >= 0 (err <= 0) */
|
|
Packit Service |
2e9770 |
abs_y[0] = mpc_imagref (op)[0];
|
|
Packit Service |
2e9770 |
if (mpfr_signbit (mpc_imagref (op)))
|
|
Packit Service |
2e9770 |
MPFR_CHANGE_SIGN (abs_y);
|
|
Packit Service |
2e9770 |
inex = mpfr_sqrt (mpc_imagref (z1), abs_y, MPFR_RNDN); /* error <= 1/2 ulp */
|
|
Packit Service |
2e9770 |
if (mpfr_signbit (mpc_imagref (op)))
|
|
Packit Service |
2e9770 |
MPFR_CHANGE_SIGN (mpc_imagref (z1));
|
|
Packit Service |
2e9770 |
/* If z1 has precision p, the error on z1 is 1/2*ulp(z1) = 2^(-p) so far,
|
|
Packit Service |
2e9770 |
and (1/12) * y^(3/2) <= (1/8) * y * y^(1/2) <= 2^(ey-3)*2^p*ulp(y^(1/2))
|
|
Packit Service |
2e9770 |
thus for p+ey-3 <= -1 we have (1/12) * y^(3/2) <= (1/2) * ulp(y^(1/2)),
|
|
Packit Service |
2e9770 |
and the total error is bounded by ulp(z1).
|
|
Packit Service |
2e9770 |
Note: if y^(1/2) is exactly representable, and ends with many zeroes,
|
|
Packit Service |
2e9770 |
then mpfr_can_round below will fail; however in that case we know that
|
|
Packit Service |
2e9770 |
Im(asin(1+i*y)) is away from +/-y^(1/2) wrt zero. */
|
|
Packit Service |
2e9770 |
if (inex == 0) /* enlarge im(z1) so that the inexact flag is correct */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
if (mpfr_signbit (mpc_imagref (op)))
|
|
Packit Service |
2e9770 |
mpfr_nextbelow (mpc_imagref (z1));
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
mpfr_nextabove (mpc_imagref (z1));
|
|
Packit Service |
2e9770 |
return 1;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
p = mpfr_get_prec (mpc_imagref (z1));
|
|
Packit Service |
2e9770 |
if (!mpfr_can_round (mpc_imagref(z1), p - 1, MPFR_RNDA, MPFR_RNDZ,
|
|
Packit Service |
2e9770 |
mpfr_get_prec (mpc_imagref(rop)) +
|
|
Packit Service |
2e9770 |
(MPC_RND_IM(rnd) == MPFR_RNDN)))
|
|
Packit Service |
2e9770 |
return 0;
|
|
Packit Service |
2e9770 |
return 1;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
int
|
|
Packit Service |
2e9770 |
mpc_asin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_prec_t p, p_re, p_im;
|
|
Packit Service |
2e9770 |
mpfr_rnd_t rnd_re, rnd_im;
|
|
Packit Service |
2e9770 |
mpc_t z1;
|
|
Packit Service |
2e9770 |
int inex, loop = 0;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* special values */
|
|
Packit Service |
2e9770 |
if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op)))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op)))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_set_nan (mpc_realref (rop));
|
|
Packit Service |
2e9770 |
mpfr_set_inf (mpc_imagref (rop), mpfr_signbit (mpc_imagref (op)) ? -1 : +1);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else if (mpfr_zero_p (mpc_realref (op)))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_set (mpc_realref (rop), mpc_realref (op), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
mpfr_set_nan (mpc_imagref (rop));
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_set_nan (mpc_realref (rop));
|
|
Packit Service |
2e9770 |
mpfr_set_nan (mpc_imagref (rop));
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
return 0;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op)))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int inex_re;
|
|
Packit Service |
2e9770 |
if (mpfr_inf_p (mpc_realref (op)))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int inf_im = mpfr_inf_p (mpc_imagref (op));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
inex_re = set_pi_over_2 (mpc_realref (rop),
|
|
Packit Service |
2e9770 |
(mpfr_signbit (mpc_realref (op)) ? -1 : 1), MPC_RND_RE (rnd));
|
|
Packit Service |
2e9770 |
mpfr_set_inf (mpc_imagref (rop), (mpfr_signbit (mpc_imagref (op)) ? -1 : 1));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (inf_im)
|
|
Packit Service |
2e9770 |
mpfr_div_2ui (mpc_realref (rop), mpc_realref (rop), 1, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_set_zero (mpc_realref (rop), (mpfr_signbit (mpc_realref (op)) ? -1 : 1));
|
|
Packit Service |
2e9770 |
inex_re = 0;
|
|
Packit Service |
2e9770 |
mpfr_set_inf (mpc_imagref (rop), (mpfr_signbit (mpc_imagref (op)) ? -1 : 1));
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
return MPC_INEX (inex_re, 0);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* pure real argument */
|
|
Packit Service |
2e9770 |
if (mpfr_zero_p (mpc_imagref (op)))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int inex_re;
|
|
Packit Service |
2e9770 |
int inex_im;
|
|
Packit Service |
2e9770 |
int s_im;
|
|
Packit Service |
2e9770 |
s_im = mpfr_signbit (mpc_imagref (op));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (mpfr_cmp_ui (mpc_realref (op), 1) > 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
if (s_im)
|
|
Packit Service |
2e9770 |
inex_im = -mpfr_acosh (mpc_imagref (rop), mpc_realref (op),
|
|
Packit Service |
2e9770 |
INV_RND (MPC_RND_IM (rnd)));
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
inex_im = mpfr_acosh (mpc_imagref (rop), mpc_realref (op),
|
|
Packit Service |
2e9770 |
MPC_RND_IM (rnd));
|
|
Packit Service |
2e9770 |
inex_re = set_pi_over_2 (mpc_realref (rop),
|
|
Packit Service |
2e9770 |
(mpfr_signbit (mpc_realref (op)) ? -1 : 1), MPC_RND_RE (rnd));
|
|
Packit Service |
2e9770 |
if (s_im)
|
|
Packit Service |
2e9770 |
mpc_conj (rop, rop, MPC_RNDNN);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else if (mpfr_cmp_si (mpc_realref (op), -1) < 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_t minus_op_re;
|
|
Packit Service |
2e9770 |
minus_op_re[0] = mpc_realref (op)[0];
|
|
Packit Service |
2e9770 |
MPFR_CHANGE_SIGN (minus_op_re);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (s_im)
|
|
Packit Service |
2e9770 |
inex_im = -mpfr_acosh (mpc_imagref (rop), minus_op_re,
|
|
Packit Service |
2e9770 |
INV_RND (MPC_RND_IM (rnd)));
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
inex_im = mpfr_acosh (mpc_imagref (rop), minus_op_re,
|
|
Packit Service |
2e9770 |
MPC_RND_IM (rnd));
|
|
Packit Service |
2e9770 |
inex_re = set_pi_over_2 (mpc_realref (rop),
|
|
Packit Service |
2e9770 |
(mpfr_signbit (mpc_realref (op)) ? -1 : 1), MPC_RND_RE (rnd));
|
|
Packit Service |
2e9770 |
if (s_im)
|
|
Packit Service |
2e9770 |
mpc_conj (rop, rop, MPC_RNDNN);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
inex_im = mpfr_set_ui (mpc_imagref (rop), 0, MPC_RND_IM (rnd));
|
|
Packit Service |
2e9770 |
if (s_im)
|
|
Packit Service |
2e9770 |
mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
inex_re = mpfr_asin (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd));
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
return MPC_INEX (inex_re, inex_im);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* pure imaginary argument */
|
|
Packit Service |
2e9770 |
if (mpfr_zero_p (mpc_realref (op)))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int inex_im;
|
|
Packit Service |
2e9770 |
int s;
|
|
Packit Service |
2e9770 |
s = mpfr_signbit (mpc_realref (op));
|
|
Packit Service |
2e9770 |
mpfr_set_ui (mpc_realref (rop), 0, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
if (s)
|
|
Packit Service |
2e9770 |
mpfr_neg (mpc_realref (rop), mpc_realref (rop), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
inex_im = mpfr_asinh (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
return MPC_INEX (0, inex_im);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* regular complex: asin(z) = -i*log(i*z+sqrt(1-z^2)) */
|
|
Packit Service |
2e9770 |
p_re = mpfr_get_prec (mpc_realref(rop));
|
|
Packit Service |
2e9770 |
p_im = mpfr_get_prec (mpc_imagref(rop));
|
|
Packit Service |
2e9770 |
rnd_re = MPC_RND_RE(rnd);
|
|
Packit Service |
2e9770 |
rnd_im = MPC_RND_IM(rnd);
|
|
Packit Service |
2e9770 |
p = p_re >= p_im ? p_re : p_im;
|
|
Packit Service |
2e9770 |
mpc_init2 (z1, p);
|
|
Packit Service |
2e9770 |
while (1)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_exp_t ex, ey, err;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
loop ++;
|
|
Packit Service |
2e9770 |
p += (loop <= 2) ? mpc_ceil_log2 (p) + 3 : p / 2;
|
|
Packit Service |
2e9770 |
mpfr_set_prec (mpc_realref(z1), p);
|
|
Packit Service |
2e9770 |
mpfr_set_prec (mpc_imagref(z1), p);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* try special code for 1+i*y with tiny y */
|
|
Packit Service |
2e9770 |
if (loop == 1 && mpc_asin_special (rop, op, rnd, z1))
|
|
Packit Service |
2e9770 |
break;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* z1 <- z^2 */
|
|
Packit Service |
2e9770 |
mpc_sqr (z1, op, MPC_RNDNN);
|
|
Packit Service |
2e9770 |
/* err(x) <= 1/2 ulp(x), err(y) <= 1/2 ulp(y) */
|
|
Packit Service |
2e9770 |
/* z1 <- 1-z1 */
|
|
Packit Service |
2e9770 |
ex = mpfr_get_exp (mpc_realref(z1));
|
|
Packit Service |
2e9770 |
mpfr_ui_sub (mpc_realref(z1), 1, mpc_realref(z1), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
ex = ex - mpfr_get_exp (mpc_realref(z1));
|
|
Packit Service |
2e9770 |
ex = (ex <= 0) ? 0 : ex;
|
|
Packit Service |
2e9770 |
/* err(x) <= 2^ex * ulp(x) */
|
|
Packit Service |
2e9770 |
ex = ex + mpfr_get_exp (mpc_realref(z1)) - p;
|
|
Packit Service |
2e9770 |
/* err(x) <= 2^ex */
|
|
Packit Service |
2e9770 |
ey = mpfr_get_exp (mpc_imagref(z1)) - p - 1;
|
|
Packit Service |
2e9770 |
/* err(y) <= 2^ey */
|
|
Packit Service |
2e9770 |
ex = (ex >= ey) ? ex : ey; /* err(x), err(y) <= 2^ex, i.e., the norm
|
|
Packit Service |
2e9770 |
of the error is bounded by |h|<=2^(ex+1/2) */
|
|
Packit Service |
2e9770 |
/* z1 <- sqrt(z1): if z1 = z + h, then sqrt(z1) = sqrt(z) + h/2/sqrt(t) */
|
|
Packit Service |
2e9770 |
ey = mpfr_get_exp (mpc_realref(z1)) >= mpfr_get_exp (mpc_imagref(z1))
|
|
Packit Service |
2e9770 |
? mpfr_get_exp (mpc_realref(z1)) : mpfr_get_exp (mpc_imagref(z1));
|
|
Packit Service |
2e9770 |
/* we have |z1| >= 2^(ey-1) thus 1/|z1| <= 2^(1-ey) */
|
|
Packit Service |
2e9770 |
mpc_sqrt (z1, z1, MPC_RNDNN);
|
|
Packit Service |
2e9770 |
ex = (2 * ex + 1) - 2 - (ey - 1); /* |h^2/4/|t| <= 2^ex */
|
|
Packit Service |
2e9770 |
ex = (ex + 1) / 2; /* ceil(ex/2) */
|
|
Packit Service |
2e9770 |
/* express ex in terms of ulp(z1) */
|
|
Packit Service |
2e9770 |
ey = mpfr_get_exp (mpc_realref(z1)) <= mpfr_get_exp (mpc_imagref(z1))
|
|
Packit Service |
2e9770 |
? mpfr_get_exp (mpc_realref(z1)) : mpfr_get_exp (mpc_imagref(z1));
|
|
Packit Service |
2e9770 |
ex = ex - ey + p;
|
|
Packit Service |
2e9770 |
/* take into account the rounding error in the mpc_sqrt call */
|
|
Packit Service |
2e9770 |
err = (ex <= 0) ? 1 : ex + 1;
|
|
Packit Service |
2e9770 |
/* err(x) <= 2^err * ulp(x), err(y) <= 2^err * ulp(y) */
|
|
Packit Service |
2e9770 |
/* z1 <- i*z + z1 */
|
|
Packit Service |
2e9770 |
ex = mpfr_get_exp (mpc_realref(z1));
|
|
Packit Service |
2e9770 |
ey = mpfr_get_exp (mpc_imagref(z1));
|
|
Packit Service |
2e9770 |
mpfr_sub (mpc_realref(z1), mpc_realref(z1), mpc_imagref(op), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
mpfr_add (mpc_imagref(z1), mpc_imagref(z1), mpc_realref(op), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
if (mpfr_cmp_ui (mpc_realref(z1), 0) == 0 || mpfr_cmp_ui (mpc_imagref(z1), 0) == 0)
|
|
Packit Service |
2e9770 |
continue;
|
|
Packit Service |
2e9770 |
ex -= mpfr_get_exp (mpc_realref(z1)); /* cancellation in x */
|
|
Packit Service |
2e9770 |
ey -= mpfr_get_exp (mpc_imagref(z1)); /* cancellation in y */
|
|
Packit Service |
2e9770 |
ex = (ex >= ey) ? ex : ey; /* maximum cancellation */
|
|
Packit Service |
2e9770 |
err += ex;
|
|
Packit Service |
2e9770 |
err = (err <= 0) ? 1 : err + 1; /* rounding error in sub/add */
|
|
Packit Service |
2e9770 |
/* z1 <- log(z1): if z1 = z + h, then log(z1) = log(z) + h/t with
|
|
Packit Service |
2e9770 |
|t| >= min(|z1|,|z|) */
|
|
Packit Service |
2e9770 |
ex = mpfr_get_exp (mpc_realref(z1));
|
|
Packit Service |
2e9770 |
ey = mpfr_get_exp (mpc_imagref(z1));
|
|
Packit Service |
2e9770 |
ex = (ex >= ey) ? ex : ey;
|
|
Packit Service |
2e9770 |
err += ex - p; /* revert to absolute error <= 2^err */
|
|
Packit Service |
2e9770 |
mpc_log (z1, z1, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
err -= ex - 1; /* 1/|t| <= 1/|z| <= 2^(1-ex) */
|
|
Packit Service |
2e9770 |
/* express err in terms of ulp(z1) */
|
|
Packit Service |
2e9770 |
ey = mpfr_get_exp (mpc_realref(z1)) <= mpfr_get_exp (mpc_imagref(z1))
|
|
Packit Service |
2e9770 |
? mpfr_get_exp (mpc_realref(z1)) : mpfr_get_exp (mpc_imagref(z1));
|
|
Packit Service |
2e9770 |
err = err - ey + p;
|
|
Packit Service |
2e9770 |
/* take into account the rounding error in the mpc_log call */
|
|
Packit Service |
2e9770 |
err = (err <= 0) ? 1 : err + 1;
|
|
Packit Service |
2e9770 |
/* z1 <- -i*z1 */
|
|
Packit Service |
2e9770 |
mpfr_swap (mpc_realref(z1), mpc_imagref(z1));
|
|
Packit Service |
2e9770 |
mpfr_neg (mpc_imagref(z1), mpc_imagref(z1), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
if (mpfr_can_round (mpc_realref(z1), p - err, MPFR_RNDN, MPFR_RNDZ,
|
|
Packit Service |
2e9770 |
p_re + (rnd_re == MPFR_RNDN)) &&
|
|
Packit Service |
2e9770 |
mpfr_can_round (mpc_imagref(z1), p - err, MPFR_RNDN, MPFR_RNDZ,
|
|
Packit Service |
2e9770 |
p_im + (rnd_im == MPFR_RNDN)))
|
|
Packit Service |
2e9770 |
break;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
inex = mpc_set (rop, z1, rnd);
|
|
Packit Service |
2e9770 |
mpc_clear (z1);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
return inex;
|
|
Packit Service |
2e9770 |
}
|