Convert log10f() to use __kernel_log(), which is more accurate and simpler.

This commit is contained in:
David Schultz 2011-03-07 03:12:08 +00:00
parent acda0929b2
commit 97a539be1f
Notes: svn2git 2020-12-20 02:59:44 +00:00
svn path=/head/; revision=219361

View File

@ -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;
}