C++ code to compute the inverse of a linearly-interpolated 1D table with values in [0, 1] x [0, 1] (note: the table size must be a power of two)

Inverse.cpp

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <vector>

// inverting piecewise linear functions
float InverseLinear(float yStart, float yEnd, float y)
{
    return (y - yStart) / (yEnd - yStart);
}

float InverseLinearInterval(float xStart, float xEnd,
                            float yStart, float yEnd,
                            float y)
{
    float scale = xEnd - xStart;
    float offset = xStart;

    return InverseLinear(yStart, yEnd, y) * scale + offset;
}

// inverts a function defined with values in [0, 1]x[0, 1]
std::vector<float> Inverse(const std::vector<float>& cdf, int sizeLog2)
{
    int size = (1 << sizeLog2);
    std::vector<float> qf(size, 0);

    for (int i = 0; i < size; ++i) {
        float u = (float)i / (size - 1); // in [0, 1]
        int index = size / 2;

        // binary search the value
        for (int j = sizeLog2 - 2; j >= 0; --j) {
            int offset = 1 << j;

            if (u <= cdf[index]) {
                index-= offset;
            } else {
                index+= offset;
            }
        }

        // invert linear segment
        if (u <= cdf[index]) {
            float xStart = (float)(index - 1) / (size - 1);
            float xEnd   = (float)(index    ) / (size - 1);

            qf[i] = InverseLinearInterval(xStart, xEnd, cdf[index - 1], cdf[index], u);
        } else {
            float xStart = (float)(index + 1) / (size - 1);
            float xEnd   = (float)(index    ) / (size - 1);

            qf[i] = InverseLinearInterval(xStart, xEnd, cdf[index + 1], cdf[index], u);
        }
    }

    return qf;
}

#if 0 // for experimenting
// dumps a plot file for GNUPLOT 
// e.g., > plot './linear' u 1:2 w lines, '' u 1:3 w lines, '' u 1:4 w lines
void print(const char *filename, const std::vector<float>& cdf, int sizeLog2)
{
    FILE *pf = fopen(filename, "w");
    int size = 1 << sizeLog2;
    std::vector<float> qf = Inverse(cdf, sizeLog2);

    for (int i = 0; i < size; ++i) {
        float u = (float)i / (size - 1);

        fprintf(pf, "%f %f %f %f \n", u, cdf[i], qf[i], u);
    }

    fclose(pf);
}


int main(int argc, char **argv)
{
    int sizeLog2 = 8;
    int size = 1 << sizeLog2;
    std::vector<float> cdf(size);

    // linear
    for (int i = 0; i < size; ++i) {
        float u = (float)i / (size - 1);

        cdf[i] = u;
    }
    print("linear", cdf, sizeLog2);

    // square
    for (int i = 0; i < size; ++i) {
        float u = (float)i / (size - 1);

        cdf[i] = u * u;
    }
    print("square", cdf, sizeLog2);

    // quintic
    for (int i = 0; i < size; ++i) {
        float u = (float)i / (size - 1);

        cdf[i] = u * u * u * u * u;
    }
    print("quintic", cdf, sizeLog2);
    
    // sqrt
    for (int i = 0; i < size; ++i) {
        float u = (float)i / (size - 1);

        cdf[i] = sqrtf(u);
    }
    print("sqrt", cdf, sizeLog2);
}
#endif