src/algorithms/trig_pi.rs

   1 use crate::traits::{Context, ConvertTo, Float};
   2
   3 mod consts {
   4     #![allow(clippy::excessive_precision)]
   5
   6     /// coefficients of taylor series for `sin(pi * x)` centered at `0`
   7     /// generated using:
   8     /// ```maxima,text
   9     /// fpprec:50$
  10     /// sinpi: bfloat(taylor(sin(%pi*x),x,0,19))$
  11     /// for i: 1 step 2 thru 19 do
  12     ///     printf(true, "pub(crate) const SINPI_KERNEL_TAYLOR_~d: f64 = ~a;~%", i, ssubst("e", "b", string(coeff(sinpi, x, i))))$
  13     /// ```
  14     pub(crate) const SINPI_KERNEL_TAYLOR_1: f64 =
  15         3.1415926535897932384626433832795028841971693993751e0;
  16     pub(crate) const SINPI_KERNEL_TAYLOR_3: f64 =
  17         -5.1677127800499700292460525111835658670375480943142e0;
  18     pub(crate) const SINPI_KERNEL_TAYLOR_5: f64 =
  19         2.550164039877345443856177583695296720669172555234e0;
  20     pub(crate) const SINPI_KERNEL_TAYLOR_7: f64 =
  21         -5.9926452932079207688773938354604004601536358636814e-1;
  22     pub(crate) const SINPI_KERNEL_TAYLOR_9: f64 =
  23         8.2145886611128228798802365523698344807837460797753e-2;
  24     pub(crate) const SINPI_KERNEL_TAYLOR_11: f64 =
  25         -7.370430945714350777259089957290781501211638236021e-3;
  26     pub(crate) const SINPI_KERNEL_TAYLOR_13: f64 =
  27         4.6630280576761256442062891447027174382819981361599e-4;
  28     pub(crate) const SINPI_KERNEL_TAYLOR_15: f64 =
  29         -2.1915353447830215827384652057094188859248708765956e-5;
  30     pub(crate) const SINPI_KERNEL_TAYLOR_17: f64 =
  31         7.9520540014755127847832068624575890327682459384282e-7;
  32     pub(crate) const SINPI_KERNEL_TAYLOR_19: f64 =
  33         -2.2948428997269873110203872385571587856074785581088e-8;
  34
  35     /// coefficients of taylor series for `cos(pi * x)` centered at `0`
  36     /// generated using:
  37     /// ```maxima,text
  38     /// fpprec:50$
  39     /// cospi: bfloat(taylor(cos(%pi*x),x,0,18))$
  40     /// for i: 0 step 2 thru 18 do
  41     ///     printf(true, "pub(crate) const COSPI_KERNEL_TAYLOR_~d: f64 = ~a;~%", i, ssubst("e", "b", string(coeff(cospi, x, i))))$
  42     /// ```
  43     pub(crate) const COSPI_KERNEL_TAYLOR_0: f64 = 1.0e0;
  44     pub(crate) const COSPI_KERNEL_TAYLOR_2: f64 =
  45         -4.9348022005446793094172454999380755676568497036204e0;
  46     pub(crate) const COSPI_KERNEL_TAYLOR_4: f64 =
  47         4.0587121264167682181850138620293796354053160696952e0;
  48     pub(crate) const COSPI_KERNEL_TAYLOR_6: f64 =
  49         -1.3352627688545894958753047828505831928711354556681e0;
  50     pub(crate) const COSPI_KERNEL_TAYLOR_8: f64 =
  51         2.3533063035889320454187935277546542154506893530856e-1;
  52     pub(crate) const COSPI_KERNEL_TAYLOR_10: f64 =
  53         -2.5806891390014060012598294252898849657186441048147e-2;
  54     pub(crate) const COSPI_KERNEL_TAYLOR_12: f64 =
  55         1.9295743094039230479033455636859576401684718150003e-3;
  56     pub(crate) const COSPI_KERNEL_TAYLOR_14: f64 =
  57         -1.0463810492484570711801672835223932761029733149091e-4;
  58     pub(crate) const COSPI_KERNEL_TAYLOR_16: f64 =
  59         4.3030695870329470072978237149669233008960901556009e-6;
  60     pub(crate) const COSPI_KERNEL_TAYLOR_18: f64 =
  61         -1.387895246221377211446808750399309343777037849978e-7;
  62 }
  63
  64 /// computes `sin(pi * x)` for `-0.25 <= x <= 0.25`
  65 /// not guaranteed to give correct sign for zero result
  66 /// has an error of up to 2ULP
  67 pub fn sin_pi_kernel_f16<Ctx: Context>(ctx: Ctx, x: Ctx::VecF16) -> Ctx::VecF16 {
  68     let x_sq = x * x;
  69     let mut v: Ctx::VecF16 = ctx.make(consts::SINPI_KERNEL_TAYLOR_5.to());
  70     v = v.mul_add_fast(x_sq, ctx.make(consts::SINPI_KERNEL_TAYLOR_3.to()));
  71     v = v.mul_add_fast(x_sq, ctx.make(consts::SINPI_KERNEL_TAYLOR_1.to()));
  72     v * x
  73 }
  74
  75 /// computes `cos(pi * x)` for `-0.25 <= x <= 0.25`
  76 /// has an error of up to 2ULP
  77 pub fn cos_pi_kernel_f16<Ctx: Context>(ctx: Ctx, x: Ctx::VecF16) -> Ctx::VecF16 {
  78     let x_sq = x * x;
  79     let mut v: Ctx::VecF16 = ctx.make(consts::COSPI_KERNEL_TAYLOR_4.to());
  80     v = v.mul_add_fast(x_sq, ctx.make(consts::COSPI_KERNEL_TAYLOR_2.to()));
  81     v.mul_add_fast(x_sq, ctx.make(consts::COSPI_KERNEL_TAYLOR_0.to()))
  82 }
  83
  84 /// computes `(sin(pi * x), cos(pi * x))`
  85 /// not guaranteed to give correct sign for zero results
  86 /// has an error of up to 2ULP
  87 pub fn sin_cos_pi_f16<Ctx: Context>(_ctx: Ctx, _x: Ctx::VecF16) -> (Ctx::VecF16, Ctx::VecF16) {
  88     todo!()
  89 }
  90
  91 #[cfg(test)]
  92 mod tests {
  93     use super::*;
  94     use crate::{f16::F16, scalar::Scalar};
  95     use std::f64;
  96
  97     struct CheckUlpCallbackArg<F, I> {
  98         distance_in_ulp: I,
  99         expected: F,
 100         result: F,
 101     }
 102
 103     #[track_caller]
 104     fn check_ulp_f16(
 105         x: F16,
 106         is_ok: impl Fn(CheckUlpCallbackArg<F16, i16>) -> bool,
 107         fn_f16: impl Fn(F16) -> F16,
 108         fn_f64: impl Fn(f64) -> f64,
 109     ) {
 110         let x_f64: f64 = x.to();
 111         let expected_f64 = fn_f64(x_f64);
 112         let expected: F16 = expected_f64.to();
 113         let result = fn_f16(x);
 114         if result == expected {
 115             return;
 116         }
 117         let distance_in_ulp = (expected.to_bits() as i16).wrapping_sub(result.to_bits() as i16);
 118         if is_ok(CheckUlpCallbackArg {
 119             distance_in_ulp,
 120             expected,
 121             result,
 122         }) {
 123             return;
 124         }
 125         panic!(
 126             "error is too big: \
 127                 x = {x:?} {x_bits:#X}, \
 128                 result = {result:?} {result_bits:#X}, \
 129                 expected = {expected:?} {expected_bits:#X}",
 130             x = x,
 131             x_bits = x.to_bits(),
 132             result = result,
 133             result_bits = result.to_bits(),
 134             expected = expected,
 135             expected_bits = expected.to_bits(),
 136         );
 137     }
 138
 139     #[test]
 140     #[cfg_attr(
 141         not(feature = "f16"),
 142         should_panic(expected = "f16 feature is not enabled")
 143     )]
 144     fn test_sin_pi_kernel_f16() {
 145         let check = |x| {
 146             check_ulp_f16(
 147                 x,
 148                 |arg| arg.distance_in_ulp <= if arg.expected == 0.to() { 0 } else { 2 },
 149                 |x| sin_pi_kernel_f16(Scalar, x),
 150                 |x| (f64::consts::PI * x).sin(),
 151             )
 152         };
 153         let quarter = F16::to_bits(0.25f32.to());
 154         for bits in (0..=quarter).rev() {
 155             check(F16::from_bits(bits));
 156             check(-F16::from_bits(bits));
 157         }
 158     }
 159
 160     #[test]
 161     #[cfg_attr(
 162         not(feature = "f16"),
 163         should_panic(expected = "f16 feature is not enabled")
 164     )]
 165     fn test_cos_pi_kernel_f16() {
 166         let check = |x| {
 167             check_ulp_f16(
 168                 x,
 169                 |arg| arg.distance_in_ulp <= 2 && arg.result <= 1.to(),
 170                 |x| cos_pi_kernel_f16(Scalar, x),
 171                 |x| (f64::consts::PI * x).cos(),
 172             )
 173         };
 174         let quarter = F16::to_bits(0.25f32.to());
 175         for bits in (0..=quarter).rev() {
 176             check(F16::from_bits(bits));
 177             check(-F16::from_bits(bits));
 178         }
 179     }
 180 }