src/algorithms/trig_pi.rs

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