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

[SwayamInSync]

Update

On more testing I realized the SLEEF trignometric operations are working correct, the error we noticed is the cause of approximating input 1e100 to floating point representation. I'll try to address the confusion happened in previous messages along with my misinterpretation.

Case-1: Input - 1e100

import numpy as np
import numpy_quaddtype as npq
from mpmath import mp

mp.dps = 100

inp = 1e100 # same as doing float(1e100)

print(np.sin(np.float64(inp)))
print(np.sin(npq.QuadPrecision(inp)))
print(mp.sin(inp))

# output
-0.3806377310050287
-0.3806377310050286660718709233321073
-0.3806377310050286660718709233321073032212675381628153852856491834993278316102538371158891733634754651

In Python, 1e100 is a float64 literal, which rounds 10^100 to the nearest representable value: exactly 10000000000000000159028911097599180468360808563945281389781327557747838772170381060813469985856815104 . NumPy's np.sin(np.float64(1e100)) and SLEEF-based np.sin(npq.QuadPrecision(1e100)) (which promotes the same float64 value to quad) both compute sin of this rounded input, yielding approximately -0.3806377310050286660718709233321073.

One thing to note here is mp.dps it stands for decimal places of precision Internally, mpmath stores numbers in binary, but dps tells it how many decimal digits you want to be accurate to.

The binary precision (mp.prec) is automatically derived from dps using

$$ \text{prec} \approx \log_2(10) * dps \approx 3.321928 * dps $$

So technically we are till now seeing mpmath results in 332 binary precision which is much more than for a 128-bit floating point (113). We can directly set mp.prec = 113 but I didn't yet just to make consistency in the discussion with earlier messages.

Tl;dr what input is getting fed to function is the float64 conversion of 1e100 and for SLEEF it is further constructing float128 from float64

Case-2: Input - "1e100"

This was the case, on which this issue was opened

import numpy as np
import numpy_quaddtype as npq
from mpmath import mp

mp.dps = 100
inp = "1e100"

print(np.sin(np.float64(inp)))
print(np.sin(npq.QuadPrecision(inp)))
print(mp.sin(inp))

# output
-0.3806377310050287
0.6798164619258695537326146384745714
-0.3723761236612766882620866955531642957196678835674347023644153882967192240437564411887366004162030232

Here we can see all 3 versions are giving different results and this is all because internally how they convert from string to float.

• np.float64 gave the same value as earlier was expected because remember Python's float and np.float64 are both float64 representation, so even input is string or not gets represented to same value
• QuadPrecision here is constructing float128 bit value directly from the string (here is no intermediate float64 => float128) which is indeed more precise and producing the corresponding result in float128 binary precision value
• mp.sin is also constructing the float representative from string but This is not 128bit float, as the mp.dps value is much higher resulting in more precise and numerically correct value (which also matches with wolframalpha)

But if we make mpmath use the 113 binary precision or 34 dps then

# mp.dps = 100
mp.prec = 113
inp = "1e100"

print(np.sin(np.float64(inp)))
print(np.sin(npq.QuadPrecision(inp)))
print(mp.sin(inp))

# Output
-0.3806377310050287
0.6798164619258695537326146384745714
0.679816461925869553732614638474571

It exactly matches with the QuadPrecision value.

Final Tl;dr the implementation seems correct, it is more just the nuances of output's floating point precision and the type of input

[SwayamInSync]

Let me know if this make sense we can close this issue

[juntyr]

This makes sense - the issue can be closed once we have tests for all trig functions

[juntyr]

Thanks for investing @SwayamInSync