Method exp2_est

Source/ThirdParty/basis_universal/encoder/cppspmd_math.h:186–257

http://www.ganssle.com/item/approximations-c-code-exponentiation-log.htm

// http://www.ganssle.com/item/approximations-c-code-exponentiation-log.htm
inline vfloat spmd_kernel::exp2_est(vfloat arg)
{
	SPMD_BEGIN_CALL

	const vfloat P00 = +7.2152891521493f;
	const vfloat P01 = +0.0576900723731f;
	const vfloat Q00 = +20.8189237930062f;
	const vfloat Q01 = +1.0f;
	const vfloat sqrt2 = 1.4142135623730950488f; // sqrt(2) for scaling

	vfloat result = 0.0f;

	// Return 0 if arg is too large.
	// We're not introducing inf/nan's into calculations, or risk doing so by returning huge default values.
	SPMD_IF(abs(arg) > 126.0f)
	{
		spmd_return();
	}
	SPMD_END_IF

	// 2**(int(a))
	vfloat two_int_a;

	// set to 1 by reduce_expb
	vint adjustment;

	// 0 if arg is +; 1 if negative
	vint negative = 0;

	// If the input is negative, invert it. At the end we'll take the reciprocal, since n**(-1) = 1/(n**x).
	SPMD_SIF(arg < 0.0f)
	{
		store(arg, -arg);
		store(negative, 1);
	}
	SPMD_SENDIF

	store_all(arg, min(arg, 126.0f));

	// reduce to [0.0, 0.5]
	reduce_expb(arg, two_int_a, adjustment);

	// The format of the polynomial is:
	//  answer=(Q(x**2) + x*P(x**2))/(Q(x**2) - x*P(x**2))
	//
	//  The following computes the polynomial in several steps:

	// Q(x**2)
	vfloat Q = vfma(Q01, (arg * arg), Q00);

	// x*P(x**2)
	vfloat x_P = arg * (vfma(P01, arg * arg, P00));

	vfloat answer = (Q + x_P) / (Q - x_P);

	// Now correct for the scaling factor of 2**(int(a))
	store_all(answer, answer * two_int_a);