src/algorithms/trig_pi.rs

   1 use crate::{
   2     f16::F16,
   3     prim::{PrimFloat, PrimSInt, PrimUInt},
   4     traits::{Compare, Context, ConvertFrom, ConvertTo, Float, Make, Select},
   5 };
   6
   7 mod consts {
   8     #![allow(clippy::excessive_precision)]
   9     #![allow(dead_code)]
  10
  11     /// coefficients of taylor series for `sin(pi * x)` centered at `0`
  12     /// generated using:
  13     /// ```maxima,text
  14     /// fpprec:50$
  15     /// sinpi: bfloat(taylor(sin(%pi*x),x,0,19))$
  16     /// for i: 1 step 2 thru 19 do
  17     ///     printf(true, "pub(crate) const SINPI_KERNEL_TAYLOR_~d: f64 = ~a;~%", i, ssubst("e", "b", string(coeff(sinpi, x, i))))$
  18     /// ```
  19     pub(crate) const SINPI_KERNEL_TAYLOR_1: f64 =
  20         3.1415926535897932384626433832795028841971693993751e0;
  21     pub(crate) const SINPI_KERNEL_TAYLOR_3: f64 =
  22         -5.1677127800499700292460525111835658670375480943142e0;
  23     pub(crate) const SINPI_KERNEL_TAYLOR_5: f64 =
  24         2.550164039877345443856177583695296720669172555234e0;
  25     pub(crate) const SINPI_KERNEL_TAYLOR_7: f64 =
  26         -5.9926452932079207688773938354604004601536358636814e-1;
  27     pub(crate) const SINPI_KERNEL_TAYLOR_9: f64 =
  28         8.2145886611128228798802365523698344807837460797753e-2;
  29     pub(crate) const SINPI_KERNEL_TAYLOR_11: f64 =
  30         -7.370430945714350777259089957290781501211638236021e-3;
  31     pub(crate) const SINPI_KERNEL_TAYLOR_13: f64 =
  32         4.6630280576761256442062891447027174382819981361599e-4;
  33     pub(crate) const SINPI_KERNEL_TAYLOR_15: f64 =
  34         -2.1915353447830215827384652057094188859248708765956e-5;
  35     pub(crate) const SINPI_KERNEL_TAYLOR_17: f64 =
  36         7.9520540014755127847832068624575890327682459384282e-7;
  37     pub(crate) const SINPI_KERNEL_TAYLOR_19: f64 =
  38         -2.2948428997269873110203872385571587856074785581088e-8;
  39
  40     /// coefficients of taylor series for `cos(pi * x)` centered at `0`
  41     /// generated using:
  42     /// ```maxima,text
  43     /// fpprec:50$
  44     /// cospi: bfloat(taylor(cos(%pi*x),x,0,18))$
  45     /// for i: 0 step 2 thru 18 do
  46     ///     printf(true, "pub(crate) const COSPI_KERNEL_TAYLOR_~d: f64 = ~a;~%", i, ssubst("e", "b", string(coeff(cospi, x, i))))$
  47     /// ```
  48     pub(crate) const COSPI_KERNEL_TAYLOR_0: f64 = 1.0e0;
  49     pub(crate) const COSPI_KERNEL_TAYLOR_2: f64 =
  50         -4.9348022005446793094172454999380755676568497036204e0;
  51     pub(crate) const COSPI_KERNEL_TAYLOR_4: f64 =
  52         4.0587121264167682181850138620293796354053160696952e0;
  53     pub(crate) const COSPI_KERNEL_TAYLOR_6: f64 =
  54         -1.3352627688545894958753047828505831928711354556681e0;
  55     pub(crate) const COSPI_KERNEL_TAYLOR_8: f64 =
  56         2.3533063035889320454187935277546542154506893530856e-1;
  57     pub(crate) const COSPI_KERNEL_TAYLOR_10: f64 =
  58         -2.5806891390014060012598294252898849657186441048147e-2;
  59     pub(crate) const COSPI_KERNEL_TAYLOR_12: f64 =
  60         1.9295743094039230479033455636859576401684718150003e-3;
  61     pub(crate) const COSPI_KERNEL_TAYLOR_14: f64 =
  62         -1.0463810492484570711801672835223932761029733149091e-4;
  63     pub(crate) const COSPI_KERNEL_TAYLOR_16: f64 =
  64         4.3030695870329470072978237149669233008960901556009e-6;
  65     pub(crate) const COSPI_KERNEL_TAYLOR_18: f64 =
  66         -1.387895246221377211446808750399309343777037849978e-7;
  67 }
  68
  69 /// computes `sin(pi * x)` for `-0.25 <= x <= 0.25`
  70 /// not guaranteed to give correct sign for zero result
  71 /// has an error of up to 2ULP
  72 pub fn sin_pi_kernel_f16<Ctx: Context>(ctx: Ctx, x: Ctx::VecF16) -> Ctx::VecF16 {
  73     let x_sq = x * x;
  74     let mut v: Ctx::VecF16 = ctx.make(consts::SINPI_KERNEL_TAYLOR_5.to());
  75     v = v.mul_add_fast(x_sq, ctx.make(consts::SINPI_KERNEL_TAYLOR_3.to()));
  76     v = v.mul_add_fast(x_sq, ctx.make(consts::SINPI_KERNEL_TAYLOR_1.to()));
  77     v * x
  78 }
  79
  80 /// computes `cos(pi * x)` for `-0.25 <= x <= 0.25`
  81 /// has an error of up to 2ULP
  82 pub fn cos_pi_kernel_f16<Ctx: Context>(ctx: Ctx, x: Ctx::VecF16) -> Ctx::VecF16 {
  83     let x_sq = x * x;
  84     let mut v: Ctx::VecF16 = ctx.make(consts::COSPI_KERNEL_TAYLOR_4.to());
  85     v = v.mul_add_fast(x_sq, ctx.make(consts::COSPI_KERNEL_TAYLOR_2.to()));
  86     v.mul_add_fast(x_sq, ctx.make(consts::COSPI_KERNEL_TAYLOR_0.to()))
  87 }
  88
  89 /// computes `sin(pi * x)` for `-0.25 <= x <= 0.25`
  90 /// not guaranteed to give correct sign for zero result
  91 /// has an error of up to 2ULP
  92 pub fn sin_pi_kernel_f32<Ctx: Context>(ctx: Ctx, x: Ctx::VecF32) -> Ctx::VecF32 {
  93     let x_sq = x * x;
  94     let mut v: Ctx::VecF32 = ctx.make(consts::SINPI_KERNEL_TAYLOR_9.to());
  95     v = v.mul_add_fast(x_sq, ctx.make(consts::SINPI_KERNEL_TAYLOR_7.to()));
  96     v = v.mul_add_fast(x_sq, ctx.make(consts::SINPI_KERNEL_TAYLOR_5.to()));
  97     v = v.mul_add_fast(x_sq, ctx.make(consts::SINPI_KERNEL_TAYLOR_3.to()));
  98     v = v.mul_add_fast(x_sq, ctx.make(consts::SINPI_KERNEL_TAYLOR_1.to()));
  99     v * x
 100 }
 101
 102 /// computes `cos(pi * x)` for `-0.25 <= x <= 0.25`
 103 /// has an error of up to 2ULP
 104 pub fn cos_pi_kernel_f32<Ctx: Context>(ctx: Ctx, x: Ctx::VecF32) -> Ctx::VecF32 {
 105     let x_sq = x * x;
 106     let mut v: Ctx::VecF32 = ctx.make(consts::COSPI_KERNEL_TAYLOR_8.to());
 107     v = v.mul_add_fast(x_sq, ctx.make(consts::COSPI_KERNEL_TAYLOR_6.to()));
 108     v = v.mul_add_fast(x_sq, ctx.make(consts::COSPI_KERNEL_TAYLOR_4.to()));
 109     v = v.mul_add_fast(x_sq, ctx.make(consts::COSPI_KERNEL_TAYLOR_2.to()));
 110     v.mul_add_fast(x_sq, ctx.make(consts::COSPI_KERNEL_TAYLOR_0.to()))
 111 }
 112
 113 /// computes `sin(pi * x)` for `-0.25 <= x <= 0.25`
 114 /// not guaranteed to give correct sign for zero result
 115 /// has an error of up to 2ULP
 116 pub fn sin_pi_kernel_f64<Ctx: Context>(ctx: Ctx, x: Ctx::VecF64) -> Ctx::VecF64 {
 117     let x_sq = x * x;
 118     let mut v: Ctx::VecF64 = ctx.make(consts::SINPI_KERNEL_TAYLOR_15.to());
 119     v = v.mul_add_fast(x_sq, ctx.make(consts::SINPI_KERNEL_TAYLOR_13.to()));
 120     v = v.mul_add_fast(x_sq, ctx.make(consts::SINPI_KERNEL_TAYLOR_11.to()));
 121     v = v.mul_add_fast(x_sq, ctx.make(consts::SINPI_KERNEL_TAYLOR_9.to()));
 122     v = v.mul_add_fast(x_sq, ctx.make(consts::SINPI_KERNEL_TAYLOR_7.to()));
 123     v = v.mul_add_fast(x_sq, ctx.make(consts::SINPI_KERNEL_TAYLOR_5.to()));
 124     v = v.mul_add_fast(x_sq, ctx.make(consts::SINPI_KERNEL_TAYLOR_3.to()));
 125     v = v.mul_add_fast(x_sq, ctx.make(consts::SINPI_KERNEL_TAYLOR_1.to()));
 126     v * x
 127 }
 128
 129 /// computes `cos(pi * x)` for `-0.25 <= x <= 0.25`
 130 /// has an error of up to 2ULP
 131 pub fn cos_pi_kernel_f64<Ctx: Context>(ctx: Ctx, x: Ctx::VecF64) -> Ctx::VecF64 {
 132     let x_sq = x * x;
 133     let mut v: Ctx::VecF64 = ctx.make(consts::COSPI_KERNEL_TAYLOR_16.to());
 134     v = v.mul_add_fast(x_sq, ctx.make(consts::COSPI_KERNEL_TAYLOR_14.to()));
 135     v = v.mul_add_fast(x_sq, ctx.make(consts::COSPI_KERNEL_TAYLOR_12.to()));
 136     v = v.mul_add_fast(x_sq, ctx.make(consts::COSPI_KERNEL_TAYLOR_10.to()));
 137     v = v.mul_add_fast(x_sq, ctx.make(consts::COSPI_KERNEL_TAYLOR_8.to()));
 138     v = v.mul_add_fast(x_sq, ctx.make(consts::COSPI_KERNEL_TAYLOR_6.to()));
 139     v = v.mul_add_fast(x_sq, ctx.make(consts::COSPI_KERNEL_TAYLOR_4.to()));
 140     v = v.mul_add_fast(x_sq, ctx.make(consts::COSPI_KERNEL_TAYLOR_2.to()));
 141     v.mul_add_fast(x_sq, ctx.make(consts::COSPI_KERNEL_TAYLOR_0.to()))
 142 }
 143
 144 /// computes `(sin(pi * x), cos(pi * x))`
 145 /// not guaranteed to give correct sign for zero results
 146 /// inherits error from `sin_pi_kernel` and `cos_pi_kernel`
 147 pub fn sin_cos_pi_impl<
 148     Ctx: Context,
 149     VecF: Float<PrimFloat = PrimF> + Make<Context = Ctx>,
 150     PrimF: PrimFloat<BitsType = PrimU>,
 151     PrimU: PrimUInt,
 152     SinPiKernel: FnOnce(Ctx, VecF) -> VecF,
 153     CosPiKernel: FnOnce(Ctx, VecF) -> VecF,
 154 >(
 155     ctx: Ctx,
 156     x: VecF,
 157     sin_pi_kernel: SinPiKernel,
 158     cos_pi_kernel: CosPiKernel,
 159 ) -> (VecF, VecF) {
 160     let two_f: VecF = ctx.make(2.0.to());
 161     let one_half: VecF = ctx.make(0.5.to());
 162     let max_contiguous_integer: VecF =
 163         ctx.make((PrimU::cvt_from(1) << (PrimF::MANTISSA_FIELD_WIDTH + 1.to())).to());
 164     // if `x` is finite and bigger than `max_contiguous_integer`, then x is an even integer
 165     let in_range = x.abs().lt(max_contiguous_integer); // use `lt` so nans are counted as out-of-range
 166     let is_finite = x.is_finite();
 167     let nan: VecF = ctx.make(f32::NAN.to());
 168     let zero_f: VecF = ctx.make(0.to());
 169     let one_f: VecF = ctx.make(1.to());
 170     let zero_i: VecF::SignedBitsType = ctx.make(0.to());
 171     let one_i: VecF::SignedBitsType = ctx.make(1.to());
 172     let two_i: VecF::SignedBitsType = ctx.make(2.to());
 173     let out_of_range_sin = is_finite.select(zero_f, nan);
 174     let out_of_range_cos = is_finite.select(one_f, nan);
 175     let xi = (x * two_f).round();
 176     let xk = x - xi * one_half;
 177     let sk = sin_pi_kernel(ctx, xk);
 178     let ck = cos_pi_kernel(ctx, xk);
 179     let xi = VecF::SignedBitsType::cvt_from(xi);
 180     let bit_0_clear = (xi & one_i).eq(zero_i);
 181     let st = bit_0_clear.select(sk, ck);
 182     let ct = bit_0_clear.select(ck, sk);
 183     let s = (xi & two_i).eq(zero_i).select(st, -st);
 184     let c = ((xi + one_i) & two_i).eq(zero_i).select(ct, -ct);
 185     (
 186         in_range.select(s, out_of_range_sin),
 187         in_range.select(c, out_of_range_cos),
 188     )
 189 }
 190
 191 /// computes `(sin(pi * x), cos(pi * x))`
 192 /// not guaranteed to give correct sign for zero results
 193 /// has an error of up to 2ULP
 194 pub fn sin_cos_pi_f16<Ctx: Context>(ctx: Ctx, x: Ctx::VecF16) -> (Ctx::VecF16, Ctx::VecF16) {
 195     sin_cos_pi_impl(ctx, x, sin_pi_kernel_f16, cos_pi_kernel_f16)
 196 }
 197
 198 /// computes `sin(pi * x)`
 199 /// not guaranteed to give correct sign for zero results
 200 /// has an error of up to 2ULP
 201 pub fn sin_pi_f16<Ctx: Context>(ctx: Ctx, x: Ctx::VecF16) -> Ctx::VecF16 {
 202     sin_cos_pi_f16(ctx, x).0
 203 }
 204
 205 /// computes `cos(pi * x)`
 206 /// not guaranteed to give correct sign for zero results
 207 /// has an error of up to 2ULP
 208 pub fn cos_pi_f16<Ctx: Context>(ctx: Ctx, x: Ctx::VecF16) -> Ctx::VecF16 {
 209     sin_cos_pi_f16(ctx, x).1
 210 }
 211
 212 /// computes `(sin(pi * x), cos(pi * x))`
 213 /// not guaranteed to give correct sign for zero results
 214 /// has an error of up to 2ULP
 215 pub fn sin_cos_pi_f32<Ctx: Context>(ctx: Ctx, x: Ctx::VecF32) -> (Ctx::VecF32, Ctx::VecF32) {
 216     sin_cos_pi_impl(ctx, x, sin_pi_kernel_f32, cos_pi_kernel_f32)
 217 }
 218
 219 /// computes `sin(pi * x)`
 220 /// not guaranteed to give correct sign for zero results
 221 /// has an error of up to 2ULP
 222 pub fn sin_pi_f32<Ctx: Context>(ctx: Ctx, x: Ctx::VecF32) -> Ctx::VecF32 {
 223     sin_cos_pi_f32(ctx, x).0
 224 }
 225
 226 /// computes `cos(pi * x)`
 227 /// not guaranteed to give correct sign for zero results
 228 /// has an error of up to 2ULP
 229 pub fn cos_pi_f32<Ctx: Context>(ctx: Ctx, x: Ctx::VecF32) -> Ctx::VecF32 {
 230     sin_cos_pi_f32(ctx, x).1
 231 }
 232
 233 #[cfg(test)]
 234 mod tests {
 235     use super::*;
 236     use crate::{
 237         f16::F16,
 238         scalar::{Scalar, Value},
 239     };
 240     use std::f64;
 241
 242     struct CheckUlpCallbackArg<F, I> {
 243         distance_in_ulp: I,
 244         x: F,
 245         expected: F,
 246         result: F,
 247     }
 248
 249     #[track_caller]
 250     fn check_ulp<T: PrimFloat>(
 251         x: T,
 252         is_ok: impl Fn(CheckUlpCallbackArg<T, u64>) -> bool,
 253         fn_f16: impl Fn(T) -> T,
 254         fn_reference: impl Fn(f64) -> f64,
 255     ) {
 256         let x_f64: f64 = x.to();
 257         let expected_f64 = fn_reference(x_f64);
 258         let expected: T = expected_f64.to();
 259         let result = fn_f16(x);
 260         if result == expected {
 261             return;
 262         }
 263         if result.is_nan() && expected.is_nan() {
 264             return;
 265         }
 266         let expected_bits: i64 = expected.to_bits().to();
 267         let result_bits: i64 = result.to_bits().to();
 268         let distance_in_ulp = (expected_bits - result_bits).unsigned_abs();
 269         if !result.is_nan()
 270             && !expected.is_nan()
 271             && is_ok(CheckUlpCallbackArg {
 272                 distance_in_ulp,
 273                 x,
 274                 expected,
 275                 result,
 276             })
 277         {
 278             return;
 279         }
 280         panic!(
 281             "error is too big: \
 282                 x = {x:?} {x_bits:#X}, \
 283                 result = {result:?} {result_bits:#X}, \
 284                 expected = {expected:?} {expected_bits:#X}, \
 285                 distance_in_ulp = {distance_in_ulp}",
 286             x = x,
 287             x_bits = x.to_bits(),
 288             result = result,
 289             result_bits = result.to_bits(),
 290             expected = expected,
 291             expected_bits = expected.to_bits(),
 292             distance_in_ulp = distance_in_ulp,
 293         );
 294     }
 295
 296     #[test]
 297     #[cfg_attr(
 298         not(feature = "f16"),
 299         should_panic(expected = "f16 feature is not enabled")
 300     )]
 301     fn test_sin_pi_kernel_f16() {
 302         let check = |x| {
 303             check_ulp(
 304                 x,
 305                 |arg| arg.distance_in_ulp <= if arg.expected == 0.to() { 0 } else { 2 },
 306                 |x| sin_pi_kernel_f16(Scalar, Value(x)).0,
 307                 |x| (f64::consts::PI * x).sin(),
 308             )
 309         };
 310         let quarter = F16::to_bits(0.25f32.to());
 311         for bits in (0..=quarter).rev() {
 312             check(F16::from_bits(bits));
 313             check(-F16::from_bits(bits));
 314         }
 315     }
 316
 317     #[test]
 318     #[cfg_attr(
 319         not(feature = "f16"),
 320         should_panic(expected = "f16 feature is not enabled")
 321     )]
 322     fn test_cos_pi_kernel_f16() {
 323         let check = |x| {
 324             check_ulp(
 325                 x,
 326                 |arg| arg.distance_in_ulp <= 2 && arg.result <= 1.to(),
 327                 |x| cos_pi_kernel_f16(Scalar, Value(x)).0,
 328                 |x| (f64::consts::PI * x).cos(),
 329             )
 330         };
 331         let quarter = F16::to_bits(0.25f32.to());
 332         for bits in (0..=quarter).rev() {
 333             check(F16::from_bits(bits));
 334             check(-F16::from_bits(bits));
 335         }
 336     }
 337
 338     #[test]
 339     #[cfg(feature = "full_tests")]
 340     fn test_sin_pi_kernel_f32() {
 341         let check = |x| {
 342             check_ulp(
 343                 x,
 344                 |arg| arg.distance_in_ulp <= if arg.expected == 0. { 0 } else { 2 },
 345                 |x| sin_pi_kernel_f32(Scalar, Value(x)).0,
 346                 |x| (f64::consts::PI * x).sin(),
 347             )
 348         };
 349         let quarter = 0.25f32.to_bits();
 350         for bits in (0..=quarter).rev() {
 351             check(f32::from_bits(bits));
 352             check(-f32::from_bits(bits));
 353         }
 354     }
 355
 356     #[test]
 357     #[cfg(feature = "full_tests")]
 358     fn test_cos_pi_kernel_f32() {
 359         let check = |x| {
 360             check_ulp(
 361                 x,
 362                 |arg| arg.distance_in_ulp <= 2 && arg.result <= 1.,
 363                 |x| cos_pi_kernel_f32(Scalar, Value(x)).0,
 364                 |x| (f64::consts::PI * x).cos(),
 365             )
 366         };
 367         let quarter = 0.25f32.to_bits();
 368         for bits in (0..=quarter).rev() {
 369             check(f32::from_bits(bits));
 370             check(-f32::from_bits(bits));
 371         }
 372     }
 373
 374     #[test]
 375     #[cfg(feature = "full_tests")]
 376     fn test_sin_pi_kernel_f64() {
 377         let check = |x| {
 378             check_ulp(
 379                 x,
 380                 sin_cos_pi_check_ulp_callback,
 381                 |x| sin_pi_kernel_f64(Scalar, Value(x)).0,
 382                 |x| reference_sin_cos_pi_f64(x).0,
 383             )
 384         };
 385         let quarter = 0.25f32.to_bits();
 386         for bits in (0..=quarter).rev().step_by(1 << 5) {
 387             check(f32::from_bits(bits) as f64);
 388             check(-f32::from_bits(bits) as f64);
 389         }
 390     }
 391
 392     #[test]
 393     #[cfg(feature = "full_tests")]
 394     fn test_cos_pi_kernel_f64() {
 395         let check = |x| {
 396             check_ulp(
 397                 x,
 398                 sin_cos_pi_check_ulp_callback,
 399                 |x| cos_pi_kernel_f64(Scalar, Value(x)).0,
 400                 |x| reference_sin_cos_pi_f64(x).1,
 401             )
 402         };
 403         let quarter = 0.25f32.to_bits();
 404         for bits in (0..=quarter).rev().step_by(1 << 5) {
 405             check(f32::from_bits(bits) as f64);
 406             check(-f32::from_bits(bits) as f64);
 407         }
 408     }
 409
 410     fn sin_cos_pi_check_ulp_callback<F: PrimFloat>(arg: CheckUlpCallbackArg<F, u64>) -> bool {
 411         if arg.x % 0.5.to() == 0.0.to() {
 412             arg.distance_in_ulp == 0
 413         } else {
 414             arg.distance_in_ulp <= 2 && arg.result.abs() <= 1.to()
 415         }
 416     }
 417
 418     #[test]
 419     #[cfg_attr(
 420         not(feature = "f16"),
 421         should_panic(expected = "f16 feature is not enabled")
 422     )]
 423     fn test_sin_pi_f16() {
 424         for bits in 0..=u16::MAX {
 425             check_ulp(
 426                 F16::from_bits(bits),
 427                 sin_cos_pi_check_ulp_callback,
 428                 |x| sin_pi_f16(Scalar, Value(x)).0,
 429                 |x| (f64::consts::PI * x).sin(),
 430             );
 431         }
 432     }
 433
 434     #[test]
 435     #[cfg_attr(
 436         not(feature = "f16"),
 437         should_panic(expected = "f16 feature is not enabled")
 438     )]
 439     fn test_cos_pi_f16() {
 440         for bits in 0..=u16::MAX {
 441             check_ulp(
 442                 F16::from_bits(bits),
 443                 sin_cos_pi_check_ulp_callback,
 444                 |x| cos_pi_f16(Scalar, Value(x)).0,
 445                 |x| (f64::consts::PI * x).cos(),
 446             );
 447         }
 448     }
 449
 450     fn reference_sin_cos_pi_f32(mut v: f64) -> (f64, f64) {
 451         if !v.is_finite() {
 452             return (f64::NAN, f64::NAN);
 453         }
 454         v %= 2.0;
 455         if v >= 1.0 {
 456             v -= 2.0;
 457         } else if v <= -1.0 {
 458             v += 2.0;
 459         }
 460         v *= 2.0;
 461         let part = v.round() as i32;
 462         v -= part as f64;
 463         v *= f64::consts::PI / 2.0;
 464         let (sin, cos) = v.sin_cos();
 465         match part {
 466             0 => (sin, cos),
 467             1 => (cos, -sin),
 468             2 => (-sin, -cos),
 469             -2 => (-sin, -cos),
 470             -1 => (-cos, sin),
 471             _ => panic!("not implemented: part={}", part),
 472         }
 473     }
 474
 475     fn reference_sin_cos_pi_f64(mut v: f64) -> (f64, f64) {
 476         use az::Cast;
 477         use rug::{float::Constant, Float};
 478         if !v.is_finite() {
 479             return (f64::NAN, f64::NAN);
 480         }
 481         v %= 2.0;
 482         if v >= 1.0 {
 483             v -= 2.0;
 484         } else if v <= -1.0 {
 485             v += 2.0;
 486         }
 487         v *= 2.0;
 488         let part = v.round() as i32;
 489         v -= part as f64;
 490         let precision = 100;
 491         let mut v = Float::with_val(precision, v);
 492         let pi = Float::with_val(precision, Constant::Pi);
 493         let pi_2 = pi / 2;
 494         v *= &pi_2;
 495         let cos = pi_2; // just a temp var, value is ignored
 496         let (sin, cos) = v.sin_cos(cos);
 497         let sin: f64 = sin.cast();
 498         let cos: f64 = cos.cast();
 499         match part {
 500             0 => (sin, cos),
 501             1 => (cos, -sin),
 502             2 => (-sin, -cos),
 503             -2 => (-sin, -cos),
 504             -1 => (-cos, sin),
 505             _ => panic!("not implemented: part={}", part),
 506         }
 507     }
 508
 509     macro_rules! test_reference_sin_cos_pi_test_cases {
 510         ($case:expr, $ty:ident) => {
 511             $case($ty::NAN, $ty::NAN, $ty::NAN);
 512             $case($ty::INFINITY, $ty::NAN, $ty::NAN);
 513             $case(-$ty::INFINITY, $ty::NAN, $ty::NAN);
 514             $case(-4., 0., 1.);
 515             $case(
 516                 -3.875,
 517                 0.38268343236508977172845998403039886676134456248563,
 518                 0.92387953251128675612818318939678828682241662586364,
 519             );
 520             $case(
 521                 -3.75,
 522                 0.70710678118654752440084436210484903928483593768847,
 523                 0.70710678118654752440084436210484903928483593768847,
 524             );
 525             $case(
 526                 -3.625,
 527                 0.92387953251128675612818318939678828682241662586364,
 528                 0.38268343236508977172845998403039886676134456248563,
 529             );
 530             $case(-3.5, 1., -0.);
 531             $case(
 532                 -3.375,
 533                 0.92387953251128675612818318939678828682241662586364,
 534                 -0.38268343236508977172845998403039886676134456248563,
 535             );
 536             $case(
 537                 -3.25,
 538                 0.70710678118654752440084436210484903928483593768847,
 539                 -0.70710678118654752440084436210484903928483593768847,
 540             );
 541             $case(
 542                 -3.125,
 543                 0.38268343236508977172845998403039886676134456248563,
 544                 -0.92387953251128675612818318939678828682241662586364,
 545             );
 546             $case(-3., -0., -1.);
 547             $case(
 548                 -2.875,
 549                 -0.38268343236508977172845998403039886676134456248563,
 550                 -0.92387953251128675612818318939678828682241662586364,
 551             );
 552             $case(
 553                 -2.75,
 554                 -0.70710678118654752440084436210484903928483593768847,
 555                 -0.70710678118654752440084436210484903928483593768847,
 556             );
 557             $case(
 558                 -2.625,
 559                 -0.92387953251128675612818318939678828682241662586364,
 560                 -0.38268343236508977172845998403039886676134456248563,
 561             );
 562             $case(-2.5, -1., 0.);
 563             $case(
 564                 -2.375,
 565                 -0.92387953251128675612818318939678828682241662586364,
 566                 0.38268343236508977172845998403039886676134456248563,
 567             );
 568             $case(
 569                 -2.25,
 570                 -0.70710678118654752440084436210484903928483593768847,
 571                 0.70710678118654752440084436210484903928483593768847,
 572             );
 573             $case(
 574                 -2.125,
 575                 -0.38268343236508977172845998403039886676134456248563,
 576                 0.92387953251128675612818318939678828682241662586364,
 577             );
 578             $case(-2., 0., 1.);
 579             $case(
 580                 -1.875,
 581                 0.38268343236508977172845998403039886676134456248563,
 582                 0.92387953251128675612818318939678828682241662586364,
 583             );
 584             $case(
 585                 -1.75,
 586                 0.70710678118654752440084436210484903928483593768847,
 587                 0.70710678118654752440084436210484903928483593768847,
 588             );
 589             $case(
 590                 -1.625,
 591                 0.92387953251128675612818318939678828682241662586364,
 592                 0.38268343236508977172845998403039886676134456248563,
 593             );
 594             $case(-1.5, 1., -0.);
 595             $case(
 596                 -1.375,
 597                 0.92387953251128675612818318939678828682241662586364,
 598                 -0.38268343236508977172845998403039886676134456248563,
 599             );
 600             $case(
 601                 -1.25,
 602                 0.70710678118654752440084436210484903928483593768847,
 603                 -0.70710678118654752440084436210484903928483593768847,
 604             );
 605             $case(
 606                 -1.125,
 607                 0.38268343236508977172845998403039886676134456248563,
 608                 -0.92387953251128675612818318939678828682241662586364,
 609             );
 610             $case(-1., -0., -1.);
 611             $case(
 612                 -0.875,
 613                 -0.38268343236508977172845998403039886676134456248563,
 614                 -0.92387953251128675612818318939678828682241662586364,
 615             );
 616             $case(
 617                 -0.75,
 618                 -0.70710678118654752440084436210484903928483593768847,
 619                 -0.70710678118654752440084436210484903928483593768847,
 620             );
 621             $case(
 622                 -0.625,
 623                 -0.92387953251128675612818318939678828682241662586364,
 624                 -0.38268343236508977172845998403039886676134456248563,
 625             );
 626             $case(-0.5, -1., 0.);
 627             $case(
 628                 -0.375,
 629                 -0.92387953251128675612818318939678828682241662586364,
 630                 0.38268343236508977172845998403039886676134456248563,
 631             );
 632             $case(
 633                 -0.25,
 634                 -0.70710678118654752440084436210484903928483593768847,
 635                 0.70710678118654752440084436210484903928483593768847,
 636             );
 637             $case(
 638                 -0.125,
 639                 -0.38268343236508977172845998403039886676134456248563,
 640                 0.92387953251128675612818318939678828682241662586364,
 641             );
 642             $case(0., 0., 1.);
 643             $case(
 644                 0.125,
 645                 0.38268343236508977172845998403039886676134456248563,
 646                 0.92387953251128675612818318939678828682241662586364,
 647             );
 648             $case(
 649                 0.25,
 650                 0.70710678118654752440084436210484903928483593768847,
 651                 0.70710678118654752440084436210484903928483593768847,
 652             );
 653             $case(
 654                 0.375,
 655                 0.92387953251128675612818318939678828682241662586364,
 656                 0.38268343236508977172845998403039886676134456248563,
 657             );
 658             $case(0.5, 1., 0.);
 659             $case(
 660                 0.625,
 661                 0.92387953251128675612818318939678828682241662586364,
 662                 -0.38268343236508977172845998403039886676134456248563,
 663             );
 664             $case(
 665                 0.75,
 666                 0.70710678118654752440084436210484903928483593768847,
 667                 -0.70710678118654752440084436210484903928483593768847,
 668             );
 669             $case(
 670                 0.875,
 671                 0.38268343236508977172845998403039886676134456248563,
 672                 -0.92387953251128675612818318939678828682241662586364,
 673             );
 674             $case(1., 0., -1.);
 675             $case(
 676                 1.125,
 677                 -0.38268343236508977172845998403039886676134456248563,
 678                 -0.92387953251128675612818318939678828682241662586364,
 679             );
 680             $case(
 681                 1.25,
 682                 -0.70710678118654752440084436210484903928483593768847,
 683                 -0.70710678118654752440084436210484903928483593768847,
 684             );
 685             $case(
 686                 1.375,
 687                 -0.92387953251128675612818318939678828682241662586364,
 688                 -0.38268343236508977172845998403039886676134456248563,
 689             );
 690             $case(1.5, -1., -0.);
 691             $case(
 692                 1.625,
 693                 -0.92387953251128675612818318939678828682241662586364,
 694                 0.38268343236508977172845998403039886676134456248563,
 695             );
 696             $case(
 697                 1.75,
 698                 -0.70710678118654752440084436210484903928483593768847,
 699                 0.70710678118654752440084436210484903928483593768847,
 700             );
 701             $case(
 702                 1.875,
 703                 -0.38268343236508977172845998403039886676134456248563,
 704                 0.92387953251128675612818318939678828682241662586364,
 705             );
 706             $case(2., -0., 1.);
 707             $case(
 708                 2.125,
 709                 0.38268343236508977172845998403039886676134456248563,
 710                 0.92387953251128675612818318939678828682241662586364,
 711             );
 712             $case(
 713                 2.25,
 714                 0.70710678118654752440084436210484903928483593768847,
 715                 0.70710678118654752440084436210484903928483593768847,
 716             );
 717             $case(
 718                 2.375,
 719                 0.92387953251128675612818318939678828682241662586364,
 720                 0.38268343236508977172845998403039886676134456248563,
 721             );
 722             $case(2.5, 1., 0.);
 723             $case(
 724                 2.625,
 725                 0.92387953251128675612818318939678828682241662586364,
 726                 -0.38268343236508977172845998403039886676134456248563,
 727             );
 728             $case(
 729                 2.75,
 730                 0.70710678118654752440084436210484903928483593768847,
 731                 -0.70710678118654752440084436210484903928483593768847,
 732             );
 733             $case(
 734                 2.875,
 735                 0.38268343236508977172845998403039886676134456248563,
 736                 -0.92387953251128675612818318939678828682241662586364,
 737             );
 738             $case(3., 0., -1.);
 739             $case(
 740                 3.125,
 741                 -0.38268343236508977172845998403039886676134456248563,
 742                 -0.92387953251128675612818318939678828682241662586364,
 743             );
 744             $case(
 745                 3.25,
 746                 -0.70710678118654752440084436210484903928483593768847,
 747                 -0.70710678118654752440084436210484903928483593768847,
 748             );
 749             $case(
 750                 3.375,
 751                 -0.92387953251128675612818318939678828682241662586364,
 752                 -0.38268343236508977172845998403039886676134456248563,
 753             );
 754             $case(3.5, -1., -0.);
 755             $case(
 756                 3.625,
 757                 -0.92387953251128675612818318939678828682241662586364,
 758                 0.38268343236508977172845998403039886676134456248563,
 759             );
 760             $case(
 761                 3.75,
 762                 -0.70710678118654752440084436210484903928483593768847,
 763                 0.70710678118654752440084436210484903928483593768847,
 764             );
 765             $case(
 766                 3.875,
 767                 -0.38268343236508977172845998403039886676134456248563,
 768                 0.92387953251128675612818318939678828682241662586364,
 769             );
 770             $case(4., -0., 1.);
 771         };
 772     }
 773
 774     #[test]
 775     fn test_reference_sin_cos_pi_f32() {
 776         fn approx_same(a: f32, b: f32) -> bool {
 777             if a.is_finite() && b.is_finite() {
 778                 (a - b).abs() < 1e-6
 779             } else {
 780                 a == b || (a.is_nan() && b.is_nan())
 781             }
 782         }
 783         #[track_caller]
 784         fn case(x: f32, expected_sin: f32, expected_cos: f32) {
 785             let (ref_sin, ref_cos) = reference_sin_cos_pi_f32(x as f64);
 786             assert!(
 787                 approx_same(ref_sin as f32, expected_sin)
 788                     && approx_same(ref_cos as f32, expected_cos),
 789                 "case failed: x={x}, expected_sin={expected_sin}, expected_cos={expected_cos}, ref_sin={ref_sin}, ref_cos={ref_cos}",
 790                 x=x,
 791                 expected_sin=expected_sin,
 792                 expected_cos=expected_cos,
 793                 ref_sin=ref_sin,
 794                 ref_cos=ref_cos,
 795             );
 796         }
 797         test_reference_sin_cos_pi_test_cases!(case, f32);
 798     }
 799
 800     #[test]
 801     fn test_reference_sin_cos_pi_f64() {
 802         fn same(a: f64, b: f64) -> bool {
 803             if a.is_finite() && b.is_finite() {
 804                 a == b
 805             } else {
 806                 a == b || (a.is_nan() && b.is_nan())
 807             }
 808         }
 809         #[track_caller]
 810         fn case(x: f64, expected_sin: f64, expected_cos: f64) {
 811             let (ref_sin, ref_cos) = reference_sin_cos_pi_f64(x);
 812             assert!(
 813                 same(ref_sin, expected_sin) && same(ref_cos, expected_cos),
 814                 "case failed: x={x}, expected_sin={expected_sin}, expected_cos={expected_cos}, ref_sin={ref_sin}, ref_cos={ref_cos}",
 815                 x=x,
 816                 expected_sin=expected_sin,
 817                 expected_cos=expected_cos,
 818                 ref_sin=ref_sin,
 819                 ref_cos=ref_cos,
 820             );
 821         }
 822         test_reference_sin_cos_pi_test_cases!(case, f64);
 823     }
 824
 825     #[test]
 826     #[cfg(feature = "full_tests")]
 827     fn test_sin_pi_f32() {
 828         for bits in 0..=u32::MAX {
 829             check_ulp(
 830                 f32::from_bits(bits),
 831                 sin_cos_pi_check_ulp_callback,
 832                 |x| sin_pi_f32(Scalar, Value(x)).0,
 833                 |x| reference_sin_cos_pi_f32(x).0,
 834             );
 835         }
 836     }
 837
 838     #[test]
 839     #[cfg(feature = "full_tests")]
 840     fn test_cos_pi_f32() {
 841         for bits in 0..=u32::MAX {
 842             check_ulp(
 843                 f32::from_bits(bits),
 844                 sin_cos_pi_check_ulp_callback,
 845                 |x| cos_pi_f32(Scalar, Value(x)).0,
 846                 |x| reference_sin_cos_pi_f32(x).1,
 847             );
 848         }
 849     }
 850 }