Maybe small improvement for accuracy for cos for regular angles

[mhvk]

@seiko2plus - Came here out of curiosity from numpy/numpy#29699 (comment), to see how things get done. For cos, the following transformation is done before any real calculation:

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

Lines 35 to 43 in 1acd145

	const V x_abs = And(abs_mask, x);
	const V x_sign = AndNot(x_abs, x);
	// Transform cosine to sine using identity: cos(x) = sin(x + π/2)
	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);
	}

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. Instead, one can move those angles to near zero and thus get optimal precision near both crossings by using that cos(x) = sin(pi/2 - |x|) , i.e. do something like,

  // Transform cosine to sine using identity cos(x) = sin(π/2 - |x|)
  // where |x| ensures optimal precision near both -pi/2 and pi/2.
  const V half_pi = Set(d, data::kHalfPi<T>);
  V x_trans = x;
  if constexpr (OP == Operation::kCos) {
    x_trans = And(abs_mask, x);
    x_trans = Sub(half_pi, x_trans);
  }
  const V x_abs = And(abs_mask, x_trans);
  const V x_sign = AndNot(x_abs, x_trans);

[seiko2plus]

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)

[mhvk]

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!