|
Packit Service |
2e9770 |
/* mpc_tan -- tangent of a complex number.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
Copyright (C) 2008-2015 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> /* for MPC_ASSERT */
|
|
Packit Service |
2e9770 |
#include <limits.h>
|
|
Packit Service |
2e9770 |
#include "mpc-impl.h"
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
int
|
|
Packit Service |
2e9770 |
mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpc_t x, y;
|
|
Packit Service |
2e9770 |
mpfr_prec_t prec;
|
|
Packit Service |
2e9770 |
mpfr_exp_t err;
|
|
Packit Service |
2e9770 |
int ok = 0;
|
|
Packit Service |
2e9770 |
int inex;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* special values */
|
|
Packit Service |
2e9770 |
if (!mpc_fin_p (op))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
if (mpfr_nan_p (mpc_realref (op)))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
if (mpfr_inf_p (mpc_imagref (op)))
|
|
Packit Service |
2e9770 |
/* tan(NaN -i*Inf) = +/-0 -i */
|
|
Packit Service |
2e9770 |
/* tan(NaN +i*Inf) = +/-0 +i */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
/* exact unless 1 is not in exponent range */
|
|
Packit Service |
2e9770 |
inex = mpc_set_si_si (rop, 0,
|
|
Packit Service |
2e9770 |
(MPFR_SIGN (mpc_imagref (op)) < 0) ? -1 : +1,
|
|
Packit Service |
2e9770 |
rnd);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
/* tan(NaN +i*y) = NaN +i*NaN, when y is finite */
|
|
Packit Service |
2e9770 |
/* tan(NaN +i*NaN) = NaN +i*NaN */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_set_nan (mpc_realref (rop));
|
|
Packit Service |
2e9770 |
mpfr_set_nan (mpc_imagref (rop));
|
|
Packit Service |
2e9770 |
inex = MPC_INEX (0, 0); /* always exact */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else if (mpfr_nan_p (mpc_imagref (op)))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
if (mpfr_cmp_ui (mpc_realref (op), 0) == 0)
|
|
Packit Service |
2e9770 |
/* tan(-0 +i*NaN) = -0 +i*NaN */
|
|
Packit Service |
2e9770 |
/* tan(+0 +i*NaN) = +0 +i*NaN */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpc_set (rop, op, rnd);
|
|
Packit Service |
2e9770 |
inex = MPC_INEX (0, 0); /* always exact */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
/* tan(x +i*NaN) = NaN +i*NaN, when x != 0 */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_set_nan (mpc_realref (rop));
|
|
Packit Service |
2e9770 |
mpfr_set_nan (mpc_imagref (rop));
|
|
Packit Service |
2e9770 |
inex = MPC_INEX (0, 0); /* always exact */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else if (mpfr_inf_p (mpc_realref (op)))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
if (mpfr_inf_p (mpc_imagref (op)))
|
|
Packit Service |
2e9770 |
/* tan(-Inf -i*Inf) = -/+0 -i */
|
|
Packit Service |
2e9770 |
/* tan(-Inf +i*Inf) = -/+0 +i */
|
|
Packit Service |
2e9770 |
/* tan(+Inf -i*Inf) = +/-0 -i */
|
|
Packit Service |
2e9770 |
/* tan(+Inf +i*Inf) = +/-0 +i */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
const int sign_re = mpfr_signbit (mpc_realref (op));
|
|
Packit Service |
2e9770 |
int inex_im;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd));
|
|
Packit Service |
2e9770 |
mpfr_setsign (mpc_realref (rop), mpc_realref (rop), sign_re, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* exact, unless 1 is not in exponent range */
|
|
Packit Service |
2e9770 |
inex_im = mpfr_set_si (mpc_imagref (rop),
|
|
Packit Service |
2e9770 |
mpfr_signbit (mpc_imagref (op)) ? -1 : +1,
|
|
Packit Service |
2e9770 |
MPC_RND_IM (rnd));
|
|
Packit Service |
2e9770 |
inex = MPC_INEX (0, inex_im);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
/* tan(-Inf +i*y) = tan(+Inf +i*y) = NaN +i*NaN, when y is
|
|
Packit Service |
2e9770 |
finite */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_set_nan (mpc_realref (rop));
|
|
Packit Service |
2e9770 |
mpfr_set_nan (mpc_imagref (rop));
|
|
Packit Service |
2e9770 |
inex = MPC_INEX (0, 0); /* always exact */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
/* tan(x -i*Inf) = +0*sin(x)*cos(x) -i, when x is finite */
|
|
Packit Service |
2e9770 |
/* tan(x +i*Inf) = +0*sin(x)*cos(x) +i, when x is finite */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_t c;
|
|
Packit Service |
2e9770 |
mpfr_t s;
|
|
Packit Service |
2e9770 |
int inex_im;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_init (c);
|
|
Packit Service |
2e9770 |
mpfr_init (s);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_sin_cos (s, c, mpc_realref (op), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd));
|
|
Packit Service |
2e9770 |
mpfr_setsign (mpc_realref (rop), mpc_realref (rop),
|
|
Packit Service |
2e9770 |
mpfr_signbit (c) != mpfr_signbit (s), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
/* exact, unless 1 is not in exponent range */
|
|
Packit Service |
2e9770 |
inex_im = mpfr_set_si (mpc_imagref (rop),
|
|
Packit Service |
2e9770 |
(mpfr_signbit (mpc_imagref (op)) ? -1 : +1),
|
|
Packit Service |
2e9770 |
MPC_RND_IM (rnd));
|
|
Packit Service |
2e9770 |
inex = MPC_INEX (0, inex_im);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_clear (s);
|
|
Packit Service |
2e9770 |
mpfr_clear (c);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
return inex;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (mpfr_zero_p (mpc_realref (op)))
|
|
Packit Service |
2e9770 |
/* tan(-0 -i*y) = -0 +i*tanh(y), when y is finite. */
|
|
Packit Service |
2e9770 |
/* tan(+0 +i*y) = +0 +i*tanh(y), when y is finite. */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int inex_im;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_set (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd));
|
|
Packit Service |
2e9770 |
inex_im = mpfr_tanh (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 |
if (mpfr_zero_p (mpc_imagref (op)))
|
|
Packit Service |
2e9770 |
/* tan(x -i*0) = tan(x) -i*0, when x is finite. */
|
|
Packit Service |
2e9770 |
/* tan(x +i*0) = tan(x) +i*0, when x is finite. */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int inex_re;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
inex_re = mpfr_tan (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd));
|
|
Packit Service |
2e9770 |
mpfr_set (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
return MPC_INEX (inex_re, 0);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* ordinary (non-zero) numbers */
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* tan(op) = sin(op) / cos(op).
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
We use the following algorithm with rounding away from 0 for all
|
|
Packit Service |
2e9770 |
operations, and working precision w:
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
(1) x = A(sin(op))
|
|
Packit Service |
2e9770 |
(2) y = A(cos(op))
|
|
Packit Service |
2e9770 |
(3) z = A(x/y)
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
the error on Im(z) is at most 81 ulp,
|
|
Packit Service |
2e9770 |
the error on Re(z) is at most
|
|
Packit Service |
2e9770 |
7 ulp if k < 2,
|
|
Packit Service |
2e9770 |
8 ulp if k = 2,
|
|
Packit Service |
2e9770 |
else 5+k ulp, where
|
|
Packit Service |
2e9770 |
k = Exp(Re(x))+Exp(Re(y))-2min{Exp(Re(y)), Exp(Im(y))}-Exp(Re(x/y))
|
|
Packit Service |
2e9770 |
see proof in algorithms.tex.
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
prec = MPC_MAX_PREC(rop);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpc_init2 (x, 2);
|
|
Packit Service |
2e9770 |
mpc_init2 (y, 2);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
err = 7;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
do
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_exp_t k, exr, eyr, eyi, ezr;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
ok = 0;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* FIXME: prevent addition overflow */
|
|
Packit Service |
2e9770 |
prec += mpc_ceil_log2 (prec) + err;
|
|
Packit Service |
2e9770 |
mpc_set_prec (x, prec);
|
|
Packit Service |
2e9770 |
mpc_set_prec (y, prec);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* rounding away from zero: except in the cases x=0 or y=0 (processed
|
|
Packit Service |
2e9770 |
above), sin x and cos y are never exact, so rounding away from 0 is
|
|
Packit Service |
2e9770 |
rounding towards 0 and adding one ulp to the absolute value */
|
|
Packit Service |
2e9770 |
mpc_sin_cos (x, y, op, MPC_RNDZZ, MPC_RNDZZ);
|
|
Packit Service |
2e9770 |
MPFR_ADD_ONE_ULP (mpc_realref (x));
|
|
Packit Service |
2e9770 |
MPFR_ADD_ONE_ULP (mpc_imagref (x));
|
|
Packit Service |
2e9770 |
MPFR_ADD_ONE_ULP (mpc_realref (y));
|
|
Packit Service |
2e9770 |
MPFR_ADD_ONE_ULP (mpc_imagref (y));
|
|
Packit Service |
2e9770 |
MPC_ASSERT (mpfr_zero_p (mpc_realref (x)) == 0);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if ( mpfr_inf_p (mpc_realref (x)) || mpfr_inf_p (mpc_imagref (x))
|
|
Packit Service |
2e9770 |
|| mpfr_inf_p (mpc_realref (y)) || mpfr_inf_p (mpc_imagref (y))) {
|
|
Packit Service |
2e9770 |
/* If the real or imaginary part of x is infinite, it means that
|
|
Packit Service |
2e9770 |
Im(op) was large, in which case the result is
|
|
Packit Service |
2e9770 |
sign(tan(Re(op)))*0 + sign(Im(op))*I,
|
|
Packit Service |
2e9770 |
where sign(tan(Re(op))) = sign(Re(x))*sign(Re(y)). */
|
|
Packit Service |
2e9770 |
int inex_re, inex_im;
|
|
Packit Service |
2e9770 |
mpfr_set_ui (mpc_realref (rop), 0, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
if (mpfr_sgn (mpc_realref (x)) * mpfr_sgn (mpc_realref (y)) < 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_neg (mpc_realref (rop), mpc_realref (rop), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
inex_re = 1;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
inex_re = -1; /* +0 is rounded down */
|
|
Packit Service |
2e9770 |
if (mpfr_sgn (mpc_imagref (op)) > 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_set_ui (mpc_imagref (rop), 1, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
inex_im = 1;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_set_si (mpc_imagref (rop), -1, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
inex_im = -1;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
/* if rounding is toward zero, fix the imaginary part */
|
|
Packit Service |
2e9770 |
if (MPC_IS_LIKE_RNDZ(MPC_RND_IM(rnd), MPFR_SIGNBIT(mpc_imagref (rop))))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_nexttoward (mpc_imagref (rop), mpc_realref (rop));
|
|
Packit Service |
2e9770 |
inex_im = -inex_im;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
if (mpfr_zero_p (mpc_realref (rop)))
|
|
Packit Service |
2e9770 |
inex_re = mpc_fix_zero (mpc_realref (rop), MPC_RND_RE(rnd));
|
|
Packit Service |
2e9770 |
if (mpfr_zero_p (mpc_imagref (rop)))
|
|
Packit Service |
2e9770 |
inex_im = mpc_fix_zero (mpc_imagref (rop), MPC_RND_IM(rnd));
|
|
Packit Service |
2e9770 |
inex = MPC_INEX(inex_re, inex_im);
|
|
Packit Service |
2e9770 |
goto end;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
exr = mpfr_get_exp (mpc_realref (x));
|
|
Packit Service |
2e9770 |
eyr = mpfr_get_exp (mpc_realref (y));
|
|
Packit Service |
2e9770 |
eyi = mpfr_get_exp (mpc_imagref (y));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* some parts of the quotient may be exact */
|
|
Packit Service |
2e9770 |
inex = mpc_div (x, x, y, MPC_RNDZZ);
|
|
Packit Service |
2e9770 |
/* OP is no pure real nor pure imaginary, so in theory the real and
|
|
Packit Service |
2e9770 |
imaginary parts of its tangent cannot be null. However due to
|
|
Packit Service |
2e9770 |
rounding errors this might happen. Consider for example
|
|
Packit Service |
2e9770 |
tan(1+14*I) = 1.26e-10 + 1.00*I. For small precision sin(op) and
|
|
Packit Service |
2e9770 |
cos(op) differ only by a factor I, thus after mpc_div x = I and
|
|
Packit Service |
2e9770 |
its real part is zero. */
|
|
Packit Service |
2e9770 |
if (mpfr_zero_p (mpc_realref (x)) || mpfr_zero_p (mpc_imagref (x)))
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
err = prec; /* double precision */
|
|
Packit Service |
2e9770 |
continue;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
if (MPC_INEX_RE (inex))
|
|
Packit Service |
2e9770 |
MPFR_ADD_ONE_ULP (mpc_realref (x));
|
|
Packit Service |
2e9770 |
if (MPC_INEX_IM (inex))
|
|
Packit Service |
2e9770 |
MPFR_ADD_ONE_ULP (mpc_imagref (x));
|
|
Packit Service |
2e9770 |
MPC_ASSERT (mpfr_zero_p (mpc_realref (x)) == 0);
|
|
Packit Service |
2e9770 |
ezr = mpfr_get_exp (mpc_realref (x));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* FIXME: compute
|
|
Packit Service |
2e9770 |
k = Exp(Re(x))+Exp(Re(y))-2min{Exp(Re(y)), Exp(Im(y))}-Exp(Re(x/y))
|
|
Packit Service |
2e9770 |
avoiding overflow */
|
|
Packit Service |
2e9770 |
k = exr - ezr + MPC_MAX(-eyr, eyr - 2 * eyi);
|
|
Packit Service |
2e9770 |
err = k < 2 ? 7 : (k == 2 ? 8 : (5 + k));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* Can the real part be rounded? */
|
|
Packit Service |
2e9770 |
ok = (!mpfr_number_p (mpc_realref (x)))
|
|
Packit Service |
2e9770 |
|| mpfr_can_round (mpc_realref(x), prec - err, MPFR_RNDN, MPFR_RNDZ,
|
|
Packit Service |
2e9770 |
MPC_PREC_RE(rop) + (MPC_RND_RE(rnd) == MPFR_RNDN));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (ok)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
/* Can the imaginary part be rounded? */
|
|
Packit Service |
2e9770 |
ok = (!mpfr_number_p (mpc_imagref (x)))
|
|
Packit Service |
2e9770 |
|| mpfr_can_round (mpc_imagref(x), prec - 6, MPFR_RNDN, MPFR_RNDZ,
|
|
Packit Service |
2e9770 |
MPC_PREC_IM(rop) + (MPC_RND_IM(rnd) == MPFR_RNDN));
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
while (ok == 0);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
inex = mpc_set (rop, x, rnd);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
end:
|
|
Packit Service |
2e9770 |
mpc_clear (x);
|
|
Packit Service |
2e9770 |
mpc_clear (y);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
return inex;
|
|
Packit Service |
2e9770 |
}
|