fmaf.c   [plain text]

/*
 *	fmaf.c
 *
 *		by Ian Ollmann
 *
 *	Copyright (c) 2007, Apple Inc. All Rights Reserved.
 *
 *	C implementation of C99 fmaf() function.
 */
 
#include <math.h>
#include <stdint.h>

float fmaf( float a, float b, float c )
{
	double product = (double) a * (double) b;		//exact
	double dc = (double) c;							//exact

	#warning fmaf not completely correct
	// Simply adding C here is incorrect about 1 in a billion times.
	// While the double precision add here is correctly rounded,
	// we take a second rounding on conversion to float on return
	// which may cause us to be off by very slightly over half an ulp
	// in round to nearest. 
	double sum = product + dc;

	// ideally, we should test here and patch up the result.
	// I think the problem only occurs in round to nearest for 
	// exact half way cases in product with a non-zero c.
	// Presumably, we could check to see if the difference between
	// (float) sum and sum is a power of two (the right exact power 
	// of two) and c is non-zero, and it rounded the wrong way, then 
	// we might tweak the answer by an ulp using something like nextafter.
	// Happily denormals are not a problem during this check.
    //
	// Alternatively, if we figure out the problem of correctly rounded
	// 3-way adds, the product could be broken into 2 floats, and we 
	// could do a 3-way add of prodHi, prodLo and c. Crlibm has a function 
	// that might do the job (DoRenormalize3), bu Im thinking that it doesnt.
	//
	// Finally, to be completely right, we'd have to detect rounding mode.
	// The half way cases are different in other rounding modes.
	
	return (float) sum;
}