libm.h File Reference

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

Include dependency graph for libm.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Macros
#define	_STATIC_ASSERT(expr)   extern char (*__static_assert__(void)) [(expr) ? 1 : -1]
#define	isnan   _isnan
#define	isnanf   _isnanf
#define	signbit(x)
#define	FORCE_EVAL(x)
#define	hidden
#define	WANT_ROUNDING   1
#define	WANT_SNAN   0
#define	issignalingf_inline(x)   0
#define	issignaling_inline(x)   0
#define	TOINT_INTRINSICS   0
#define	predict_true(x)   (x)
#define	predict_false(x)   (x)
#define	EXTRACT_WORDS(hi, lo, d)
#define	GET_HIGH_WORD(hi, d)
#define	GET_LOW_WORD(lo, d)
#define	INSERT_WORDS(d, hi, lo)
#define	SET_HIGH_WORD(d, hi)    INSERT_WORDS(d, hi, (uint32_t)asuint64(d))
#define	SET_LOW_WORD(d, lo)    INSERT_WORDS(d, asuint64(d)>>32, lo)
#define	GET_FLOAT_WORD(w, d)
#define	SET_FLOAT_WORD(d, w)
#define	FP_ILOGB0   (-0x7fffffff - 1)
#define	FP_ILOGBNAN   0x7fffffff
#define	isfinite(x)   (_finite((double)(x)) != 0)

Functions
	_STATIC_ASSERT (sizeof(float)==4)
	_STATIC_ASSERT (sizeof(double)==8)
	_STATIC_ASSERT (sizeof(long double)==8)
_Check_return_ int __cdecl	_isnanf (_In_ float _X)
static __inline uint32_t	__FLOAT_BITS (float __f)
static __inline uint64_t	__DOUBLE_BITS (double __f)
static void	fp_force_evalf (float x)
static void	fp_force_eval (double x)
static float	eval_as_float (float x)
static double	eval_as_double (double x)
static float	fp_barrierf (float x)
static double	fp_barrier (double x)
static uint32_t	asuint (float f)
static float	asfloat (uint32_t i)
static uint64_t	asuint64 (double f)
static double	asdouble (uint64_t i)
_Check_return_ int __cdecl	_finite (_In_ double _X)
int	ilogb (double)
double	logb (double)
double	scalbn (double, int)
hidden int	__rem_pio2_large (double *, double *, int, int, int)
hidden int	__rem_pio2 (double, double *)
hidden double	__tan (double, double, int)
hidden double	__math_xflow (uint32_t, double)
hidden double	__math_uflow (uint32_t)
hidden double	__math_oflow (uint32_t)
hidden double	__math_divzero (uint32_t)
hidden double	__math_invalid (double)

Macro Definition Documentation

_STATIC_ASSERT

#define _STATIC_ASSERT

(

expr

)

extern char (*__static_assert__(void)) [(expr) ? 1 : -1]

Definition at line 19 of file libm.h.

EXTRACT_WORDS

#define EXTRACT_WORDS	(		hi,
			lo,
			d
	)

Value:

do {                                              \
  uint64_t __u = asuint64(d);                     \
  (hi) = __u >> 32;                               \
  (lo) = (uint32_t)__u;                           \
} while (0)

Definition at line 149 of file libm.h.

FORCE_EVAL

#define FORCE_EVAL

(

x

)

Value:

    do {                    \
    if (sizeof(x) == sizeof(float)) {         \
        fp_force_evalf(x);                    \
    } else {                                  \
        fp_force_eval(x);                     \
    }                                         \
} while(0)

Definition at line 64 of file libm.h.

FP_ILOGB0

#define FP_ILOGB0   (-0x7fffffff - 1)

Definition at line 188 of file libm.h.

FP_ILOGBNAN

#define FP_ILOGBNAN   0x7fffffff

Definition at line 191 of file libm.h.

GET_FLOAT_WORD

#define GET_FLOAT_WORD	(		w,
			d
	)

Value:

do {                                              \
  (w) = asuint(d);                                \
} while (0)

Definition at line 177 of file libm.h.

GET_HIGH_WORD

#define GET_HIGH_WORD	(		hi,
			d
	)

Value:

do {                                              \
  (hi) = asuint64(d) >> 32;                       \
} while (0)

Definition at line 156 of file libm.h.

GET_LOW_WORD

#define GET_LOW_WORD	(		lo,
			d
	)

Value:

do {                                              \
  (lo) = (uint32_t)asuint64(d);                   \
} while (0)

Definition at line 161 of file libm.h.

hidden

uint32_t hidden

Definition at line 73 of file libm.h.

Referenced by ewmh_get_window_state(), fsctl_set_xattr(), id(), TOOLBAR_SetRows(), and ui_seamless_begin().

INSERT_WORDS

#define INSERT_WORDS	(		d,
			hi,
			lo
	)

Value:

do {                                              \
  (d) = asdouble(((uint64_t)(hi)<<32) | (uint32_t)(lo)); \
} while (0)

Definition at line 166 of file libm.h.

isfinite

#define isfinite

(

x

)

(_finite((double)(x)) != 0)

Definition at line 197 of file libm.h.

isnan

#define isnan   _isnan

Definition at line 28 of file libm.h.

isnanf

#define isnanf   _isnanf

Definition at line 29 of file libm.h.

issignaling_inline

#define issignaling_inline

(

x

)

0

Definition at line 81 of file libm.h.

issignalingf_inline

#define issignalingf_inline

(

x

)

0

Definition at line 80 of file libm.h.

predict_false

#define predict_false

(

x

)

(x)

Definition at line 93 of file libm.h.

predict_true

#define predict_true

(

x

)

(x)

Definition at line 92 of file libm.h.

SET_FLOAT_WORD

#define SET_FLOAT_WORD	(		d,
			w
	)

Value:

do {                                              \
  (d) = asfloat(w);                               \
} while (0)

Definition at line 182 of file libm.h.

SET_HIGH_WORD

#define SET_HIGH_WORD	(		d,
			hi
	)		INSERT_WORDS(d, hi, (uint32_t)asuint64(d))

Definition at line 171 of file libm.h.

SET_LOW_WORD

#define SET_LOW_WORD	(		d,
			lo
	)		INSERT_WORDS(d, asuint64(d)>>32, lo)

Definition at line 174 of file libm.h.

signbit

#define signbit

(

x

)

Value:

    ((sizeof(x) == sizeof(float)) ? (int)(__FLOAT_BITS(x)>>31) : \
                                    (int)(__DOUBLE_BITS(x)>>63))

Definition at line 46 of file libm.h.

TOINT_INTRINSICS

#define TOINT_INTRINSICS   0

Definition at line 84 of file libm.h.

WANT_ROUNDING

#define WANT_ROUNDING   1

Definition at line 76 of file libm.h.

WANT_SNAN

#define WANT_SNAN   0

Definition at line 78 of file libm.h.

Function Documentation

__DOUBLE_BITS()

static __inline uint64_t __DOUBLE_BITS ( double  __f )

static

Definition at line 39 of file libm.h.

{
    union {double __f; uint64_t __i;} __u;
    __u.__f = __f;
    return __u.__i;
}

__FLOAT_BITS()

static __inline uint32_t __FLOAT_BITS ( float  __f )

static

Definition at line 32 of file libm.h.

{
    union {float __f; uint32_t __i;} __u;
    __u.__f = __f;
    return __u.__i;
}

__math_divzero()

hidden double __math_divzero

(

uint32_t

sign

)

Definition at line 13 of file __math_divzero.c.

{
    return fp_barrier(sign ? -1.0 : 1.0) / 0.0;
}

Referenced by log().

__math_invalid()

hidden double __math_invalid

(

double

x

)

Definition at line 13 of file __math_invalid.c.

{
    return (x - x) / (x - x);
}

Referenced by log(), pow(), and sqrt().

__math_oflow()

hidden double __math_oflow

(

uint32_t

sign

)

Definition at line 13 of file __math_oflow.c.

{
    return __math_xflow(sign, 0x1p769);
}

Referenced by exp(), exp_inline(), and pow().

__math_uflow()

hidden double __math_uflow

(

uint32_t

sign

)

Definition at line 13 of file __math_uflow.c.

{
    return __math_xflow(sign, 0x1p-767);
}

Referenced by exp(), exp_inline(), and pow().

__math_xflow()

hidden double __math_xflow	(	uint32_t	sign,
		double	y
	)

Definition at line 13 of file __math_xflow.c.

{
    return eval_as_double(fp_barrier(sign ? -y : y) * y);
}

Referenced by __math_oflow(), and __math_uflow().

__rem_pio2()

hidden int __rem_pio2	(	double	x,
		double *	y
	)

Definition at line 59 of file __rem_pio2.c.

