src/algorithms/trig_pi.rs

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