Approximations of cosPi and sinPi, C++ (gcc)

Language: Layout:

// I hereby release this code into the public domain.

#include <cstdio>
#include <cmath>
#include <cfenv>
#include <limits>

// Computes cos(pi*x) for |x|<1/4.
static inline double k_cosPi(double x)
{
    double x2=x*x;
    // Coefficients computed via minimizing L2-norm of cos(pi*x)-(1+P(x^2)).
    return 1.0+
        x2*(-4.934802200544679276620401+
        x2*(+4.058712126416754257530665+
        x2*(-1.335262768852578716599486+
        x2*(+2.353306302218198462156966e-1+
        x2*(-2.580688640321055216598382e-2+
        x2*(+1.929474803148289697385952e-3+
        x2*(-1.036104942062675021343889e-4)))))));
}

// Computes sin(pi*x) for |x|<1/4.
static inline double k_sinPi(double x)
{
    double x2=x*x;
    // Coefficients computed via minimizing L2-norm of sin(pi*x)-x*(pi+P(x^2)).
    return 2.0*x+
        x*(1.141592653589793238462643+ // pi split in 2 parts, for extra accuracy.
        x2*(-5.167712780049953463142196+
        x2*(+2.55016403987302778905671+
        x2*(-5.992645289176800629240315e-1+
        x2*(+8.214586862112711268943473e-2+
        x2*(-7.370015454628586641275853e-3+
        x2*(+4.615085719100285309744452e-4)))))));
}

// Returns approximation of cos(pi*x).
// cosPi(n)=-1^n
// cosPi(n+1/2)=+0
// |cosPi(x)|<=1
// Estimated error bound: 1.5 ulp (NOT TESTED).
double cosPi(double x)
{
    if(x<0.0) return cosPi(-x);                        // cosPi(-x) is the same as cosPi(x).
    if(!(x<1.0e+16)) return 1.0+(x-x);                 // All doubles > 10^16 are even integers. Also handles inf and NaN.
    if(x<1e-50) return 1.0-x;                          // Avoid underflow. Handle inexact.
    long long n=(long long)(1.0+4.0*x)>>1;             // 2*x rounded to nearest integer; avoids raising INEXACT for (half)integers.
    x-=0.5*double(n);                                  // Exact, no rounding.
    switch(n&3)
    {
        case 0: return +k_cosPi(x);
        case 1: return -k_sinPi(x);
        case 2: return -k_cosPi(x);
        case 3: return +k_sinPi(x);
    }
    return 0.0/0.0;
}

// Returns approximation of sin(pi*x).
// sinPi(n)=0, same sign as n.
// sinPi(n+1/2)=-1^(n+1)
// |sinPi(x)|<=1
// Estimated error bound: 1.5 ulp (NOT TESTED).
double sinPi(double x)
{
    if(x<0.0) return -sinPi(-x);                       // sinPi(-x) is the same as -sinPi(x).
    if(!(x<1.0e+16)) return x-x;                       // All doubles > 10^16 are even integers. Also handles inf and NaN.
    if(x<1e-50) return 3.141592653589793238462643*x;   // Handle underflow and inexact.
    long long n=(long long)(1.0+4.0*x)>>1;             // 2*x rounded to nearest integer; avoids raising INEXACT for (half)integers.
    x-=0.5*double(n);                                  // Exact, no rounding.
    switch(n&3)
    {
        case 0: return +k_sinPi(x);
        case 1: return +k_cosPi(x);
        case 2: return -k_sinPi(x);
        case 3: return -k_cosPi(x);
    }
    return 0.0/0.0;
}

static void print_exceptions(int e)
{
    struct {int a;const char *b;} s[]={
        {FE_DIVBYZERO,"FE_DIVBYZERO"},
        {FE_INEXACT,"FE_INEXACT"},
        {FE_INVALID,"FE_INVALID"},
        {FE_OVERFLOW,"FE_OVERFLOW"},
        {FE_UNDERFLOW,"FE_UNDERFLOW"}
    };
    for(int i=0,n=sizeof(s)/sizeof(s[0]);i<n;++i)
    {
        if(e&s[i].a) std::printf("%s ",s[i].b);
    }
}

