Fast IEEE754 float logarithm approximation in C

fast_log.h

/**
	Fast logarithm approximation (around 1.5) with Taylor series.
	
	sage: ln(x).taylor(x, 1.5, 5).horner(x)
	(((x*(0.02633744855967078*x - 0.24691358024691357) + 0.9876543209876543)*x - 2.2222222222222223)*x + 3.333333333333333)*x - 1.8778682252251688
*/
static inline float fast_logf_taylor( float x )
{
	return fmaf( fmaf( fmaf( fmaf( fmaf(0.02633744855967078f, x, -0.24691358024691357f), x, 0.9876543209876543f), x, -2.2222222222222223f), x, 3.333333333333333f), x, -1.8778682252251688f);
};

static inline float fast_log2f_taylor( float x )
{
	return fast_logf_taylor( x ) * 1.44269504089f;
};

/**
	\brief Computes log2(x) based on Taylor series and floating point
	bitwise manipulation.

	\author Jacek Wieczorek
	
	Theory:
		Let 'x' denote a 32-bit floating point number.
		Thus, it can be written as:
		x = s * m * 2^e
		
		Where:
			- s is a sign bit assumed to be one
			- m is a floating point mantissa in range [1;2)
			- e is an integer exponent in range [-128;127]
		
		Floating point mantissa and exponent can be expressed as:
		m = 1 + mi / mimax
		e = ei - 127

		Where:
			- mi is integer value of mantissa bits
			- mimax is max value mantissa can take
			- ei is integer value of exponent bits
		
		Therefore, we can rewrite x as:
		x = (1 + mi / mimax) * 2^(ei - 127)

		And it's base 2 logarithm:
		lg(x) = ei - 127 + lg(1 + mi/mimax)

		Since mimax is constant, the division can be replaced with multiplication
		by reciprocal:
		lg(x) = ei - 127 + lg(1 + mi * (1/mimax))

		Since the logarithm's on the right argument ranges from 1 to 2 it
		can be accurately approximated with Taylor series.
		
		Denormal numbers require special treatment, however, since they are
		located near zero, the logarithm's output value changes abruptly in that region
		anyway, therefore the resulting inaccuracy is negligible.

*/
static inline float fast_log2f( float x )
{
	union { float xf; uint32_t xi; };

	// Extract exponent and mantissa
	xf = x;
	int ei = ( ( xi >> 23 ) & 0xff ) - 127; // Exponent
	int mi = xi & 0x7fffff;                 // Mantissa
	
	// Real mantissa value
	float mf = 1.f + mi * ( 1.f / 0x7fffff );
	
	// Denormal numbers (optional)
	// if ( ei == -127 ) mf = mi * ( 1.f / 0x7fffff );

	return ei + fast_log2f_taylor( mf );
}



/**
	Fast natural logarithm
*/
static inline float fast_logf( float x )
{
	static const float c = std::log( 2 );
	return c * fast_log2f( x );
}