Sine and Cosine Accuracy on AMD64 and Pentium 4

Let's consider the accuracy of sine and cosine. I've run tests as
follows, using a program provided at the end of this message.


On the Opteron, using GCC 4.0.0 release, the command lines produce these
outputs:


-lm -O3 -march=k8 -funsafe-math-optimizations -mfpmath=387


    generates:
            fsincos


    cumulative accuracy: 60.830074998557684 (binary)
                                                  18.311677213055471 (decimal)


-lm -O3 -march=k8 -mfpmath=387


    generates:
            call sin
            call cos


    cumulative accuracy: 49.415037499278846 (binary)
                                                  14.875408524143376 (decimal)


-lm -O3 -march=k8 -funsafe-math-optimizations


    generates:
            call sin
            call cos


    cumulative accuracy: 47.476438043942984 (binary)
                                                  14.291831938509427 (decimal)


-lm -O3 -march=k8


    generates:
            call sin
            call cos


    cumulative accuracy: 47.476438043942984 (binary)
                                                  14.291831938509427 (decimal)


The default for Opteron is -mfpmath=sse; as has been discussed in other
threads, this may not be a good choice. I also note that using
-funsafe-math-optimizations (and thus the combined fsincos instruction)
*increases* accuracy.


On the Pentium4, using the same version of GCC, I get:


-lm -O3 -march=pentium4 -funsafe-math-optimizations


    cumulative accuracy: 63.000000000000000 (binary)
                                                  18.964889726830815 (decimal)


-lm -O3 -march=pentium4


    cumulative accuracy: 49.299560281858909 (binary)
                                                  14.840646417884166 (decimal)


-lm -O3 -march=pentium4 -funsafe-math-optimizations -mfpmath=sse


    cumulative accuracy: 47.476438043942984 (binary)
                                                  14.291831938509427 (decimal)


The program used is below. I'm very open to suggestions about this
program, which is a subset of a larger accuracy benchmark I'm writing
(Subtilis).


#include <fenv.h>
#pragma STDC FENV_ACCESS ON
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdbool.h>
#include <string.h>


static bool verbose = false;
#define PI 3.14159265358979323846


// Test floating point accuracy
inline double binary_accuracy(double x)
{
        return -(log(fabs(x)) / log(2.0));
}


inline double decimal_accuracy(double x)
{
        return -(log(fabs(x)) / log(10.0));
}


// accuracy of trigonometric functions
void trigtest()
{
        static const double range = PI; // * 2.0;
        static const double incr = PI / 100.0;


        if (verbose)
              printf(" x diff accuracy\n");


        double final = 1.0;
        double x;


        for (x = -range; x <= range; x += incr)
        {
                double s1 = sin(x);
                double c1 = cos(x);
                double one = s1 * s1 + c1 * c1;
                double diff = one - 1.0;
                final *= one;


                double accuracy1 = binary_accuracy(diff);


                if (verbose)
                        printf("%20.15f %14g %20.15f\n",x,diff,accuracy1);
        }


        final -= 1.0;


        printf("\ncumulative accuracy: %20.15f (binary)\n",
                      binary_accuracy(final));


        printf(" %20.15f (decimal)\n",
                      decimal_accuracy(final));
}


// Entry point
int main(int argc, char ** argv)
{
        int i;


        // do we have verbose output?
        if (argc > 1)
        {
                for (i = 1; i < argc; ++i)
                {
                        if (!strcmp(argv[i],"-v"))
                        {
                                verbose = true;
                                break;
                        }
                }
        }




        // run tests
        trigtest();


        // done
        return 0;
}


...Scott
Coyote Gulch Productions
http://www.coyotegulch.com

Return to the comp.compilers page.
Search the comp.compilers archives again.