Fix degree trig functions for large arguments

[Kayvan-Zahiri]

Fixes scipy/scipy#20723.

sindg, cosdg, tandg and cotdg bail out with return 0.0 above 1e14, so
sindg(2e14) gives 0.0 instead of -sin(20 deg).

Replaced the bail-out with an exact std::fmod(x, 360.0) reduction. 360 is
exactly representable and the remainder is exact, so the reduced argument is the
true one.

• sindg(2e14) -> -0.3420201433256687
• bit-for-bit unchanged for |x| <= 1e14
• correct for every finite argument up to 1e300, measured 0.64 ulp
• +-inf still returns 0.0 with SF_ERROR_NO_RESULT, unchanged; NaN still NaN

Test added to tests/xsf_tests/test_trig.cpp: 12 argument pairs across all four
functions, both signs, asserting the reduced and unreduced arguments agree.

The scipy-side tests need this merged and the submodule bumped first, so I will
open that separately.

Note the thread is not settled on whether to fix this or document it. I went with
the reduction because it turns out to be exact, but happy to be told otherwise.

AI disclosure: an AI assistant helped investigate and draft this. I compiled the
headers standalone and ran the tests, and can explain the change.

[j-bowhay]

Some of the comments are useful but there is a lot of llm waffle, please reduce this

[j-bowhay]

Could you also comment on how to decided on this approach Vs the ideas discussed in the original issue

[Kayvan-Zahiri]

Comments trimmed.

I took @fancidev's fmod suggestion over the degree-to-radian helper: converting first multiplies by pi/180 before reducing and loses accuracy, e.g. sin(radians(2e14)) gives -0.34204442 against -0.34202014 reducing first. 360 is exactly representable and pi is not, so degrees get an exact reduction for free.

@nickodell is right that large x carries few meaningful digits, but that is the input rather than the reduction. fmod adds no error, so this returns the correctly rounded result for the double it was handed, as libm's sin does. That is also why I did not just document the limit.

[nickodell]

fmod exact?

Replaced the bail-out with an exact std::fmod(x, 360.0) reduction. 360 is exactly representable and the remainder is exact, so the reduced argument is the true one.

I was fairly surprised by the claim that fmod is exact, so I looked into it further. Here's what the IEEE754 spec has to say about the remainder operation.

5.1 Arithmetic
An implementation shall provide the add, subtract, multiply, divide, and remainder operations for any two operands of
the same format, for each supported format; it should also provide the operations for operands of differing formats.
The destination format (regardless of the rounding precision control of 4.3) shall be at least as wide as the wider
operand’s format. All results shall be rounded as specified in Section 4.
When y ≠ 0, the remainder r = x REM y is defined regardless of the rounding mode by the mathematical
relation r = x − y × n, where n is the integer nearest the exact value x/y; whenever |n−x/y| =1/2, then n is even.
Thus, the remainder is always exact. If r = 0, its sign shall be that of x. Precision control (4.3) shall not apply to
the remainder operation.

Surprisingly, even though the operation x / y is inexact, and the result is not even always representable as a floating point number, the operation std::remainder(x, y) is.

The C++ spec is a little less specific about the behavior of std::fmod(), but I think most implementations implement it in substantially the same way. (Minus that fmod()'s range is [0, y) and remainder()'s range is [-y/2, y/2].)

The thing that is more interesting than whether x = std::fmod(x, 360.0); is correct is what this implies about the rest of the function.

If you accept this argument that fmod is exact, it implies that a huge portion of sindg() is unneeded. For example, it computes z = y - 16 * floor(y / 16) using the following formula:

    /* strip high bits of integer part to prevent integer overflow */
    z = std::ldexp(y, -4);
    z = std::floor(z);        /* integer part of y/8 */
    z = y - std::ldexp(z, 4); /* y - 16 * (y/16) */

This is designed to avoid directly doing a division. But if we could do std::fmod(y, 16.0) and have it be exact, this could save 3 lines of code.

Then it does bit-twiddling on j to get an integer 0-7.

j = j & 07

But you can recognize this as j = j % 8. I suspect many simplifications like this are possible.

Conceptual issue

Previously, I said this.

I'm not sure there's an accurate way to do this for large x. If x is large, then x is not very precise in an absolute sense. If x is about 1e14, then x % 360 can be accurate, at most, to 5 decimal digits of precision.

