#include <math.h>
#include <heport.h>
#include <hemath.h>

#ifndef _GLIBCXX_CMATH
#define NPOLY 4

// Check and verification is in /home/cc/qq/qqsin

float sin22p5(const float x) //  0.0 <= x <= 45.0 degrees expressed as 0...1.0
{
    static float apolys[NPOLY] ={
        -0.0000705737583431, // apolys[0]
        0.7868158711878633, // apolys[1]
        -0.0057807310387922, // apolys[2]
        -0.0739438871780882 // apolys[3]
    };
    float y;
    uint i;

    y = apolys[NPOLY - 1];
    for (i = NPOLY - 1; i;) {
        --i;
        y *= x;
        y += apolys[i];
    }

    return y;
}

float cos22p5(const float x) //  0.0 <= x <= 45.0 Grad, ausgedrueckt als 0...1.0
{
    static float apolyc[NPOLY] = // Berechnet und getestet durch qq/qqsin.cpp
    {
        0.9998076596624137, // apolyc[0]
        0.0041579158400429, // apolyc[1]
        -0.3276732635797221, // apolyc[2]
        0.0306285609237086 // apolyc[3]
    };
    float y;
    uint i;

    y = apolyc[NPOLY - 1];
    for (i = NPOLY - 1; i;) {
        --i;
        y *= x;
        y += apolyc[i];
    }

    return y;
}






// This sine routine splits the circle in 8 parts, each quadrant in two parts

float sin(const float x) // x as radians
{
    // requires for hr1 = sin(-10.0) ; 4519 cycles =282 usec
    float w, r;
    u32 u;
    bool neg;
    u8 segment;
    float y;

    neg = false;
    r = x * 0.15915494309189533576888376337251F; // 1.0/(2.0*cPI) ;
    if (r < 0.0F) {
        neg = true;
        r = -r;
    }
    u = (u32) r;
    w = r - u; // 0...360 degrees = 0.0.99999
    w = fmulpow2(w, 3); // w = w * 8.0F ;
    segment = (u8) (u32) w; // 0...7

    switch (segment) {
        case 0: y = sin22p5(w);
            break;
        case 1: y = cos22p5(2.0F - w);
            break;
        case 2: y = cos22p5(w - 2.0F);
            break;
        case 3: y = sin22p5(4.0F - w);
            break;
        case 4: y = -sin22p5(w - 4.0F);
            break;
        case 5: y = -cos22p5(6.0F - w);
            break;
        case 6: y = -cos22p5(w - 6.0F);
            break;
        default: y = -sin22p5(8.0F - w);
            break;
    }
    if (neg) return -y;
    return y;
}

FUNCT float cos(const float x) // x as radians
{
    return sin((const float) (x + 1.5707963267948966192313216916398F)); // PI / 2.0
}
#endif