int main()
{
    //std::printf("Testing cosPi:\n");
    //std::printf("    x       |cosPi(x)-cos(pi*x)| in ulp\n");
    double e;
    e=0.0;
    int N=5000000;
    for(int i=0;i<N;i+=1)
    {
        double x=double(i)/double(N);
        double a=cosPi(x);
        long double b=cosl(3.141592653589793238462643383279502884197169399375L*x);
        if(std::fabs(x-0.5)<0.25) b=sinl(3.141592653589793238462643383279502884197169399375L*(0.5-x));
        double c=std::fabs(a-b);
        double d=std::abs(b-std::nextafter(double(b),2.0*double(b)));
        if(d==0.0) d=std::nextafter(0.0,1.0);
        if(c/d>e) e=c/d;
        // if(i%1000==0) std::printf("%5.3f      %5.2f\n",x,c/d);
    }
    std::printf("Max. measured error cosPi: %.2f ulp\n",e);
    //std::printf("Testing sinPi:\n");
    //std::printf("    x       |sinPi(x)-sin(pi*x)| in ulp\n");
    e=0.0;
    for(int i=0;i<N;i+=1)
    {
        double x=double(i)/double(N);
        double a=sinPi(x);
        long double b=sinl(3.141592653589793238462643383279502884197169399375L*x);
        if(std::fabs(x-0.5)<0.25) b=cosl(3.141592653589793238462643383279502884197169399375L*(0.5-x));
        if(x>0.75) b=sinl(3.141592653589793238462643383279502884197169399375L*(1.0-x));
        double c=std::fabs(a-b);
        double d=std::abs(b-std::nextafter(double(b),2.0*double(b)));
        if(d==0.0) d=std::nextafter(0.0,1.0);
        if(c/d>e) e=c/d;
        // if(i%1000==0) std::printf("%5.3f      %5.2f\n",x,c/d);
    }
    std::printf("Max. measured error sinPi: %.2f ulp\n",e);

if(0){
    for(double x=-70.0;x<713.0;x+=0.0037)
        std::printf("%s","!"+(std::fabs(cosPi(x)-cos(3.14159265358979323846*x))<0.000000000001));
    for(double x=-70.0;x<713.0;x+=0.0037)
       std::printf("%s","!"+(std::fabs(sinPi(x)-sin(3.14159265358979323846*x))<0.000000000001));
    }

double x[]={+0.0,-0.0,+0.7,-0.7,+10.0,-10.0,+73.5,-73.5,+1.5,-1.5,+1e-170,-1e-170,+1e-310,-1e-310,1e20,-1e20,1.0/0.0,-1.0/0.0,0.0/0.0,std::numeric_limits<double>::signaling_NaN()};
    int K=sizeof(x)/sizeof(x[0]);
    for(int i=0;i<K;++i)
    {
        std::feclearexcept(FE_ALL_EXCEPT);
        double y=cosPi(x[i]);
        int e=std::fetestexcept(FE_ALL_EXCEPT);
        if(i<K-1) std::printf("cosPi(%g)=%.17g ",x[i],y);
        else std::printf("cosPi(sNaN)=%g ",y);
        print_exceptions(e);
        std::printf("\n");
    }
    for(int i=0;i<K;++i)
    {
        std::feclearexcept(FE_ALL_EXCEPT);
        double y=sinPi(x[i]);
        int e=std::fetestexcept(FE_ALL_EXCEPT);
        if(i<K-1) std::printf("sinPi(%g)=%.17g ",x[i],y);
        else std::printf("sinPi(sNaN)=%g ",y);
        print_exceptions(e);
        std::printf("\n");
    }
    return 0;
}

Show compiler warnings [ + ] Compiler args [ + ] Show input
Compilation time: 0.45 sec, absolute running time: 4.63 sec, cpu time: 5.96 sec, memory peak: 3 Mb, absolute service time: 5,09 sec