fma_ieee754_soft.cc

template <class T
//cpp11 , typename std::enable_if<std::is_floating_point<T>::value>::type* = nullptr	
>
T fma(const T& x, const T& y, const T& z)
{
	const T cs = TwoProduct<T>::split_value;

	T s = z, w = T(0), h, q, r, x1, x2, y1, y2;

	// can be a loop where x and y are values at index and w, s are accumulators. fmav.
	{
//!# TwoProduct(x,y,h,r).
	q  = x;
//!# split x into x1,x2.
	r  = cs * q;
	x2 = r - q;
	x1 = r - x2;
	x2 = q - x1;
	r  = y;
//!# h=x*y.
	h  = q * r;
//!# split y into y1,y2.
	q  = cs * r;
	y2 = q - r;
	y1 = q - y2;
	y2 = r - y1;
//!# r=x2*y2-(((h-x1*y1) - x2*y1) - x1*y2
	q  = x1 * y1;
	q  = h - q;
	y1 = y1 * x2;
	q  = q - y1;
	x1 = x1 * y2;
	q  = q - x1;
	x2 = x2 * y2;
	r  = x2 - q;
//!# (w,q)=TwoSum(w,h).
	x1 = w + h;
	x2 = x1 - w;
	y1 = x1 - x2;
	y2 = h - x2;
	q  = w - y1;
	q  = q + y2;
	w  = x1;
//!# s=s+(q+r).
	q = q + r;
	s = s + q;
	}
			
	return w + s;
}