/*
 *   4. For negative x, since (G is gamma function)
 *		-x*G(-x)*G(x) = pi/sin(pi*x),
 * 	we have
 * 		G(x) = pi/(sin(pi*x)*(-x)*G(-x))
 *	since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
 *	Hence, for x<0, signgam = sign(sin(pi*x)) and 
 *		lgamma(x) = log(|Gamma(x)|)
 *			  = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
 *	Note: one should avoid compute pi*(-x) directly in the 
 *	      computation of sin(pi*(-x)).
 */

/* sin(x)
 * Return sine function of x.
 *
 * kernel function:
 *  k_sin        ... sine function on [-pi/4,pi/4]
 *  k_cos        ... cos  function on [-pi/4,pi/4]
 *  __ieee754_rem_pio2  ... argument reduction routine
 *
 * Method.
 *      Let S and C denote the sin and cos respectively on 
 *  [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 
 *  in [-pi/4 , +pi/4], and let n = k mod 4.
 *  We have
 *
 *          n        sin(x)      cos(x)
 *     ----------------------------------------------------------
 *          0          S       C
 *          1          C      -S
 *          2         -S      -C
 *          3         -C       S
 *     ----------------------------------------------------------
 */

#include <stdint.h>
#include <math.h>

//#define MLA mlaf
//#define C2V(x) (x)
//#include "estrin.h"
#define C2V(x) (x)
#define POLY2(x, c1, c0) ((x)*C2V(c1) + C2V(c0))
#define POLY3(x, c2, c1, c0) ((C2V(c0) + (x)*(C2V(c1) + (x)*C2V(c2))))
#define POLY4(x, x2, c3, c2, c1, c0) ((C2V(c0) + (x)*C2V(c1)) + (x2)*(C2V(c2) + (x)*C2V(c3)))
#define POLY5(x, x2, x4, c4, c3, c2, c1, c0) ((x4)*C2V(c4) + POLY4(x, x2, c3, c2, c1, c0))
#define POLY6(x, x2, x4, c5, c4, c3, c2, c1, c0) ((x4)*POLY2(x, c5, c4) + POLY4(x, x2, c3, c2, c1, c0))
#define POLY7(x, x2, x4, c6, c5, c4, c3, c2, c1, c0) ((x4)*POLY3(x, c6, c5, c4) + POLY4(x, x2, c3, c2, c1, c0))
#define POLY8(x, x2, x4, c7, c6, c5, c4, c3, c2, c1, c0) ((x4)*POLY4(x, x2, c7, c6, c5, c4) + POLY4(x, x2, c3, c2, c1, c0))

static __attribute__((always_inline)) inline float mlaf(float x, float y, float z) { return x * y + z; }

static const float
C1  = -5.0000000000e-01,
C2  =  4.1666667908e-02, /* 0x3d2aaaab */
C3  = -1.3888889225e-03, /* 0xbab60b61 */
C4  =  2.4801587642e-05, /* 0x37d00d01 */
C5  = -2.7557314297e-07, /* 0xb493f27c */
C6  =  2.0875723372e-09, /* 0x310f74f6 */
C7  = -1.1359647598e-11; /* 0xad47d74e */

// Differs from libc cosf on [0, pi/2] at only the point 0x3fc460f8 (~1.5342093) by 0.0000001005828f
static inline float k_cosf(float x)
{
    float z,z2,z4,r;
    z  = x*x;
    z2 = z*z;
    z4 = z2*z2;

    r = POLY6(z, z2, z4, C7, C6, C5, C4, C3, C2);
    r = z*r;

    return 1.0f + (z*r + C1*z);
}

static const float
S1 = -0x1.555556p-03f,
S2 =  0x1.111111p-07f,
S3 = -0x1.a01906p-13f,
S4 =  0x1.7196aap-19f,
S5 = -0x1.9db448p-26f;
                                                                                                                                                                          