{
    union {double f; uint64_t i;} u = {x};
    double_t z,w,t,r,fn;
    double tx[3],ty[2];
    uint32_t ix;
    int sign, n, ex, ey, i;
 
    sign = u.i>>63;
    ix = u.i>>32 & 0x7fffffff;
    if (ix <= 0x400f6a7a) {  /* |x| ~<= 5pi/4 */
        if ((ix & 0xfffff) == 0x921fb)  /* |x| ~= pi/2 or 2pi/2 */
            goto medium;  /* cancellation -- use medium case */
        if (ix <= 0x4002d97c) {  /* |x| ~<= 3pi/4 */
            if (!sign) {
                z = x - pio2_1;  /* one round good to 85 bits */
                y[0] = z - pio2_1t;
                y[1] = (z-y[0]) - pio2_1t;
                return 1;
            } else {
                z = x + pio2_1;
                y[0] = z + pio2_1t;
                y[1] = (z-y[0]) + pio2_1t;
                return -1;
            }
        } else {
            if (!sign) {
                z = x - 2*pio2_1;
                y[0] = z - 2*pio2_1t;
                y[1] = (z-y[0]) - 2*pio2_1t;
                return 2;
            } else {
                z = x + 2*pio2_1;
                y[0] = z + 2*pio2_1t;
                y[1] = (z-y[0]) + 2*pio2_1t;
                return -2;
            }
        }
    }
    if (ix <= 0x401c463b) {  /* |x| ~<= 9pi/4 */
        if (ix <= 0x4015fdbc) {  /* |x| ~<= 7pi/4 */
            if (ix == 0x4012d97c)  /* |x| ~= 3pi/2 */
                goto medium;
            if (!sign) {
                z = x - 3*pio2_1;
                y[0] = z - 3*pio2_1t;
                y[1] = (z-y[0]) - 3*pio2_1t;
                return 3;
            } else {
                z = x + 3*pio2_1;
                y[0] = z + 3*pio2_1t;
                y[1] = (z-y[0]) + 3*pio2_1t;
                return -3;
            }
        } else {
            if (ix == 0x401921fb)  /* |x| ~= 4pi/2 */
                goto medium;
            if (!sign) {
                z = x - 4*pio2_1;
                y[0] = z - 4*pio2_1t;
                y[1] = (z-y[0]) - 4*pio2_1t;
                return 4;
            } else {
                z = x + 4*pio2_1;
                y[0] = z + 4*pio2_1t;
                y[1] = (z-y[0]) + 4*pio2_1t;
                return -4;
            }
        }
    }
    if (ix < 0x413921fb) {  /* |x| ~< 2^20*(pi/2), medium size */
medium:
        /* rint(x/(pi/2)) */
        fn = (double_t)x*invpio2 + toint - toint;
        n = (int32_t)fn;
        r = x - fn*pio2_1;
        w = fn*pio2_1t;  /* 1st round, good to 85 bits */
        /* Matters with directed rounding. */
        if (predict_false(r - w < -pio4)) {
            n--;
            fn--;
            r = x - fn*pio2_1;
            w = fn*pio2_1t;
        } else if (predict_false(r - w > pio4)) {
            n++;
            fn++;
            r = x - fn*pio2_1;
            w = fn*pio2_1t;
        }
        y[0] = r - w;
        u.f = y[0];
        ey = u.i>>52 & 0x7ff;
        ex = ix>>20;
        if (ex - ey > 16) { /* 2nd round, good to 118 bits */
            t = r;
            w = fn*pio2_2;
            r = t - w;
            w = fn*pio2_2t - ((t-r)-w);
            y[0] = r - w;
            u.f = y[0];
            ey = u.i>>52 & 0x7ff;
            if (ex - ey > 49) {  /* 3rd round, good to 151 bits, covers all cases */
                t = r;
                w = fn*pio2_3;
                r = t - w;
                w = fn*pio2_3t - ((t-r)-w);
                y[0] = r - w;
            }
        }
        y[1] = (r - y[0]) - w;
        return n;
    }
    /*
     * all other (large) arguments
     */
    if (ix >= 0x7ff00000) {  /* x is inf or NaN */
        y[0] = y[1] = x - x;
        return 0;
    }
    /* set z = scalbn(|x|,-ilogb(x)+23) */
    u.f = x;
    u.i &= (uint64_t)-1>>12;
    u.i |= (uint64_t)(0x3ff + 23)<<52;
    z = u.f;
    for (i=0; i < 2; i++) {
        tx[i] = (double)(int32_t)z;
        z     = (z-tx[i])*0x1p24;
    }
    tx[i] = z;
    /* skip zero terms, first term is non-zero */
    while (tx[i] == 0.0)
        i--;
    n = __rem_pio2_large(tx,ty,(int)(ix>>20)-(0x3ff+23),i+1,1);
    if (sign) {
        y[0] = -ty[0];
        y[1] = -ty[1];
        return -n;
    }
    y[0] = ty[0];
    y[1] = ty[1];
    return n;
}

