Different results between c++ and c# sin function with large values

Question

I came across such strange behavior of the Math.Sin function in C#, when I use large numbers; for example:

C#: .Net 4.7.2: Math.Sin(6.2831853071795856E+45) = 6.2831853071795856E+45

C++: sin(6.2831853071795856E+45) = -0.089650623841643268

Any ideas how to get the same results as C++?

C# sample:

double value = 6.2831853071795856E+45;    
Console.WriteLine(Math.Sin(value));

C++ sample:

double x = 6.2831853071795856E+45;
double result;
    
result = sin(x);
cout << "sin(x) = " << result << endl;

Answer 1

f(x) = sin(x) oscillates between plus and minus 1 with period 2 * pi, i.e. sin(x) = sin(N * 2 * pi + x).

Forget about programming language and forget about sin() for a moment and think about what the value in scientific notation 1.23E+45 means. It could represent any value from 1225000...(41 zeros) to 1234999...(42 nines), i.e. 1E43 possible integers (or an infinite number of values if fractions are allowed).

I'm playing fast and loose with strict definitions here by speaking in decimal terms and equating significant digits with precision. Find links here to learn more: https://en.wikipedia.org/wiki/Precision. I also hope I've got my exponent arithmetic right, but the implication should be clear regardless.

Back to sin(), C# and C++: since the pure mathematical function output repeats every 2 * pi (~6.28) and double precision floats (IEEE 754 64 bit doubles https://en.wikipedia.org/wiki/IEEE_754) in both C# and C++ have only 15-17 significant decimal digits, no output from the sin() function can possibly be meaningful because the input values are "missing" (imprecise by) at least 30 significant/relevant decimal digits.

Answer 2

Both answers are very wrong—but the question you're asking is likely to get you into trouble.

The true value of sin(6.2831853071795856 × 10⁴⁵) is approximately 0.09683996046341126; this approximation differs from the true value by less than 10⁻¹⁶ parts of the true value. (I computed this approximation using Sollya with 165 bits of intermediate precision.)

However, you won't get this answer by asking a C# function with the signature public static double Sin (double a), or a C++ function with the signature double sin(double). Why? 6.2831853071795856 × 10⁴⁵ is not an IEEE 754 binary64, or ‘double’, floating-point number, so at best you will learn what the sin of a nearby floating-point number is. The nearest floating-point number, and what you will usually get by typing 6.2831853071795856E+45 into a program, is 6283185307179585571582855233194226059181031424, which differs from 6.2831853071795856 × 10⁴⁵ by 28417144766805773940818968576 ≈ 2.84 × 10²⁸.

The IEEE 754 binary64 floating-point number 6283185307179585571582855233194226059181031424 is a good approximation to 6.2831853071795856 × 10⁴⁵ in relative error (it differs by less than 10⁻¹⁷ parts of the true value), but the absolute error ~2.84 × 10²⁸ is far beyond the period 2𝜋 of sin (and nowhere near an integral multiple of 𝜋). So the answer you get by asking a double function will bear no resemblance to the question your source code seems to ask: where you write sin(6.2831853071795856E+45), instead of sin(6283185307179585600000000000000000000000000000) at best you will get sin(6283185307179585571582855233194226059181031424), which is about 0.8248163906169679 (again, plus or minus 10⁻¹⁶ parts of the true value).

It's not the fault of floating-point that this goes wrong. Floating-point arithmetic can do a good job of keeping relative error small—and a good math library can easily use binary64 floating-point to compute a good answer to the question you asked. The error in your input to sin could just as well have come from a small measurement error, if your ruler doesn't have many more than 10⁴⁵ gradations. The error could have come from some kind of approximation error, say by using a truncated series to evaluate whatever function gave you sin's input, no matter what kind of arithmetic you used to compute that input. The problem is that you asked to evaluate the periodic function (sin) at a point where a small relative error corresponds to an absolute error far beyond the period of the function.

---

So if you find yourself trying to answer the question of what the sin of 6.2831853071795856 × 10⁴⁵ is, you're probably doing something wrong—and naively using double floating-point math library routines is not going to help you to answer your question. But compounding this problem, your C# and C++ implementations both fail to return anything near the true value of sin(6283185307179585571582855233194226059181031424):

• The C# documentation for Math.Sin advertises that there may be machine-dependent restrictions on the domain.

  Most likely, you are on an Intel CPU, and your C# implementation simply executes the Intel x87 fsin instruction, which according to the Intel manual is restricted to inputs in the domain [−2⁶³, 2⁶³], whereas yours is beyond 2¹⁵². Inputs outside this domain are, as you observed, returned verbatim, even though they are totally nonsensical values for the sine function.

  A quick and dirty way to get the input into a range that will definitely work is to write:

  Math.Sin(Math.IEEERemainder(6.2831853071795856E+45, 2*Math.PI))
  

  That way, you aren't misusing the Math.Sin library routine, so the answer will at least lie in [−1,1] as a sine should. And you can arrange to get the same or nearby result in C/C++ with sin(fmod(6.2831853071795856E+45, 2*M_PI)). But you'll probably get a result near 0.35680453559729486, which is also wrong—see below on argument reduction.

