Floating Point Error Accumulation Demo

floating_point_err_accum.cpp

#include <stdio.h>
#include <cmath>

void 
//__attribute__ ((optimize("-O0")))
CoordsAccumulative(float xMin, float xMax, float yMin, float yMax, int N, float* result)
{
    float PerY = (yMax - yMin) / N;
    float PerX = N * PerY;;

    float ixPerX = 0.0f;

    for (int i = 0; i < N; i++)
    {
    float iyPerY = 0.0f;
        for (int j = 0; j < N; j++)
        {
            result[i * N + j] *= ixPerX + iyPerY;
            iyPerY += PerY;
            //printf("%f %f\r\n", ixPerX, iyPerY);
        }
        ixPerX += PerX;
    }
}

void
//__attribute__ ((optimize("-O0")))
CoordsDirect(float xMin, float xMax, float yMin, float yMax, int N, float* result)
{
    float PerY = (yMax - yMin) / N;
    float PerX = N * PerY;;

    float ixPerX = 0.0f;
    float iyPerY = 0.0f;

    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < N; j++)
        {
            result[i * N + j] *= i * PerX + j * PerY;
            //printf("%f %f\r\n", i * PerX, j * PerY);
        }
    }
}

int main()
{
    int N = 1024;
    float xMin = -1e-5;
    float xMax = 1e-5;
    float yMin = -1e-5;
    float yMax = 1e-5;

    float* accum_result = new float[N*N];
    float* direct_result = new float[N*N];

    for (int i = 0; i < N * N; i++)
    {
        direct_result[i] = accum_result[i] = (float)rand() / RAND_MAX * 1e5;
    }

    CoordsAccumulative(xMin, xMax, yMin, yMax, N, accum_result);
    CoordsDirect(xMin, xMax, yMin, yMax, N, direct_result);

    //return 0;
    for (int j = 0; j < N; j++)
    {
        int i = j * N;
        //printf("%.8f %.8f", accum_result[i], direct_result[i]);
        printf("\t%.8f\r\n", fabs(accum_result[i] - direct_result[i]));
    }
}