Referenced by tan().

__rem_pio2_large()

hidden int __rem_pio2_large	(	double *	x,
		double *	y,
		int	e0,
		int	nx,
		int	prec
	)

Definition at line 283 of file __rem_pio2_large.c.

{
    int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
    double z,fw,f[20],fq[20],q[20];
 
    /* initialize jk*/
    jk = init_jk[prec];
    jp = jk;
 
    /* determine jx,jv,q0, note that 3>q0 */
    jx = nx-1;
    jv = (e0-3)/24;  if(jv<0) jv=0;
    q0 = e0-24*(jv+1);
 
    /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
    j = jv-jx; m = jx+jk;
    for (i=0; i<=m; i++,j++)
        f[i] = j<0 ? 0.0 : (double)ipio2[j];
 
    /* compute q[0],q[1],...q[jk] */
    for (i=0; i<=jk; i++) {
        for (j=0,fw=0.0; j<=jx; j++)
            fw += x[j]*f[jx+i-j];
        q[i] = fw;
    }
 
    jz = jk;
recompute:
    /* distill q[] into iq[] reversingly */
    for (i=0,j=jz,z=q[jz]; j>0; i++,j--) {
        fw    = (double)(int32_t)(0x1p-24*z);
        iq[i] = (int32_t)(z - 0x1p24*fw);
        z     = q[j-1]+fw;
    }
 
    /* compute n */
    z  = scalbn(z,q0);       /* actual value of z */
    z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */
    n  = (int32_t)z;
    z -= (double)n;
    ih = 0;
    if (q0 > 0) {  /* need iq[jz-1] to determine n */
        i  = iq[jz-1]>>(24-q0); n += i;
        iq[jz-1] -= i<<(24-q0);
        ih = iq[jz-1]>>(23-q0);
    }
    else if (q0 == 0) ih = iq[jz-1]>>23;
    else if (z >= 0.5) ih = 2;
 
    if (ih > 0) {  /* q > 0.5 */
        n += 1; carry = 0;
        for (i=0; i<jz; i++) {  /* compute 1-q */
            j = iq[i];
            if (carry == 0) {
                if (j != 0) {
                    carry = 1;
                    iq[i] = 0x1000000 - j;
                }
            } else
                iq[i] = 0xffffff - j;
        }
        if (q0 > 0) {  /* rare case: chance is 1 in 12 */
            switch(q0) {
            case 1:
                iq[jz-1] &= 0x7fffff; break;
            case 2:
                iq[jz-1] &= 0x3fffff; break;
            }
        }
        if (ih == 2) {
            z = 1.0 - z;
            if (carry != 0)
                z -= scalbn(1.0,q0);
        }
    }
 
    /* check if recomputation is needed */
    if (z == 0.0) {
        j = 0;
        for (i=jz-1; i>=jk; i--) j |= iq[i];
        if (j == 0) {  /* need recomputation */
            for (k=1; iq[jk-k]==0; k++);  /* k = no. of terms needed */
 
            for (i=jz+1; i<=jz+k; i++) {  /* add q[jz+1] to q[jz+k] */
                f[jx+i] = (double)ipio2[jv+i];
                for (j=0,fw=0.0; j<=jx; j++)
                    fw += x[j]*f[jx+i-j];
                q[i] = fw;
            }
            jz += k;
            goto recompute;
        }
    }
 
    /* chop off zero terms */
    if (z == 0.0) {
        jz -= 1;
        q0 -= 24;
        while (iq[jz] == 0) {
            jz--;
            q0 -= 24;
        }
    } else { /* break z into 24-bit if necessary */
        z = scalbn(z,-q0);
        if (z >= 0x1p24) {
            fw = (double)(int32_t)(0x1p-24*z);
            iq[jz] = (int32_t)(z - 0x1p24*fw);
            jz += 1;
            q0 += 24;
            iq[jz] = (int32_t)fw;
        } else
            iq[jz] = (int32_t)z;
    }
 
    /* convert integer "bit" chunk to floating-point value */
    fw = scalbn(1.0,q0);
    for (i=jz; i>=0; i--) {
        q[i] = fw*(double)iq[i];
        fw *= 0x1p-24;
    }
 
    /* compute PIo2[0,...,jp]*q[jz,...,0] */
    for(i=jz; i>=0; i--) {
        for (fw=0.0,k=0; k<=jp && k<=jz-i; k++)
            fw += PIo2[k]*q[i+k];
        fq[jz-i] = fw;
    }
 
    /* compress fq[] into y[] */
    switch(prec) {
    case 0:
        fw = 0.0;
        for (i=jz; i>=0; i--)
            fw += fq[i];
        y[0] = ih==0 ? fw : -fw;
        break;
    case 1:
    case 2:
        fw = 0.0;
        for (i=jz; i>=0; i--)
            fw += fq[i];
        // TODO: drop excess precision here once double_t is used
        fw = (double)fw;
        y[0] = ih==0 ? fw : -fw;
        fw = fq[0]-fw;
        for (i=1; i<=jz; i++)
            fw += fq[i];
        y[1] = ih==0 ? fw : -fw;
        break;
    case 3:  /* painful */
        for (i=jz; i>0; i--) {
            fw      = fq[i-1]+fq[i];
            fq[i]  += fq[i-1]-fw;
            fq[i-1] = fw;
        }
        for (i=jz; i>1; i--) {
            fw      = fq[i-1]+fq[i];
            fq[i]  += fq[i-1]-fw;
            fq[i-1] = fw;
        }
        for (fw=0.0,i=jz; i>=2; i--)
            fw += fq[i];
        if (ih==0) {
            y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
        } else {
            y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
        }
    }
    return n&7;
}

