#include <float.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
/* returns the correctly rounded to float result of
the mathematical expression a*b+c.
Assumes round-to-nearest, a≥0, b≥0, c≥0, all finite. */
float myfma(float a, float b, float c)
{
double p1 = a * (double) b;
double p2 = c;
double s1 = p1 + p2;
if (p1 < p2)
{
double tmp = p1;
p1 = p2;
p2 = tmp;
}
double r1 = s1 - p1 - p2; /* fma = s1 + r1 */
float f1 = s1;
if (r1 == 0) return f1;
double t = s1 - f1;
if (t == 0) return f1; /* definitely not a halfway point */
double dir = copysign(1.0 / 0.0, t);
float f2 = nextafterf(f1, dir);
double ulp = f2 - (double) f1;
#if 0
printf("f1:%a\nf2:%a\np1:%a\np2:%a\ns1:%a\nt :%a\nr1:%a\nul:%a\n", f1, f2, p1, p2, s1, t, r1, ulp);
#endif
return f2;
return f1;
double r = t - 0.5 * ulp + r1;
if (r > 0 && t > 0 || r < 0 && t < 0)
return f2;
if (r > 0 && t < 0 || r < 0 && t > 0)
return f1;
/* Round to even. */
return f1 + 0.5 * ulp;
}
int main()
{
while (1)
{
int a
= (rand() & 0xFFFFF) << (rand() & 7); int b
= (rand() & 0xFFFFF) << (rand() & 7); long long c
= (long long)(rand() & 0xFFFFF) << (rand() & 31);
float truefma = (long long) a * b + c;
float r = myfma(a, b, c);
if (r != truefma)
{
printf("bad: %d %d %lld\ntr:%a\nmy:%a\n\n", a
, b
, c
, truefma
, r
); }
}
}
I2luY2x1ZGUgPGZsb2F0Lmg+CiNpbmNsdWRlIDxtYXRoLmg+CiNpbmNsdWRlIDxzdGRsaWIuaD4KI2luY2x1ZGUgPHN0ZGlvLmg+CgovKiByZXR1cm5zIHRoZSBjb3JyZWN0bHkgcm91bmRlZCB0byBmbG9hdCByZXN1bHQgb2YgCiAgIHRoZSBtYXRoZW1hdGljYWwgZXhwcmVzc2lvbiBhKmIrYy4gCiAgIEFzc3VtZXMgcm91bmQtdG8tbmVhcmVzdCwgYeKJpTAsIGLiiaUwLCBj4omlMCwgYWxsIGZpbml0ZS4gKi8KZmxvYXQgbXlmbWEoZmxvYXQgYSwgZmxvYXQgYiwgZmxvYXQgYykgCnsKICBkb3VibGUgcDEgPSBhICogKGRvdWJsZSkgYjsgCiAgZG91YmxlIHAyID0gYzsKICBkb3VibGUgczEgPSBwMSArIHAyOwogIGlmIChwMSA8IHAyKQogICAgewogICAgICBkb3VibGUgdG1wID0gcDE7CiAgICAgIHAxID0gcDI7CiAgICAgIHAyID0gdG1wOwogICAgfQogIGRvdWJsZSByMSA9IHMxIC0gcDEgLSBwMjsgLyogZm1hID0gczEgKyByMSAqLwoKICBmbG9hdCBmMSA9IHMxOyAgICAgICAgICAgIAoKICBpZiAocjEgPT0gMCkgcmV0dXJuIGYxOwoKICBkb3VibGUgdCA9IHMxIC0gZjE7CgogIGlmICh0ID09IDApIHJldHVybiBmMTsgLyogZGVmaW5pdGVseSBub3QgYSBoYWxmd2F5IHBvaW50ICovCgogIGRvdWJsZSBkaXIgPSBjb3B5c2lnbigxLjAgLyAwLjAsIHQpOwogIGZsb2F0IGYyID0gbmV4dGFmdGVyZihmMSwgZGlyKTsKICBkb3VibGUgdWxwID0gZjIgLSAoZG91YmxlKSBmMTsKICAKI2lmIDAKICBwcmludGYoImYxOiVhXG5mMjolYVxucDE6JWFcbnAyOiVhXG5zMTolYVxudCA6JWFcbnIxOiVhXG51bDolYVxuIiwgCiAgICAgICAgICBmMSwgICAgZjIsICAgIHAxLCAgICBwMiwgICAgczEsICAgIHQsICAgICByMSwgICAgdWxwKTsKI2VuZGlmCiAgCiAgaWYgKGZhYnModCkgPiBmYWJzKDAuNzUgKiB1bHApKQogICAgcmV0dXJuIGYyOwoKICBpZiAoZmFicyh0KSA8IGZhYnMoMC4yNSAqIHVscCkpCiAgICByZXR1cm4gZjE7CgogIGRvdWJsZSByID0gdCAtIDAuNSAqIHVscCArIHIxOwoKICBpZiAociA+IDAgJiYgdCA+IDAgfHwgciA8IDAgJiYgdCA8IDApCiAgICByZXR1cm4gZjI7CgogIGlmIChyID4gMCAmJiB0IDwgMCB8fCByIDwgMCAmJiB0ID4gMCkKICAgIHJldHVybiBmMTsKCiAgLyogUm91bmQgdG8gZXZlbi4gKi8KICByZXR1cm4gZjEgKyAwLjUgKiB1bHA7Cn0KICAKaW50IG1haW4oKQp7CiAgd2hpbGUgKDEpCiAgICB7CiAgICAgIGludCBhID0gKHJhbmQoKSAmIDB4RkZGRkYpIDw8IChyYW5kKCkgJiA3KTsKICAgICAgaW50IGIgPSAocmFuZCgpICYgMHhGRkZGRikgPDwgKHJhbmQoKSAmIDcpOwogICAgICBsb25nIGxvbmcgYyA9IChsb25nIGxvbmcpKHJhbmQoKSAmIDB4RkZGRkYpIDw8IChyYW5kKCkgJiAzMSk7CgogICAgICBmbG9hdCB0cnVlZm1hID0gKGxvbmcgbG9uZykgYSAqIGIgKyBjOwogICAgICBmbG9hdCByID0gbXlmbWEoYSwgYiwgYyk7CiAgICAgIGlmIChyICE9IHRydWVmbWEpCiAgICAgICAgewogICAgICAgICAgcHJpbnRmKCJiYWQ6ICVkICVkICVsbGRcbnRyOiVhXG5teTolYVxuXG4iLCBhLCBiLCBjLCB0cnVlZm1hLCByKTsKICAgICAgICB9CiAgICB9Cn0K