1 /* mpfr_acosu -- acosu(x) = acos(x)*u/(2*pi)
2 mpfr_acospi -- acospi(x) = acos(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 acos(x)*u/(2*pi) */
28 int
29 mpfr_acosu (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 compared, inexact;
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 inexact));
43
44 /* Singular cases */
45 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
46 {
47 if (MPFR_IS_NAN (x) || MPFR_IS_INF (x))
48 {
49 MPFR_SET_NAN (y);
50 MPFR_RET_NAN;
51 }
52 else /* necessarily x=0 */
53 {
54 MPFR_ASSERTD(MPFR_IS_ZERO(x));
55 /* acos(0)=Pi/2 thus acosu(0)=u/4 */
56 return mpfr_set_ui_2exp (y, u, -2, rnd_mode);
57 }
58 }
59
60 compared = mpfr_cmpabs_ui (x, 1);
61 if (compared > 0)
62 {
63 /* acosu(x) = NaN for |x| > 1, included for u=0, since NaN*0 = NaN */
64 MPFR_SET_NAN (y);
65 MPFR_RET_NAN;
66 }
67
68 if (u == 0) /* return +0 since acos(x)>=0 */
69 {
70 MPFR_SET_ZERO (y);
71 MPFR_SET_POS (y);
72 MPFR_RET (0);
73 }
74
75 if (compared == 0)
76 {
77 /* |x| = 1: acosu(1,u) = +0, acosu(-1,u)=u/2 */
78 if (MPFR_SIGN(x) > 0) /* IEEE-754 2019: acosPi(1) = +0 */
79 return mpfr_set_ui (y, 0, rnd_mode);
80 else
81 return mpfr_set_ui_2exp (y, u, -1, rnd_mode);
82 }
83
84 /* acos(1/2) = pi/6 and acos(-1/2) = pi/3, thus in these cases acos(x,u)
85 is exact when u is a multiple of 3 */
86 if (mpfr_cmp_si_2exp (x, MPFR_SIGN(x), -1) == 0 && (u % 3) == 0)
87 return mpfr_set_si_2exp (y, u / 3, MPFR_IS_NEG (x) ? 0 : -1, rnd_mode);
88
89 prec = MPFR_PREC (y);
90
91 MPFR_SAVE_EXPO_MARK (expo);
92
93 /* For |x|<0.5, we have acos(x) = pi/2 - x*r(x) with |r(x)| < 1.05
94 thus acosu(x,u) = u/4*(1 - x*s(x)) with 0 <= s(x) < 1.
95 If EXP(x) <= -prec-3, then |u/4*x*s(x)| < u/4*2^(-prec-3) < ulp(u/4)/8
96 <= ulp(RN(u/4))/4, thus the result will be u/4, nextbelow(u/4) or
97 nextabove(u/4).
98 Warning: when u/4 is a power of two, the difference between u/4 and
99 nextbelow(u/4) is only 1/4*ulp(u/4).
100 We also require x < 2^-64, so that in the case u/4 is not exact,
101 the contribution of x*s(x) is smaller compared to the last bit of u. */
102 expx = MPFR_GET_EXP(x);
103 if (expx <= -64 && expx <= - (mpfr_exp_t) prec - 3)
104 {
105 prec = (MPFR_PREC(y) <= 63) ? 65 : MPFR_PREC(y) + 2;
106 /* now prec > 64 and prec > MPFR_PREC(y)+1 */
107 mpfr_init2 (tmp, prec);
108 inexact = mpfr_set_ui (tmp, u, MPFR_RNDN); /* exact since prec >= 64 */
109 MPFR_ASSERTD(inexact == 0);
110 /* for x>0, we have acos(x) < pi/2; for x<0, we have acos(x) > pi/2 */
111 if (MPFR_IS_POS(x))
112 mpfr_nextbelow (tmp);
113 else
114 mpfr_nextabove (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 /* since prec >= PREC(y)+2, the rounding of tmp is correct */
120 inexact = mpfr_div_2ui (y, tmp, 2, rnd_mode);
121 mpfr_clear (tmp);
122 goto end;
123 }
124
125 prec += MPFR_INT_CEIL_LOG2(prec) + 10;
126
127 mpfr_init2 (tmp, prec);
128 mpfr_init2 (pi, prec);
129
130 MPFR_ZIV_INIT (loop, prec);
131 for (;;)
132 {
133 /* In the error analysis below, each thetax denotes a variable such that
134 |thetax| <= 2^-prec */
135 mpfr_acos (tmp, x, MPFR_RNDN);
136 /* tmp = acos(x) * (1 + theta1) */
137 mpfr_const_pi (pi, MPFR_RNDN);
138 /* pi = Pi * (1 + theta2) */
139 mpfr_div (tmp, tmp, pi, MPFR_RNDN);
140 /* tmp = acos(x)/Pi * (1 + theta3)^3 */
141 mpfr_mul_ui (tmp, tmp, u, MPFR_RNDN);
142 /* tmp = acos(x)*u/Pi * (1 + theta4)^4 */
143 mpfr_div_2ui (tmp, tmp, 1, MPFR_RNDN); /* exact */
144 /* tmp = acos(x)*u/(2*Pi) * (1 + theta4)^4 */
145 /* since |(1 + theta4)^4 - 1| <= 8*|theta4| for prec >= 2,
146 the relative error is less than 2^(3-prec) */
147 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - 3,
148 MPFR_PREC (y), rnd_mode)))
149 break;
150 MPFR_ZIV_NEXT (loop, prec);
151 mpfr_set_prec (tmp, prec);
152 mpfr_set_prec (pi, prec);
153 }
154 MPFR_ZIV_FREE (loop);
155
156 inexact = mpfr_set (y, tmp, rnd_mode);
157 mpfr_clear (tmp);
158 mpfr_clear (pi);
159
160 end:
161 MPFR_SAVE_EXPO_FREE (expo);
162 return mpfr_check_range (y, inexact, rnd_mode);
163 }
164
165 int
166 mpfr_acospi (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
167 {
168 return mpfr_acosu (y, x, 2, rnd_mode);
169 }