fp-round.cpp

// Fixed version of code from
// http://blog.frama-c.com/index.php?post/2013/05/02/nearbyintf1#pr

#include <stdio.h>
#include <limits.h>
#include <math.h>

// Original function
// float myround(float f)
// {
//   if (f >= 0x1.0p23) return f;
//   return (float) (unsigned int) (f + 0.49999997f);
// }

float myround(float f)
{
    if (f <= -0x1.0p23 || f >= 0x1.0p23)
        return f;
    if (-0.5 < f && f < 0.5f)
        return 0.0f;
    if (f > 0.0f)
        return (float) (int) (f + 0.49999997f);
    else
        return (float) (int) (f - 0.49999997f);
}

void expect_round(float f, float expected_rounded)
{
    float rounded = myround(f);
    if(rounded != expected_rounded)
    {
        printf("FAIL myround(%1.9e) = %1.9e but expected %1.9e\n", f, rounded, expected_rounded);
    }

    float std_round = roundf(f);
    if(rounded != std_round)
    {
        printf("FAIL myround(%1.9e) = %1.9e but round is %1.9e", f, rounded, std_round);
    }
}

void check(float f, float expected_rounded)
{
    expect_round(f, expected_rounded);
    expect_round(-f, -expected_rounded);
}

void check(float f_is_int)
{
    check(f_is_int, f_is_int);
}

int main(int argc, char const *argv[])
{
    check(0.0f);
    check(1.0f);
    check(0.5f, 1.0f);
    check(1.5f, 2.0f);
    check(101.5f, 102.0f);

    // Very large value that would overflow int
    check(float(INT_MAX) + 1000.0f);

    // Predecessor of 0.5/1.5. Mathematically correct result will be exactly
    // between 1.0/2.0 and it's predecessor.
    check(0.49999997f, 0.0f);
    check(nexttowardf(1.5f, 0.0f), 1.0f);

    // Value whose ULP (distance to last/prev float) is exactly 1.0
    check(8388609.0f);
    check(8388610.0f);

    return 0;
}