Referenced by __rem_pio2().

__tan()

hidden double __tan	(	double	x,
		double	y,
		int	odd
	)

Definition at line 76 of file __tan.c.

{
    double_t z, r, v, w, s, a;
    double w0, a0;
    uint32_t hx;
    int big, sign;
 
    GET_HIGH_WORD(hx,x);
    big = (hx&0x7fffffff) >= 0x3FE59428; /* |x| >= 0.6744 */
    if (big) {
        sign = hx>>31;
        if (sign) {
            x = -x;
            y = -y;
        }
        x = (pio4 - x) + (pio4lo - y);
        y = 0.0;
    }
    z = x * x;
    w = z * z;
    /*
     * Break x^5*(T[1]+x^2*T[2]+...) into
     * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
     * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
     */
    r = T[1] + w*(T[3] + w*(T[5] + w*(T[7] + w*(T[9] + w*T[11]))));
    v = z*(T[2] + w*(T[4] + w*(T[6] + w*(T[8] + w*(T[10] + w*T[12])))));
    s = z * x;
    r = y + z*(s*(r + v) + y) + s*T[0];
    w = x + r;
    if (big) {
        s = 1 - 2*odd;
        v = s - 2.0 * (x + (r - w*w/(w + s)));
        return sign ? -v : v;
    }
    if (!odd)
        return w;
    /* -1.0/(x+r) has up to 2ulp error, so compute it accurately */
    w0 = w;
    SET_LOW_WORD(w0, 0);
    v = r - (w0 - x);       /* w0+v = r+x */
    a0 = a = -1.0 / w;
    SET_LOW_WORD(a0, 0);
    return a0 + a*(1.0 + a0*w0 + a0*v);
}

Referenced by tan().

_finite()

_Check_return_ int __cdecl _finite

(

_In_ double

_X

)

Definition at line 14 of file _finite.c.

{
    union { double f; uint64_t ui64; } u = { _X };
    uint64_t exp = u.ui64 & 0x7FF0000000000000;
    return (exp != 0x7FF0000000000000);
}

_isnanf()

_Check_return_ int __cdecl _isnanf

(

_In_ float

_X

)

Definition at line 13 of file _isnanf.c.

{
    union { float f; uint32_t ui32; } u = { _X };
    return (u.ui32 & ~0x80000000) > 0x7F800000;
}

_STATIC_ASSERT() [1/3]

_STATIC_ASSERT

(

sizeof(double)

= =8

)

_STATIC_ASSERT() [2/3]

_STATIC_ASSERT

(

sizeof(float)

= =4

)

