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

问题:

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

回答1:

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:

• BigNum
• W3B.Sine

回答2:

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();
}

回答3:

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.

回答4:

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.

回答5:

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.

回答6:

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);
        }
    }