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

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.