/* Approximate sine on [-PI/4,+PI/4]. Maximum ulp error = 0.64721
 *  Returns -0.0f for an argument of -0.0f
 *  Polynomial approximation based on unpublished work by T. Myklebust
 *  Note: S3 & S4 replaced with versions from a Remez Exchange minimax polynomial computation,
 *        and S5 added. S1 & S2 remain the same. Now this has error of at most 1ulp on [-PI/2,PI/2].
 *
 *        - bofh
 */
static inline float k_sinf(float x)
{
    float z,z2,z4,r;
    z  = x*x;
    z2 = z*z;
    z4 = z2*z2;

    r = POLY6(z, z2, z4, S5, S4, S3, S2, S1, 0.0f);
    return x+x*r;
}

static const float pi = 3.14159265358979323846f;

typedef union _f32 {
    float f;
    uint32_t i;
} _f32;

void sincosf_pi(float x, float *cosx, float *sinx)
{
    unsigned int sgnx = (x < 0.0f);
    unsigned int n0 = lrintf(fabsf(x));  /* n0 = int(|x| mod 1.0) */
	float y, s, c, sy, cy;
	int n = 0;
    _f32 tmp;

	y = fabsf(x);
    y = (y - (float)n0);

    /* argument reduction */
	//x  *= 0.5f;

    /* spurious inexact if odd int */
    //y = 2.0f*(x - lrintf(x));  /* y = |x| mod 2.0 */
    n = (int)(y*4.0f);
    n = ((n+1)>>1);
    y -= n*0.5f;

    s =  k_sinf(pi*y);
    c =  k_cosf(pi*y);
	sy = ((n&1) ?  c : s);
	cy = ((n&1) ? -s : c);

    sy = (sgnx ? -sy : sy);

	*sinx = sy;
	*cosx = cy;

handle_signs:
    n0 <<= 31;
    //tmp.f = *cosx; tmp.i ^= n0; *cosx = tmp.f; 
    //tmp.f = *sinx; tmp.i ^= n0; *sinx = tmp.f; 
    //if (n0&1) { *sinx = -*sinx; *cosx = -*cosx; }
    //tmp.f = cy; tmp.i ^= n0; *cosx = tmp.f; 
    //tmp.f = sy; tmp.i ^= n0; *sinx = tmp.f; 
    tmp.f = cy; tmp.i ^= n0; cy = tmp.f; 
    tmp.f = sy; tmp.i ^= n0; sy = tmp.f; 
	*sinx = sy;
	*cosx = cy;
}

float cosf_pi(float x)
{
	unsigned int n = 0;
	float y;
    _f32 c;

	x = fabsf(x);

    /* argument reduction */
    n = lrintf(x);
    y = (x - (float)n);  /* y = |x| mod 1.0 */

    c.f = k_cosf(pi*y);
    /* both gcc and clang keep stubbornly compiling this to a branch, so we need to handhold their hands. sigh.
    if (n&1) c = -c;*/
    c.i ^= (n << 31);
    return c.f;
}

float sinf_pi(float x)
{
    int n = lrintf(x);
	float y;
    _f32 s;

    y = (x - (float)n);  /* y = |x| mod 1.0 */

    s.f = k_sinf(pi*y);
    /* both gcc and clang keep stubbornly compiling this to a branch, so we need to handhold their hands. sigh.
    if (n&1) s = -s;*/
    s.i ^= (n << 31);
    return s.f;
}

#ifdef TEST

#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <sys/time.h>

static float rel_err_diff(float x1, float x2)
{
  if(x1 == 0.0f && x2 == 0.0f) {
    return 0.0f;
  } else if (x1 == 0.0f) {
    return fabsf(x2);
  } else if(x1 <= 1e38f && x2 <= 1e38f && (x1 + x2 != 0.0f)) {
    //return fabsf((x1-x2)/(x1+x2));
    return fabsf(x1-x2);
  } else {
    return 1.0f;
  }
}