I continue to think this is true: in some sense if someone is asking what sindg(10.0**100) is, they are in big trouble the moment their argument to sindg is rounded. This issue is largely unfixable. If they are trying to get sin of numbers exactly representable as a float, then they're fine, but my feeling is that this is a pretty unusual case.

However, I do agree that defining sindg(x) to be a best-effort attempt to compute the sin is more mathematically elegant than just returning zero.

Test script?

correct for every finite argument up to 1e300, measured 0.64 ulp

How'd you test this? Can you provide your testing script? Is this ULP error figure based on RMSE, max error, something else?

[Kayvan-Zahiri]

Good catch, and thanks for going to the spec on this.

The 0.64 figure doesn't hold up. I didn't keep the script that produced it, and re-measuring doesn't give me that number however I slice it, so I'm withdrawing it.

What I can reproduce is the part that matters. Over 4000 arguments between 1e14 and 1e300, both signs, sindg(x) returns exactly the same bits as sindg(fmod(x, 360)), for all four functions. That is what the test in this PR checks, and it is really the whole argument: if the reduction is exact, then how accurate the function is at huge x is just how accurate it is between 0 and 360. Big arguments never needed an accuracy claim of their own.

Over that range I measure a worst case of 1.43 ulp and an average of 0.11, against mpmath at 60 digits, skipping points right beside a zero where the comparison doesn't mean much. That is the existing code, which this patch doesn't change.

Script below. On the simplifications you spotted, I think you're right, but I'd rather keep this PR to the fix and open a follow-up for them.

Script

"""Accuracy check for xsf's degree trig after the fmod reduction. Needs mpmath.

Build the probe first:

    c++ -std=c++17 -O2 -I include probe.cpp -o probe

where probe.cpp reads doubles as hex from stdin and prints
sindg/cosdg/tandg/cotdg as hex, so bit patterns round-trip.

It answers two separate questions:

A. Is the reduction exact?  sindg(x) must be bit-identical to
   sindg(fmod(x, 360)).  If that holds, accuracy at large x IS accuracy on
   [0, 360), so large arguments need no separate error budget.

B. How accurate is the kernel on [0, 360)?  This is the pre-existing cephes
   polynomial and is unchanged by the patch; it is here as the baseline that
   claim A inherits.
"""
import math, random, struct, subprocess, sys
from mpmath import mp, mpf, sin as msin, cos as mcos, pi as mpi

mp.dps = 60
PROBE = sys.argv[1] if len(sys.argv) > 1 else "./probe"
N = 4000

b = lambda x: struct.unpack("<Q", struct.pack("<d", x))[0]
d = lambda u: struct.unpack("<d", struct.pack("<Q", u))[0]

def run(xs):
    out = subprocess.run([PROBE], input="\n".join(f"{b(x):x}" for x in xs),
                         capture_output=True, text=True).stdout.split()
    return [[d(int(h, 16)) for h in out[i*4:i*4+4]] for i in range(len(xs))]

random.seed(20723)
NAMES = ("sindg", "cosdg", "tandg", "cotdg")

big = []
while len(big) < N:
    x = math.exp(random.uniform(math.log(1e14), math.log(1e300)))
    if math.isfinite(x):
        big.append(random.choice([1, -1]) * x)
gb = run(big)
gr = run([math.fmod(x, 360.0) for x in big])
print(f"A. reduction exactness, {N} args in [1e14, 1e300], both signs")
for i, n in enumerate(NAMES):
    ok = sum(b(p[i]) == b(q[i]) or (math.isnan(p[i]) and math.isnan(q[i]))
             for p, q in zip(gb, gr))
    print(f"   {n:6} bit-identical to reduced argument: {ok}/{N}")

xs = [random.uniform(0, 360) for _ in range(N)]
res = run(xs)
print(f"\nB. kernel accuracy on [0,360) vs mpmath, {mp.dps} digits")
for i, n in enumerate(("sindg", "cosdg")):
    e = []
    for x, r in zip(xs, res):
        ref = (msin if n == "sindg" else mcos)(mpi * mpf(x) / 180)
        if abs(float(ref)) < 1e-8:      # near a zero, a relative ulp says nothing
            continue
        e.append(0.0 if r[i] == float(ref)
                 else abs(float(mpf(r[i]) - ref)) / math.ulp(float(ref)))
    e.sort()
    rms = math.sqrt(sum(v * v for v in e) / len(e))
    print(f"   {n:6} n={len(e)}  max {e[-1]:.2f}  mean {sum(e)/len(e):.2f}  "
          f"rms {rms:.2f}  p99 {e[int(.99*len(e))]:.2f}  >1ulp {sum(v>1 for v in e)}")

