1 /* mpfr_asinu -- asinu(x) = asin(x)*u/(2*pi) 2 mpfr_asinpi -- asinpi(x) = asin(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 asin(x)*u/(2*pi) */ 28 int 29 mpfr_asinu (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 int compared, inexact; 34 MPFR_SAVE_EXPO_DECL (expo); 35 MPFR_ZIV_DECL (loop); 36 37 MPFR_LOG_FUNC 38 (("x[%Pu]=%.*Rg u=%lu rnd=%d", mpfr_get_prec(x), mpfr_log_prec, x, u, 39 rnd_mode), 40 ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, 41 inexact)); 42 43 /* Singular cases */ 44 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 45 { 46 if (MPFR_IS_NAN (x) || MPFR_IS_INF (x)) 47 { 48 MPFR_SET_NAN (y); 49 MPFR_RET_NAN; 50 } 51 else /* necessarily x=0 */ 52 { 53 MPFR_ASSERTD(MPFR_IS_ZERO(x)); 54 /* asin(0)=0 with same sign, even when u=0 to ensure 55 asinu(-x,u) = -asinu(x,u) */ 56 MPFR_SET_ZERO (y); 57 MPFR_SET_SAME_SIGN (y, x); 58 MPFR_RET (0); /* exact result */ 59 } 60 } 61 62 compared = mpfr_cmpabs_ui (x, 1); 63 if (compared > 0) 64 { 65 /* asinu(x) = NaN for |x| > 1, included for u=0, since NaN*0 = NaN */ 66 MPFR_SET_NAN (y); 67 MPFR_RET_NAN; 68 } 69 70 if (u == 0) /* return 0 with sign of x */ 71 { 72 MPFR_SET_ZERO (y); 73 MPFR_SET_POS (y); 74 MPFR_RET (0); 75 } 76 77 if (compared == 0) 78 { 79 /* |x| = 1: asinu(1,u) = u/4, asinu(-1,u)=-u/4 */ 80 /* we can't use mpfr_set_si_2exp with -u since -u might not be 81 representable as long */ 82 if (MPFR_SIGN(x) > 0) 83 return mpfr_set_ui_2exp (y, u, -2, rnd_mode); 84 else 85 { 86 inexact = mpfr_set_ui_2exp (y, u, -2, MPFR_INVERT_RND(rnd_mode)); 87 MPFR_CHANGE_SIGN(y); 88 return -inexact; 89 } 90 } 91 92 /* asin(+/-1/2) = +/-pi/6, thus asin(+/-1/2,u) = +/-u/12 is exact when u is 93 a multiple of 3 */ 94 if (mpfr_cmp_si_2exp (x, MPFR_SIGN(x), -1) == 0 && (u % 3) == 0) 95 { 96 long v = u / 3; 97 if (MPFR_IS_NEG (x)) 98 v = -v; 99 return mpfr_set_si_2exp (y, v, -2, rnd_mode); 100 } 101 102 prec = MPFR_PREC (y); 103 104 MPFR_SAVE_EXPO_MARK (expo); 105 106 prec += MPFR_INT_CEIL_LOG2(prec) + 10; 107 108 mpfr_init2 (tmp, prec); 109 mpfr_init2 (pi, prec); 110 111 MPFR_ZIV_INIT (loop, prec); 112 for (;;) 113 { 114 /* In the error analysis below, each thetax denotes a variable such that 115 |thetax| <= 2^(1-prec) */ 116 mpfr_asin (tmp, x, MPFR_RNDA); 117 /* tmp = asin(x) * (1 + theta1), and tmp cannot be zero since we rounded 118 away from zero, and the case x=0 was treated before */ 119 /* first multiply by u to avoid underflow issues */ 120 mpfr_mul_ui (tmp, tmp, u, MPFR_RNDA); 121 /* tmp = asin(x)*u * (1 + theta2)^2, and |tmp| >= 0.5*2^emin */ 122 mpfr_const_pi (pi, MPFR_RNDZ); /* round toward zero since we we will 123 divide by pi, to round tmp away */ 124 /* pi = Pi * (1 + theta3) */ 125 mpfr_div (tmp, tmp, pi, MPFR_RNDA); 126 /* tmp = asin(x)*u/Pi * (1 + theta4)^4, with |tmp| > 0 */ 127 /* since we rounded away from 0, if we get 0.5*2^emin here, it means 128 |asinu(x,u)| < 0.25*2^emin (pi is not exact) thus we have underflow */ 129 if (MPFR_EXP(tmp) == __gmpfr_emin) 130 { 131 /* mpfr_underflow rounds away for RNDN */ 132 mpfr_clear (tmp); 133 mpfr_clear (pi); 134 MPFR_SAVE_EXPO_FREE (expo); 135 return mpfr_underflow (y, 136 (rnd_mode == MPFR_RNDN) ? MPFR_RNDZ : rnd_mode, 1); 137 } 138 mpfr_div_2ui (tmp, tmp, 1, MPFR_RNDA); /* exact */ 139 /* tmp = asin(x)*u/(2*Pi) * (1 + theta4)^4 */ 140 /* since |(1 + theta4)^4 - 1| <= 8*|theta4| for prec >= 3, 141 the relative error is less than 2^(4-prec) */ 142 MPFR_ASSERTD(!MPFR_IS_ZERO(tmp)); 143 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - 4, 144 MPFR_PREC (y), rnd_mode))) 145 break; 146 MPFR_ZIV_NEXT (loop, prec); 147 mpfr_set_prec (tmp, prec); 148 mpfr_set_prec (pi, prec); 149 } 150 MPFR_ZIV_FREE (loop); 151 152 inexact = mpfr_set (y, tmp, rnd_mode); 153 mpfr_clear (tmp); 154 mpfr_clear (pi); 155 156 MPFR_SAVE_EXPO_FREE (expo); 157 return mpfr_check_range (y, inexact, rnd_mode); 158 } 159 160 int 161 mpfr_asinpi (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 162 { 163 return mpfr_asinu (y, x, 2, rnd_mode); 164 } 165