/* specfunc/sincos_pi.c

 * 

 * Copyright (C) 2017 Gerard Jungman, Konrad Griessinger (konradg@gmx.net)

 * 

 * This program is free software; you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation; either version 3 of the License, or (at

 * your option) any later version.

 * 

 * This program is distributed in the hope that it will be useful, but

 * WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

 * General Public License for more details.

 * 

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

 */

/* routines for computing sin(pi*x) and cos(pi*x), respectively, with argument reduction */

#include <math.h>

#include <gsl/gsl_math.h>

#include <gsl/gsl_sf_trig.h>

#include <gsl/gsl_errno.h>

#include <gsl/gsl_sf_sincos_pi.h>

/* Any double precision number bigger than this is automatically an even integer. */

#define TWOBIG (2.0 / GSL_DBL_EPSILON)

/* routine computing sin(pi*x) valid for |x| <= 0.25 using a Taylor expansion around the origin and otherwise a rational approximation from the reference below. Spot-checked to give around 2e-16 relative accuracy. */

/* I. Koren and O. Zinaty. Evaluating elementary functions in a numerical

coprocessor based on rational approximations. IEEE Transactions on

Computers, Vol.39, No.8, August 1990, pp 1030-1037. */

/*

static int

sin_pi_koren(const double x, gsl_sf_result *result)

{

  result->val = 0.0;

  result->err = 0.0;

  if (16.0*fabs(x) < 1.0) {

    const double y = M_PI * x;

    const double a = y*y;

    result->val = y*(1.0 - a*(1.0 - a*(1.0 - a*(1.0 - a*(1.0 - a/110.0)/72.0)/42.0)/20.0)/6.0);

  }

  else {

    const double a0 = 1805490264.690988571178600370234394843221;

    const double a1 = -164384678.227499837726129612587952660511;

    const double a2 =    3664210.647581261810227924465160827365;

    const double a3 =     -28904.140246461781357223741935980097;

    const double a4 =         76.568981088717405810132543523682;

    const double b0 = 2298821602.638922662086487520330827251172;

    const double b1 =   27037050.118894436776624866648235591988;

    const double b2 =     155791.388546947693206469423979505671;

    const double b3 =        540.567501261284024767779280700089;

    const double t = 16.0*x*x;

    result->val = 4.0*x*(((( a4*t + a3 )*t + a2 )*t + a1 )*t + a0)/(((( t + b3 )*t + b2 )*t + b1 )*t + b0);

  }

  result->err = GSL_DBL_EPSILON*fabs(result->val);

  return GSL_SUCCESS;

}

*/

/* routine computing cos(pi*x) valid for |x| <= 0.25 using a Taylor expansion around the origin and otherwise a rational approximation from the reference below. Spot-checked to give around 2e-16 relative accuracy. */

/* I. Koren and O. Zinaty. Evaluating elementary functions in a numerical

coprocessor based on rational approximations. IEEE Transactions on

Computers, Vol.39, No.8, August 1990, pp 1030-1037. */

/*

static int

cos_pi_koren(const double x, gsl_sf_result *result)

{

  result->val = 0.0;

  result->err = 0.0;

  if (20.0*fabs(x) < 1.0) {

    const double y = M_PI * x;

    const double a = y*y;

    result->val = 1.0 - 0.5*a*(1.0 - a*(1.0 - a*(1.0 - a*(1.0 - a/90.0)/56.0)/30.0)/12.0);

  }

  else {

    const double a0 = 1090157078.174871420428849017262549038606;

    const double a1 = -321324810.993150712401352959397648541681;

    const double a2 =   12787876.849523878944051885325593878177;

    const double a3 =    -150026.206045948110568310887166405972;

    const double a4 =        538.333564203182661664319151379451;

    const double b0 = 1090157078.174871420428867295670039506886;

    const double b1 =   14907035.776643879767410969509628406502;

    const double b2 =     101855.811943661368302608146695082218;

    const double b3 =        429.772865107391823245671264489311;

    const double t = 16.0*x*x;

    result->val = (((( a4*t + a3 )*t + a2 )*t + a1 )*t + a0)/(((( t + b3 )*t + b2 )*t + b1 )*t + b0);

  }

  result->err = GSL_DBL_EPSILON*fabs(result->val);

  return GSL_SUCCESS;

}

*/

/* routine computing sin(pi*x) using a Taylor expansion around the origin and otherwise the library function. */

static int

sin_pi_taylor(const double x, gsl_sf_result *result)

