#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <stdint.h>

#define RANGE 256
uint32_t log2_lut[RANGE]; 

uint16_t fast_gamma_16b(uint32_t* l2_lut, uint32_t g_q24, int n) {
    if (n == 0) return 0;
    
    uint32_t p_q24 = ((uint64_t)l2_lut[n] * g_q24) >> 24;
    int32_t c_fixed = -(int32_t)p_q24; 
    int i = c_fixed >> 24;
    uint32_t f_31 = (c_fixed & 0xFFFFFF) << 7;

    const uint32_t c[5] = {3676939, 6937762, 16216082, 23252085, 16777278};
    uint64_t res = c[0];
    for (int k = 1; k < 5; k++) {
        res = ((res * f_31) >> 32) + c[k];
    }

    int shift = 8 - i;
    uint32_t result = 0;
    if (shift < 32) {
        uint32_t round = (shift == 0) ? 0 : (1u << (shift - 1));
        result = (res + round) >> shift;
    }
    return (result > 65535) ? 65535 : (uint16_t)result;
}

void calc_log2_lut(uint32_t* lut, int range) {
    lut[0] = 0;
    for (int n = 1; n < range; n++) {
        lut[n] = (uint32_t)(-log2(n / ((range - 1) * 1.0)) * 16777216.0 + 0.5);
    }
}

int main() {
    const float gamma = 2.2f;
    calc_log2_lut(log2_lut, RANGE);

    int max_err = 0;
    uint32_t g_q24 = (uint32_t)(gamma * 16777216.0);

    for (int n = 0; n < RANGE; n++) {
        uint16_t fast_int = fast_gamma_16b(log2_lut, g_q24, n);
        uint16_t standard = (uint16_t)(pow((double)n / (RANGE - 1), gamma) * 65535.0 + 0.5);
        int err = (int)fast_int - (int)standard;
        if (abs(err) > abs(max_err)) max_err = err;
        
        if (err != 0 || n == 0 || n == RANGE - 1) {
            printf("%4d | %6d | %6d | %d\n", n, standard, fast_int, err);
        }
    }
    printf(">> Max Error: %d\n", max_err);
    return 0;
}