s_ilogb.c revision 1.3 
1/* @(#)s_ilogb.c 5.1 93/09/24 */
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13#if defined(LIBM_SCCS) && !defined(lint)
14static char rcsid[] = "$NetBSD: s_ilogb.c,v 1.9 1995/05/10 20:47:28 jtc Exp $";
15#endif
16
17/* ilogb(double x)
18 * return the binary exponent of non-zero x
19 * ilogb(0) = 0x80000001
20 * ilogb(inf/NaN) = 0x7fffffff (no signal is raised)
21 */
22
23#include <machine/cdefs.h>
24#include <float.h>
25#include <math.h>
26
27#include "math_private.h"
28
29int
30ilogb(double x)
31{
32	int32_t hx,lx,ix;
33
34	GET_HIGH_WORD(hx,x);
35	hx &= 0x7fffffff;
36	if(hx<0x00100000) {
37	    GET_LOW_WORD(lx,x);
38	    if((hx|lx)==0)
39		return 0x80000001;	/* ilogb(0) = 0x80000001 */
40	    else			/* subnormal x */
41		if(hx==0) {
42		    for (ix = -1043; lx>0; lx<<=1) ix -=1;
43		} else {
44		    for (ix = -1022,hx<<=11; hx>0; hx<<=1) ix -=1;
45		}
46	    return ix;
47	}
48	else if (hx<0x7ff00000) return (hx>>20)-1023;
49	else return 0x7fffffff;
50}
51
52#if LDBL_MANT_DIG == 53
53#ifdef __weak_alias
54__weak_alias(ilogbl, ilogb);
55#endif /* __weak_alias */
56#endif /* LDBL_MANT_DIG == 53 */
57