/* e_asinf.c -- float version of e_asin.c.
 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
 */

/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 */

#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: e_asinf.c,v 1.5 1995/05/12 04:57:25 jtc Exp $";
#endif

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

static const float 
one =  1.0000000000e+00, /* 0x3F800000 */
pio2_hi =  1.5707962513e+00, /* 0x3fc90fda */
pio2_lo =  7.5497894159e-08; /* 0x33a22168 */

#if 0
/* coefficient for R(x^2) */
static const float 
pS0 =  1.6666667163e-01, /* 0x3e2aaaab */
pS1 = -3.2556581497e-01, /* 0xbea6b090 */
pS2 =  2.0121252537e-01, /* 0x3e4e0aa8 */
pS3 = -4.0055535734e-02, /* 0xbd241146 */
pS4 =  7.9153501429e-04, /* 0x3a4f7f04 */
pS5 =  3.4793309169e-05, /* 0x3811ef08 */
qS1 = -2.4033949375e+00, /* 0xc019d139 */
qS2 =  2.0209457874e+00, /* 0x4001572d */
qS3 = -6.8828397989e-01, /* 0xbf303361 */
qS4 =  7.7038154006e-02; /* 0x3d9dc62e */

float R2(float t)
{	
    float p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
	float q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
	float w = p/q;
    return w;
}
#endif

float R2(float t)
{	
    const float
    c0=0.33333841285043691065,
    c1=0.14985117853050373729,
    c2=0.091268980381668241081,
    c3=0.047723382122511172410,
    c4=0.084263405223358602512;
	float w = t*(c0 + t*(c1 + t*(c2 + t*(c3 + t*c4))));
    return 0.5f*w;
}

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

float fast_asinf(float x)
{
	float t,w,s;
    _f32 ix;

    ix.f = fabsf(x);
	if(ix.i==0x3f800000) {
		/* asin(1)=+-pi/2 with inexact */
	    return x*(pio2_hi+pio2_lo);
	} else if(ix.i> 0x3f800000) {	/* |x|>= 1 */
	    return 0.0f/0.0f;		/* asin(|x|>1) is NaN */   
	} else if (ix.i<0x3f000000) {	/* |x|<0.5 */
	    if(ix.i<0x32000000) {		/* if |x| < 2**-27 */
		    return x;/* return x with inexact if x!=0*/
	    } 
		t = x*x;
		w = R2(t);
		return x+x*w;
	}

	/* 1> |x|>= 0.5 */
	t = 0.5f*(1.0f-fabsf(x));
	w = R2(t);
	s = sqrtf(t);
	t = pio2_hi-(2.0f*(s+s*w)-pio2_lo);
	if(x<0) { t = -t; }
    return t;
}

#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, nerrs = 0;
    float s, c, diff, maxdiff = 0.0f;

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

    nerrs = 0;
    maxdiff = 0.0f;
    f.i = 0x3f000000;
    printf("start: %.7f\n", f.f-0.5f);
    for (i = 0; i < (1 << 23U); i++) {
        f.i = 0x3f000000|i;
        f.f -= 0.5f;
        g1.f = asinf(f.f);
        c = fast_asinf(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: asinf(%.7g): %.7g != fast_asinf(%.7g): %.7g (diff: %.7e / %a\n", f.f, g1.f, f.f, c, diff, diff);
        }
    }
    printf("end: %.7f\n", f.f);
    printf("Total errors: %u, largest error: cos: %.7e\n", nerrs, maxdiff);

    nerrs = 0;
    maxdiff = 0.0f;
    f.i = 0x3f000000;
    printf("start: %.7f\n", f.f);
    for (i = 0; i < (1 << 23U); i++) {
        f.i = 0x3f000000|i;
        g1.f = asinf(f.f);
        c = fast_asinf(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: asinf(%.7g): %.7g != fast_asinf(%.7g): %.7g (diff: %.7e / %a\n", f.f, g1.f, f.f, c, diff, diff);
        }
    }
    printf("end: %.7f\n", f.f);
    printf("Total errors: %u, largest error: cos: %.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 = asinf(f.f);
        c = fast_asinf(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: asinf(%.7g): %.7g != fast_asinf(%.7g): %.7g (diff: %.7e / %a\n", f.f, g1.f, f.f, c, diff, diff);
        }
    }
    printf("end: %.7f\n", f.f);
    printf("Total errors: %u, largest error: cos: %.7e\n", nerrs, maxdiff);

    return 0;
}

#endif