• The C++ implementation of sin that you are using, however, is simply broken; there is no such restriction on the domain in the C++ standard, and with widely available high-quality software to compute argument reduction modulo 𝜋 and to compute sin on the reduced domain, there's no excuse for screwing this up (even if this is not a good question to ask!).

  I don't know what the bug is just from eyeballing the output, but most likely it is in the argument reduction step: since sin(𝑥 + 2𝜋) = sin(𝑥) and sin(−𝑥) = −sin(𝑥), if you want to compute sin(𝑥) for an arbitrary real number it suffices to compute sin(𝑦) where 𝑦 = 𝑥 + 2𝜋𝑘 lies in [−𝜋,𝜋], for some integer 𝑘. Argument reduction is the task of computing 𝑦 given 𝑥.

  Typical x87-based implementations of sin use the fldpi instruction to load an approximation to 𝜋 in binary80 format with 64 bits of precision, and then use fprem1 to reduce modulo that approximation to 𝜋. This approximation is not very good: internally, the Intel architecture approximates 𝜋 by 0x0.c90fdaa22168c234cp+2 = 3.1415926535897932384585988507819109827323700301349163055419921875 with 66 bits of precision, and fldpi then rounds it to the binary80 floating-point number 0x0.c90fdaa22168c235p+2 = 3.14159265358979323851280895940618620443274267017841339111328125 with only 64 bits of precision.

  In contrast, typical math libraries, such as the venerable fdlibm, usually use an approximation with well over 100 bits of precision for argument reduction modulo 𝜋, which is why fdlibm derivatives are able to compute sin(6283185307179585571582855233194226059181031424) quite accurately.

  However, the obvious x87 computation with fldpi/fprem1/fsin gives about −0.8053589558881794, and the same with the x87 unit set to binary64 arithmetic (53-bit precision) instead of binary80 arithmetic (64-bit precision), or just using a binary64 approximation to 𝜋 in the first place, gives about 0.35680453559729486. So evidently your C++ math library is doing something else to give the wrong answer to a bad question!

---

Judging by the number you fed in, I would guess you may have been trying to see what happens when you try to evaluate the sin of a large multiple of 2𝜋: 6.2831853071795856 × 10⁴⁵ has a small relative error from 2𝜋 × 10⁴⁵. Of course, such a sin is always zero, but perhaps you more generally want to compute the function sin(2𝜋𝑡) where 𝑡 is a floating-point number.

The standard for floating-point arithmetic, IEEE 754-2019, recommends (but does not mandate) operations sinPi, cosPi, tanPi, etc., with sinPi(𝑡) = sin(𝜋⋅𝑡). If your math library supported them, you could use sinPi(2*t) to get an approximation to sin(2𝜋𝑡). In this case, since 𝑡 is an integer (and indeed for any binary64 floating-point numbers at least 2⁵² in magnitude, which are all integers), you would get exactly 0.

Unfortunately, .NET does not yet have these functions (as of 2021-02-04), nor does the C or C++ standard math library include them, although you can easily find example code for sinPi and cosPi floating around.

Of course, your question still has the problem that evaluating even the sinPi function at very inputsis asking for trouble, because even a small relative error (say 10⁻¹⁷) still means an absolute error far beyond the function's period, 2. But the problem won't be magnified by bad argument reduction modulo a transcendental number: computing the remainder after dividing a floating-point number by 2 is easy to do exactly.

Answer 3

