C++ SET_FLOAT_WORD函数代码示例

float
__ieee754_log10f(float x)
{
	float f,hi,lo,y,z;
	int32_t i,k,hx;

	GET_FLOAT_WORD(hx,x);

        k=0;
        if (hx < 0x00800000) {                  /* x < 2**-126  */
            if ((hx&0x7fffffff)==0)
                return -two25/zero;             /* log(+-0)=-inf */
            if (hx<0) return (x-x)/zero;        /* log(-#) = NaN */
            k -= 25; x *= two25; /* subnormal number, scale up x */
	    GET_FLOAT_WORD(hx,x);
        }
	if (hx >= 0x7f800000) return x+x;
	k += (hx>>23)-127;
	hx &= 0x007fffff;
	i = (hx+(0x4afb0d))&0x800000;
	SET_FLOAT_WORD(x,hx|(i^0x3f800000));	/* normalize x or x/2 */
	k += (i>>23);
	y = (float)k;
	f = __kernel_logf(x);
	x = x - (float)1.0;
	GET_FLOAT_WORD(hx,x);
	SET_FLOAT_WORD(hi,hx&0xfffff000);
	lo = x - hi;
	z = y*log10_2lo + (x+f)*ivln10lo + (lo+f)*ivln10hi + hi*ivln10hi;
	return  z+y*log10_2hi;
}

float
cbrtf(float x)
{
	float r,s,t;
	int32_t hx;
	uint32_t sign;
	uint32_t high;

	GET_FLOAT_WORD(hx,x);
	sign=hx&0x80000000; 		/* sign= sign(x) */
	hx  ^=sign;
	if(hx>=0x7f800000) return(x+x); /* cbrt(NaN,INF) is itself */
	if(hx==0)
	    return(x);		/* cbrt(0) is itself */

	SET_FLOAT_WORD(x,hx);	/* x <- |x| */
    /* rough cbrt to 5 bits */
	if(hx<0x00800000) 		/* subnormal number */
	  {SET_FLOAT_WORD(t,0x4b800000); /* set t= 2**24 */
	   t*=x; GET_FLOAT_WORD(high,t); SET_FLOAT_WORD(t,high/3+B2);
	  }
	else
	  SET_FLOAT_WORD(t,hx/3+B1);


    /* new cbrt to 23 bits */
	r=t*t/x;
	s=C+r*t;
	t*=G+F/(s+E+D/s);

    /* retore the sign bit */
	GET_FLOAT_WORD(high,t);
	SET_FLOAT_WORD(t,high|sign);
	return(t);
}

DLLEXPORT float
cbrtf(float x)
{
	double r,T;
	float t;
	int32_t hx;
	u_int32_t sign;
	u_int32_t high;

	GET_FLOAT_WORD(hx,x);
	sign=hx&0x80000000; 		/* sign= sign(x) */
	hx  ^=sign;
	if(hx>=0x7f800000) return(x+x); /* cbrt(NaN,INF) is itself */

    /* rough cbrt to 5 bits */
	if(hx<0x00800000) { 		/* zero or subnormal? */
	    if(hx==0)
		return(x);		/* cbrt(+-0) is itself */
	    SET_FLOAT_WORD(t,0x4b800000); /* set t= 2**24 */
#if _TMS320C6X
            float __fmpy_by_0x1p24(float);
            x = __fmpy_by_0x1p24(x);
#endif
	    t*=x;

	    GET_FLOAT_WORD(high,t);
	    SET_FLOAT_WORD(t,sign|((high&0x7fffffff)/3+B2));
	} else
	    SET_FLOAT_WORD(t,sign|(hx/3+B1));

    /*
     * First step Newton iteration (solving t*t-x/t == 0) to 16 bits.  In
     * double precision so that its terms can be arranged for efficiency
     * without causing overflow or underflow.
     */
	T=t;
	r=T*T*T;
	T=T*((double)x+x+r)/(x+r+r);

    /*
     * Second step Newton iteration to 47 bits.  In double precision for
     * efficiency and accuracy.
     */
	r=T*T*T;
	T=T*((double)x+x+r)/(x+r+r);

#if _TMS320C6X
	if(hx<0x00800000) { 		/* zero or subnormal? */
            T /= 0x1p8;
        }
#endif

    /* rounding to 24 bits is perfect in round-to-nearest mode */
	return(T);
}