static __attribute__((always_inline)) inline void disable_denormals(void)
{
    unsigned int a = 0x9f80;
    __asm__ __volatile__ ("ldmxcsr %0" : "+m" (a));
}

int main(void) {
    _f32 f, g1, g2;
    unsigned int i, j, nerrs = 0;
    float s, c, diff, maxdiff = 0.0f, maxdiff2 = 0.0f;

    disable_denormals();
    setvbuf(stdout, NULL, _IONBF, 0);

    nerrs = 0;
    maxdiff = 0.0f;
    f.i = 0x3f800000;
    printf("start: %.7f\n", f.f-1.0f);
    for (i = 0; i < (1 << 23U); i++) {
        f.i = 0x3f800000|i;
        f.f -= 1.0f;
        g1.f = cosf(pi*f.f);
        sincosf_pi(f.f, &g2.f, &s);
        if ((diff = rel_err_diff(g1.f, g2.f)) >=  0x1.0p-23f) {
            nerrs++;
            if (diff > fabsf(maxdiff)) maxdiff = diff;
            if (!(i & 0xffff))
                printf("Error: cosf(%.7g): %.7g != k_cosf(%.7g): %.7g (diff: %.7e / 0x%08x)\n", f.f, g1.f, f.f, g2.f, diff, (g1.i ^ g2.i));
        }
    }
    printf("end: %.7f\n", f.f);
    printf("Total errors, cosf: %u, largest error: %.7e\n", nerrs, maxdiff); nerrs = 0;

    nerrs = 0;
    maxdiff = 0.0f;
    f.i = 0x3f800000;
    printf("start: %.7f\n", f.f-1.0f);
    for (i = 0; i < (1 << 23U); i++) {
        f.i = 0x3f800000|i;
        f.f -= 1.0f;
        g1.f = sinf(pi*f.f);
        sincosf_pi(f.f, &s, &g2.f);
        if ((diff = rel_err_diff(g1.f, g2.f)) >=  0x1.0p-23f) {
            nerrs++;
            if (diff > fabsf(maxdiff)) maxdiff = diff;
            if (!(i & 0xffff))
                printf("Error: sinf(%.7g): %.7g != k_sinf(%.7g): %.7g (diff: %.7e / 0x%08x)\n", f.f, g1.f, f.f, g2.f, diff, (g1.i ^ g2.i));
        }
    }
    printf("end: %.7f\n", f.f);
    printf("Total errors, sinf: %u, largest error: %.7e\n", nerrs, maxdiff); nerrs = 0;
#if 0
    printf("start: %.7f\n", f.f);
    for (i = 0; i < (1 << 23U); i++) {
        f.i = 0x3f800000|i;
        g1.f = cosf(pi*f.f);
        sincosf_pi(f.f, &s, &g2.f);
        if ((diff = rel_err_diff(g1.f, g2.f)) >= 0x1.0p-23f) {
            nerrs++;
            if (diff > fabsf(maxdiff)) maxdiff = diff;
            if (!(i & 0xffff))
                printf("Error: cosf(%.7g): %.7g != k_cosf(%.7g): %.7g (diff: %.7e / 0x%08x)\n", f.f, g1.f, f.f, g2.f, diff, (g1.i ^ g2.i));
        }
    }
    printf("end: %.7f\n", f.f);
    printf("Total errors: %u, largest error: %.7e\n", nerrs, maxdiff); nerrs = 0;
    printf("start: %.7f\n", f.f);
    for (i = 0; i < (1 << 23U); i++) {
        f.i = 0x40000000|i;
        g1.f = cosf(pi*f.f);
        sincosf_pi(f.f, &s, &g2.f);
        if ((diff = rel_err_diff(g1.f, g2.f)) >= 0x1.0p-23f) {
            nerrs++;
            if (diff > fabsf(maxdiff)) maxdiff = diff;
            if (!(i & 0xffff) && ((g1.i ^ g2.i) > 3))
                printf("Error: cosf(%.7g): %.7g != k_cosf(%.7g): %.7g (diff: %.7e / 0x%08x)\n", f.f, g1.f, f.f, g2.f, diff, (g1.i ^ g2.i));
        }
    }
    printf("end: %.7f\n", f.f);
    printf("Total errors: %u, largest error: %.7e\n", nerrs, maxdiff);
    nerrs = 0;
    maxdiff = 0.0f;
    f.i = 0x3f800000;
    printf("start: %.7f\n", f.f-1.0f);
    for (i = 0; i < (1 << 23U); i++) {
        f.i = 0x3f800000|i;
        f.f -= 1.0f;
        g1.f = sinf(pi*f.f);
        sincosf_pi(f.f, &g2.f, &c);
        if ((diff = rel_err_diff(g1.f, g2.f)) >= 0x1.0p-23f) {
            nerrs++;
            if (diff > fabsf(maxdiff)) maxdiff = diff;
            if (!(i & 0xffff))
                printf("Error: sinf(%.7g): %.7g != k_sinf(%.7g): %.7g (diff: %.7e / 0x%08x)\n", f.f, g1.f, f.f, g2.f, diff, (g1.i ^ g2.i));
        }
    }
    printf("end: %.7f\n", f.f);
    printf("Total errors: %u, largest error: %.7e\n", nerrs, maxdiff); nerrs = 0;
    printf("start: %.7f\n", f.f);
    for (i = 0; i < (1 << 23U); i++) {
        f.i = 0x3f800000|i;
        g1.f = sinf(pi*f.f);
        sincosf_pi(f.f, &g2.f, &c);
        if ((diff = rel_err_diff(g1.f, g2.f)) >= 0x1.0p-23f) {
            nerrs++;
            if (diff > fabsf(maxdiff)) maxdiff = diff;
            if (!(i & 0xffff))
                printf("Error: sinf(%.7g): %.7g != k_sinf(%.7g): %.7g (diff: %.7e / 0x%08x)\n", f.f, g1.f, f.f, g2.f, diff, (g1.i ^ g2.i));
        }
    }
    printf("end: %.7f\n", f.f);
    printf("Total errors: %u, largest error: %.7e\n", nerrs, maxdiff); nerrs = 0;
    printf("start: %.7f\n", f.f);
    for (i = 0; i < (1 << 23U); i++) {
        f.i = 0x40000000|i;
        g1.f = sinf(pi*f.f);
        sincosf_pi(f.f, &g2.f, &c);
        if ((diff = rel_err_diff(g1.f, g2.f)) >= 100000) {
            nerrs++;
            if (diff > fabsf(maxdiff)) maxdiff = diff;
            if (!(i & 0xffff))
                printf("Error: sinf(%.7g): %.7g != k_sinf(%.7g): %.7g (diff: %.7e / 0x%08x)\n", f.f, g1.f, f.f, g2.f, diff, (g1.i ^ g2.i));
        }
    }
    printf("end: %.7f\n", f.f);
    printf("Total errors: %u, largest error: %.7e\n", nerrs, maxdiff);
#endif

#if 0
    nerrs = 0;
    for (j = 120; j < 255; j++) {
        printf("start: %.7f\n", f.f);
        for (i = 0; i < (1 << 23U); i++) {
            f.i = (j << 23)|i;
            g1.f = cosf(pi*f.f);
            sincosf_pi(f.f, &s, &g2.f);
            if ((diff = rel_err_diff(g1.f, g2.f)) >= 2e-7f) {
                nerrs++;
                if (diff > fabsf(maxdiff)) maxdiff = diff;
                if (!(i & 0xffff))
                    printf("Error: cosf(%.7g): %.7g != k_cosf(%.7g): %.7g (diff: %.7e / 0x%08x)\n", f.f, g1.f, f.f, g2.f, diff, (g1.i ^ g2.i));
            }
        }
        printf("end: %.7f\n", f.f);
    }
    printf("Total errors (cos_pi): %u, largest error: %.7e\n", nerrs, maxdiff);
    nerrs = 0;
    for (j = 120; j < 255; j++) {
        printf("start: %.7f\n", f.f);
        for (i = 0; i < (1 << 23U); i++) {
            f.i = (j << 23)|i;
            g1.f = sinf(pi*f.f);
            sincosf_pi(f.f, &g2.f, &c);
            if ((diff = rel_err_diff(g1.f, g2.f)) >= 2e-7f) {
                nerrs++;
                if (diff > fabsf(maxdiff)) maxdiff = diff;
                if (!(i & 0xffff))
                    printf("Error: sinf(%.7g): %.7g != k_sinf(%.7g): %.7g (diff: %.7e / 0x%08x)\n", f.f, g1.f, f.f, g2.f, diff, (g1.i ^ g2.i));
            }
        }
        printf("end: %.7f\n", f.f);
    }
    printf("Total errors (sin_pi): %u, largest error: %.7e\n", nerrs, maxdiff);
#endif

    nerrs = 0;
    maxdiff = 0.0f;
    for (j = 120; j < 151; j++) {
    //f.i = 0x3f800000;
    //printf("start: %.7f\n", f.f-1.0f);
    f.i = (j<<23);
        printf("start: %.7f (j: %u)\n", f.f, j);
    for (i = 0; i < (1 << 23U); i++) {
        //f.i = 0x3f800000|i;
        //f.f -= 1.0f;
        f.i = (j<<23)|i;
        //g1.f = cosf(pi*f.f);
        g1.f = (float)cos(M_PI*(double)f.f);
        c = cosf_pi(f.f);
        if ((diff = rel_err_diff(c, g1.f)) >=  0x1.0p-22f) {
            nerrs++;
            if (diff > fabsf(maxdiff)) maxdiff = diff;
            if (!(i & 0xffff))
                printf("Error: cosf(%.7g): %.7g != k_cosf(%.7g): %.7g (diff: %.7e / %a\n", f.f, c, f.f, g1.f, diff, diff);
        }
    }
    printf("end: %.7f\n", f.f);
    }
    printf("Total errors: %u, largest error: cos: %.7e\n", nerrs, maxdiff);

    nerrs = 0;
    maxdiff = 0.0f;
    maxdiff2 = 0.0f;
    for (j = 120; j < 151; j++) {
    //f.i = 0x3f800000;
    //printf("start: %.7f\n", f.f-1.0f);
    f.i = (j<<23);
        printf("start: %.7f (j: %u)\n", f.f, j);
    for (i = 0; i < (1 << 23U); i++) {
        //f.i = 0x3f800000|i;
        //f.f -= 1.0f;
        f.i = (j<<23)|i;
        //g1.f = cosf(pi*f.f);
        //g2.f = sinf(pi*f.f);
        g1.f = (float)cos(M_PI*(double)f.f);
        g2.f = (float)sin(M_PI*(double)f.f);
        sincosf_pi(f.f, &s, &c);
        if ((diff = rel_err_diff(c, g1.f)) >=  0x1.0p-22f) {
            nerrs++;
            if (diff > fabsf(maxdiff)) maxdiff = diff;
            if (!(i & 0xffff))
                printf("Error: cosf(%.7g): %.7g != k_cosf(%.7g): %.7g (diff: %.7e / %a\n", f.f, c, f.f, g1.f, diff, diff);
        }
        if ((diff = rel_err_diff(s, g2.f)) >=  0x1.0p-22f) {
            nerrs++;
            if (diff > fabsf(maxdiff)) maxdiff2 = diff;
            if (!(i & 0xffff))
                printf("Error: sinf(%.7g): %.7g != k_sinf(%.7g): %.7g (diff: %.7e / %a\n", f.f, s, f.f, g2.f, diff, diff);
        }
    }
    printf("end: %.7f\n", f.f);
    }
    printf("Total errors: %u, largest error: cos: %.7e, sin: %.7e\n", nerrs, maxdiff, maxdiff2);

    return 0;
}

#endif