BUG: Give low-accuracy f64 cosine its own minimax fit

Cosine reduces around (N+0.5)pi, so its polynomial needs its own
coefficient set rather than reusing the sine fit to bitexact with SVML.
low-inl.h now selects cosine-tuned c3..c15 for the kCos path.
Also raise the low-accuracy cosine large-argument threshold to
0x1.921fb2200366fp+24 (largest double with trunc(|x|/pi) <= 8388606)
so it stays on the fast path up to the real boundary instead of
falling back to extended precision at 2^24.

Author: seiko2plus

npsr/trig/inl.h

@@ -49,7 +49,7 @@ namespace npsr::HWY_NAMESPACE::trig {
  * - Double: |x| > 2^24 (16,777,216 - where 53-bit mantissa loses precision)
  */
 template <Operation OP, typename Prec, typename V>
-NPSR_INTRIN V Trig(Prec &prec, V x) {
+NPSR_INTRIN V Trig(Prec& prec, V x) {
   using namespace hwy::HWY_NAMESPACE;
   constexpr bool kIsSingle = std::is_same_v<TFromV<V>, float>;
   const DFromV<V> d;
@@ -81,8 +81,14 @@ NPSR_INTRIN V Trig(Prec &prec, V x) {
     // Thresholds chosen based on when standard reduction loses accuracy:
     // - Float: 10,000 is conservative but ensures < 1 ULP error
     // - Double: 2^24 is where mantissa can't represent x and x-2π distinctly
+    // - Double low-accuracy cosine:  stays on its fast path
+    //   until trunc(|x|/π) > 8388606; this is the largest double below that
+    //   boundary (~8388607π)
+    constexpr bool kIsLowCos =
+        Prec::kLowAccuracy && OP == trig::Operation::kCos;
+    constexpr double kLargeF64 = kIsLowCos ? 0x1.921fb2200366fp+24 : 16777216.0;
     auto has_large_arg =
-        And(Gt(Abs(x), Set(d, kIsSingle ? 10000.0f : 16777216.0)), is_finite);
+        And(Gt(Abs(x), Set(d, kIsSingle ? 10000.0f : kLargeF64)), is_finite);
 
     // Extended precision is expensive, only use when necessary
     if (HWY_UNLIKELY(!AllFalse(d, has_large_arg))) {
@@ -121,7 +127,7 @@ namespace npsr::HWY_NAMESPACE {
  * ```
  */
 template <typename Prec, typename V>
-NPSR_INTRIN V Sin(Prec &prec, V x) {
+NPSR_INTRIN V Sin(Prec& prec, V x) {
   return trig::Trig<trig::Operation::kSin>(prec, x);
 }
 
@@ -143,7 +149,7 @@ NPSR_INTRIN V Sin(Prec &prec, V x) {
  * ```
  */
 template <typename Prec, typename V>
-NPSR_INTRIN V Cos(Prec &prec, V x) {
+NPSR_INTRIN V Cos(Prec& prec, V x) {
   return trig::Trig<trig::Operation::kCos>(prec, x);
 }
 

npsr/trig/low-inl.h

@@ -54,15 +54,17 @@ NPSR_INTRIN V PolyLow(V r, V r2) {
 template <Operation OP, typename V, HWY_IF_F64(TFromV<V>)>
 NPSR_INTRIN V PolyLow(V r, V r2) {
   using namespace hn;
-
   const DFromV<V> d;
-  const V c15 = Set(d, -0x1.9f1517e9f65fp-41);
-  const V c13 = Set(d, 0x1.60e6bee01d83ep-33);
-  const V c11 = Set(d, -0x1.ae6355aaa4a53p-26);
-  const V c9 = Set(d, 0x1.71de3806add1ap-19);
-  const V c7 = Set(d, -0x1.a01a019a659ddp-13);
-  const V c5 = Set(d, 0x1.111111110a573p-7);
-  const V c3 = Set(d, -0x1.55555555554a8p-3);
+  // Sine and cosine use distinct minimax fits of sin(r): the cosine path
+  // reduces around (N + 0.5)*pi, so its coefficient set differs 
+  constexpr bool kCos = OP == Operation::kCos;
+  const V c15 = Set(d, kCos ? -0x1.9f0d60811aac8p-41 : -0x1.9f1517e9f65fp-41);
+  const V c13 = Set(d, kCos ? 0x1.60e6857a2f220p-33 : 0x1.60e6bee01d83ep-33);
+  const V c11 = Set(d, kCos ? -0x1.ae63546002231p-26 : -0x1.ae6355aaa4a53p-26);
+  const V c9 = Set(d, kCos ? 0x1.71de38030fea0p-19 : 0x1.71de3806add1ap-19);
+  const V c7 = Set(d, kCos ? -0x1.a01a019a5b87bp-13 : -0x1.a01a019a659ddp-13);
+  const V c5 = Set(d, kCos ? 0x1.111111110a4a8p-7 : 0x1.111111110a573p-7);
+  const V c3 = Set(d, kCos ? -0x1.55555555554a7p-3 : -0x1.55555555554a8p-3);
   V poly = MulAdd(c15, r2, c13);
   poly = MulAdd(r2, poly, c11);
   poly = MulAdd(r2, poly, c9);