1 /* e_log10f.c -- float version of e_log10.c.
2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
6 * ====================================================
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
9 * Developed at SunPro, a Sun Microsystems, Inc. business.
10 * Permission to use, copy, modify, and distribute this
11 * software is freely granted, provided that this notice
13 * ====================================================
17 static char rcsid
[] = "$FreeBSD: src/lib/msun/src/e_log10f.c,v 1.13 2011/10/16 05:36:23 das Exp $";
21 * Float version of e_log10.c. See the latter for most comments.
25 #include "math_private.h"
29 two25
= 3.3554432000e+07, /* 0x4c000000 */
30 ivln10hi
= 4.3432617188e-01, /* 0x3ede6000 */
31 ivln10lo
= -3.1689971365e-05, /* 0xb804ead9 */
32 log10_2hi
= 3.0102920532e-01, /* 0x3e9a2080 */
33 log10_2lo
= 7.9034151668e-07; /* 0x355427db */
35 static const float zero
= 0.0;
36 static const volatile float vzero
__attribute__ ((__section__(".rodata"))) = 0.0;
39 __ieee754_log10f(float x
)
41 float f
,hfsq
,hi
,lo
,r
,y
;
47 if (hx
< 0x00800000) { /* x < 2**-126 */
48 if ((hx
&0x7fffffff)==0)
49 return -two25
/vzero
; /* log(+-0)=-inf */
50 if (hx
<0) return (x
-x
)/zero
; /* log(-#) = NaN */
51 k
-= 25; x
*= two25
; /* subnormal number, scale up x */
54 if (hx
>= 0x7f800000) return x
+x
;
56 return zero
; /* log(1) = +0 */
59 i
= (hx
+(0x4afb0d))&0x800000;
60 SET_FLOAT_WORD(x
,hx
|(i
^0x3f800000)); /* normalize x or x/2 */
64 hfsq
= (float)0.5*f
*f
;
67 /* See e_log2f.c and e_log2.c for details. */
68 if (sizeof(float_t
) > sizeof(float))
69 return (r
- hfsq
+ f
) * ((float_t
)ivln10lo
+ ivln10hi
) +
70 y
* ((float_t
)log10_2lo
+ log10_2hi
);
72 GET_FLOAT_WORD(hx
,hi
);
73 SET_FLOAT_WORD(hi
,hx
&0xfffff000);
74 lo
= (f
- hi
) - hfsq
+ r
;
75 return y
*log10_2lo
+ (lo
+hi
)*ivln10lo
+ lo
*ivln10hi
+ hi
*ivln10hi
+