#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;
}

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CiNpbmNsdWRlIDxzdGRsaWIuaD4KI2luY2x1ZGUgPHN0ZGludC5oPgoKI2RlZmluZSBSQU5HRSAyNTYKdWludDMyX3QgbG9nMl9sdXRbUkFOR0VdOyAKCnVpbnQxNl90IGZhc3RfZ2FtbWFfMTZiKHVpbnQzMl90KiBsMl9sdXQsIHVpbnQzMl90IGdfcTI0LCBpbnQgbikgewogICAgaWYgKG4gPT0gMCkgcmV0dXJuIDA7CiAgICAKICAgIHVpbnQzMl90IHBfcTI0ID0gKCh1aW50NjRfdClsMl9sdXRbbl0gKiBnX3EyNCkgPj4gMjQ7CiAgICBpbnQzMl90IGNfZml4ZWQgPSAtKGludDMyX3QpcF9xMjQ7IAogICAgaW50IGkgPSBjX2ZpeGVkID4+IDI0OwogICAgdWludDMyX3QgZl8zMSA9IChjX2ZpeGVkICYgMHhGRkZGRkYpIDw8IDc7CgogICAgY29uc3QgdWludDMyX3QgY1s1XSA9IHszNjc2OTM5LCA2OTM3NzYyLCAxNjIxNjA4MiwgMjMyNTIwODUsIDE2Nzc3Mjc4fTsKICAgIHVpbnQ2NF90IHJlcyA9IGNbMF07CiAgICBmb3IgKGludCBrID0gMTsgayA8IDU7IGsrKykgewogICAgICAgIHJlcyA9ICgocmVzICogZl8zMSkgPj4gMzIpICsgY1trXTsKICAgIH0KCiAgICBpbnQgc2hpZnQgPSA4IC0gaTsKICAgIHVpbnQzMl90IHJlc3VsdCA9IDA7CiAgICBpZiAoc2hpZnQgPCAzMikgewogICAgICAgIHVpbnQzMl90IHJvdW5kID0gKHNoaWZ0ID09IDApID8gMCA6ICgxdSA8PCAoc2hpZnQgLSAxKSk7CiAgICAgICAgcmVzdWx0ID0gKHJlcyArIHJvdW5kKSA+PiBzaGlmdDsKICAgIH0KICAgIHJldHVybiAocmVzdWx0ID4gNjU1MzUpID8gNjU1MzUgOiAodWludDE2X3QpcmVzdWx0Owp9Cgp2b2lkIGNhbGNfbG9nMl9sdXQodWludDMyX3QqIGx1dCwgaW50IHJhbmdlKSB7CiAgICBsdXRbMF0gPSAwOwogICAgZm9yIChpbnQgbiA9IDE7IG4gPCByYW5nZTsgbisrKSB7CiAgICAgICAgbHV0W25dID0gKHVpbnQzMl90KSgtbG9nMihuIC8gKChyYW5nZSAtIDEpICogMS4wKSkgKiAxNjc3NzIxNi4wICsgMC41KTsKICAgIH0KfQoKaW50IG1haW4oKSB7CiAgICBjb25zdCBmbG9hdCBnYW1tYSA9IDIuMmY7CiAgICBjYWxjX2xvZzJfbHV0KGxvZzJfbHV0LCBSQU5HRSk7CgogICAgaW50IG1heF9lcnIgPSAwOwogICAgdWludDMyX3QgZ19xMjQgPSAodWludDMyX3QpKGdhbW1hICogMTY3NzcyMTYuMCk7CgogICAgZm9yIChpbnQgbiA9IDA7IG4gPCBSQU5HRTsgbisrKSB7CiAgICAgICAgdWludDE2X3QgZmFzdF9pbnQgPSBmYXN0X2dhbW1hXzE2Yihsb2cyX2x1dCwgZ19xMjQsIG4pOwogICAgICAgIHVpbnQxNl90IHN0YW5kYXJkID0gKHVpbnQxNl90KShwb3coKGRvdWJsZSluIC8gKFJBTkdFIC0gMSksIGdhbW1hKSAqIDY1NTM1LjAgKyAwLjUpOwogICAgICAgIGludCBlcnIgPSAoaW50KWZhc3RfaW50IC0gKGludClzdGFuZGFyZDsKICAgICAgICBpZiAoYWJzKGVycikgPiBhYnMobWF4X2VycikpIG1heF9lcnIgPSBlcnI7CiAgICAgICAgCiAgICAgICAgaWYgKGVyciAhPSAwIHx8IG4gPT0gMCB8fCBuID09IFJBTkdFIC0gMSkgewogICAgICAgICAgICBwcmludGYoIiU0ZCB8ICU2ZCB8ICU2ZCB8ICVkXG4iLCBuLCBzdGFuZGFyZCwgZmFzdF9pbnQsIGVycik7CiAgICAgICAgfQogICAgfQogICAgcHJpbnRmKCI+PiBNYXggRXJyb3I6ICVkXG4iLCBtYXhfZXJyKTsKICAgIHJldHVybiAwOwp9