probe.cpp

// Prints sindg/cosdg/tandg/cotdg for doubles read from stdin as hex, so the
// exact bit pattern round-trips and nothing is lost to decimal formatting.
#include <cstdio>
#include <cstring>
#include <cstdint>
#include "xsf/cephes/sindg.h"
#include "xsf/cephes/tandg.h"
int main() {
    char line[128];
    while (fgets(line, sizeof line, stdin)) {
        uint64_t bits; if (sscanf(line, "%llx", (unsigned long long*)&bits) != 1) continue;
        double x; memcpy(&x, &bits, 8);
        double s = xsf::cephes::sindg(x), c = xsf::cephes::cosdg(x);
        double t = xsf::cephes::tandg(x), o = xsf::cephes::cotdg(x);
        uint64_t bs, bc, bt, bo;
        memcpy(&bs,&s,8); memcpy(&bc,&c,8); memcpy(&bt,&t,8); memcpy(&bo,&o,8);
        printf("%llx %llx %llx %llx\n", (unsigned long long)bs, (unsigned long long)bc,
               (unsigned long long)bt, (unsigned long long)bo);
    }
    return 0;
}

[nickodell]

What I can reproduce is the part that matters. Over 4000 arguments between 1e14 and 1e300, both signs, sindg(x) returns exactly the same bits as sindg(fmod(x, 360)), for all four functions.

Sorry for taking a few days to get back to you.

My first reaction to this is that 4000 seems a bit low. If you think about the Pentium FDIV bug, that was a 1 in 9 billion bug. I'm not expecting you to test that far - we don't have Intel's resources. But if I 100x the size of this test, it takes 16 seconds to run. This is not especially relevant to your claims about the mean of errors, (these hardly change at larger sizes) but it's quite important to your claims about the worst case. There are relatively few cases where we can just crank a handle, pour in more compute, and get better software out the other side, but randomized testing against a reliable benchmark is one of them.

LLMs will steer you wrong here, because they always seem to suggest parameters for randomized testing that are too small. Two days ago, I wrote a fuzzer for a parser at work. The LLM ran it, and said, "I tried 30K iterations and different seeds, and there are no findings. This fuzzer is probably a dead end." I ran the same piece of software unchanged for 30M iterations, and found a security bug.

So my advice would be to either make the test way way bigger, or to identify a way to test more efficiently, or both. Hypothesis is a useful library here, because it'll pick lots of numbers that trigger corners cases, and you get something like the log uniform distribution you're using ( x = math.exp(random.uniform(math.log(1e14), math.log(1e300))) ) out of the box. The SciPy testing suite has a number of good examples of how to use hypothesis, which you can find here.

Over 4000 arguments between 1e14 and 1e300, both signs, sindg(x) returns exactly the same bits as sindg(fmod(x, 360)), for all four functions. That is what the test in this PR checks, and it is really the whole argument: if the reduction is exact, then how accurate the function is at huge x is just how accurate it is between 0 and 360. Big arguments never needed an accuracy claim of their own.

Should I understand this to mean that you get 0/4000 for all of the tests under A shown below?

$ python3 probe.py
A. reduction exactness, 4000 args in [1e14, 1e300], both signs
   sindg  bit-identical to reduced argument: 58/4000
   cosdg  bit-identical to reduced argument: 0/4000
   tandg  bit-identical to reduced argument: 58/4000
   cotdg  bit-identical to reduced argument: 0/4000

B. kernel accuracy on [0,360) vs mpmath, 60 digits
   sindg  n=4000  max 1.43  mean 0.12  rms 0.29  p99 0.93  >1ulp 16
   cosdg  n=4000  max 1.32  mean 0.11  rms 0.28  p99 0.92  >1ulp 16

Because that is not what I get on my computer. (Linux, gcc, AMD Ryzen 9 5900X.)