{

  result->val = 0.0;

  result->err = 0.0;

  if (16.0*fabs(x) < 1.0) {

    const double y = M_PI * x;

    const double a = y*y;

    result->val = y*(1.0 - a*(1.0 - a*(1.0 - a*(1.0 - a*(1.0 - a/110.0)/72.0)/42.0)/20.0)/6.0);

  }

  else {

    result->val = sin(M_PI*x);

  }

  result->err = GSL_DBL_EPSILON*fabs(result->val);

  return GSL_SUCCESS;

}

/* routine computing sin(pi*x) using a Taylor expansion around the origin and otherwise the library function. */

static int

cos_pi_taylor(const double x, gsl_sf_result *result)

{

  result->val = 0.0;

  result->err = 0.0;

  if (20.0*fabs(x) < 1.0) {

    const double y = M_PI * x;

    const double a = y*y;

    result->val = 1.0 - 0.5*a*(1.0 - a*(1.0 - a*(1.0 - a*(1.0 - a/90.0)/56.0)/30.0)/12.0);

  }

  else {

    result->val = cos(M_PI*x);

  }

  result->err = GSL_DBL_EPSILON*fabs(result->val);

  return GSL_SUCCESS;

}

int

gsl_sf_sin_pi_e(const double x, gsl_sf_result *result)

{

  double intx = 0.0, fracx = 0.0;

  long q;

  int sign = 1, status;

  result->val = 0.0;

  result->err = 0.0;

  fracx = modf(x,&intx;;

  if (fracx == 0.0) return GSL_SUCCESS;

  if(fabs(intx) >= TWOBIG) return GSL_SUCCESS; /* to be sure. Actually should be covered by the line above */

  q = ( ( (intx >= LONG_MIN) && (intx <= LONG_MAX) ) ? intx : fmod(intx, 2.0) );

  sign = ( q % 2 ? -1 : 1 );

  /* int sign = 1 - 2*((int)round(fmod(fabs(intx),2.0))); */

  if (fabs(fracx) == 0.5) { /* probably unnecessary */

    if (fracx < 0.0) sign = -sign;

    result->val = ( sign != 1 ? -1.0 : 1.0 );

    return GSL_SUCCESS;

  }

  if (fabs(fracx) > 0.5) {

    sign = -sign;

    fracx = ( fracx > 0.0 ? fracx-1.0 : fracx+1.0 );

  }

  status = 0;

  if (fracx > 0.25) {

    status = cos_pi_taylor((fracx-0.5), result);

  }

  else if (fracx < -0.25) {

    status = cos_pi_taylor((fracx+0.5), result);

    sign = -sign;

  }

  else {

    status = sin_pi_taylor(fracx, result);

  }

  if (sign != 1) result->val = -result->val;

  return status;

}

int

gsl_sf_cos_pi_e(const double x, gsl_sf_result *result)

{

  double intx = 0.0, fracx = 0.0;

  long q;

  int sign = 1, status;

  result->val = 0.0;

  result->err = 0.0;

  fracx = modf(x,&intx;;

  if (fabs(fracx) == 0.5) return GSL_SUCCESS;

  if(fabs(intx) >= TWOBIG) {

    result->val = 1.0;

    return GSL_SUCCESS;

  }

  q = ( ( (intx >= LONG_MIN) && (intx <= LONG_MAX) ) ? intx : fmod(intx, 2.0) );

  sign = ( q % 2 ? -1 : 1 );

  /* int sign = 1 - 2*((int)round(fmod(fabs(intx),2.0))); */

  if (fracx == 0.0) { /* probably unnecessary */

    result->val = ( sign != 1 ? -1.0 : 1.0 );

    return GSL_SUCCESS;

  }

  if (fabs(fracx) > 0.5) {

    sign = -sign;

    fracx = ( fracx > 0.0 ? fracx-1.0 : fracx+1.0 );

  }

  status = 0;

  if (fracx > 0.25) {

    status = sin_pi_taylor((fracx-0.5), result);

    sign = -sign;

  }

  else if (fracx < -0.25) {

    status = sin_pi_taylor((fracx+0.5), result);

  }

  else {

    status = cos_pi_taylor(fracx, result);

  }

  if (sign != 1) result->val = -result->val;

  return status;

}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double

gsl_sf_sin_pi(const double x)

{

  EVAL_RESULT(gsl_sf_sin_pi_e(x, &result));

}

double

gsl_sf_cos_pi(const double x)

{

  EVAL_RESULT(gsl_sf_cos_pi_e(x, &result));

}