1 /* 2009 for Newlib:  Sun's s_ilogb.c converted to be s_logb.c.  */
2 /* @(#)s_ilogb.c 5.1 93/09/24 */
3 /*
4  * ====================================================
5  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6  *
7  * Developed at SunPro, a Sun Microsystems, Inc. business.
8  * Permission to use, copy, modify, and distribute this
9  * software is freely granted, provided that this notice
10  * is preserved.
11  * ====================================================
12  */
13 /*
14 FUNCTION
15        <<logb>>, <<logbf>>---get exponent of floating-point number
16 INDEX
17 	logb
18 INDEX
19 	logbf
20 
21 SYNOPSIS
22 	#include <math.h>
23         double logb(double <[x]>);
24         float logbf(float <[x]>);
25 
26 DESCRIPTION
27 The <<logb>> functions extract the exponent of <[x]>, as a signed integer value
28 in floating-point format.  If <[x]> is subnormal it is treated as though it were
29 normalized; thus, for positive finite <[x]>,
30 @ifnottex
31 1 <= (<[x]> * FLT_RADIX to the power (-logb(<[x]>))) < FLT_RADIX.
32 @end ifnottex
33 @tex
34 $1 \leq ( x \cdot FLT\_RADIX ^ {-logb(x)} ) < FLT\_RADIX$.
35 @end tex
36 A domain error may occur if the argument is zero.
37 In this floating-point implementation, FLT_RADIX is 2.  Which also means
38 that for finite <[x]>, <<logb>>(<[x]>) = <<floor>>(<<log2>>(<<fabs>>(<[x]>))).
39 
40 All nonzero, normal numbers can be described as
41 @ifnottex
42 <[m]> * 2**<[p]>, where 1.0 <= <[m]> < 2.0.
43 @end ifnottex
44 @tex
45 $m \cdot 2^p$, where $1.0 \leq m < 2.0$.
46 @end tex
47 The <<logb>> functions examine the argument <[x]>, and return <[p]>.
48 The <<frexp>> functions are similar to the <<logb>> functions, but
49 returning <[m]> adjusted to the interval [.5, 1) or 0, and <[p]>+1.
50 
51 RETURNS
52 @comment Formatting note:  "$@" forces a new line
53 When <[x]> is:@*
54 +inf or -inf, +inf is returned;@*
55 NaN, NaN is returned;@*
56 0, -inf is returned, and the divide-by-zero exception is raised;@*
57 otherwise, the <<logb>> functions return the signed exponent of <[x]>.
58 
59 PORTABILITY
60 ANSI C, POSIX
61 
62 SEEALSO
63 frexp, ilogb
64 */
65 
66 /* double logb(double x)
67  * return the binary exponent of non-zero x
68  * logb(0) = -inf, raise divide-by-zero floating point exception
69  * logb(+inf|-inf) = +inf (no signal is raised)
70  * logb(NaN) = NaN (no signal is raised)
71  * Per C99 recommendation, a NaN argument is returned unchanged.
72  */
73 
74 #include "fdlibm.h"
75 
76 #ifdef _NEED_FLOAT64
77 
78 __float64
logb64(__float64 x)79 logb64(__float64 x)
80 {
81 	__int32_t hx,lx,ix;
82 
83 	EXTRACT_WORDS(hx,lx,x);
84 	hx &= 0x7fffffff;		/* high |x| */
85 	if(hx<0x00100000) {		/* 0 or subnormal */
86 	    if((hx|lx)==0)  {
87 		/* arg==0:  return -inf and raise divide-by-zero exception */
88 		return -1./fabs64(x);	/* logb(0) = -inf */
89 		}
90 	    else			/* subnormal x */
91 		if(hx==0) {
92 		    for (ix = -1043; lx>0; lx<<=1) ix -=1;
93 		} else {
94 		    for (ix = -1022,hx<<=11; hx>0; hx<<=1) ix -=1;
95 		}
96 	    return (__float64) ix;
97 	}
98 	else if (hx<0x7ff00000) return (hx>>20)-1023;	/* normal # */
99 	else if (hx>0x7ff00000 || lx)  return x+x;	/* x==NaN */
100 	else  return HUGE_VAL;	/* x==inf (+ or -) */
101 }
102 
103 _MATH_ALIAS_d_d(logb)
104 
105 #endif /* _NEED_FLOAT64 */
106