From 97a539be1f02398505a3331996bed3e83a91664b Mon Sep 17 00:00:00 2001
From: David Schultz <das@FreeBSD.org>
Date: Mon, 7 Mar 2011 03:12:08 +0000
Subject: [PATCH] Convert log10f() to use __kernel_log(), which is more
 accurate and simpler.

---
 lib/msun/src/e_log10f.c | 30 +++++++++++++++++++-----------
 1 file changed, 19 insertions(+), 11 deletions(-)

diff --git a/lib/msun/src/e_log10f.c b/lib/msun/src/e_log10f.c
index a0b1618a4437..940b8319f804 100644
--- a/lib/msun/src/e_log10f.c
+++ b/lib/msun/src/e_log10f.c
@@ -1,7 +1,3 @@
-/* e_log10f.c -- float version of e_log10.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -16,12 +12,18 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
+/*
+ * Return the base 10 logarithm of x. See k_log.c for details on the algorithm.
+ */
+
 #include "math.h"
 #include "math_private.h"
+#include "k_logf.h"
 
 static const float
 two25      =  3.3554432000e+07, /* 0x4c000000 */
-ivln10     =  4.3429449201e-01, /* 0x3ede5bd9 */
+ivln10hi   =  4.3432617188e-01, /* 0x3ede6000 */
+ivln10lo   = -3.1689971365e-05, /* 0xb804ead9 */
 log10_2hi  =  3.0102920532e-01, /* 0x3e9a2080 */
 log10_2lo  =  7.9034151668e-07; /* 0x355427db */
 
@@ -30,7 +32,7 @@ static const float zero   =  0.0;
 float
 __ieee754_log10f(float x)
 {
-	float y,z;
+	float f,hi,lo,y,z;
 	int32_t i,k,hx;
 
 	GET_FLOAT_WORD(hx,x);
@@ -45,10 +47,16 @@ __ieee754_log10f(float x)
         }
 	if (hx >= 0x7f800000) return x+x;
 	k += (hx>>23)-127;
-	i  = ((u_int32_t)k&0x80000000)>>31;
-        hx = (hx&0x007fffff)|((0x7f-i)<<23);
-        y  = (float)(k+i);
-	SET_FLOAT_WORD(x,hx);
-	z  = y*log10_2lo + ivln10*__ieee754_logf(x);
+	hx &= 0x007fffff;
+	i = (hx+(0x4afb0d))&0x800000;
+	SET_FLOAT_WORD(x,hx|(i^0x3f800000));	/* normalize x or x/2 */
+	k += (i>>23);
+	y = (float)k;
+	f = __kernel_logf(x);
+	x = x - (float)1.0;
+	GET_FLOAT_WORD(hx,x);
+	SET_FLOAT_WORD(hi,hx&0xfffff000);
+	lo = x - hi;
+	z = y*log10_2lo + (x+f)*ivln10lo + (lo+f)*ivln10hi + hi*ivln10hi;
 	return  z+y*log10_2hi;
 }