1 /* mpfr_sinu -- sinu(x) = sin(2*pi*x/u)
2 mpfr_sinpi -- sinpi(x) = sin(pi*x)
3
4 Copyright 2020-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 /* References:
28 * Steve Kargl wrote sinpi and friends for FreeBSD's libm under BSD
29 license:
30 https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=218514
31 */
32
33 /* put in y the correctly rounded value of sin(2*pi*x/u) */
34 int
mpfr_sinu(mpfr_ptr y,mpfr_srcptr x,unsigned long u,mpfr_rnd_t rnd_mode)35 mpfr_sinu (mpfr_ptr y, mpfr_srcptr x, unsigned long u, mpfr_rnd_t rnd_mode)
36 {
37 mpfr_srcptr xp;
38 mpfr_prec_t precy, prec;
39 mpfr_exp_t expx, expt, err;
40 mpfr_t t, xr;
41 int inexact = 0, nloops = 0, underflow = 0;
42 MPFR_ZIV_DECL (loop);
43 MPFR_SAVE_EXPO_DECL (expo);
44
45 MPFR_LOG_FUNC (
46 ("x[%Pd]=%.*Rg u=%lu rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, u,
47 rnd_mode),
48 ("y[%Pd]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y,
49 inexact));
50
51 if (u == 0 || MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
52 {
53 /* for u=0, return NaN */
54 if (u == 0 || MPFR_IS_NAN (x) || MPFR_IS_INF (x))
55 {
56 MPFR_SET_NAN (y);
57 MPFR_RET_NAN;
58 }
59 else /* x is zero */
60 {
61 MPFR_ASSERTD (MPFR_IS_ZERO (x));
62 MPFR_SET_ZERO (y);
63 MPFR_SET_SAME_SIGN (y, x);
64 MPFR_RET (0);
65 }
66 }
67
68 MPFR_SAVE_EXPO_MARK (expo);
69
70 /* Range reduction. We do not need to reduce the argument if it is
71 already reduced (|x| < u).
72 Note that the case |x| = u is better in the "else" branch as it
73 will give xr = 0. */
74 if (mpfr_cmpabs_ui (x, u) < 0)
75 {
76 xp = x;
77 }
78 else
79 {
80 mpfr_exp_t p = MPFR_GET_PREC (x) - MPFR_GET_EXP (x);
81 int inex;
82
83 /* Let's compute xr = x mod u, with signbit(xr) = signbit(x), which
84 may be important when x is a multiple of u, in which case xr = 0
85 (but this property is actually not needed in the code below).
86 The precision of xr is chosen to ensure that x mod u is exactly
87 representable in xr, e.g., the maximum size of u + the length of
88 the fractional part of x. Note that since |x| >= u in this branch,
89 the additional memory amount will not be more than the one of x.
90 */
91 mpfr_init2 (xr, sizeof (unsigned long) * CHAR_BIT + (p < 0 ? 0 : p));
92 MPFR_DBGRES (inex = mpfr_fmod_ui (xr, x, u, MPFR_RNDN)); /* exact */
93 MPFR_ASSERTD (inex == 0);
94 if (MPFR_IS_ZERO (xr))
95 {
96 mpfr_clear (xr);
97 MPFR_SAVE_EXPO_FREE (expo);
98 MPFR_SET_ZERO (y);
99 MPFR_SET_SAME_SIGN (y, x);
100 MPFR_RET (0);
101 }
102 xp = xr;
103 }
104
105 /* now |xp/u| < 1 */
106
107 precy = MPFR_GET_PREC (y);
108 expx = MPFR_GET_EXP (xp);
109 /* For x large, since argument reduction is expensive, we want to avoid
110 any failure in Ziv's strategy, thus we take into account expx too. */
111 prec = precy + MAX(expx, MPFR_INT_CEIL_LOG2 (precy)) + 8;
112 MPFR_ASSERTD(prec >= 2);
113 mpfr_init2 (t, prec);
114 MPFR_ZIV_INIT (loop, prec);
115 for (;;)
116 {
117 nloops ++;
118 /* In the error analysis below, xp stands for x.
119 We first compute an approximation t of 2*pi*x/u, then call sin(t).
120 If t = 2*pi*x/u + s, then |sin(t) - sin(2*pi*x/u)| <= |s|. */
121 mpfr_set_prec (t, prec);
122 mpfr_const_pi (t, MPFR_RNDN); /* t = pi * (1 + theta1) where
123 |theta1| <= 2^-prec */
124 mpfr_mul_2ui (t, t, 1, MPFR_RNDN); /* t = 2*pi * (1 + theta1) */
125 mpfr_mul (t, t, xp, MPFR_RNDN); /* t = 2*pi*x * (1 + theta2)^2 where
126 |theta2| <= 2^-prec */
127 mpfr_div_ui (t, t, u, MPFR_RNDN); /* t = 2*pi*x/u * (1 + theta3)^3 where
128 |theta3| <= 2^-prec */
129 /* if t is zero here, it means the division by u underflows, then
130 sin(t) also underflows, since |sin(x)| <= |x|. */
131 if (MPFR_UNLIKELY (MPFR_IS_ZERO (t)))
132 {
133 inexact = mpfr_underflow (y, rnd_mode, MPFR_SIGN(t));
134 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_INEXACT
135 | MPFR_FLAGS_UNDERFLOW);
136 underflow = 1;
137 goto end;
138 }
139 /* since prec >= 2, |(1 + theta3)^3 - 1| <= 4*theta3 <= 2^(2-prec) */
140 expt = MPFR_GET_EXP (t);
141 /* we have |s| <= 2^(expt + 2 - prec) */
142 mpfr_sin (t, t, MPFR_RNDA);
143 /* t cannot be zero here, since we excluded t=0 before, which is the
144 only exact case where sin(t)=0, and we round away from zero */
145 err = expt + 2 - prec;
146 expt = MPFR_GET_EXP (t); /* new exponent of t */
147 /* the total error is bounded by 2^err + ulp(t) = 2^err + 2^(expt-prec)
148 thus if err <= expt-prec, it is bounded by 2^(expt-prec+1),
149 otherwise it is bounded by 2^(err+1). */
150 err = (err <= expt - prec) ? expt - prec + 1 : err + 1;
151 /* normalize err for mpfr_can_round */
152 err = expt - err;
153 if (MPFR_CAN_ROUND (t, err, precy, rnd_mode))
154 break;
155 /* Check exact cases only after the first level of Ziv' strategy, to
156 avoid slowing down the average case. Exact cases are:
157 (a) 2*pi*x/u is a multiple of pi/2, i.e., x/u is a multiple of 1/4
158 (b) 2*pi*x/u is +/-pi/6 modulo pi, i.e., x/u = +/-1/12 mod 1/2 */
159 if (nloops == 1)
160 {
161 /* detect case (a) */
162 inexact = mpfr_div_ui (t, xp, u, MPFR_RNDA);
163 mpfr_mul_2ui (t, t, 2, MPFR_RNDA);
164 if (inexact == 0 && mpfr_integer_p (t))
165 {
166 if (!mpfr_odd_p (t))
167 /* t is even: we have a multiple of pi, thus sinu = 0,
168 for the sign, we follow IEEE 754-2019: sinPi(+n) is +0
169 and sinPi(-n) is -0 for positive integers n, so that the
170 function is odd. */
171 mpfr_set_zero (y, MPFR_IS_POS(t) ? +1 : -1);
172 else /* t is odd */
173 {
174 inexact = mpfr_sub_ui (t, t, 1, MPFR_RNDZ);
175 MPFR_ASSERTD(inexact == 0);
176 inexact = mpfr_div_2ui (t, t, 1, MPFR_RNDZ);
177 MPFR_ASSERTD(inexact == 0);
178 if (MPFR_IS_ZERO (t) || !mpfr_odd_p (t))
179 /* case pi/2: sinu = 1 */
180 mpfr_set_ui (y, 1, MPFR_RNDZ);
181 else
182 mpfr_set_si (y, -1, MPFR_RNDZ);
183 }
184 goto end;
185 }
186 /* detect case (b): this can only occur if u is divisible by 3 */
187 if ((u % 3) == 0)
188 {
189 inexact = mpfr_div_ui (t, xp, u / 3, MPFR_RNDZ);
190 /* t should be +/-1/4 mod 3/2 */
191 mpfr_mul_2ui (t, t, 2, MPFR_RNDZ);
192 /* t should be +/-1 mod 6, i.e., in {1,5,7,11} mod 12:
193 t = 1 mod 6: case pi/6: return 1/2
194 t = 5 mod 6: case 5pi/6: return 1/2
195 t = 7 mod 6: case 7pi/6: return -1/2
196 t = 11 mod 6: case 11pi/6: return -1/2 */
197 if (inexact == 0 && mpfr_integer_p (t))
198 {
199 mpz_t z;
200 unsigned long mod12;
201 mpz_init (z);
202 inexact = mpfr_get_z (z, t, MPFR_RNDZ);
203 MPFR_ASSERTN(inexact == 0);
204 mod12 = mpz_fdiv_ui (z, 12);
205 mpz_clear (z);
206 if (mod12 == 1 || mod12 == 5)
207 {
208 mpfr_set_ui_2exp (y, 1, -1, MPFR_RNDZ);
209 goto end;
210 }
211 else if (mod12 == 7 || mod12 == 11)
212 {
213 mpfr_set_si_2exp (y, -1, -1, MPFR_RNDZ);
214 goto end;
215 }
216 }
217 }
218 }
219 MPFR_ZIV_NEXT (loop, prec);
220 }
221 MPFR_ZIV_FREE (loop);
222
223 inexact = mpfr_set (y, t, rnd_mode);
224
225 end:
226 mpfr_clear (t);
227 if (xp != x)
228 {
229 MPFR_ASSERTD (xp == xr);
230 mpfr_clear (xr);
231 }
232 MPFR_SAVE_EXPO_FREE (expo);
233 return underflow ? inexact : mpfr_check_range (y, inexact, rnd_mode);
234 }
235
236 int
mpfr_sinpi(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode)237 mpfr_sinpi (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
238 {
239 return mpfr_sinu (y, x, 2, rnd_mode);
240 }
241