arcsine unction demo in c

#include "math.h"
#include "errno.h"

double arcsine();

double asin(x)
double x;
{
        return arcsine(x,0);
}

double acos(x)
double x;
{
        return arcsine(x,1);
}

#define P1 -0.27368494524164255994e+2
#define P2 +0.57208227877891731407e+2
#define P3 -0.39688862997504877339e+2
#define P4 +0.10152522233806463645e+2
#define P5 -0.69674573447350646411
#define Q0 -0.16421096714498560795e+3
#define Q1 +0.41714430248260412556e+3
#define Q2 -0.38186303361750149284e+3
#define Q3 +0.15095270841030604719e+3
#define Q4 -0.23823859153670238830e+2

#define P(g) ((((P5*g P4)*g P3)*g P2)*g P1)
#define Q(g) (((((g Q4)*g Q3)*g Q2)*g Q1)*g Q0)

double arcsine(x,flg)
double x;
{
        double y, g, r;
        register int i;
        extern int errno;
        static double a[2] = { 0.0, 0.78539816339744830962 };
        static double b[2] = { 1.57079632679489661923, 0.78539816339744830962 };

        y = fabs(x);
        i = flg;
        if (y < 2.3e-10)
                r = y;
        else {
                if (y > 0.5) {
                        i = 1-i;
                        if (y > 1.0) {
                                errno = EDOM;
                                return 0.0;
                        }
                        g = (0.5-y)+0.5;
                        g = ldexp(g,-1);
                        y = sqrt(g);
                        y = -(y+y);
                } else
                        g = y*y;
                r = y + y*
                                ((P(g)*g)
                                /Q(g));
        }
        if (flg) {
                if (x < 0.0)
                        r = (b[i] + r) + b[i];
                else
                        r = (a[i] - r) + a[i];
        } else {
                r = (a[i] + r) + a[i];
                if (x < 0.0)
                        r = -r;
        }
        return r;
}

Related Links :