/* e_powf.c -- float version of e_pow.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.
 * ====================================================
 */

#include <ansi_parms.h>

#if _EWL_FLOATING_POINT && _EWL_USES_SUN_SP_MATH_LIB

#include <fdlibm.h>
#include <fenv.h>

#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: e_powf.c,v 1.8 1997/10/09 11:29:42 lukem Exp $");
#endif


#ifdef __STDC__
static const f32_t
#else
static f32_t
#endif
bp[] = {1.0F, 1.5F,},
dp_h[] = { 0.0F, 5.84960938e-01F,}, /* 0x3f15c000 */
dp_l[] = { 0.0F, 1.56322085e-06F,}, /* 0x35d1cfdc */
zero    =  0.0F,
one	=  1.0F,
two	=  2.0F,
two24	=  16777216.0F,	/* 0x4b800000 */
hugeval	=  1.0e30F,
tiny    =  1.0e-30F,
	/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
L1  =  6.0000002384e-01F, /* 0x3f19999a */
L2  =  4.2857143283e-01F, /* 0x3edb6db7 */
L3  =  3.3333334327e-01F, /* 0x3eaaaaab */
L4  =  2.7272811532e-01F, /* 0x3e8ba305 */
L5  =  2.3066075146e-01F, /* 0x3e6c3255 */
L6  =  2.0697501302e-01F, /* 0x3e53f142 */
P1   =  1.6666667163e-01F, /* 0x3e2aaaab */
P2   = -2.7777778450e-03F, /* 0xbb360b61 */
P3   =  6.6137559770e-05F, /* 0x388ab355 */
P4   = -1.6533901999e-06F, /* 0xb5ddea0e */
P5   =  4.1381369442e-08F, /* 0x3331bb4c */
lg2  =  6.9314718246e-01F, /* 0x3f317218 */
lg2_h  =  6.93145752e-01F, /* 0x3f317200 */
lg2_l  =  1.42860654e-06F, /* 0x35bfbe8c */
ovt =  4.2995665694e-08F, /* -(128-log2(ovfl+.5ulp)) */
cp    =  9.6179670095e-01F, /* 0x3f76384f =2/(3ln2) */
cp_h  =  9.6179199219e-01F, /* 0x3f763800 =head of cp */
cp_l  =  4.7017383622e-06F, /* 0x369dc3a0 =tail of cp_h */
ivln2    =  1.4426950216e+00F, /* 0x3fb8aa3b =1/ln2 */
ivln2_h  =  1.4426879883e+00F, /* 0x3fb8aa00 =16b 1/ln2*/
ivln2_l  =  7.0526075433e-06F; /* 0x36eca570 =1/ln2 tail*/

MISRA_EXCEPTION_RULE_10_1()
MISRA_EXCEPTION_RULE_12_7()
#ifdef __STDC__
	f32_t __ieee754_powf(f32_t x, f32_t y)
#else
	f32_t __ieee754_powf(x,y)
	f32_t x, y;
