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