/* reduce_sin.c  sin(-x) = -sin(x)  sin(n*2Pi+x) = sin(x)  */
/*               sin(x>Pi) = sin(x-2Pi) = -sin(2Pi-x)      */
/*               sin(x>Pi/2) = sin(Pi-x)                   */
/*               sin(x>Pi/4) = cos(Pi/2-x)                 */
/*                                                         */
/*               cos(-x) = cos(x)  cos(n*2Pi+x) = cos(x)   */
/*               cos(x>Pi) = cos(2Pi-x)                    */
/*               cos(x>Pi/2) = -cos(Pi-x)                  */
/*               cos(x>Pi/4) = sin(Pi/2-x)                 */

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

static double Pi = 3.1415926535897932384626433832795028841971 ;
static double twoPi, halfPi, fourthPi;

static void sin_reduce(double x); /* sin */
static void cos_reduce(double x); /* cos */


int main(int argc, char * argv[])
{
  double x;

  twoPi = 2.0*Pi;
  halfPi = Pi/2.0;
  fourthPi = Pi/4.0;

  printf("reduce_sin.c fourthPi=%f \n",fourthPi);
  x = twoPi+0.1;
  sin_reduce(x);
  x = twoPi-0.1;
  sin_reduce(x);
  x = Pi+0.1;
  sin_reduce(x);
  x = Pi-0.1;
  sin_reduce(x);
  x = halfPi+0.1;
  sin_reduce(x);
  x = halfPi-0.1;
  sin_reduce(x);
  x = fourthPi+0.1;
  sin_reduce(x);
  x = fourthPi-0.1;
  sin_reduce(x);
  x = twoPi+Pi+halfPi+fourthPi+0.1;
  sin_reduce(x);
  x = 100.0*twoPi+Pi+halfPi+fourthPi+0.1;
  sin_reduce(x);

  x = twoPi+0.1;
  cos_reduce(x);
  x = twoPi-0.1;
  cos_reduce(x);
  x = Pi+0.1;
  cos_reduce(x);
  x = Pi-0.1;
  cos_reduce(x);
  x = halfPi+0.1;
  cos_reduce(x);
  x = halfPi-0.1;
  cos_reduce(x);
  x = fourthPi+0.1;
  cos_reduce(x);
  x = fourthPi-0.1;
  cos_reduce(x);
  x = twoPi+Pi+halfPi+fourthPi+0.1;
  cos_reduce(x);
  x = 100.0*twoPi+Pi+halfPi+fourthPi+0.1;
  cos_reduce(x);

  return 0;
} /* end main */

static void sin_reduce(double x)
{
  int neg=0; /* negate final value */
  double sinx, result;
  int n;

  sinx = sin(x);
  if(x<0.0)
  {
    x = -x;
    neg = 1;
  }
  n = (int)(x/twoPi);
  if(n>0) x = x-(double)n*twoPi; /* should be higher precision */
  if(x>Pi)
  {
    x = twoPi-x;
    neg = 1-neg;
  }
  if(x>halfPi)
  {
    x = Pi-x;
  }
  if(x>fourthPi)
  {
    x = halfPi-x;
    result = cos(x);
  }
  else
  {
    result = sin(x);
  }
  if(neg) result = -result;
  printf("x=%f, sinx=%f, res=%f, err=%g \n",x,sinx,result,sinx-result);
} /* end sin_reduce */

static void cos_reduce(double x)
{
  int neg=0; /* negate final value */
  double cosx, result;
  int n;

  cosx = cos(x);
  if(x<0.0) x = -x;
  n = (int)(x/twoPi);
  if(n>0) x = x-(double)n*twoPi; /* should be higher precision */
  if(x>Pi)
  {
    x = twoPi-x;
  }
  if(x>halfPi)
  {
    x = Pi-x;
    neg = 1;
  }
  if(x>fourthPi)
  {
    x = halfPi-x;
    result = sin(x);
  }
  else
  {
    result = cos(x);
  }
  if(neg) result = -result;
  printf("x=%f, cosx=%f, res=%f, err=%g \n",x,cosx,result,cosx-result);
} /* end cos_reduce */

/* end reduce_sin.c */