#endif
{
	f32_t z,ax,z_h,z_l,p_h,p_l;
	f32_t y1,t1,t2,r,s,sn,t,u,v,w;
	int32_t i,j,k,yisint,n;
	int32_t hx,hy,is;
	uint32_t ix,iy;

#if _EWL_IEC_559_ADDITIONS
	MISRA_EXCEPTION_RULE_13_3()
	if (x == one) {
		MISRA_EXCEPTION_RULE_14_7()
		return one;
	}
	MISRA_EXCEPTION_RULE_13_3()
	MISRA_EXCEPTION_RULE_12_4()
	if ((x == -one) && isinf(y)) {
		MISRA_EXCEPTION_RULE_14_7()
		return one;
	}
#endif /* _EWL_IEC_559_ADDITIONS */
	GET_FLOAT_WORD(hx,x);
	GET_FLOAT_WORD(hy,y);
	ix = (uint32_t)hx&0x7fffffffUL;
	iy = (uint32_t)hy&0x7fffffffUL;

	#if _EWL_C99
		if ((iy == 0u) && (ix == 0u)) {	/* domain error may occur */
			MISRA_EXCEPTION_RULE_13_7()
			if ((uint_t)math_errhandling & (uint_t)MATH_ERRNO) {
				MISRA_EXCEPTION_RULE_20_5()
				errno=EDOM;
			}
			MISRA_EXCEPTION_RULE_13_7()
			if ((uint_t)math_errhandling & (uint_t)MATH_ERREXCEPT) {
				MISRA_EXCEPTION_RULE_20_5()
				feraiseexcept((int_t)FE_INVALID);
			}
		}		
		if ((ix == 0u) && (hy < 0)) {	/* domain error may occur */
			MISRA_EXCEPTION_RULE_13_7()
			if ((uint_t)math_errhandling & (uint_t)MATH_ERRNO) {
				MISRA_EXCEPTION_RULE_20_5()
				errno=EDOM;
			}
			MISRA_EXCEPTION_RULE_13_7()
			if ((uint_t)math_errhandling & (uint_t)MATH_ERREXCEPT) {
				MISRA_EXCEPTION_RULE_20_5()
				feraiseexcept((int_t)FE_INVALID);
			}
		}
	#endif

    /* y==zero: x**0 = 1 */
	if(iy==0UL) {
		MISRA_EXCEPTION_RULE_14_7()
		return one;
	}

    /* +-NaN return x+y */
	if((ix > 0x7f800000UL) ||
	   (iy > 0x7f800000UL))
	{
		#if _EWL_C99
			MISRA_EXCEPTION_RULE_13_7()
			if ((uint_t)math_errhandling & (uint_t)MATH_ERRNO) {
				MISRA_EXCEPTION_RULE_20_5()
				errno=ERANGE;
			}
		#else
		#ifdef __STDC__
			MISRA_EXCEPTION_RULE_20_5()
			errno=EDOM;
		#endif
		#endif
		MISRA_EXCEPTION_RULE_14_7()
		return x+y;
	}
    /* determine if y is an odd int when x < 0
     * yisint = 0	... y is not an integer
     * yisint = 1	... y is an odd int
     * yisint = 2	... y is an even int
     */
	yisint  = 0L;
	if(hx<0L) {
	    if(iy>=0x4b800000UL) {
	    	yisint = 2L; /* even integer y */
	    } else if(iy>=0x3f800000UL) {
			k = (int32_t)(iy>>23)-0x7fL;	   /* exponent */
			j = (int32_t)(iy>>(23-k));
			if((j<<(23-k))==iy) {
				yisint = 2L-(j&1);
			}
	    } else {}
	}

    /* special value of y */
	if (iy==0x7f800000UL) {	/* y is +-inf */
	    if (ix==0x3f800000UL) {
			MISRA_EXCEPTION_RULE_14_7()
	        return  y - y;	/* inf**+-1 is NaN */
	    } else if (ix > 0x3f800000UL) {/* (|x|>1)**+-inf = inf,0 */
			MISRA_EXCEPTION_RULE_14_7()
	        return (hy>=0L)? y: zero;
	    } else {			/* (|x|<1)**-,+inf = inf,0 */
			MISRA_EXCEPTION_RULE_14_7()
	        return (hy<0)?-y: zero;
		}
	}
	if(iy==0x3f800000UL) {	/* y is  +-1 */
	    if(hy<0L) {
			MISRA_EXCEPTION_RULE_14_7()
	    	return one/x;
	    } else {
			MISRA_EXCEPTION_RULE_14_7()
	    	return x;
		}
	}
	if(hy==0x40000000L) {
		MISRA_EXCEPTION_RULE_14_7()
		return x*x; /* y is  2 */
	}
	if(hy==0x3f000000L) {	/* y is  0.5 */
	    if(hx>=0L) {	/* x >= +0 */
			MISRA_EXCEPTION_RULE_14_7()
	    	return __ieee754_sqrtf(x);
		}
	}

	ax   = fabsf(x);
    /* special value of x */
	if((ix==0x7f800000UL)||(ix==0UL)||(ix==0x3f800000UL)){
	    z = ax;			/*x is +-0,+-inf,+-1*/
	    if(hy<0L) {z = one/z;}	/* z = (1/|x|) */
	    if(hx<0L) {
			if(((ix-0x3f800000UL)|(uint32_t)yisint)==0UL) {
			    z = (z-z)/(z-z); /* (-1)**non-int is NaN */
			} else {
				if(yisint==1L) {
			    	z = -z;		/* (x<0)**odd = -(|x|**odd) */
			    }
	    	}
	    }
		MISRA_EXCEPTION_RULE_14_7()
	    return z;
	}

    MISRA_EXCEPTION_RULE_12_7()
	n = (int32_t)((uint32_t)hx>>31)-1L;

    /* (x<0)**(non-int) is NaN */
	if((n|yisint)==0L)
	{
		#if _EWL_C99
			MISRA_EXCEPTION_RULE_13_7()
			if ((uint_t)math_errhandling & (uint_t)MATH_ERRNO) {	
				MISRA_EXCEPTION_RULE_20_5()
				errno=EDOM;
			}
			MISRA_EXCEPTION_RULE_13_7()
			if ((uint_t)math_errhandling & (uint_t)MATH_ERREXCEPT) {
				MISRA_EXCEPTION_RULE_20_5()
				feraiseexcept((int_t)FE_INVALID);	
			}
		#else	
		#ifdef __STDC__
			MISRA_EXCEPTION_RULE_20_5()
			errno=EDOM;  /* mf-- added to conform to old ANSI standard */
        #endif
		#endif
		#ifdef __CWCC__
			MISRA_EXCEPTION_RULE_14_7()
			return NAN;
		#else
			MISRA_EXCEPTION_RULE_14_7()
			return (x-x)/(x-x);
		#endif
	}

	sn = one; /* s (sign of result -ve**odd) = -1 else = 1 */
	if((n|(yisint-1))==0L) {sn = -one;}/* (-ve)**(odd int) */

    /* |y| is huge */
	if(iy>0x4d000000UL) { /* if |y| > 2**27 */
	/* over/underflow if x is not close to one */
	    if(ix<0x3f7ffff8UL) {
			MISRA_EXCEPTION_RULE_14_7()
			return (hy<0L)? (sn*hugeval*hugeval):(sn*tiny*tiny);
		}
	    if(ix>0x3f800007UL) {
			MISRA_EXCEPTION_RULE_14_7()
			return (hy>0L)? (sn*hugeval*hugeval):(sn*tiny*tiny);
		}
	/* now |1-x| is tiny <= 2**-20, suffice to compute
	   log(x) by x-x^2/2+x^3/3-x^4/4 */
	    t = ax-1;		/* t has 20 trailing zeros */
	    w = (t*t)*(0.5F-(t*(0.333333333333F-(t*0.25F))));
	    u = ivln2_h*t;	/* ivln2_h has 16 sig. bits */
	    v = (t*ivln2_l)-(w*ivln2);
	    t1 = u+v;
	    GET_FLOAT_WORD(is,t1);
	    SET_FLOAT_WORD(t1,is&0xfffff000L);
	    t2 = v-(t1-u);
	} else {
	    f32_t s2,s_h,s_l,t_h,t_l;
	    n = 0;
	/* take care subnormal number */
	    if(ix<0x00800000UL) {
			ax *= two24; n -= 24L; GET_FLOAT_UWORD(ix,ax);
		}
	    n  += (int32_t)((ix)>>23)-0x7fL;
	    j  = (int32_t)(ix&0x007fffffUL);
	/* determine interval */
	    ix = j|0x3f800000UL;		/* normalize ix */
	    if(j<=0x1cc471L) {k=0L;}	/* |x|<sqrt(3/2) */
	    else if(j<0x5db3d7L) {k=1L;}	/* |x|<sqrt(3)   */
	    else {k=0;n+=1L;ix -= 0x00800000UL;}
	    SET_FLOAT_UWORD(ax,ix);

	/* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
	    u = ax-bp[k];		/* bp[0]=1.0, bp[1]=1.5 */
	    v = one/(ax+bp[k]);
	    s = u*v;
	    s_h = s;
	    GET_FLOAT_WORD(is,s_h);
	    SET_FLOAT_WORD(s_h,is&0xfffff000UL);
	/* t_h=ax+bp[k] High */
	    is = (int32_t)(((ix>>1)&0xfffff000UL)|0x20000000UL);
	    SET_FLOAT_WORD(t_h,is+0x00400000UL+(k<<21));
	    t_l = ax - (t_h-bp[k]);
	    s_l = v*((u-(s_h*t_h))-(s_h*t_l));
	/* compute log(ax) */
	    s2 = s*s;
	    r = (s2*s2)*(L1+(s2*(L2+(s2*(L3+(s2*(L4+(s2*(L5+(s2*L6))))))))));
	    r += s_l*(s_h+s);
	    s2  = s_h*s_h;
	    t_h = 3.0F+s2+r;
	    GET_FLOAT_WORD(is,t_h);
	    SET_FLOAT_WORD(t_h,is&0xfffff000UL);
	    t_l = r-((t_h-3.0F)-s2);
	/* u+v = s*(1+...) */
	    u = s_h*t_h;
	    v = (s_l*t_h)+(t_l*s);
	/* 2/(3log2)*(s+...) */
	    p_h = u+v;
	    GET_FLOAT_WORD(is,p_h);
	    SET_FLOAT_WORD(p_h,is&0xfffff000UL);
	    p_l = v-(p_h-u);
	    z_h = cp_h*p_h;		/* cp_h+cp_l = 2/(3*log2) */
	    z_l = (cp_l*p_h)+(p_l*cp)+dp_l[k];
	/* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
	    t = (f32_t)n;
	    t1 = (((z_h+z_l)+dp_h[k])+t);
	    GET_FLOAT_WORD(is,t1);
	    SET_FLOAT_WORD(t1,is&0xfffff000UL);
	    t2 = z_l-(((t1-t)-dp_h[k])-z_h);
	}

    /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
	GET_FLOAT_WORD(is,y);
	SET_FLOAT_WORD(y1,is&0xfffff000UL);
	p_l = ((y-y1)*t1)+(y*t2);
	p_h = y1*t1;
	z = p_l+p_h;
	GET_FLOAT_WORD(j,z);
	if (j>0x43000000L) {				/* if z > 128 */
		MISRA_EXCEPTION_RULE_14_7()
	    return sn*hugeval*hugeval;			/* overflow */
	} else if (j==0x43000000L) {			/* if z == 128 */
	    if((p_l+ovt)>(z-p_h)) {				/* overflow */
			MISRA_EXCEPTION_RULE_14_7()
			return sn*hugeval*hugeval;
		}
	} else if ((j&0x7fffffffUL)>0x43160000UL) {		/* z <= -150 */
		MISRA_EXCEPTION_RULE_14_7()
	    return sn*tiny*tiny;			/* underflow */
	} else if (j==(int32_t)0xc3160000L){			/* z == -150 */
	    if(p_l<=(z-p_h)) {
			MISRA_EXCEPTION_RULE_14_7()
	    	return sn*tiny*tiny;		/* underflow */
		}
	}
    /*
     * compute 2**(p_h+p_l)
     */
	i = j&0x7fffffffUL;
	MISRA_EXCEPTION_RULE_10_3()
	k = (i>>23)-0x7fUL;
	n = 0;
	MISRA_EXCEPTION_RULE_10_3()
	if(i>0x3f000000L) {		/* if |z| > 0.5, set n = [z+0.5] */
	    n = j+(int32_t)(0x00800000UL>>(k+1L));
	    k = ((n&0x7fffffffL)>>23)-0x7fL;	/* new k for n */
	    SET_FLOAT_WORD(t,n&~(0x007fffffL>>k));
	    n = ((n&0x007fffffL)|0x00800000L)>>(23L-k);
	    if(j<0) {n = -n;}
	    p_h -= t;
	}
	t = p_l+p_h;
	GET_FLOAT_WORD(is,t);
	SET_FLOAT_WORD(t,is&0xffff8000UL);
	u = t*lg2_h;
	v = ((p_l-(t-p_h))*lg2)+(t*lg2_l);
	z = u+v;
	w = v-(z-u);
	t  = z*z;
	t1 = z - (t*(P1+(t*(P2+(t*(P3+(t*(P4+(t*P5)))))))));
	r  = ((z*t1)/(t1-two))-(w+(z*w));
	z  = one-(r-z);
	GET_FLOAT_WORD(j,z);
	MISRA_EXCEPTION_RULE_12_7()
	j += (n<<23);
	MISRA_EXCEPTION_RULE_12_7()
	MISRA_QUIET_LINT_02()
	if((j>>23)<=0L) {
		z = scalbnf(z,n);	/* subnormal output */
	} else {
		SET_FLOAT_WORD(z,j);
	}
	return sn*z;
}
#endif /* _EWL_FLOATING_POINT  */
