/* mpf_sin_cos.c  reduce and use series                     */
/* sin(x) = series_sin(xc)   -pi/4 < xc < pi/4              */
/* x - x^3/3! + x^5/5! - x^7/7! ...  use rfact[i] = 1/i!    */
/* cos(x) = series_cos(xc)   -pi/4 < xc < pi/4              */
/* 1 - x^2/2! + x^4/4! - x^6/6! ...  use rfact[i] = 1/i!    */
/* needs less that n=digits/2  for 0.78^n/n! < 10^(-digits) */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mpf_math.h" /* includes gmp */

static int n;  /* need  0.78^n/n! <  10^(-digits) */
static mpf_t zero;
static mpf_t twopi;
static mpf_t pi;
static mpf_t halfpi;
static mpf_t fourthpi;
static mpf_t temp;
static mpf_t sum;
static mpf_t xp;
static mpf_t xr;
static mpf_t xc2;
static mpf_t *rfact; /* reciprocal factorial */
static init = 1;
static int debug = 0;


void series_sin(mpf_t sum, mpf_t xc) /* |xc|<=0.78 */
{
  int i;

  if(debug) gmp_printf("series_sin xc=%50.45Ff \n", xc);
  mpf_mul(xc2,xc,xc);
  n = ((digits+7)/8)*4+1;
  mpf_set_ui(sum,0);
  for(i=n; i>1; i=i-4)
  {
    mpf_add(sum, sum, rfact[i]);
    mpf_mul(sum, sum, xc2);
    mpf_sub(sum, sum, rfact[i-2]);
    mpf_mul(sum, sum, xc2);
  }
  mpf_add_ui(sum,sum,1); /* rfact[1] */
  mpf_mul(sum,sum,xc); /* sum is returned */
} /* end series_sin */

void series_cos(mpf_t sum, mpf_t xc) /* |xc|<=0.78 */
{
  int i;
  int neg = 1;

  if(debug) gmp_printf("series_cos xc=%50.45Ff \n", xc);
  mpf_mul(xc2, xc, xc);
  mpf_set_si(sum,0);
  n = ((digits+7)/8)*4;
  for(i=n; i>1; i=i-4)
  {
    mpf_add(sum,sum,rfact[i]);
    mpf_mul(sum,sum,xc2);
    mpf_sub(sum,sum,rfact[i-2]);
    mpf_mul(sum,sum,xc2);
  }
  mpf_add_ui(sum,sum,1); /* sum is returned */
} /* end series_cos */

void mpf_sin(mpf_t result, mpf_t x)
{
  int i, j;
  int neg=0; /* negate result */

  if(init==1) /* one time only */
  {
    mpf_set_default_prec(digits*3.32); /* in bits */
    fflush(stdout);
    mpf_init(twopi);
    mpf_twopi(twopi); /* cache constants locally */
    mpf_init(pi);
    mpf_pi(pi);
    mpf_init(halfpi);
    mpf_halfpi(halfpi);
    mpf_init(fourthpi);
    mpf_fourthpi(fourthpi);
    rfact = (mpf_t *)malloc((digits+1)*sizeof(mpf_t));
    for(i=0; i<=digits; i++) mpf_init(rfact[i]);
    mpf_rfact(rfact);

    mpf_init(temp); /* init all mpf variables */
    mpf_init(sum);
    mpf_init(xp);
    mpf_init(xr);
    mpf_init(xc2);
    mpf_init_set_si(zero,0);
    init = 0;
  }

  mpf_set(xr,x);
  /*               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)                 */
  if(mpf_cmp(xr,zero)<0)
  {
    mpf_abs(xr,xr);
    neg = 1;
  }

  if(mpf_cmp(xr, twopi)>=0)
  {
    mpf_div(temp,xr,twopi);
    j = mpf_get_si(temp); /* number of 2Pi's */
    mpf_set_si(temp,j);
    mpf_mul(temp,temp,twopi);
    mpf_sub(xr,xr,temp); /* should use higher precision */
  } 

  if(mpf_cmp(xr, pi)>=0)
  {
    mpf_sub(xr,twopi,xr);
    neg = 1-neg;
  } 

  if(mpf_cmp(xr, halfpi)>=0)
  {
    mpf_sub(xr,pi,xr);
  } 

  if(mpf_cmp(xr,fourthpi)>=0)
  {
    mpf_sub(xr,halfpi,xr);
    series_cos(result,xr);
    if(neg) mpf_neg(result,result);
    return;
  } 
  else
  {
    series_sin(result,xr);
    if(neg) mpf_neg(result,result);
    return;
  }
} /* end mpf_sin */

void mpf_cos(mpf_t result, mpf_t x)
{
  int i, j;
  int neg=0; /* negate result */

  if(init==1) /* one time only */
  {
    mpf_set_default_prec(digits*3.32); /* in bits */
    fflush(stdout);
    mpf_init(twopi);
    mpf_twopi(twopi); /* cache constants locally */
    mpf_init(pi);
    mpf_pi(pi);
    mpf_init(halfpi);
    mpf_halfpi(halfpi);
    mpf_init(fourthpi);
    mpf_fourthpi(fourthpi);
    rfact = (mpf_t *)malloc((digits+1)*sizeof(mpf_t));
    for(i=0; i<=digits; i++) mpf_init(rfact[i]);
    mpf_rfact(rfact);

    mpf_init(temp); /* init all mpf variables */
    mpf_init(sum);
    mpf_init(xp);
    mpf_init(xr);
    mpf_init(xc2);
    mpf_init_set_si(zero,0);
    init = 0;
  }

  gmp_printf("mpf_cos x=%50.45Ff \n", x);
  mpf_set(xr,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)                 */
  if(mpf_cmp(xr,zero)<0)
  {
    mpf_abs(xr,xr);
  }

  if(mpf_cmp(xr, twopi)>=0)
  {
    mpf_div(temp,xr,twopi);
    j = mpf_get_si(temp); /* number of 2Pi's */
    mpf_set_si(temp,j);
    mpf_mul(temp,temp,twopi);
    mpf_sub(xr,xr,temp); /* should use higher precision */
  } 

  if(mpf_cmp(xr, pi)>=0)
  {
    mpf_sub(xr,twopi,xr);
  } 

  if(mpf_cmp(xr, halfpi)>=0)
  {
    mpf_sub(xr,pi,xr);
    neg = 1-neg;
  } 

  if(mpf_cmp(xr,fourthpi)>=0)
  {
    mpf_sub(xr,halfpi,xr);
    series_sin(result,xr);
    if(neg) mpf_neg(result,result);
    return;
  } 
  else
  {
    series_cos(result,xr);
    if(neg) mpf_neg(result,result);
    return;
  }
} /* end mpf_cos */

/* end mpf_sin_cos.c */