Fix degree trig functions for large arguments

include/xsf/cephes/sindg.h

@@ -37,7 +37,7 @@
  * ERROR MESSAGES:
  *
  *   message           condition        value returned
- * sindg total loss   x > 1.0e14 (IEEE)     0.0
+ * sindg no result     x = +-INFINITY        0.0
  *
  */
 /*							cosdg.c
@@ -98,8 +98,6 @@ namespace cephes {
                                      -2.48015872936186303776E-5, 1.38888888888806666760E-3,  -4.16666666666666348141E-2,
                                      4.99999999999999999798E-1};
 
-        constexpr double sindg_lossth = 1.0e14;
-
     } // namespace detail
 
     XSF_HOST_DEVICE inline double sindg(double x) {
@@ -113,11 +111,16 @@ namespace cephes {
             sign = -1;
         }
 
-        if (x > detail::sindg_lossth) {
+        if (std::isinf(x)) {
             set_error("sindg", SF_ERROR_NO_RESULT, NULL);
             return (0.0);
         }
 
+        /* Reduce modulo a full turn. Exact: 360 is representable and fmod is
+         * exact. The octant index below is taken modulo 8, so dropping whole
+         * turns does not change the result. */
+        x = std::fmod(x, 360.0);
+
         y = std::floor(x / 45.0); /* integer part of x/M_PI_4 */
 
         /* strip high bits of integer part to prevent integer overflow */
@@ -163,11 +166,16 @@ namespace cephes {
         if (x < 0)
             x = -x;
 
-        if (x > detail::sindg_lossth) {
+        if (std::isinf(x)) {
             set_error("cosdg", SF_ERROR_NO_RESULT, NULL);
             return (0.0);
         }
 
+        /* Reduce modulo a full turn. Exact: 360 is representable and fmod is
+         * exact. The octant index below is taken modulo 8, so dropping whole
+         * turns does not change the result. */
+        x = std::fmod(x, 360.0);
+
         y = std::floor(x / 45.0);
         z = std::ldexp(y, -4);
         z = std::floor(z);        /* integer part of y/8 */

include/xsf/cephes/tandg.h

@@ -35,7 +35,7 @@
  * ERROR MESSAGES:
  *
  *   message         condition          value returned
- * tandg total loss   x > 1.0e14 (IEEE)     0.0
+ * tandg no result    x = +-INFINITY        0.0
  * tandg singularity  x = 180 k  +  90     INFINITY
  */
 /*							cotdg.c
@@ -64,7 +64,7 @@
  * ERROR MESSAGES:
  *
  *   message         condition          value returned
- * cotdg total loss   x > 1.0e14 (IEEE)     0.0
+ * cotdg no result    x = +-INFINITY        0.0
  * cotdg singularity  x = 180 k            INFINITY
  */
 
@@ -82,7 +82,6 @@ namespace xsf {
 namespace cephes {
 
     namespace detail {
-        constexpr double tandg_lossth = 1.0e14;
 
         XSF_HOST_DEVICE inline double tancot(double xx, int cotflg) {
             double x;
@@ -97,11 +96,16 @@ namespace cephes {
                 sign = 1;
             }
 
-            if (x > detail::tandg_lossth) {
-                set_error("tandg", SF_ERROR_NO_RESULT, NULL);
+            if (std::isinf(x)) {
+                set_error((cotflg ? "cotdg" : "tandg"), SF_ERROR_NO_RESULT, NULL);
                 return 0.0;
             }
 
+            /* Reduce modulo a full turn. Exact: 360 is representable and fmod
+             * is exact. A turn is an even number of half turns, so the parity
+             * of k below is unchanged. */
+            x = std::fmod(x, 360.0);
+
             /* modulo 180 */
             double k = std::floor(x / 180.0);
             x = x - 180.0 * k;

tests/xsf_tests/test_trig.cpp

@@ -99,3 +99,37 @@ TEST_CASE("cotdg IEEE infinity sign scipy/20731", "[cotdg][xsf_tests]") {
         }
     }
 }
+
+/* These returned 0.0 for |x| > 1e14 because `x - 45 * floor(x / 45)` loses
+ * accuracy there. An exact fmod by 360 first removes the limit. scipy/scipy#20723
+ */
+TEST_CASE("degree trig large arguments scipy/20723", "[sindg][cosdg][tandg][cotdg][xsf_tests]") {
+    /* {x, x mod 360}. The reduction is exact, so the two must agree exactly. */
+    std::vector<std::pair<double, double>> cases{
+        {1e14 + 1.0, 281.0}, /* just past the old 1e14 cutoff */
+        {2e14, 200.0},       /* the example in the issue */
+        {1e15, 280.0},
+        {9007199254740991.0, 31.0}, /* 2**53 - 1, the last odd integer */
+        {1e16, 280.0},              /* beyond 2**53 */
+        {1e17, 280.0},
+        {1e18, 280.0},
+        {1e300, 0.0},
+        {360000000000030.0, 30.0},
+        {360000000000090.0, 90.0}, /* pole of tandg */
+        {360000000000123.0, 123.0},
+        {360000000000180.0, 180.0}, /* pole of cotdg */
+    };
+    for (auto [x, reduced] : cases) {
+        for (double sign : {1.0, -1.0}) {
+            double y = sign * x;
+            double y_reduced = sign * reduced;
+            CAPTURE(x, reduced, sign);
+            REQUIRE(xsf::sindg(y) == xsf::sindg(y_reduced));
+            REQUIRE(xsf::cosdg(y) == xsf::cosdg(y_reduced));
+            REQUIRE(xsf::tandg(y) == xsf::tandg(y_reduced));
+            REQUIRE(xsf::cotdg(y) == xsf::cotdg(y_reduced));
+        }
+    }
+    /* sindg(2e14) used to be 0.0; sindg(2e14 mod 360) = sindg(200) */
+    REQUIRE(std::abs(xsf::sindg(2e14) + 0.3420201433256687) < 1e-15);
+}