Maybe small improvement for accuracy for cos for regular angles #6

Activity

1. seiko2plus commented on Dec 10, 2025

   

   seiko2plus

   on Dec 10, 2025

   Member

   More actions

   Came here out of curiosity from numpy/numpy#29699 (comment),

   Welcome!

   For cos, the following transformation is done before any real calculation:

   Yes, but to be clear, this function is only for single-precision with ~1 ULP error and is valid for argument ranges |x| <= 10000.0. Larger argument reduction is handled by the slow extended prec path:
   

   numpy-simd-routines/npsr/trig/extended-inl.h

   Lines 15 to 18 in 1acd145

   	template <Operation OP, class V>
   	NPSR_INTRIN V Extended(V x) {
   	using namespace hn;
   	namespace data = ::npsr::trig::data;

   This would seem to loose some precision for the common case where angle input will be between -pi and pi, by moving those angles to pi/2 to 3*pi/2.

   It wouldn’t matter, since it will eventually be rounded to calculate the N' (adjusted multiples of π). x_trans is not used to calculate the remainder; we use |x|, so we need to adjust N for cosine. Can we still use cos(x) = sin(π/2 − |x|)? Yes, but it will require an extra sign-flip operation. Here is an example:

   diff --git a/npsr/trig/high-inl.h b/npsr/trig/high-inl.h
   index 68cdd9f..fd570c1 100644
   --- a/npsr/trig/high-inl.h
   +++ b/npsr/trig/high-inl.h
   @@ -39,7 +39,7 @@ NPSR_INTRIN V High(V x) {
      const V half_pi = Set(d, data::kHalfPi<T>);
      V x_trans = x_abs;
      if constexpr (OP == Operation::kCos) {
   -    x_trans = Add(x_abs, half_pi);
   +    x_trans = Sub(x_abs, half_pi);
      }
      // check zero input/subnormal for cosine (cos(~0) = 1)
      const auto is_cos_near_zero = Eq(x_trans, half_pi);
   @@ -51,11 +51,11 @@ NPSR_INTRIN V High(V x) {
    
      // Adjust quotient for cosine (accounts for π/2 phase shift)
      if constexpr (OP == Operation::kCos) {
   -    // For cosine, we computed N = round((x + π/2)/π) but need N' for x:
   -    //   N = round((x + π/2)/π) = round(x/π + 0.5)
   -    // This is often 1 more than round(x/π), so we subtract 0.5:
   -    //   N' = N - 0.5
   -    n = Sub(n, Set(d, 0.5f));
   +    // For cosine, we computed N = round((π/2 - |x|)/π) but need N' for x:
   +    //   N = round((π/2 - |x|)/π) = round(x/π - 0.5)
   +    // This is often 1 less than round(x/π), so we add 0.5:
   +    //   N' = N + 0.5
   +    n = Add(n, Set(d, 0.5f));
      }
      auto WideCal = [](const VW &nh, const VW &xh_abs) -> VW {
        const DFromV<VW> dw;
   @@ -85,6 +85,10 @@ NPSR_INTRIN V High(V x) {
      poly = Xor(poly,
                 BitCast(d, ShiftLeft<sizeof(T) * 8 - 1>(BitCast(du, n_biased))));
      if constexpr (OP == Operation::kCos) {
   +    // We computed sin(x - π/2) = -cos(x).
   +    // We need to negate the result to get cos(x).
   +    // Cosine is even, so we do NOT apply x_sign.
   +    poly = Xor(poly, Set(d, -0.0f));
        poly = IfThenElse(is_cos_near_zero, Set(d, 1.0f), poly);
      } else {
        // Restore original sign for sine (odd function)
   
2. mhvk commented on Dec 10, 2025

   

   mhvk

   on Dec 10, 2025

   Author

   More actions

   Yes, agreed it would cost another sign flip or absolute value operation (at least if one wants to keep the precision symmetric). Not quite sure whether your suggestion would work, since now x_trans can be negative. If that is OK, then in my suggestion there would be no need to take the second absolute value (or to calculate the sign?).

   Anyway, this probably all is in the noise!