float
__ieee754_log2f(float x)
{
	float f,hfsq,hi,lo,r,y;
	int32_t i,k,hx;

	GET_FLOAT_WORD(hx,x);

	k=0;
	if (hx < 0x00800000) {			/* x < 2**-126  */
	    if ((hx&0x7fffffff)==0)
		return -two25/vzero;		/* log(+-0)=-inf */
	    if (hx<0) return (x-x)/zero;	/* log(-#) = NaN */
	    k -= 25; x *= two25; /* subnormal number, scale up x */
	    GET_FLOAT_WORD(hx,x);
	}
	if (hx >= 0x7f800000) return x+x;
	if (hx == 0x3f800000)
	    return zero;			/* log(1) = +0 */
	k += (hx>>23)-127;
	hx &= 0x007fffff;
	i = (hx+(0x4afb0d))&0x800000;
	SET_FLOAT_WORD(x,hx|(i^0x3f800000));	/* normalize x or x/2 */
	k += (i>>23);
	y = (float)k;
	f = x - (float)1.0;
	hfsq = (float)0.5*f*f;
	r = k_log1pf(f);

	/*
	 * We no longer need to avoid falling into the multi-precision
	 * calculations due to compiler bugs breaking Dekker's theorem.
	 * Keep avoiding this as an optimization.  See e_log2.c for more
	 * details (some details are here only because the optimization
	 * is not yet available in double precision).
	 *
	 * Another compiler bug turned up.  With gcc on i386,
	 * (ivln2lo + ivln2hi) would be evaluated in float precision
	 * despite runtime evaluations using double precision.  So we
	 * must cast one of its terms to float_t.  This makes the whole
	 * expression have type float_t, so return is forced to waste
	 * time clobbering its extra precision.
	 */
	if (sizeof(float_t) > sizeof(float))
		return (r - hfsq + f) * ((float_t)ivln2lo + ivln2hi) + y;

	hi = f - hfsq;
	GET_FLOAT_WORD(hx,hi);
	SET_FLOAT_WORD(hi,hx&0xfffff000);
	lo = (f - hi) - hfsq + r;
	return (lo+hi)*ivln2lo + lo*ivln2hi + hi*ivln2hi + y;
}

float nextafterf(float x, float y)
{
	int32_t hx, hy, ix, iy;

	GET_FLOAT_WORD(hx, x);
	GET_FLOAT_WORD(hy, y);
	ix = hx & 0x7fffffff;		/* |x| */
	iy = hy & 0x7fffffff;		/* |y| */

	/* x is nan or y is nan? */
	if ((ix > 0x7f800000) || (iy > 0x7f800000))
		return x + y;

	if (x == y)
		return y;

	if (ix == 0) { /* x == 0? */
/* glibc 2.4 does not seem to set underflow? */
/*		float u; */
		/* return +-minsubnormal */
		SET_FLOAT_WORD(x, (hy & 0x80000000) | 1);
/*		u = x * x;     raise underflow flag */
/*		math_force_eval(u); */
		return x;
	}

	if (hx >= 0) { /* x > 0 */
		if (hx > hy) { /* x > y: x -= ulp */
			hx -= 1;
		} else { /* x < y: x += ulp */
			hx += 1;
		}
	} else { /* x < 0 */
		if (hy >= 0 || hx > hy) { /* x < y: x -= ulp */
			hx -= 1;
		} else { /* x > y: x += ulp */
			hx += 1;
		}
	}
	hy = hx & 0x7f800000;
	if (hy >= 0x7f800000) {
		x = x + x; /* overflow */
		return x; /* overflow */
	}
	if (hy < 0x00800000) {
		float u = x * x; /* underflow */
		math_force_eval(u); /* raise underflow flag */
	}
	SET_FLOAT_WORD(x, hx);
	return x;
}

float npy_nextafterf(float x, float y)
{
    volatile float t;
    npy_int32 hx, hy, ix, iy;

    GET_FLOAT_WORD(hx, x);
    GET_FLOAT_WORD(hy, y);
    ix = hx & 0x7fffffff;       /* |x| */
    iy = hy & 0x7fffffff;       /* |y| */

    if ((ix > 0x7f800000) ||    /* x is nan */
        (iy > 0x7f800000))      /* y is nan */
        return x + y;
    if (x == y)
        return y;               /* x=y, return y */
    if (ix == 0) {              /* x == 0 */
        SET_FLOAT_WORD(x, (hy & 0x80000000) | 1); /* return +-minsubnormal */
        t = x * x;
        if (t == x)
            return t;
        else
            return x;           /* raise underflow flag */
    }
    if (hx >= 0) {              /* x > 0 */
        if (hx > hy) {          /* x > y, x -= ulp */
            hx -= 1;
        } else {                /* x < y, x += ulp */
            hx += 1;
        }
    } else {                    /* x < 0 */
        if (hy >= 0 || hx > hy) {       /* x < y, x -= ulp */
            hx -= 1;
        } else {                /* x > y, x += ulp */
            hx += 1;
        }
    }
    hy = hx & 0x7f800000;
    if (hy >= 0x7f800000)
        return x + x;           /* overflow  */
    if (hy < 0x00800000) {      /* underflow */
        t = x * x;
        if (t != x) {           /* raise underflow flag */
            SET_FLOAT_WORD(y, hx);
            return y;
        }
    }
    SET_FLOAT_WORD(x, hx);
    return x;
}

