/* Phase 2c: fp8 e4m3 decode + 128x128 block-scale dequant + matmul, validated vs Python.
 * Build: gcc -O2 -o phase2c_fp8 phase2c_fp8.c -lm
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>

#define R 2048
#define C 6144
#define SR 16   /* scale rows = ceil(R/128) */
#define SC 48   /* scale cols = ceil(C/128) */

/* fp8 e4m3fn: 1 sign, 4 exp (bias 7), 3 mant. exp==15&&mant==7 = NaN. */
static float e4m3_to_f32(uint8_t b) {
    int sign = (b >> 7) & 1, exp = (b >> 3) & 0xF, mant = b & 0x7;
    float v;
    if (exp == 0xF && mant == 0x7) return sign ? -NAN : NAN;
    if (exp == 0) v = ldexpf((float)mant / 8.0f, -6);          /* subnormal */
    else          v = ldexpf(1.0f + (float)mant / 8.0f, exp - 7);/* normal */
    return sign ? -v : v;
}

static void* load(const char* path, long bytes) {
    FILE* f = fopen(path, "rb"); if (!f) { perror(path); exit(1); }
    void* p = malloc(bytes);
    if (fread(p, 1, bytes, f) != (size_t)bytes) { fprintf(stderr, "short %s\n", path); exit(1); }
    fclose(f); return p;
}

int main(void) {
    const char* D = "ref"; char p[512];
    snprintf(p, sizeof p, "%s/expert_fp8.u8", D);     uint8_t* w8   = load(p, (long)R * C);
    snprintf(p, sizeof p, "%s/expert_scale.f32", D);  float*   sc   = load(p, (long)SR * SC * 4);
    snprintf(p, sizeof p, "%s/expert_dequant.f32", D);float*   dref = load(p, (long)R * C * 4);
    snprintf(p, sizeof p, "%s/matmul_x.f32", D);      float*   x    = load(p, (long)C * 4);
    snprintf(p, sizeof p, "%s/matmul_y.f32", D);      float*   yref = load(p, (long)R * 4);

    /* dequant + validate against Python dequant */
    float* w = malloc((long)R * C * sizeof(float));
    double dmax = 0;
    for (int i = 0; i < R; i++)
        for (int j = 0; j < C; j++) {
            float v = e4m3_to_f32(w8[(long)i * C + j]) * sc[(i / 128) * SC + (j / 128)];
            w[(long)i * C + j] = v;
            double d = fabs((double)v - dref[(long)i * C + j]);
            if (d > dmax) dmax = d;
        }
    printf("dequant max abs diff vs Python: %.3e\n", dmax);

    /* matmul y = w @ x, validate */
    float* y = malloc(R * sizeof(float));
    double ymax = 0;
    for (int i = 0; i < R; i++) {
        double acc = 0;
        const float* wr = w + (long)i * C;
        for (int j = 0; j < C; j++) acc += (double)wr[j] * x[j];
        y[i] = (float)acc;
        double d = fabs(acc - yref[i]);
        if (d > ymax) ymax = d;
    }
    printf("matmul  y[0..8]: "); for (int i = 0; i < 8; i++) printf("%.4f ", y[i]);
    printf("\nmatmul max abs diff vs Python: %.3e\n", ymax);
    int pass = (dmax < 1e-4) && (ymax < 1e-3);
    printf("VERDICT: %s\n", pass ? "PASS — fp8 dequant + matmul match Python" : "FAIL");
    return pass ? 0 : 1;
}