1/* mpfr_tanu -- tanu(x) = tan(2*pi*x/u) 2 mpfr_tanpi -- tanpi(x) = tan(pi*x) 3 4Copyright 2020-2023 Free Software Foundation, Inc. 5Contributed by the AriC and Caramba projects, INRIA. 6 7This file is part of the GNU MPFR Library. 8 9The GNU MPFR Library is free software; you can redistribute it and/or modify 10it under the terms of the GNU Lesser General Public License as published by 11the Free Software Foundation; either version 3 of the License, or (at your 12option) any later version. 13 14The GNU MPFR Library is distributed in the hope that it will be useful, but 15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 16or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 17License for more details. 18 19You should have received a copy of the GNU Lesser General Public License 20along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 21https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2251 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 23 24#define MPFR_NEED_LONGLONG_H 25#include "mpfr-impl.h" 26 27/* put in y the correctly rounded value of tan(2*pi*x/u) */ 28int 29mpfr_tanu (mpfr_ptr y, mpfr_srcptr x, unsigned long u, mpfr_rnd_t rnd_mode) 30{ 31 mpfr_srcptr xp; 32 mpfr_prec_t precy, prec; 33 mpfr_exp_t expx, expt, err; 34 mpfr_t t, xr; 35 int inexact = 0, nloops = 0, underflow = 0; 36 MPFR_ZIV_DECL (loop); 37 MPFR_SAVE_EXPO_DECL (expo); 38 39 MPFR_LOG_FUNC ( 40 ("x[%Pd]=%.*Rg u=%lu rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, u, 41 rnd_mode), 42 ("y[%Pd]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, 43 inexact)); 44 45 if (u == 0 || MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 46 { 47 /* for u=0, return NaN */ 48 if (u == 0 || MPFR_IS_NAN (x) || MPFR_IS_INF (x)) 49 { 50 MPFR_SET_NAN (y); 51 MPFR_RET_NAN; 52 } 53 else /* x is zero */ 54 { 55 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 56 MPFR_SET_ZERO (y); 57 MPFR_SET_SAME_SIGN (y, x); 58 MPFR_RET (0); 59 } 60 } 61 62 MPFR_SAVE_EXPO_MARK (expo); 63 64 /* Range reduction. We do not need to reduce the argument if it is 65 already reduced (|x| < u). 66 Note that the case |x| = u is better in the "else" branch as it 67 will give xr = 0. */ 68 if (mpfr_cmpabs_ui (x, u) < 0) 69 { 70 xp = x; 71 } 72 else 73 { 74 mpfr_exp_t p = MPFR_GET_PREC (x) - MPFR_GET_EXP (x); 75 int inex; 76 77 /* Let's compute xr = x mod u, with signbit(xr) = signbit(x), which 78 may be important when x is a multiple of u, in which case xr = 0 79 (but this property is actually not needed in the code below). 80 The precision of xr is chosen to ensure that x mod u is exactly 81 representable in xr, e.g., the maximum size of u + the length of 82 the fractional part of x. Note that since |x| >= u in this branch, 83 the additional memory amount will not be more than the one of x. 84 Note that due to the rules on the special values, we needed to 85 consider a period of u instead of u/2. */ 86 mpfr_init2 (xr, sizeof (unsigned long) * CHAR_BIT + (p < 0 ? 0 : p)); 87 MPFR_DBGRES (inex = mpfr_fmod_ui (xr, x, u, MPFR_RNDN)); /* exact */ 88 MPFR_ASSERTD (inex == 0); 89 if (MPFR_IS_ZERO (xr)) 90 { 91 mpfr_clear (xr); 92 MPFR_SAVE_EXPO_FREE (expo); 93 MPFR_SET_ZERO (y); 94 MPFR_SET_SAME_SIGN (y, x); 95 MPFR_RET (0); 96 } 97 xp = xr; 98 } 99 100 /* now |xp/u| < 1 */ 101 102 precy = MPFR_GET_PREC (y); 103 expx = MPFR_GET_EXP (xp); 104 /* For x large, since argument reduction is expensive, we want to avoid 105 any failure in Ziv's strategy, thus we take into account expx too. */ 106 prec = precy + MAX(expx,MPFR_INT_CEIL_LOG2(precy)) + 8; 107 MPFR_ASSERTD(prec >= 2); 108 mpfr_init2 (t, prec); 109 MPFR_ZIV_INIT (loop, prec); 110 for (;;) 111 { 112 int inex; 113 nloops ++; 114 /* In the error analysis below, xp stands for x. 115 We first compute an approximation t of 2*pi*x/u, then call tan(t). 116 If t = 2*pi*x/u + s, then 117 |tan(t) - tan(2*pi*x/u)| = |s| * (1 + tan(v)^2) where v is in the 118 interval [t, t+s]. If we ensure that |t| >= |2*pi*x/u|, since tan() is 119 increasing, we can bound tan(v)^2 by tan(t)^2. */ 120 mpfr_set_prec (t, prec); 121 mpfr_const_pi (t, MPFR_RNDU); /* t = pi * (1 + theta1) where 122 |theta1| <= 2^(1-prec) */ 123 mpfr_mul_2ui (t, t, 1, MPFR_RNDN); /* t = 2*pi * (1 + theta1) */ 124 mpfr_mul (t, t, xp, MPFR_RNDA); /* t = 2*pi*x * (1 + theta2)^2 where 125 |theta2| <= 2^(1-prec) */ 126 inex = mpfr_div_ui (t, t, u, MPFR_RNDN); 127 /* t = 2*pi*x/u * (1 + theta3)^3 where |theta3| <= 2^(1-prec) */ 128 /* if t is zero here, it means the division by u underflows, then 129 tan(t) also underflows, since |tan(x)| <= |x|. */ 130 if (MPFR_UNLIKELY (MPFR_IS_ZERO (t))) 131 { 132 inexact = mpfr_underflow (y, rnd_mode, MPFR_SIGN(t)); 133 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_INEXACT 134 | MPFR_FLAGS_UNDERFLOW); 135 underflow = 1; 136 goto end; 137 } 138 /* emulate mpfr_div_ui (t, t, u, MPFR_RNDA) above, so that t is rounded 139 away from zero */ 140 if (MPFR_SIGN(t) > 0 && inex < 0) 141 mpfr_nextabove (t); 142 else if (MPFR_SIGN(t) < 0 && inex > 0) 143 mpfr_nextbelow (t); 144 expt = MPFR_GET_EXP (t); 145 /* since prec >= 3, |(1 + theta3)^3 - 1| <= 4*theta3 <= 2^(3-prec) 146 thus |s| = |t - 2*pi*x/u| <= |t| * 2^(3-prec) */ 147 mpfr_tan (t, t, MPFR_RNDA); 148 { 149 /* compute an upper bound for 1+tan(t)^2 */ 150 mpfr_t z; 151 mpfr_init2 (z, 64); 152 mpfr_sqr (z, t, MPFR_RNDU); 153 mpfr_add_ui (z, z, 1, MPFR_RNDU); 154 expt += MPFR_GET_EXP (z); 155 /* now |t - tan(2*pi*x/u)| <= ulp(t) + 2^(expt + 3 - prec) */ 156 mpfr_clear (z); 157 } 158 /* t cannot be zero here, since we excluded t=0 before, which is the 159 only exact case where tan(t)=0, and we round away from zero */ 160 err = expt + 3 - prec; 161 expt = MPFR_GET_EXP (t); /* new exponent of t */ 162 /* the total error is bounded by 2^err + ulp(t) = 2^err + 2^(expt-prec) 163 thus if err <= expt-prec, it is bounded by 2^(expt-prec+1), 164 otherwise it is bounded by 2^(err+1). */ 165 err = (err <= expt - prec) ? expt - prec + 1 : err + 1; 166 /* normalize err for mpfr_can_round */ 167 err = expt - err; 168 if (MPFR_CAN_ROUND (t, err, precy, rnd_mode)) 169 break; 170 /* Check exact cases only after the first level of Ziv' strategy, to 171 avoid slowing down the average case. Exact cases are when 2*pi*x/u 172 is a multiple of pi/4, i.e., x/u a multiple of 1/8: 173 (a) x/u = {0,1/2} mod 1: return +0 or -0 174 (b) x/u = {1/4,3/4} mod 1: return +Inf or -Inf 175 (c) x/u = {1/8,3/8,5/8,7/8} mod 1: return 1 or -1 */ 176 if (nloops == 1) 177 { 178 inexact = mpfr_div_ui (t, xp, u, MPFR_RNDA); 179 mpfr_mul_2ui (t, t, 3, MPFR_RNDA); 180 if (inexact == 0 && mpfr_integer_p (t)) 181 { 182 mpz_t z; 183 unsigned long mod8; 184 mpz_init (z); 185 inexact = mpfr_get_z (z, t, MPFR_RNDZ); 186 MPFR_ASSERTN(inexact == 0); 187 mod8 = mpz_fdiv_ui (z, 8); 188 mpz_clear (z); 189 if (mod8 == 0 || mod8 == 4) /* case (a) */ 190 mpfr_set_zero (y, ((mod8 == 0) ? +1 : -1) * MPFR_SIGN (x)); 191 else if (mod8 == 2 || mod8 == 6) /* case (b) */ 192 { 193 mpfr_set_inf (y, (mod8 == 2) ? +1 : -1); 194 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_DIVBY0); 195 } 196 else /* case (c) */ 197 { 198 if (mod8 == 1 || mod8 == 5) 199 mpfr_set_ui (y, 1, rnd_mode); 200 else 201 mpfr_set_si (y, -1, rnd_mode); 202 } 203 goto end; 204 } 205 } 206 MPFR_ZIV_NEXT (loop, prec); 207 } 208 MPFR_ZIV_FREE (loop); 209 210 inexact = mpfr_set (y, t, rnd_mode); 211 212 end: 213 mpfr_clear (t); 214 if (xp != x) 215 { 216 MPFR_ASSERTD (xp == xr); 217 mpfr_clear (xr); 218 } 219 MPFR_SAVE_EXPO_FREE (expo); 220 return underflow ? inexact : mpfr_check_range (y, inexact, rnd_mode); 221} 222 223int 224mpfr_tanpi (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 225{ 226 return mpfr_tanu (y, x, 2, rnd_mode); 227} 228