float __nexttowardf(float x, long double y)
{
	int32_t hx,ix,iy;
	uint32_t hy,ly,esy;

	GET_FLOAT_WORD(hx,x);
	GET_LDOUBLE_WORDS(esy,hy,ly,y);
	ix = hx&0x7fffffff;		/* |x| */
	iy = esy&0x7fff;		/* |y| */

	if((ix>0x7f800000) ||			/* x is nan */
	   (iy>=0x7fff&&((hy|ly)!=0)))		/* y is nan */
	   return x+y;
	if((long double) x==y) return y;	/* x=y, return y */
	if(ix==0) {				/* x == 0 */
	    float u;
	    SET_FLOAT_WORD(x,((esy&0x8000)<<16)|1);/* return +-minsub*/
	    u = math_opt_barrier (x);
	    u = u * u;
	    math_force_eval (u);		/* raise underflow flag */
	    return x;
	}
	if(hx>=0) {				/* x > 0 */
	    if(x > y) {				/* x -= ulp */
		hx -= 1;
	    } else {				/* x < y, x += ulp */
		hx += 1;
	    }
	} else {				/* x < 0 */
	    if(x < y) {				/* x -= ulp */
		hx -= 1;
	    } else {				/* x > y, x += ulp */
		hx += 1;
	    }
	}
	hy = hx&0x7f800000;
	if(hy>=0x7f800000) {
	  float u = x+x;			/* overflow  */
	  math_force_eval (u);
	  __set_errno (ERANGE);
	}
	if(hy<0x00800000) {
	    float u = x*x;			/* underflow */
	    math_force_eval (u);		/* raise underflow flag */
	    __set_errno (ERANGE);
	}
	SET_FLOAT_WORD(x,hx);
	return x;
}

float nextafterf(float x, float y)
{
	int32_t hx,hy,ix,iy;

	GET_FLOAT_WORD(hx,x);
	GET_FLOAT_WORD(hy,y);
	ix = hx&0x7fffffff;		/* |x| */
	iy = hy&0x7fffffff;		/* |y| */

	if((ix>0x7f800000) ||   /* x is nan */
	   (iy>0x7f800000))     /* y is nan */
    return x+y;
	if(x==y) return y;		/* x=y, return y */
	if(ix==0) {				/* x == 0 */
    //float u;
    SET_FLOAT_WORD(x,(hy&0x80000000)|1);/* return +-minsubnormal */
    //u = math_opt_barrier (x);
    //u = u*u;
    //math_force_eval (u);		/* raise underflow flag */
    return x;
	}
	if(hx>=0) {				/* x > 0 */
    if(hx>hy) {				/* x > y, x -= ulp */
      hx -= 1;
    } else {				/* x < y, x += ulp */
      hx += 1;
    }
	} else {				/* x < 0 */
    if(hy>=0||hx>hy){			/* x < y, x -= ulp */
      hx -= 1;
    } else {				/* x > y, x += ulp */
      hx += 1;
    }
	}
	hy = hx&0x7f800000;
	if(hy>=0x7f800000) {
	  x = x+x;	/* overflow  */
	  //if (FLT_EVAL_METHOD != 0)
	  //  asm ("" : "+m"(x));
	  return x;	/* overflow  */
	}
	//if(hy<0x00800000) {
  //  float u = x*x;			/* underflow */
  //  math_force_eval (u);		/* raise underflow flag */
	//}
	SET_FLOAT_WORD(x,hx);
	return x;
}

float
__kernel_cosf(float x, float y)
{
	float a,hz,z,r,qx;
	int32_t ix;
	GET_FLOAT_WORD(ix,x);
	ix &= 0x7fffffff;			/* ix = |x|'s high word*/
	if(ix<0x32000000) {			/* if x < 2**27 */
	    if(((int)x)==0) return one;		/* generate inexact */
	}
	z  = x*x;
	r  = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))));
	if(ix < 0x3e99999a) 			/* if |x| < 0.3 */ 
	    return one - ((float)0.5*z - (z*r - x*y));
	else {
	    if(ix > 0x3f480000) {		/* x > 0.78125 */
		qx = (float)0.28125;
	    } else {
	        SET_FLOAT_WORD(qx,ix-0x01000000);	/* x/4 */
	    }
	    hz = (float)0.5*z-qx;
	    a  = one-qx;
	    return a - (hz - (z*r-x*y));
	}
}

	float nanf(const char *unused)
{
	float x;

	SET_FLOAT_WORD(x,0x7fc00000);
	return x;
}

