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 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 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