#include "vlog.hpp"
#include <cmath>
#include <cstdio>
#include <cstdint>
#include <vector>

// Checks the 0x1802a24c0 transcription against the reference natural logarithm.
//
// The kernel is a vectorized single-precision ln(x), so the reference is
// std::log (double) evaluated on the same float inputs. The reported metric is
// the max relative error over well-conditioned points (|ln(x)| > 1e-6); for
// inputs where ln(x) ~ 0 (x ~ 1) an absolute error is reported instead.
// A float minimax polynomial reaches ~1-2 ulp, so a 1e-6 relative gate is the
// right tolerance (a double FFT-style 1e-9 gate would be unreachable for float).

static double rel_err(double a, double b) {
    double denom = std::abs(b) > 1e-6 ? std::abs(b) : 1.0;
    return std::abs(a - b) / denom;
}

int main() {
    const uint32_t N = 8192;
    std::vector<float> src(N), dst(N);

    // Sweep a wide dynamic range plus a fine neighbourhood around 1.0.
    uint32_t i = 0;
    for (int k = -30; k <= 30 && i < N; ++k) {
        double v = std::pow(10.0, double(k) * 0.5);
        src[i++] = static_cast<float>(v);
    }
    for (int k = -60; k <= 60 && i < N; ++k) {
        src[i++] = static_cast<float>(1.0 + double(k) * 1e-3);
    }
    for (int k = 0; k < 1024 && i < N; ++k) {
        src[i++] = static_cast<float>(double(k + 1) / 1024.0);
    }
    while (i < N) src[i++] = static_cast<float>(i);

    vlog::log_f32(src.data(), dst.data(), N);

    double max_rel = 0.0, max_abs = 0.0;
    uint32_t rel_i = 0, abs_i = 0;
    for (uint32_t j = 0; j < N; ++j) {
        double ref  = std::log(static_cast<double>(src[j]));
        double mine = static_cast<double>(dst[j]);
        double e    = std::abs(mine - ref);
        if (e > max_abs) { max_abs = e; abs_i = j; }
        double r = rel_err(mine, ref);
        if (r > max_rel) { max_rel = r; rel_i = j; }
    }

    std::printf("n = %u\n", N);
    std::printf("max relative error = %.6e  (at src=%.9g, got %.12g, ref %.12g)\n",
                max_rel, src[rel_i], dst[rel_i], std::log(static_cast<double>(src[rel_i])));
    std::printf("max absolute error = %.6e  (at src=%.9g)\n", max_abs, src[abs_i]);

    if (max_rel < 1e-6) {
        std::printf("ALL OK\n");
        return 0;
    }
    std::printf("FAILED\n");
    return 1;
}