_STATIC_ASSERT() [3/3]

_STATIC_ASSERT

(

sizeof(long double)

= =8

)

asdouble()

static double asdouble ( uint64_t  i )

inlinestatic

Definition at line 142 of file libm.h.

{
    union {uint64_t _i; double _f;} u;
    u._i = i;
    return u._f;
}

Referenced by exp(), exp_inline(), log(), log_inline(), pow(), specialcase(), and sqrt().

asfloat()

static float asfloat ( uint32_t  i )

inlinestatic

Definition at line 128 of file libm.h.

{
    union {uint32_t _i; float _f;} u;
    u._i = i;
    return u._f;
}

asuint()

static uint32_t asuint ( float  f )

inlinestatic

Definition at line 121 of file libm.h.

{
    union {float _f; uint32_t _i;} u;
    u._f = f;
    return u._i;
}

Referenced by _logbf().

asuint64()

static uint64_t asuint64 ( double  f )

inlinestatic

Definition at line 135 of file libm.h.

{
    union {double _f; uint64_t _i;} u;
    u._f = f;
    return u._i;
}

Referenced by _logb(), exp(), exp_inline(), log(), pow(), sqrt(), top12(), top16(), and zeroinfnan().

eval_as_double()

static double eval_as_double ( double  x )

inlinestatic

Definition at line 102 of file libm.h.

{
    double y = x;
    return y;
}

Referenced by __math_xflow(), exp(), exp_inline(), log(), specialcase(), and sqrt().

eval_as_float()

static float eval_as_float ( float  x )

inlinestatic

Definition at line 96 of file libm.h.

{
    float y = x;
    return y;
}

fp_barrier()

static double fp_barrier ( double  x )

inlinestatic

Definition at line 114 of file libm.h.

{
    volatile double y = x;
    return y;
}

Referenced by __math_divzero(), __math_xflow(), pow(), and specialcase().

fp_barrierf()

static float fp_barrierf ( float  x )

inlinestatic

Definition at line 108 of file libm.h.

{
    volatile float y = x;
    return y;
}

Referenced by ilogb().

fp_force_eval()

static void fp_force_eval ( double  x )

inlinestatic

Definition at line 57 of file libm.h.

{
    volatile double y;
    y = x;
    (void)y;
}

Referenced by specialcase().

fp_force_evalf()

static void fp_force_evalf ( float  x )

inlinestatic

Definition at line 50 of file libm.h.

{
    volatile float y;
    y = x;
    (void)y;
}

ilogb()

int ilogb

(

double

x

)

Definition at line 14 of file ilogb.c.

{
    #pragma STDC FENV_ACCESS ON
    union {double f; uint64_t i;} u = {x};
    uint64_t i = u.i;
    int e = i>>52 & 0x7ff;
 
    if (!e) {
        i <<= 12;
        if (i == 0) {
            FORCE_EVAL(fp_barrierf(0.0f)/0.0f);
            return FP_ILOGB0;
        }
        /* subnormal x */
        for (e = -0x3ff; i>>63 == 0; e--, i<<=1);
        return e;
    }
    if (e == 0x7ff) {
        FORCE_EVAL(fp_barrierf(0.0f)/0.0f);
        return i<<12 ? FP_ILOGBNAN : INT_MAX;
    }
    return e - 0x3ff;
}

logb()

double logb

(

double

x

)

Definition at line 243 of file mingw_math.h.

  {
    double res = 0.0;
    __asm__ __volatile__ ("fxtract\n\t"
      "fstp %%st" : "=t" (res) : "0" (x));
    return res;
  }

scalbn()

double scalbn	(	double	x,
		int	n
	)

Definition at line 14 of file scalbn.c.

{
    union {double f; uint64_t i;} u;
    double_t y = x;
 
    if (n > 1023) {
        y *= 0x1p1023;
        n -= 1023;
        if (n > 1023) {
            y *= 0x1p1023;
            n -= 1023;
            if (n > 1023)
                n = 1023;
        }
    } else if (n < -1022) {
        /* make sure final n < -53 to avoid double
           rounding in the subnormal range */
        y *= 0x1p-1022 * 0x1p53;
        n += 1022 - 53;
        if (n < -1022) {
            y *= 0x1p-1022 * 0x1p53;
            n += 1022 - 53;
            if (n < -1022)
                n = -1022;
        }
    }
    u.i = (uint64_t)(0x3ff+n)<<52;
    x = y * u.f;
    return x;
}