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[%Pu]=%.*Rg u=%lu rnd=%d", mpfr_get_prec(x), mpfr_log_prec, x, u, 40 rnd_mode), 41 ("y[%Pu]=%.*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 } 170