fast sine approximations

fast_sin.cpp

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

// P Found using maxima
//
// y(x) := 4 * x * (%pi-x) / (%pi^2) ; 
// z(x) := (1-p)*y(x) + p * y(x)^2;
// e(x) := z(x) - sin(x);
// solve( diff( integrate( e(x)^2, x, 0, %pi/2 ), p ) = 0, p ),numer;
// 
// [p = .2248391013559941]

double fast_sin1(double x)
{
  const double A = 4.0/(M_PI*M_PI);
  const double P =  0.2248391013559941;
  double y = A* x * ( M_PI - x );

  
  return y* ( (1-P)  + y * P );
}

// P and Q found using maxima
//
// y(x) := 4 * x * (%pi-x) / (%pi^2) ; 
// zz(x) := (1-p-q)*y(x) + p * y(x)^2 + q * y(x)^3
// ee(x) := zz(x) - sin(x)
// solve( [ integrate( diff(ee(x)^2, p ), x, 0, %pi/2 ) = 0, integrate( diff(ee(x)^2,q), x, 0, %pi/2 ) = 0 ] , [p,q] ),numer;
//
// [[p = .1952403377008734, q = .01915214119105392]]

double fast_sin2(double x)
{
  const double A = 4.0/(M_PI*M_PI);
  const double P = 0.1952403377008734;
  const double Q = 0.01915214119105392;

  double y = A* x * ( M_PI - x );

  
  return y*( (1-P-Q) + y*( P + y * Q ) ) ;
}


int main()
{
  printf( "x        sin(x)   fs1(x)   fs2(x)   e1(x)    e2(x)\n" );
  for(unsigned int i=0; i<100; ++i)
  {
    double x = i*(M_PI/(100-1));
    printf( "%04.6f %04.6f %04.6f %04.6f %04.6f %04.6f\n", x, sin(x), fast_sin1(x), fast_sin2(x), sin(x)-fast_sin1(x), sin(x)-fast_sin2(x) );
  }
}