/*
* fmaf.c
*
* by Ian Ollmann
*
* Copyright (c) 2007, Apple Inc. All Rights Reserved.
*
* C implementation of C99 fmaf() function.
*/
#include
#include
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;
}