1 /* mpfr_atanu -- atanu(x) = atan(x)*u/(2*pi)
2 mpfr_atanpi -- atanpi(x) = atan(x)/pi
3
4 Copyright 2021-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 atan(x)*u/(2*pi) */
28 int
29 mpfr_atanu (mpfr_ptr y, mpfr_srcptr x, unsigned long u, mpfr_rnd_t rnd_mode)
30 {
31 mpfr_t tmp, pi;
32 mpfr_prec_t prec;
33 mpfr_exp_t expx;
34 int inex;
35 MPFR_SAVE_EXPO_DECL (expo);
36 MPFR_ZIV_DECL (loop);
37
38 MPFR_LOG_FUNC
39 (("x[%Pd]=%.*Rg u=%lu rnd=%d", mpfr_get_prec(x), mpfr_log_prec, x, u,
40 rnd_mode),
41 ("y[%Pd]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y,
42 inex));
43
44 /* Singular cases */
45 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
46 {
47 if (MPFR_IS_NAN (x))
48 {
49 MPFR_SET_NAN (y);
50 MPFR_RET_NAN;
51 }
52 else if (MPFR_IS_INF (x))
53 {
54 /* atanu(+Inf,u) = u/4, atanu(-Inf,u) = -u/4 */
55 if (MPFR_IS_POS (x))
56 return mpfr_set_ui_2exp (y, u, -2, rnd_mode);
57 else
58 {
59 inex = mpfr_set_ui_2exp (y, u, -2, MPFR_INVERT_RND (rnd_mode));
60 MPFR_CHANGE_SIGN (y);
61 return -inex;
62 }
63 }
64 else /* necessarily x=0 */
65 {
66 MPFR_ASSERTD(MPFR_IS_ZERO(x));
67 /* atan(0)=0 with same sign, even when u=0 to ensure
68 atanu(-x,u) = -atanu(x,u) */
69 MPFR_SET_ZERO (y);
70 MPFR_SET_SAME_SIGN (y, x);
71 MPFR_RET (0); /* exact result */
72 }
73 }
74
75 if (u == 0) /* return 0 with sign of x, which is coherent with case x=0 */
76 {
77 MPFR_SET_ZERO (y);
78 MPFR_SET_SAME_SIGN (y, x);
79 MPFR_RET (0);
80 }
81
82 if (mpfr_cmpabs_ui (x, 1) == 0)
83 {
84 /* |x| = 1: atanu(1,u) = u/8, atanu(-1,u)=-u/8 */
85 /* we can't use mpfr_set_si_2exp with -u since -u might not be
86 representable as long */
87 if (MPFR_SIGN(x) > 0)
88 return mpfr_set_ui_2exp (y, u, -3, rnd_mode);
89 else
90 {
91 inex = mpfr_set_ui_2exp (y, u, -3, MPFR_INVERT_RND(rnd_mode));
92 MPFR_CHANGE_SIGN(y);
93 return -inex;
94 }
95 }
96
97 /* For x>=1, we have pi/2-1/x < atan(x) < pi/2, thus
98 u/4-u/(2*pi*x) < atanu(x,u) < u/4, and the relative difference between
99 atanu(x,u) and u/4 is less than 2/(pi*x) < 1/x <= 2^(1-EXP(x)).
100 If the relative difference is <= 2^(-prec-2), then the difference
101 between atanu(x,u) and u/4 is <= 1/4*ulp(u/4) <= 1/2*ulp(RN(u/4)).
102 We also require x >= 2^64, which implies x > 2*u/pi, so that
103 (u-1)/4 < u/4-u/(2*pi*x) < u/4. */
104 expx = MPFR_GET_EXP(x);
105 if (expx >= 65 && expx - 1 >= MPFR_PREC(y) + 2)
106 {
107 prec = (MPFR_PREC(y) <= 63) ? 65 : MPFR_PREC(y) + 2;
108 /* now prec > 64 and prec > MPFR_PREC(y)+1 */
109 mpfr_init2 (tmp, prec);
110 /* since expx >= 65, we have emax >= 65, thus u is representable here,
111 and we don't need to work in an extended exponent range */
112 inex = mpfr_set_ui (tmp, u, MPFR_RNDN); /* exact since prec >= 64 */
113 MPFR_ASSERTD(inex == 0);
114 mpfr_nextbelow (tmp);
115 /* Since prec >= 65, the last significant bit of tmp is 1, and since
116 prec > PREC(y), tmp is not representable in the target precision,
117 which ensures we will get a correct ternary value below. */
118 MPFR_ASSERTD(mpfr_min_prec(tmp) > MPFR_PREC(y));
119 if (MPFR_SIGN(x) < 0)
120 MPFR_CHANGE_SIGN(tmp);
121 /* since prec >= PREC(y)+2, the rounding of tmp is correct */
122 inex = mpfr_div_2ui (y, tmp, 2, rnd_mode);
123 mpfr_clear (tmp);
124 return inex;
125 }
126
127 prec = MPFR_PREC (y);
128
129 MPFR_SAVE_EXPO_MARK (expo);
130
131 prec += MPFR_INT_CEIL_LOG2(prec) + 10;
132
133 mpfr_init2 (tmp, prec);
134 mpfr_init2 (pi, prec);
135
136 MPFR_ZIV_INIT (loop, prec);
137 for (;;)
138 {
139 /* In the error analysis below, each thetax denotes a variable such that
140 |thetax| <= 2^(1-prec) */
141 mpfr_atan (tmp, x, MPFR_RNDA);
142 /* tmp = atan(x) * (1 + theta1), and tmp cannot be zero since we rounded
143 away from zero, and the case x=0 was treated before */
144 /* first multiply by u to avoid underflow issues */
145 mpfr_mul_ui (tmp, tmp, u, MPFR_RNDA);
146 /* tmp = atan(x)*u * (1 + theta2)^2, and |tmp| >= 0.5*2^emin */
147 mpfr_const_pi (pi, MPFR_RNDZ); /* round toward zero since we we will
148 divide by pi, to round tmp away */
149 /* pi = Pi * (1 + theta3) */
150 mpfr_div (tmp, tmp, pi, MPFR_RNDA);
151 /* tmp = atan(x)*u/Pi * (1 + theta4)^4, with |tmp| > 0 */
152 /* since we rounded away from 0, if we get 0.5*2^emin here, it means
153 |atanu(x,u)| < 0.25*2^emin (pi is not exact) thus we have underflow */
154 if (MPFR_EXP(tmp) == __gmpfr_emin)
155 {
156 /* mpfr_underflow rounds away for RNDN */
157 mpfr_clear (tmp);
158 mpfr_clear (pi);
159 MPFR_SAVE_EXPO_FREE (expo);
160 return mpfr_underflow (y,
161 (rnd_mode == MPFR_RNDN) ? MPFR_RNDZ : rnd_mode, 1);
162 }
163 mpfr_div_2ui (tmp, tmp, 1, MPFR_RNDA); /* exact */
164 /* tmp = atan(x)*u/(2*Pi) * (1 + theta4)^4 */
165 /* since |(1 + theta4)^4 - 1| <= 8*|theta4| for prec >= 3,
166 the relative error is less than 2^(4-prec) */
167 MPFR_ASSERTD(!MPFR_IS_ZERO(tmp));
168 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - 4,
169 MPFR_PREC (y), rnd_mode)))
170 break;
171 MPFR_ZIV_NEXT (loop, prec);
172 mpfr_set_prec (tmp, prec);
173 mpfr_set_prec (pi, prec);
174 }
175 MPFR_ZIV_FREE (loop);
176
177 inex = mpfr_set (y, tmp, rnd_mode);
178 mpfr_clear (tmp);
179 mpfr_clear (pi);
180
181 MPFR_SAVE_EXPO_FREE (expo);
182 return mpfr_check_range (y, inex, rnd_mode);
183 }
184
185 int
186 mpfr_atanpi (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
187 {
188 return mpfr_atanu (y, x, 2, rnd_mode);
189 }