double.Cos and double.SinCos can produce very different results

[stephentoub]

Repro:

double d = 1e18;
double cos1 = double.Cos(d);
(_, double cos2) = double.SinCos(d);
Console.WriteLine(cos1);
Console.WriteLine(cos2);

outputs:

0.11837199021871073
-0.9929693207404051

on .NET 8

[ghost]

Tagging subscribers to this area: @dotnet/area-system-numerics
See info in area-owners.md if you want to be subscribed.

Issue Details

---

Repro:

double d = 1e18;
double cos1 = double.Cos(d);
(_, double cos2) = double.SinCos(d);
Console.WriteLine(cos1);
Console.WriteLine(cos2);

outputs:

0.11837199021871073
-0.9929693207404051

on .NET 8

Author:	stephentoub
Assignees:	-
Labels:	area-System.Numerics, untriaged
Milestone:	-

[tannergooding]

SinCos is returning sin for both values. Looks to repro for debug and release

[tannergooding]

Looks to be Windows specific as well, but repros back to .NET 6 when SinCos was introduced

[tannergooding]

Looks to be a bug in MSVC for float_control(precise, off) (that is /fp:fast) which is used in the ecall wrappers. It repros in a standalone C++ application.

We're generating:

mov         qword ptr [rsp+8],rbx  
push        rdi  
sub         rsp,20h  
movaps      xmm1,xmm0  
mov         rdi,r8  
xorps       xmm0,xmm0  
mov         rbx,rdx  
movsd       xmm0,xmm1  
call        __libm_sse2_sincos_ (07FF705521F90h)  
movsd       mmword ptr [rbx],xmm0  
mov         rbx,qword ptr [rsp+30h]  
shufpd      xmm0,xmm0,1  
movsd       mmword ptr [rdi],xmm0  
add         rsp,20h  
pop         rdi  
ret  

So it looks like the expectation is that __libm_sse2_sincos takes in value in xmm0 (these matches the expecting calling convention) and then returns the result in xmm0[0] and xmm0[1]. However, sin is in both elements instead.

The API is internally calling __vdecl_sin2 which dispatches to __sse4_sin2@@16. The bug looks to be here as it is not setting the results correctly, I think it's because sin2 is meant to do 2xSin operations and it needs to be calling a different API for SinCos, but I'm working to confirm

This also explains why it reproduces on .NET 6, as we are dynamically invoking the C runtime and the underlying implementation has changed in the past couple years (it appears we are using SVML in some places now).

[tannergooding]

There is indeed a __vdecl_sincos2 that should have been the target instead, looks to be an explicit bug in the __libm_sse2_sincos_ helper.

This might date back to ~2019 given the announcement for using MSVC here: https://devblogs.microsoft.com/cppblog/msvc-backend-updates-in-visual-studio-2019-preview-2/

It's unclear when the optimization started occurring to use the internal helper, but it looks like it wasn't always being done, but since it's a dynamic call into the CRT, it happened somewhat implicitly as part of a tooling upgrade.

[tannergooding]

Just copy/pasting my analysis from the PR

It looks like __libm_sse2_sincos is trying to rely on the fact that sin(x + PI / 2) == cos(x) and so adjusts the second input by that much. This is correct for infinitely precise arithmetic, but for floating-point means that cos is less accurate than it should be for regular values and will not work at all for large values due to the increasing gap between which numbers are representable (it breaks down entirely around 1e16, but starts showing issues far earlier).

Given that its just calling into SVML and SVML provides a dedicated __vdecl_sincos2, that should really be getting used instead.

We missed this because we were testing normal input ranges and not testing anything quite so large, effectively relying on the fact that the underlying libm implementation is expected to be largely correct, and that we simply want to minimally validate the bindings are accurate and that it matches any IEEE 754 required behavior (such as sin(+0) == +0 and sin(-0) == -0, with no exception)

[tannergooding]

The MSVC bug has been filed here: https://developercommunity.visualstudio.com/t/MSVCs-sincos-implementation-is-incorrec/10582378