Open
Description
I was reading an interesting article about posits (thanks @simonbyrne!) and came across Kahan's algorithm for calculating discriminant which avoids catastrophic cancellation and looks cheap-ish on hardware with FMA. At some stage we should probably see whether it's a reasonable performance tradeoff for 2x2 _det
and related computations:
// computes: ad-bc within +/- 3/2 ulp of exact
double discriminant(double a, double b, double c, double d)
{
double w = b*c;
double e = fma(-b,c,w);
double f = fma(a,d,-w);
return f+e;
}
Further analysis: http://www.ams.org/journals/mcom/2013-82-284/S0025-5718-2013-02679-8/home.html