A sample showing an example of finding limits of floating point accumulation

fp_accumulation_limits.c

/*------------------------------------------------------------------------------------------------------------------
 * A sample that demonstrates 32-bit floating point precision
 *
 *  'compute_upper_bound_f32' finds such floating point number for given 'x' that
 *      upper_bound + x == upper_bound
 *
 *  'compute_lower_bound_f32' finds such floating point number for given 'x' that
 *      x + lower_bound == x
 *----------------------------------------------------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

typedef float f32;
typedef unsigned int u32;

typedef union _f32_to_u32
{
    struct
    {
        u32 m : 23;
        u32 e :  8;
        u32 s :  1;
    }   as_bit;
    f32 as_f32;
    u32 as_u32;
} f32_to_u32;

static f32 next_f32(f32 x)
{
    f32_to_u32 cvt;
    cvt.as_f32  = x;
    cvt.as_u32 += 1;
    return cvt.as_f32;
}

static f32 prev_f32(f32 x)
{
    f32_to_u32 cvt;
    cvt.as_f32  = x;
    cvt.as_u32 -= 1;
    return cvt.as_f32;
}

static f32 compute_upper_bound_f32(f32 x)
{
    f32_to_u32 cvt;
    cvt.as_f32 = x;
    cvt.as_bit.e += (cvt.as_bit.m != 0);
    cvt.as_bit.m  =  0;
    cvt.as_bit.e += 24;
    return cvt.as_f32;
}

static f32 compute_lower_bound_f32(f32 x)
{
    f32_to_u32 cvt;
    cvt.as_f32 = x;
    cvt.as_bit.m  = 0x7fffff;
    cvt.as_bit.e -= 25;
    return cvt.as_f32;
}

static u32 num_passed_sums(f32_to_u32 sum, f32 x)
{
    u32 passed = 0;
    sum.as_bit.m = 0;
    for (u32 i = 0; i < (1 << 23); ++i, sum.as_bit.m += 1)
    {
        passed += (sum.as_f32 + x != sum.as_f32);
    }
    return passed;
}

static u32 test(f32 upper_bound, f32 x)
{
    f32_to_u32 cvt;
    cvt.as_f32 = upper_bound;

    cvt.as_bit.m = 0x7fffff;
    f32 ub = cvt.as_f32;

    cvt.as_bit.m = 0;
    f32 lb = cvt.as_f32;

    printf("\t Checking range [%.9f, %.9f] ... ", ub, lb);
    u32 num_passed = num_passed_sums(cvt, x);
    printf("%u test accumulations passed out of %u \n", num_passed, 1 << 23);

    cvt.as_bit.m  = 0x7fffff;
    cvt.as_bit.e -= 1;
    ub = cvt.as_f32;

    cvt.as_bit.m  = 0;
    lb = cvt.as_f32;

    printf("\t Checking range [%.9f, %.9f] ... ", ub, lb);
    u32 num_failed = (1 << 23) - num_passed_sums(cvt, x);
    printf("%u test accumulations failed out of %u \n", num_failed, 1 << 23);

    if (num_passed != 1 << 23 && num_failed != 1 << 23 && (x + upper_bound == upper_bound))
    {
        printf("[Succeed] ");
        return 1;
    }
    else
    {
        printf("[Failed] \n");
        return 0;
    }
}

static void test_upper_bound(f32 x)
{
    printf("Trying to find Sup(%.9f) ... \n ", x);
    f32 upper_bound = compute_upper_bound_f32(x);

    if (test(upper_bound, x))
    {
        printf("Sup(%.9f) = %.9f \n", x, upper_bound);
    }
}

static void test_lower_bound(f32 x)
{
    printf("Trying to find Inf(%.8f) ... \n", x);
    f32 lower_bound = compute_lower_bound_f32(x);

    if (test(x, lower_bound))
    {
        printf("Inf(%.9f) = %.9f \n", x, lower_bound);
    }
}

int main()
{
    test_upper_bound(0.03334f);
    test_upper_bound(0.01667f);

    test_upper_bound(prev_f32(prev_f32(0.0625f)));
    test_upper_bound(prev_f32(0.0625f));
    test_upper_bound(0.0625f);
    test_upper_bound(next_f32(0.0625f));
    test_upper_bound(next_f32(next_f32(0.0625f)));


    test_lower_bound(prev_f32(prev_f32(1048576.0f)));
    test_lower_bound(prev_f32(1048576.0f));
    test_lower_bound(1048576.0f);
    test_lower_bound(next_f32(1048576.0f));
    test_lower_bound(next_f32(next_f32(1048576.0f)));

    srand(time(0));

    for (u32 i = 0; i < 12; ++i)
    {
        f32 random_f32 = 1024.0f * (static_cast<f32>(rand()) / static_cast<f32>(RAND_MAX));
        test_upper_bound(random_f32);
        test_lower_bound(random_f32);
    }
    return 0;
}