sincostest.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#ifndef M_PI
#define M_PI 3.14159265358979323846264338327950288419716939937511
#endif

#ifndef DEBUG
#define DEBUG 0
#endif

#ifndef NBTESTS
#define NBTESTS 1000000000
#endif

#ifndef DOSRAND
#define DOSRAND 0
#endif

double
random_double (void)
  // Returns random real between 0 and 1.
{
  return (double)(rand()) / ((double)RAND_MAX+1) / ((double)RAND_MAX+1) / ((double)RAND_MAX+1) / ((double)RAND_MAX+1)
    + (double)(rand()) / ((double)RAND_MAX+1) / ((double)RAND_MAX+1) / ((double)RAND_MAX+1)
    + (double)(rand()) / ((double)RAND_MAX+1) / ((double)RAND_MAX+1)
    + (double)(rand()) / ((double)RAND_MAX+1);
}

int
main (void)
{
  const int nbtests = NBTESTS;
  double errabsmax = -1., err2absmax = -1., errabssum = 0., err2abssum = 0.;
  int errposcnt = 0, err2poscnt = 0;
#if DOSRAND
  srand(time(NULL));
#endif
  for ( int i=0 ; i<nbtests ; i++ ) {
    double x = random_double() * M_PI;
    double u = cos(x);
    double v = sin(x);
    double errabs = fabs(1 - (u*u + v*v));
    double renorm = sqrt(u*u + v*v);
    double u2 = u/renorm;
    double v2 = v/renorm;
    double err2abs = fabs(1 - (u2*u2 + v2*v2));
    errabssum += errabs;
    if ( errabs > errabsmax )
      errabsmax = errabs;
    if ( errabs > 0. )
      errposcnt++;
    err2abssum += err2abs;
    if ( err2abs > err2absmax )
      err2absmax = err2abs;
    if ( err2abs > 0. )
      err2poscnt++;
#if DEBUG
    fprintf (stderr, "Test %d:\n", i);
    fprintf (stderr, "angle = %e = %a\n", x, x);
    fprintf (stderr, "Without normalizing:\n");
    fprintf (stderr, "cos = %e = %a\n", u, u);
    fprintf (stderr, "sin = %e = %a\n", v, v);
    fprintf (stderr, "error = %e = %a\n", errabs, errabs);
    fprintf (stderr, "With normalizing:\n");
    fprintf (stderr, "cos = %e = %a\n", u2, u2);
    fprintf (stderr, "sin = %e = %a\n", v2, v2);
    fprintf (stderr, "error = %e = %a\n", err2abs, err2abs);
    fprintf (stderr, "\n");
#endif
  }
  printf ("After %d tests:\n", nbtests);
  printf ("Without normalizing: error sum = %e, max error = %e, pos count = %d\n",
	  errabssum, errabsmax, errposcnt);
  printf ("With    normalizing: error sum = %e, max error = %e, pos count = %d\n",
	  err2abssum, err2absmax, err2poscnt);
}

Context: https://twitter.com/gro_tsen/status/1754812805553782951