1*05a0b428SJohn Marino /* @(#)s_logb.c 5.1 93/09/24 */
2*05a0b428SJohn Marino /*
3*05a0b428SJohn Marino * ====================================================
4*05a0b428SJohn Marino * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5*05a0b428SJohn Marino *
6*05a0b428SJohn Marino * Developed at SunPro, a Sun Microsystems, Inc. business.
7*05a0b428SJohn Marino * Permission to use, copy, modify, and distribute this
8*05a0b428SJohn Marino * software is freely granted, provided that this notice
9*05a0b428SJohn Marino * is preserved.
10*05a0b428SJohn Marino * ====================================================
11*05a0b428SJohn Marino */
12*05a0b428SJohn Marino
13*05a0b428SJohn Marino /*
14*05a0b428SJohn Marino * double logb(x)
15*05a0b428SJohn Marino * IEEE 754 logb. Included to pass IEEE test suite. Not recommend.
16*05a0b428SJohn Marino * Use ilogb instead.
17*05a0b428SJohn Marino */
18*05a0b428SJohn Marino
19*05a0b428SJohn Marino #include "math.h"
20*05a0b428SJohn Marino #include "math_private.h"
21*05a0b428SJohn Marino
22*05a0b428SJohn Marino double
logb(double x)23*05a0b428SJohn Marino logb(double x)
24*05a0b428SJohn Marino {
25*05a0b428SJohn Marino int32_t lx,ix;
26*05a0b428SJohn Marino EXTRACT_WORDS(ix,lx,x);
27*05a0b428SJohn Marino ix &= 0x7fffffff; /* high |x| */
28*05a0b428SJohn Marino if((ix|lx)==0) return -1.0/fabs(x);
29*05a0b428SJohn Marino if(ix>=0x7ff00000) return x*x;
30*05a0b428SJohn Marino if((ix>>=20)==0) /* IEEE 754 logb */
31*05a0b428SJohn Marino return -1022.0;
32*05a0b428SJohn Marino else
33*05a0b428SJohn Marino return (double) (ix-1023);
34*05a0b428SJohn Marino }
35