Accuracy of Math.Sin() and Math.Cos() in C#

Question

I am terribly annoyed by the inaccuracy of the intrinsic trig functions in the CLR. It is well know that

Math.Sin(Math.PI)=0.00000000000000012246063538223773

instead of 0. Something similar happens with Math.Cos(Math.PI/2).

But when I am doing a long series of calculations that on special cases evaluate to

Math.Sin(Math.PI/2+x)-Math.Cos(x)

and the result is zero for x=0.2, but not zero for x=0.1 (try it). Another issue is when the argument is a large number, the inaccuracy gets proportionally large.

So I wonder if anyone has coded some better representation of the trig functions in C# for sharing with the world. Does the CLR call some standard C math library implementing CORDIC or something similar? link:wikipedia CORDIC

How accurate do you believe the representation of pi as a double is? — Jon Skeet, Jul 14 '10 at 19:31
If you want symbolic maths, do symbolic maths. If you use floating point types, you get finite precision. — AakashM, Jul 14 '10 at 19:45
-1 for not "doing your homework", and also for thinking that `System.Math` is part of C# (hint: it's part of the .NET Framework). — John Saunders, Jul 14 '10 at 19:49
Taking a quick look at the CORDIC link, why do you think it would be an improvement? It looks like a different method of calculation, but it would have similar precision issues. — David Thornley, Jul 14 '10 at 19:55
If you are measuring the position of the Earth in miles, Math.Sin(Math.PI) as you show it above will give a position that is off by seven tenths of an inch. What application are you developing that requires such precision?! — Jeffrey L Whitledge, Jul 14 '10 at 20:06
My problems arise in multibody dynamics where quantities are 6 component vectors (and not scalar), and transformations have 6`x`6=36 elements. Often times you expect results to evaluate to zero. I guess this is a hybrid application where it uses some expression trees for the joint variable and in the end produces expressions (functions) to be integrated numerically. When things do not simplify properly it bombs. — John Alexiou, Jul 15 '10 at 13:09
@jalexiou: If you're doing dynamics problems that involve differential equations and numeric integration all in floats then the inaccuracy of a few quadrillionths in the value of pi is likely to be the least of your problems. There are many well-known error-terms-growing-too-large problems with numeric approximation algorithms for differential equations and quadrature. If you want to continue to pursue a numeric approximation solution, pick up a good undergraduate text on numerical methods. If you want to do math in exact arithmetic, get mathematica or maple or some such tool. — Eric Lippert, Jul 15 '10 at 16:14

Darin Dimitrov · Answer 1 · 2010-07-14T19:54:54.923

This has nothing to do with accuracy of trigonometric functions but more with the CLS type system. According to the documentation a double has 15-16 digits precision (which is exactly what you get) so you can't be more precise with this type. So if you want more precision you will need to create a new type that is capable of storing it.

Also notice that you should never be writing a code like this:

double d = CalcFromSomewhere();
if (d == 0)
{
    DoSomething();
}

You should do instead:

double d = CalcFromSomewhere();
double epsilon = 1e-5; // define the precision you are working with
if (Math.Abs(d) < epsilon)
{
    DoSomething();
}

score 10 · Answer 2 · answered Jul 14 '10 at 22:48

I hear you. I am terribly annoyed by the inaccuracy of division. The other day I did:

Console.WriteLine(1.0 / 3.0);

and I got 0.333333333333333, instead of the correct answer which is 0.333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333...

Perhaps now you see what the problem is. Math.Pi is not equal to pi any more than 1.0 / 3.0 is equal to one third. Both of them differ from the true value by a few hundred quadrillionths, and therefore any calculations you perform with Math.Pi or 1.0/3.0 are also going to be off by a few hundred quadrillionths, including taking the sine.

If you don't like that approximate arithmetic is approximate then don't use approximate arithmetic. Use exact arithmetic. I used to use Waterloo Maple when I needed exact arithmetic; perhaps you should buy a copy of that.

Out of curiosity, if one wants to calculate exact arithmetic in C#, would it be possible? Just don't know what one can use to do that? Not even decimals would cut it, right? — Joan Venge, Jul 14 '10 at 23:47
@Joan: Yes, we can do it easily for integers, and even for rational numbers (just store numerator/denominator as large integers), and throw rules into our library for whatever specific real numbers we want: pi, e, square-roots, etc. However, a library for arbitrary-precision arithmetic on *any* imaginable real number is impossible; even assuming you had some way of storing them *(say as a formula which will give us an arbitrary digit of the number)*, it is computationally-impossible to even compare two arbitrary real numbers for equality!! — BlueRaja - Danny Pflughoeft, Jul 15 '10 at 00:53
@BlueRaja: indeed. Typically what you do in that case is manipulate the arithmetic quantities symbolically; you have a symbol for pi and a symbol for e, the same way that you have a symbol for 1, 2, 3, and then you encode all the arithmetic and trigonometric identities, like that sine of pi is zero, and so on. Symbolic math is hard. — Eric Lippert, Jul 15 '10 at 14:11

