Inaccurate trignometric functions for "very" large/small inputs

[SwayamInSync]

print(np.sin(1e100))
print(np.sin(npq.QuadPrecision("1e100")))
print(np.sin(npq.QuadPrecision("1e100", backend="longdouble")))

# output
-0.3806377310050287
0.679816461925869553732614638474571
-0.380637731005028678854529289310449

This is kind of expected, since using taylor expansion for large arguments can cause higher ULP.

A possible workaround in numpy_quaddtype is we can patch the trignometric functions to use any argument reduction strategies like Payne Hanek

@ngoldbaum what is your opinion here? We can also go in the current trivial way and mention in documentation about the working ranges but this won't be similar to NumPy's default behaviour as it seems to apply the reductions before trig functions.

[SwayamInSync]

It seems sleef also deploys a vectorized implementation of Payne Hanek reduction, but it is only available for single and double precision trig functions double Sleef_sin_u10(double a) and float Sleef_sinf_u10(float a)

[juntyr]

I am surprised by this issue since the sleef quad sine is documented as having 1.0 ULP error (independent of the input).

Is -0.38… the mathematically correct value here … ok Wolfram Alpha says it should be -0.37… (https://www.wolframalpha.com/input?i=sin%2810%5E100%29&assumption=%22TrigRD%22+-%3E+%22R%22) so sleef is indeed wrong here

[SwayamInSync]

There's an old discussion on this shibatch/sleef#7 and it seems the bound of 1.0 ULP is for the arguments within range

[WarrenWeckesser]

Drive by comment: with double precision floating point, 1e100 $\ne 10^{100}$. The exact integer value of the double precision number represented as 1e100 is 10000000000000000159028911097599180468360808563945281389781327557747838772170381060813469985856815104. One way to get this is to use fractions.Fraction(1e100):

In [6]: from fractions import Fraction

In [7]: Fraction(1e100)
Out[7]: Fraction(10000000000000000159028911097599180468360808563945281389781327557747838772170381060813469985856815104, 1)

Since we know 1e100 must be an integer, an even easier method to get the exact integer value is simply int(1e100):

In [8]: int(1e100)
Out[9]: 10000000000000000159028911097599180468360808563945281389781327557747838772170381060813469985856815104

I generally use Fraction() when the argument isn't necessarily an integer.

Use that integer as the argument of sin() in Wolfram Alpha. Alpha gives the result:

-0.38063773100502866607187092333210730322126753816281538528564918349932783161025383711588917336347546507272963047064947647620158196441772704920084...

This agrees with mpmath:

In [13]: from mpmath import mp

In [14]: mp.dps = 100

In [15]: mp.sin(1e100)
Out[15]: mpf('-0.3806377310050286660718709233321073032212675381628153852856491834993278316102538371158891733634754650736')

[ngoldbaum]

As for handling this, I think it's totally reasonable to just document that this is a problem in the SLEEF backend.

This also might be another reason to explore the quad precision type available in recent LLVM versions, which does have support in math.h and implements some trig functions.

[juntyr]

I don't think just documenting it is enough. In my library I rely on being able to use numpy_quaddtype as a drop-in with ufuncs that have the expected behaviour. While I currently maintain several wrappers to fix up some cases that @SwayamInSync and I have been fixing, I'd like to not have to. So I think we should implement some solution.

[ngoldbaum]

Fair enough.

Either way, you should document that you're either departing from SLEEF or inherit the issue from SLEEF.

If it's feasible to work around this in the dtype then I'd say go for that, especially if you already have code you can upstream.

[SwayamInSync]

If it's feasible to work around this in the dtype then I'd say go for that, especially if you already have code you can upstream.

One solution in my mind is to have a util file that implements the Payne Hanek (which sleef also uses for double and single precision) (can be non-vectorized for now but atleast should be working) and then we can monkey patch the Sleef's trig function to take the reduced value as input.

so the workflow would be like, inside umath/ops.hpp

Sleef_quad quad_sin(const Sleef_quad input)
{
	Sleef_quad reduced_arg = payne_hanek_reduction(input);
	return Sleef_sinq1_u10(reduced_arg);
}

This to all trig functions