float log10f(float x)
{
	float y, z;
	s32_t i, k, hx;

	GET_FLOAT_WORD(hx,x);

	k = 0;
	if (hx < 0x00800000)
	{
		if ((hx & 0x7fffffff) == 0)
			return -two25 / zero;
		if (hx < 0)
			return (x - x) / zero;
		k -= 25;
		x *= two25;
		GET_FLOAT_WORD(hx,x);
	}
	if (hx >= 0x7f800000)
		return x + x;
	k += (hx >> 23) - 127;
	i = ((u32_t) k & 0x80000000) >> 31;
	hx = (hx & 0x007fffff) | ((0x7f - i) << 23);
	y = (float) (k + i);
	SET_FLOAT_WORD(x,hx);
	z = y * log10_2lo + ivln10 * logf(x);
	return z + y * log10_2hi;
}

float frexpf(float x, int *eptr)
{
	s32_t hx, ix;

	GET_FLOAT_WORD(hx,x);
	ix = 0x7fffffff & hx;
	*eptr = 0;

	if (ix >= 0x7f800000 || (ix == 0))
		return x;

	if (ix < 0x00800000)
	{
		x *= two25;
		GET_FLOAT_WORD(hx,x);
		ix = hx & 0x7fffffff;
		*eptr = -25;
	}

	*eptr += (ix >> 23) - 126;
	hx = (hx & 0x807fffff) | 0x3f000000;
	SET_FLOAT_WORD(x,hx);

	return x;
}

float
log10f(float x)
{
	float y,z;
	int32_t i,k,hx;

	GET_FLOAT_WORD(hx,x);

        k=0;
        if (hx < 0x00800000) {                  /* x < 2**-126  */
            if ((hx&0x7fffffff)==0)
                return -two25/zero;             /* log(+-0)=-inf */
            if (hx<0) return (x-x)/zero;        /* log(-#) = NaN */
            k -= 25; x *= two25; /* subnormal number, scale up x */
	    GET_FLOAT_WORD(hx,x);
        }
	if (hx >= 0x7f800000) return x+x;
	k += (hx>>23)-127;
	i  = ((u_int32_t)k&0x80000000)>>31;
        hx = (hx&0x007fffff)|((0x7f-i)<<23);
        y  = (float)(k+i);
	SET_FLOAT_WORD(x,hx);
	z  = y*log10_2lo + ivln10*logf(x);
	return  z+y*log10_2hi;
}

float complex ctanhf(float complex z)
{
	float x, y;
	float t, beta, s, rho, denom;
	uint32_t hx, ix;

	x = crealf(z);
	y = cimagf(z);

	GET_FLOAT_WORD(hx, x);
	ix = hx & 0x7fffffff;

	if (ix >= 0x7f800000) {
		if (ix & 0x7fffff)
			return CMPLXF(x, (y == 0 ? y : x * y));
		SET_FLOAT_WORD(x, hx - 0x40000000);
		return CMPLXF(x, copysignf(0, isinf(y) ? y : sinf(y) * cosf(y)));
	}

	if (!isfinite(y))
		return CMPLXF(y - y, y - y);

	if (ix >= 0x41300000) { /* x >= 11 */
		float exp_mx = expf(-fabsf(x));
		return CMPLXF(copysignf(1, x), 4 * sinf(y) * cosf(y) * exp_mx * exp_mx);
	}

	t = tanf(y);
	beta = 1.0 + t * t;
	s = sinhf(x);
	rho = sqrtf(1 + s * s);
	denom = 1 + beta * s * s;
	return CMPLXF((beta * rho * s) / denom, t / denom);
}

float scalbnf(float x, int n)
{
	float scale;

	if (n > 127) {
		x *= 0x1p127f;
		n -= 127;
		if (n > 127) {
			x *= 0x1p127f;
			n -= 127;
			if (n > 127) {
				x = (float)(x * 0x1p127f);
				return x;
			}
		}
	} else if (n < -126) {
		x *= 0x1p-126f;
		n += 126;
		if (n < -126) {
			x *= 0x1p-126f;
			n += 126;
			if (n < -126) {
				x = (float)(x * 0x1p-126f);
				return x;
			}
		}
	}
	SET_FLOAT_WORD(scale, (uint32_t)(0x7f+n)<<23);
	x = (float)(x * scale);
	return x;
}