The sin of large numbers like 6.2831853071795856E+45; makes no sense. Remember that Pi is approximately 3.14 and that your floating numbers are probably IEEE 754 (see the http://floating-point-gui.de/ for more)

The computed number is probably entirely wrong.

Consider using static analysis tools like Fluctuat on your C or C++ code. Or dynamic tools like CADNA

As a rule of thumb, trigonometric functions like sin, cos, tan should not be used with large numbers (e.g. more than a few thousands times PI ℼ = 3.14 in absolute value).

Otherwise, use a bignum library like GMPlib. It will slow down calculations by a factor of thousands, but it could give you meaningful results for sin (pi * 1.e+45)

And please read something about Taylor series (they are related to trigonometric functions)

Try also to use GNUplot to "visualize" the sin function (for reasonable numbers).

Mathematically the sin of every finite real number is between -1 and 1. So Math.Sin(6.2831853071795856E+45) = 6.2831853071795856E+45 is wrong.

Read also conference papers like this one and that one.

Answer 4

As a sidenote there was a change in how sin/cos/tan are calculated .NET. It happened on 2 Jun 2016, with this commit on .NET Core. As the author wrote in floatdouble.cpp:

Sin, Cos, and Tan on AMD64 Windows were previously implemented in vm\amd64\JitHelpers_Fast.asm by calling x87 floating point code (fsin, fcos, fptan) because the CRT helpers were too slow. This is no longer the case and the CRT call is used on all platforms.

Note that the CRT is the C Language Runtime. This explains why newer versions of .NET Core give the same result as C++ when used on the same platform. They are using the same CRT other C++ programs are using.

The last version of the "old" code for those interested: 1 and 2.

It is unclear if/when .NET Framework inherited these changes. From some tests it seems that

.NET Core >= 2.0 (haven't tested previous versions): Math.Sin(6.2831853071795856E+45) == 0.824816390616968

.NET Framework 4.8 (32 and 64 bits): Math.Sin(6.2831853071795856E+45) == 6.28318530717959E+45

so it didn't inherit them.

Funny addendum

Made some checks, and Math.Sin(6.2831853071795856E+45) == 6.28318530717959E+45 is the "official" answer of the fsin opcode of x87 assembly, so it is the "official" answer of Intel (about 1980) about how much is the sin of 6.2831853071795856E+45, so if you use an Intel, you must trust that Math.Sin(6.2831853071795856E+45) == 6.28318530717959E+45, otherwise you are a traitor! 🤣🤣🤣

If you want to doublecheck:

public static class TrigAsm
{
    [DllImport("kernel32.dll", ExactSpelling = true, SetLastError = true)]
    private static extern IntPtr VirtualAlloc(IntPtr lpAddress, IntPtr dwSize, uint flAllocationType, uint flProtect);

    [DllImport("kernel32.dll", ExactSpelling = true, SetLastError = true)]
    [return: MarshalAs(UnmanagedType.Bool)]
    private static extern bool VirtualProtect(IntPtr lpAddress, IntPtr dwSize, uint flAllocationType, out uint lpflOldProtect);

    [DllImport("kernel32.dll", ExactSpelling = true, SetLastError = true)]
    [return: MarshalAs(UnmanagedType.Bool)]
    private static extern bool VirtualFree(IntPtr lpAddress, IntPtr dwSize, uint dwFreeType);

    private const uint PAGE_READWRITE = 0x04;
    private const uint PAGE_EXECUTE = 0x10;
    private const uint MEM_COMMIT = 0x1000;
    private const uint MEM_RELEASE = 0x8000;

    [SuppressUnmanagedCodeSecurity]
    [UnmanagedFunctionPointer(CallingConvention.StdCall)]
    public delegate double Double2Double(double d);

    public static readonly Double2Double Sin;

    static TrigAsm()
    {
        // Opcoes generated with https://defuse.ca/online-x86-assembler.htm
        byte[] body = Environment.Is64BitProcess ?
            new byte[] 
            { 
                0xF2, 0x0F, 0x11, 0x44, 0x24, 0x08, // movsd  QWORD PTR [rsp+0x8],xmm0
                0xDD, 0x44, 0x24, 0x08,             // fld    QWORD PTR [rsp+0x8] 
                0xD9, 0xFE,                         // fsin
                0xDD, 0x5C, 0x24, 0x08,             // fstp   QWORD PTR [rsp+0x8]
                0xF2, 0x0F, 0x10, 0x44, 0x24, 0x08, // movsd  xmm0,QWORD PTR [rsp+0x8]
                0xC3,                               // ret
            } :
            new byte[] 
            { 
                0xDD, 0x44, 0x24, 0x04, // fld    QWORD PTR [esp+0x4]
                0xD9, 0xFE,             // fsin
                0xC2, 0x08, 0x00,       // ret    0x8
            };

        IntPtr buf = IntPtr.Zero;

        try
        {
            // We VirtualAlloc body.Length bytes, with R/W access
            // Note that from what I've read, MEM_RESERVE is useless
            // if the first parameter is IntPtr.Zero
            buf = VirtualAlloc(IntPtr.Zero, (IntPtr)body.Length, MEM_COMMIT, PAGE_READWRITE);

            if (buf == IntPtr.Zero)
            {
                throw new Win32Exception();
            }

            // Copy our instructions in the buf
            Marshal.Copy(body, 0, buf, body.Length);

            // Change the access of the allocated memory from R/W to Execute
            uint oldProtection;
            bool result = VirtualProtect(buf, (IntPtr)body.Length, PAGE_EXECUTE, out oldProtection);

            if (!result)
            {
                throw new Win32Exception();
            }

            // Create a delegate to the "function"
            Sin = (Double2Double)Marshal.GetDelegateForFunctionPointer(buf, typeof(Double2Double));

            buf = IntPtr.Zero;
        }
        finally
        {
            // There was an error!
            if (buf != IntPtr.Zero)
            {
                // Free the allocated memory
                bool result = VirtualFree(buf, IntPtr.Zero, MEM_RELEASE);

                if (!result)
                {
                    throw new Win32Exception();
                }
            }
        }
    }
}

In .NET you can't have inline assembly (and the inline assembler of Visual Studio is 32 bits only). I solved it putting the assembly opcodes in a block of memory and marking the memory as executable.

Post scriptum: note that no-one uses anymore the x87 instruction set, because it is slow.