1/* mpfr_sinu -- sinu(x) = sin(2*pi*x/u) 2 mpfr_sinpi -- sinpi(x) = sin(pi*x) 3 4Copyright 2020-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/* 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) */ 34int 35mpfr_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 236int 237mpfr_sinpi (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 238{ 239 return mpfr_sinu (y, x, 2, rnd_mode); 240} 241