score 6 · Answer 3 · answered Jul 14 '10 at 19:53

This is a result of floating-point precision. You get a certain number of significant digits possible, and anything that can't be represented exactly is approximated. For example, pi is not a rational number, and so it's impossible to get an exact representation. Since you can't get an exact value of pi, you aren't going to get exact sines and cosines of numbers including pi (nor will you get exact values of sines and cosines most of the time).

The best intermediate explanation is "What Every Computer Scientist Should Know About Floating-Point Arithmetic". If you don't want to go into that, just remember that floating point numbers are usually approximations, and that floating-point calculations are like moving piles of sand on the ground: with everything you do with them, you lose a little sand and pick up a little dirt.

If you want exact representation, you'll need to find yourself a symbolic algebra system.

I understand if PI is not defined exactly due to IEEE-754 arithmetic, and I do wish I could drive a symbolic algebra system with C#, but I can't right now. — John Alexiou, Jul 14 '10 at 20:04

score 2 · Answer 4 · answered Jun 27 '12 at 04:07

I reject the idea the the errors are due to round-off. What can be done is define sin(x) as follows, using a Taylor's expansion with 6 terms:

    const double π=Math.PI;
    const double π2=Math.PI/2;
    const double π4=Math.PI/4;

    public static double Sin(double x)
    {

        if (x==0) { return 0; }
        if (x<0) { return -Sin(-x); }
        if (x>π) { return -Sin(x-π); }
        if (x>π4) { return Cos(π2-x); }

        double x2=x*x;

        return x*(x2/6*(x2/20*(x2/42*(x2/72*(x2/110*(x2/156-1)+1)-1)+1)-1)+1);
    }

    public static double Cos(double x)
    {
        if (x==0) { return 1; }
        if (x<0) { return Cos(-x); }
        if (x>π) { return -Cos(x-π); }
        if (x>π4) { return Sin(π2-x); }

        double x2=x*x;

        return x2/2*(x2/12*(x2/30*(x2/56*(x2/90*(x2/132-1)+1)-1)+1)-1)+1;
    }

Typical error is 1e-16 and worst case is 1e-11. It is worse than the CLR, but it is controllable by adding more terms. The good news is that for the special cases in the OP and for Sin(45°) the answer is exact.

Related post on which angles have exact trig. values http://math.stackexchange.com/q/176889/3301 — John Alexiou, Dec 20 '13 at 23:17
If the only problem are the special cases of the OP, one might as well write them out as if-else statements instead of going with such a solution. These methods may work when called directly with e.g. zero, but not an approximate value. If e.g. you call Sin with the approximate result of a previous operation, it won't work (e.g. because x is very small, but not exactly zero – even if it "should" be). The real problem here is the use of double and the expectation to get exact results, which no algorithm or formula in the world can fix. — enzi, Jan 25 '14 at 18:50

score 2 · Accepted Answer · answered Jul 14 '10 at 23:21

2

You need to use an arbitrary-precision decimal library. (.Net 4.0 has an arbitrary integer class, but not decimal).

A few popular ones are available:

answered Jul 14 '10 at 23:21

BlueRaja - Danny Pflughoeft

84,206
33
197
283

score 1 · Answer 6 · answered May 18 '18 at 14:14

Our current implementation of sine and cosine is

    public static double Sin(double d) {
        d = d % (2 * Math.PI); // Math.Sin calculates wrong results for values larger than 1e6
        if (d == 0 || d == Math.PI || d == -Math.PI) {
            return 0.0;
        }
        else {
            return Math.Sin(d);
        }
    }

    public static double Cos(double d) {
        d = d % (2 * Math.PI); // Math.Cos calculates wrong results for values larger than 1e6
        double multipleOfPi = d / Math.PI; // avoid calling the expensive modulo function twice
        if (multipleOfPi == 0.5 || multipleOfPi == -0.5 || multipleOfPi == 1.5 || multipleOfPi == -1.5) { 
            return 0.0;
        }
        else {
            return Math.Cos(d);
        }
    }

Try `Sin(1e16*Math.PI)` which `C#` yields `-0.375175399369948` versus the correct 0. In fact, this is the same result you get in MATLAB, Excel and any other numeric library that uses the `x87` co-processor for trig. Since the original post, I have realized the problem isn't `C#` but the underlying IEEE math platform that all modern PC use. — John Alexiou, May 18 '18 at 14:39

Accuracy of Math.Sin() and Math.Cos() in C#

6 Answers6

Linked

Related