Method log2_est

engine/dlib/src/basis/encoder/cppspmd_math.h:65–121

See https://tech.ebayinc.com/engineering/fast-approximate-logarithms-part-iii-the-formulas/

// See https://tech.ebayinc.com/engineering/fast-approximate-logarithms-part-iii-the-formulas/
inline vfloat spmd_kernel::log2_est(vfloat v)
{
	vfloat signif, fexp;

	// Just clamp to a very small value, instead of checking for invalid inputs.
	vfloat x = max(v, 2.2e-38f);

	/*
	 * Assume IEEE representation, which is sgn(1):exp(8):frac(23)
	 * representing (1+frac)*2^(exp-127).  Call 1+frac the significand
	 */

	 // get exponent
	vint ux1_i = cast_vfloat_to_vint(x);

	vint exp = VUINT_SHIFT_RIGHT(ux1_i & 0x7F800000, 23);

	// actual exponent is exp-127, will subtract 127 later

	vint ux2_i;
	vfloat ux2_f;

	vint greater = ux1_i & 0x00400000;  // true if signif > 1.5
	SPMD_SIF(greater != 0)
	{
		// signif >= 1.5 so need to divide by 2.  Accomplish this by stuffing exp = 126 which corresponds to an exponent of -1 
		store_all(ux2_i, (ux1_i & 0x007FFFFF) | 0x3f000000);

		store_all(ux2_f, cast_vint_to_vfloat(ux2_i));

		// 126 instead of 127 compensates for division by 2
		store_all(fexp, vfloat(exp - 126));    
	}
	SPMD_SELSE(greater != 0)
	{
		// get signif by stuffing exp = 127 which corresponds to an exponent of 0
		store(ux2_i, (ux1_i & 0x007FFFFF) | 0x3f800000);

		store(ux2_f, cast_vint_to_vfloat(ux2_i));

		store(fexp, vfloat(exp - 127));
	}
	SPMD_SENDIF

	store_all(signif, ux2_f);
	store_all(signif, signif - 1.0f);

	const float a = 0.1501692f, b = 3.4226132f, c = 5.0225057f, d = 4.1130283f, e = 3.4813372f;

	vfloat xm1 = signif;
	vfloat xm1sqr = xm1 * xm1;
		
	return fexp + ((a * (xm1sqr * xm1) + b * xm1sqr + c * xm1) / (xm1sqr + d * xm1 + e));
	
	// fma lowers accuracy for SSE4.1 - no idea why (compiler reordering?)
	//return fexp + ((vfma(a, (xm1sqr * xm1), vfma(b, xm1sqr, c * xm1))) / (xm1sqr + vfma(d, xm1, e)));
}