/*
 *   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>

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,r;
    z  = x*x;
    r  = z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*C7)))));
    return 1.0f + (z*r + C1*z);
}

static const float pi = 3.14159265358979323846f;

float cosf_pi(float x)
{
	float y, s, c;
	int n = 0;

	x = fabsf(x);

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

    c = k_cosf(pi*y);
    if (n&1) c = -c;
    return c;
}

#ifdef TEST

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

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

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;

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

    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);

    return 0;
}

#endif