1 /* mpfr_tanu -- tanu(x) = tan(2*pi*x/u)
2 mpfr_tanpi -- tanpi(x) = tan(pi*x)
3
4 Copyright 2020-2023 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramba projects, INRIA.
6
7 This file is part of the GNU MPFR Library.
8
9 The GNU MPFR Library is free software; you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
13
14 The GNU MPFR Library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
21 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
22 51 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) */
28 int
29 mpfr_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
223 int
224 mpfr_tanpi (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
225 {
226 return mpfr_tanu (y, x, 2, rnd_mode);
227 }