xref: /freebsd/lib/msun/src/e_log2f.c (revision f5f40dd63bc7acbb5312b26ac1ea1103c12352a6)
1 /*
2  * ====================================================
3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4  *
5  * Developed at SunPro, a Sun Microsystems, Inc. business.
6  * Permission to use, copy, modify, and distribute this
7  * software is freely granted, provided that this notice
8  * is preserved.
9  * ====================================================
10  */
11 
12 /*
13  * Float version of e_log2.c.  See the latter for most comments.
14  */
15 
16 #include "math.h"
17 #include "math_private.h"
18 #include "k_logf.h"
19 
20 static const float
21 two25      =  3.3554432000e+07, /* 0x4c000000 */
22 ivln2hi    =  1.4428710938e+00, /* 0x3fb8b000 */
23 ivln2lo    = -1.7605285393e-04; /* 0xb9389ad4 */
24 
25 static const float zero   =  0.0;
26 static volatile float vzero = 0.0;
27 
28 float
29 log2f(float x)
30 {
31 	float f,hfsq,hi,lo,r,y;
32 	int32_t i,k,hx;
33 
34 	GET_FLOAT_WORD(hx,x);
35 
36 	k=0;
37 	if (hx < 0x00800000) {			/* x < 2**-126  */
38 	    if ((hx&0x7fffffff)==0)
39 		return -two25/vzero;		/* log(+-0)=-inf */
40 	    if (hx<0) return (x-x)/zero;	/* log(-#) = NaN */
41 	    k -= 25; x *= two25; /* subnormal number, scale up x */
42 	    GET_FLOAT_WORD(hx,x);
43 	}
44 	if (hx >= 0x7f800000) return x+x;
45 	if (hx == 0x3f800000)
46 	    return zero;			/* log(1) = +0 */
47 	k += (hx>>23)-127;
48 	hx &= 0x007fffff;
49 	i = (hx+(0x4afb0d))&0x800000;
50 	SET_FLOAT_WORD(x,hx|(i^0x3f800000));	/* normalize x or x/2 */
51 	k += (i>>23);
52 	y = (float)k;
53 	f = x - (float)1.0;
54 	hfsq = (float)0.5*f*f;
55 	r = k_log1pf(f);
56 
57 	/*
58 	 * We no longer need to avoid falling into the multi-precision
59 	 * calculations due to compiler bugs breaking Dekker's theorem.
60 	 * Keep avoiding this as an optimization.  See e_log2.c for more
61 	 * details (some details are here only because the optimization
62 	 * is not yet available in double precision).
63 	 *
64 	 * Another compiler bug turned up.  With gcc on i386,
65 	 * (ivln2lo + ivln2hi) would be evaluated in float precision
66 	 * despite runtime evaluations using double precision.  So we
67 	 * must cast one of its terms to float_t.  This makes the whole
68 	 * expression have type float_t, so return is forced to waste
69 	 * time clobbering its extra precision.
70 	 */
71 	if (sizeof(float_t) > sizeof(float))
72 		return (r - hfsq + f) * ((float_t)ivln2lo + ivln2hi) + y;
73 
74 	hi = f - hfsq;
75 	GET_FLOAT_WORD(hx,hi);
76 	SET_FLOAT_WORD(hi,hx&0xfffff000);
77 	lo = (f - hi) - hfsq + r;
78 	return (lo+hi)*ivln2lo + lo*ivln2hi + hi*ivln2hi + y;
79 }
80