1 /* mpc_tan -- tangent of a complex number. 2 3 Copyright (C) 2008, 2009, 2010, 2011 INRIA 4 5 This file is part of GNU MPC. 6 7 GNU MPC is free software; you can redistribute it and/or modify it under 8 the terms of the GNU Lesser General Public License as published by the 9 Free Software Foundation; either version 3 of the License, or (at your 10 option) any later version. 11 12 GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY 13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for 15 more details. 16 17 You should have received a copy of the GNU Lesser General Public License 18 along with this program. If not, see http://www.gnu.org/licenses/ . 19 */ 20 21 #include <stdio.h> /* for MPC_ASSERT */ 22 #include <limits.h> 23 #include "mpc-impl.h" 24 25 int 26 mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) 27 { 28 mpc_t x, y; 29 mpfr_prec_t prec; 30 mpfr_exp_t err; 31 int ok = 0; 32 int inex; 33 34 /* special values */ 35 if (!mpc_fin_p (op)) 36 { 37 if (mpfr_nan_p (mpc_realref (op))) 38 { 39 if (mpfr_inf_p (mpc_imagref (op))) 40 /* tan(NaN -i*Inf) = +/-0 -i */ 41 /* tan(NaN +i*Inf) = +/-0 +i */ 42 { 43 /* exact unless 1 is not in exponent range */ 44 inex = mpc_set_si_si (rop, 0, 45 (MPFR_SIGN (mpc_imagref (op)) < 0) ? -1 : +1, 46 rnd); 47 } 48 else 49 /* tan(NaN +i*y) = NaN +i*NaN, when y is finite */ 50 /* tan(NaN +i*NaN) = NaN +i*NaN */ 51 { 52 mpfr_set_nan (mpc_realref (rop)); 53 mpfr_set_nan (mpc_imagref (rop)); 54 inex = MPC_INEX (0, 0); /* always exact */ 55 } 56 } 57 else if (mpfr_nan_p (mpc_imagref (op))) 58 { 59 if (mpfr_cmp_ui (mpc_realref (op), 0) == 0) 60 /* tan(-0 +i*NaN) = -0 +i*NaN */ 61 /* tan(+0 +i*NaN) = +0 +i*NaN */ 62 { 63 mpc_set (rop, op, rnd); 64 inex = MPC_INEX (0, 0); /* always exact */ 65 } 66 else 67 /* tan(x +i*NaN) = NaN +i*NaN, when x != 0 */ 68 { 69 mpfr_set_nan (mpc_realref (rop)); 70 mpfr_set_nan (mpc_imagref (rop)); 71 inex = MPC_INEX (0, 0); /* always exact */ 72 } 73 } 74 else if (mpfr_inf_p (mpc_realref (op))) 75 { 76 if (mpfr_inf_p (mpc_imagref (op))) 77 /* tan(-Inf -i*Inf) = -/+0 -i */ 78 /* tan(-Inf +i*Inf) = -/+0 +i */ 79 /* tan(+Inf -i*Inf) = +/-0 -i */ 80 /* tan(+Inf +i*Inf) = +/-0 +i */ 81 { 82 const int sign_re = mpfr_signbit (mpc_realref (op)); 83 int inex_im; 84 85 mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd)); 86 mpfr_setsign (mpc_realref (rop), mpc_realref (rop), sign_re, GMP_RNDN); 87 88 /* exact, unless 1 is not in exponent range */ 89 inex_im = mpfr_set_si (mpc_imagref (rop), 90 mpfr_signbit (mpc_imagref (op)) ? -1 : +1, 91 MPC_RND_IM (rnd)); 92 inex = MPC_INEX (0, inex_im); 93 } 94 else 95 /* tan(-Inf +i*y) = tan(+Inf +i*y) = NaN +i*NaN, when y is 96 finite */ 97 { 98 mpfr_set_nan (mpc_realref (rop)); 99 mpfr_set_nan (mpc_imagref (rop)); 100 inex = MPC_INEX (0, 0); /* always exact */ 101 } 102 } 103 else 104 /* tan(x -i*Inf) = +0*sin(x)*cos(x) -i, when x is finite */ 105 /* tan(x +i*Inf) = +0*sin(x)*cos(x) +i, when x is finite */ 106 { 107 mpfr_t c; 108 mpfr_t s; 109 int inex_im; 110 111 mpfr_init (c); 112 mpfr_init (s); 113 114 mpfr_sin_cos (s, c, mpc_realref (op), GMP_RNDN); 115 mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd)); 116 mpfr_setsign (mpc_realref (rop), mpc_realref (rop), 117 mpfr_signbit (c) != mpfr_signbit (s), GMP_RNDN); 118 /* exact, unless 1 is not in exponent range */ 119 inex_im = mpfr_set_si (mpc_imagref (rop), 120 (mpfr_signbit (mpc_imagref (op)) ? -1 : +1), 121 MPC_RND_IM (rnd)); 122 inex = MPC_INEX (0, inex_im); 123 124 mpfr_clear (s); 125 mpfr_clear (c); 126 } 127 128 return inex; 129 } 130 131 if (mpfr_zero_p (mpc_realref (op))) 132 /* tan(-0 -i*y) = -0 +i*tanh(y), when y is finite. */ 133 /* tan(+0 +i*y) = +0 +i*tanh(y), when y is finite. */ 134 { 135 int inex_im; 136 137 mpfr_set (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd)); 138 inex_im = mpfr_tanh (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd)); 139 140 return MPC_INEX (0, inex_im); 141 } 142 143 if (mpfr_zero_p (mpc_imagref (op))) 144 /* tan(x -i*0) = tan(x) -i*0, when x is finite. */ 145 /* tan(x +i*0) = tan(x) +i*0, when x is finite. */ 146 { 147 int inex_re; 148 149 inex_re = mpfr_tan (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd)); 150 mpfr_set (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd)); 151 152 return MPC_INEX (inex_re, 0); 153 } 154 155 /* ordinary (non-zero) numbers */ 156 157 /* tan(op) = sin(op) / cos(op). 158 159 We use the following algorithm with rounding away from 0 for all 160 operations, and working precision w: 161 162 (1) x = A(sin(op)) 163 (2) y = A(cos(op)) 164 (3) z = A(x/y) 165 166 the error on Im(z) is at most 81 ulp, 167 the error on Re(z) is at most 168 7 ulp if k < 2, 169 8 ulp if k = 2, 170 else 5+k ulp, where 171 k = Exp(Re(x))+Exp(Re(y))-2min{Exp(Re(y)), Exp(Im(y))}-Exp(Re(x/y)) 172 see proof in algorithms.tex. 173 */ 174 175 prec = MPC_MAX_PREC(rop); 176 177 mpc_init2 (x, 2); 178 mpc_init2 (y, 2); 179 180 err = 7; 181 182 do 183 { 184 mpfr_exp_t k, exr, eyr, eyi, ezr; 185 186 ok = 0; 187 188 /* FIXME: prevent addition overflow */ 189 prec += mpc_ceil_log2 (prec) + err; 190 mpc_set_prec (x, prec); 191 mpc_set_prec (y, prec); 192 193 /* rounding away from zero: except in the cases x=0 or y=0 (processed 194 above), sin x and cos y are never exact, so rounding away from 0 is 195 rounding towards 0 and adding one ulp to the absolute value */ 196 mpc_sin_cos (x, y, op, MPC_RNDZZ, MPC_RNDZZ); 197 MPFR_ADD_ONE_ULP (mpc_realref (x)); 198 MPFR_ADD_ONE_ULP (mpc_imagref (x)); 199 MPFR_ADD_ONE_ULP (mpc_realref (y)); 200 MPFR_ADD_ONE_ULP (mpc_imagref (y)); 201 MPC_ASSERT (mpfr_zero_p (mpc_realref (x)) == 0); 202 203 if ( mpfr_inf_p (mpc_realref (x)) || mpfr_inf_p (mpc_imagref (x)) 204 || mpfr_inf_p (mpc_realref (y)) || mpfr_inf_p (mpc_imagref (y))) { 205 /* If the real or imaginary part of x is infinite, it means that 206 Im(op) was large, in which case the result is 207 sign(tan(Re(op)))*0 + sign(Im(op))*I, 208 where sign(tan(Re(op))) = sign(Re(x))*sign(Re(y)). */ 209 int inex_re, inex_im; 210 mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN); 211 if (mpfr_sgn (mpc_realref (x)) * mpfr_sgn (mpc_realref (y)) < 0) 212 { 213 mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN); 214 inex_re = 1; 215 } 216 else 217 inex_re = -1; /* +0 is rounded down */ 218 if (mpfr_sgn (mpc_imagref (op)) > 0) 219 { 220 mpfr_set_ui (mpc_imagref (rop), 1, GMP_RNDN); 221 inex_im = 1; 222 } 223 else 224 { 225 mpfr_set_si (mpc_imagref (rop), -1, GMP_RNDN); 226 inex_im = -1; 227 } 228 inex = MPC_INEX(inex_re, inex_im); 229 goto end; 230 } 231 232 exr = mpfr_get_exp (mpc_realref (x)); 233 eyr = mpfr_get_exp (mpc_realref (y)); 234 eyi = mpfr_get_exp (mpc_imagref (y)); 235 236 /* some parts of the quotient may be exact */ 237 inex = mpc_div (x, x, y, MPC_RNDZZ); 238 /* OP is no pure real nor pure imaginary, so in theory the real and 239 imaginary parts of its tangent cannot be null. However due to 240 rouding errors this might happen. Consider for example 241 tan(1+14*I) = 1.26e-10 + 1.00*I. For small precision sin(op) and 242 cos(op) differ only by a factor I, thus after mpc_div x = I and 243 its real part is zero. */ 244 if (mpfr_zero_p (mpc_realref (x)) || mpfr_zero_p (mpc_imagref (x))) 245 { 246 err = prec; /* double precision */ 247 continue; 248 } 249 if (MPC_INEX_RE (inex)) 250 MPFR_ADD_ONE_ULP (mpc_realref (x)); 251 if (MPC_INEX_IM (inex)) 252 MPFR_ADD_ONE_ULP (mpc_imagref (x)); 253 MPC_ASSERT (mpfr_zero_p (mpc_realref (x)) == 0); 254 ezr = mpfr_get_exp (mpc_realref (x)); 255 256 /* FIXME: compute 257 k = Exp(Re(x))+Exp(Re(y))-2min{Exp(Re(y)), Exp(Im(y))}-Exp(Re(x/y)) 258 avoiding overflow */ 259 k = exr - ezr + MPC_MAX(-eyr, eyr - 2 * eyi); 260 err = k < 2 ? 7 : (k == 2 ? 8 : (5 + k)); 261 262 /* Can the real part be rounded? */ 263 ok = (!mpfr_number_p (mpc_realref (x))) 264 || mpfr_can_round (mpc_realref(x), prec - err, GMP_RNDN, GMP_RNDZ, 265 MPC_PREC_RE(rop) + (MPC_RND_RE(rnd) == GMP_RNDN)); 266 267 if (ok) 268 { 269 /* Can the imaginary part be rounded? */ 270 ok = (!mpfr_number_p (mpc_imagref (x))) 271 || mpfr_can_round (mpc_imagref(x), prec - 6, GMP_RNDN, GMP_RNDZ, 272 MPC_PREC_IM(rop) + (MPC_RND_IM(rnd) == GMP_RNDN)); 273 } 274 } 275 while (ok == 0); 276 277 inex = mpc_set (rop, x, rnd); 278 279 end: 280 mpc_clear (x); 281 mpc_clear (y); 282 283 return inex; 284 } 285