Sine and Cosine Accuracy

[Scott Robert Ladd]

Let's consider the accuracy of sice 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