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
mpfr_asinu(mpfr_ptr y,mpfr_srcptr x,unsigned long u,mpfr_rnd_t rnd_mode)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
mpfr_asinpi(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode)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