1/* mpfr_acosu  -- acosu(x)  = acos(x)*u/(2*pi)
2   mpfr_acospi -- acospi(x) = acos(x)/pi
3
4Copyright 2021-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 acos(x)*u/(2*pi) */
28int
29mpfr_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
165int
166mpfr_acospi (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
167{
168  return mpfr_acosu (y, x, 2, rnd_mode);
169}
170