// delaybench.cpp - Limenka fork delay-function benchmark
//
// Candidates for the per-block sequential delay (D = 60s, CPU-bound,
// low hardware variance, non-shortcuttable, verified by recomputation):
//
//  a1. M=1 schoolbook long division, hardware 128/64 DIV per step
//  a2. M=1 schoolbook long division, Barrett reciprocal estimate
//      (division-free: mul/add/sub/compare only - the uniform variant)
//  b1. M=4 Knuth D, hardware-DIV qhat
//  b2. M=4 Knuth D, Barrett reciprocal qhat (division-free)
//  c.  quadratic map x = x*x + c mod 2^128 (record only; shortcut
//      status unresolved)
//  d.  SHA256 single chain (hardware-variance control baseline)
//
// The dividend stream comes from xorshift64* seeded from the block
// input: the work binds to the block, cannot be precomputed, and the
// working set is a few words regardless of stream length.
//
// Build: g++ -O2 -march=native -o delaybench delaybench.cpp
#include <cstdint>
#include <cstdio>
#include <cstring>
#include <chrono>
#include <cstdlib>

using u64 = uint64_t;
using u128 = unsigned __int128;
using u32 = uint32_t;
using u8 = uint8_t;

static u64 prng_state;

static inline u64 prng_next() {
    u64 x = prng_state;
    x ^= x >> 12;
    x ^= x << 25;
    x ^= x >> 27;
    prng_state = x;
    return x * 0x2545F4914F6CDD1Dull;
}

// ---------------- a1. M=1, hardware 128/64 DIV ----------------
static u64 run_div1(u64 d, u64 words, bool verbose) {
    u64 r = 0;
    for (u64 i = 0; i < words; i++) {
        u64 w = prng_next();
        r = (u64)((((u128)r << 64) | w) % d);
        if (verbose)
            printf("  div1 step %llu: r=%016llx\n", (unsigned long long)i, (unsigned long long)r);
    }
    return r;
}

// high 128 bits of (u1:u0) * (v1:v0) - GCC truncates u128*u128 to the
// low 128 bits, so the schoolbook is done explicitly
static inline u64 mulhi128(u64 u1, u64 u0, u64 v1, u64 v0) {
    u64 hi0 = (u64)((u128)u0 * v0 >> 64);
    u128 t1 = (u128)u0 * v1 + hi0;
    u64 lo1 = (u64)t1, hi1 = (u64)(t1 >> 64);
    u128 t2 = (u128)u1 * v0 + lo1;
    u64 lo2 = (u64)t2, hi2 = (u64)(t2 >> 64);
    (void)lo2;
    return (u64)((u128)u1 * v1 + hi1 + hi2);
}

// ---------------- a2. M=1, Barrett reciprocal (division-free) ----------------
// v = floor((2^128 - 1) / d), precomputed once per block.
// qhat = (u * v) >> 128 for u = (r << 64 | w) < d*2^64; estimate is
// exact or 1 too small; refine with one conditional subtract.
static inline u64 barrett_step(u64 r, u64 w, u64 d, u128 v) {
    u128 u = ((u128)r << 64) | w;
    u64 q = mulhi128(r, w, (u64)(v >> 64), (u64)v);
    u128 rem = u - (u128)q * d;
    if (rem >= d) {
        rem -= d;
        q++;
    }
    (void)q;
    return (u64)rem;
}

static u64 run_div1_barrett(u64 d, u64 words, bool verbose) {
    u128 v = (((u128)0 - 1) << 64 | (u128)0xffffffffffffffffULL) / d; // floor((2^128-1)/d)
    u64 r = 0;
    for (u64 i = 0; i < words; i++) {
        u64 w = prng_next();
        r = barrett_step(r, w, d, v);
        if (verbose)
            printf("  bar1 step %llu: r=%016llx\n", (unsigned long long)i, (unsigned long long)r);
    }
    return r;
}

// ---------------- b1. M=4 Knuth D, hardware-DIV qhat ----------------
static inline void knuth_step_div(u64* u, const u64* d, u64 w) {
    u64 q = (u64)((((u128)u[3] << 64) | u[2]) / d[3]);
    u64 r[5] = {w, u[0], u[1], u[2], u[3]};
    bool neg = false;
    for (int i = 0; i < 4; i++) {
        u128 prod = (u128)q * d[i];
        u64 lo = (u64)prod, hi = (u64)(prod >> 64);
        u64 t = r[i] - lo;
        u64 b0 = (r[i] < lo);
        r[i] = t;
        // r[i+1] -= hi + b0, split to avoid the hi+b0 wrap when
        // hi == 2^64-1 and b0 == 1
        t = r[i + 1] - b0;
        u64 b1 = (r[i + 1] < b0);
        r[i + 1] = t;
        t = r[i + 1] - hi;
        b1 |= (r[i + 1] < hi);
        r[i + 1] = t;
        for (int j = i + 2; j < 5 && b1; j++) {
            t = r[j] - 1;
            b1 = (r[j] < 1);
            r[j] = t;
        }
        if (b1) neg = true;  // borrow escaped the 5-word value
    }
    // negative test: top bit of the 5-word value (two's complement).
    // The add-back carry must propagate into r[4].
    while ((r[4] >> 63) != 0) {
        u64 carry = 0;
        for (int i = 0; i < 4; i++) {
            u64 s = r[i] + d[i];
            u64 c1 = (s < r[i]);
            u64 s2 = s + carry;
            u64 c2 = (s2 < s);
            r[i] = s2;
            carry = c1 | c2;
        }
        r[4] += carry;
    }
    for (int rep = 0; rep < 4; rep++) {
        // r >= d? (equality also subtracts once)
        int ge = 1;
        for (int i = 3; i >= 0; i--) {
            if (r[i] > d[i]) break;
            if (r[i] < d[i]) { ge = 0; break; }
        }
        if (!ge) break;
        u64 borrow = 0;
        for (int i = 0; i < 4; i++) {
            u64 t = r[i] - d[i];
            u64 b1 = (r[i] < d[i]);
            u64 t2 = t - borrow;
            u64 b2 = (t < borrow);
            r[i] = t2;
            borrow = b1 | b2;
        }
    }
    u[0] = r[0]; u[1] = r[1]; u[2] = r[2]; u[3] = r[3];
}

static u64 run_knuth4_div(const u64* d, u64 words, bool verbose) {
    u64 u[4] = {0, 0, 0, 0};
    for (u64 i = 0; i < words; i++) {
        u64 w = prng_next();
        knuth_step_div(u, d, w);
        if (verbose)
            printf("  knuthdiv step %llu: u=%016llx %016llx %016llx %016llx\n",
                   (unsigned long long)i,
                   (unsigned long long)u[0], (unsigned long long)u[1],
                   (unsigned long long)u[2], (unsigned long long)u[3]);
    }
    return u[0] ^ u[1] ^ u[2] ^ u[3];
}

// ---------------- b2. M=4 Knuth D, Barrett reciprocal qhat ----------------
// v3 = floor((2^128 - 1) / d3), precomputed per block. qhat =
// ((u3:u2) * v3) >> 128, exact or off by 1; refined by two compares.
static u64 run_knuth4_barrett(const u64* d, u64 words, bool verbose) {
    u64 u[4] = {0, 0, 0, 0};
    u128 v3 = (((u128)0 - 1) << 64 | (u128)0xffffffffffffffffULL) / d[3];
    for (u64 i = 0; i < words; i++) {
        u64 w = prng_next();
        u64 q = mulhi128(u[3], u[2], (u64)(v3 >> 64), (u64)v3);
        // refine: while (u3:u2) < q*d3, q--; while >= (q+1)*d3, q++
        {
            u128 top = (u128)u[3] << 64 | u[2];
            if ((u128)q * d[3] > top) q--;
            else if (top >= ((u128)q + 1) * d[3]) q++;
        }
        u64 r[5] = {w, u[0], u[1], u[2], u[3]};
        bool neg = false;
        for (int i = 0; i < 4; i++) {
            u128 prod = (u128)q * d[i];
            u64 lo = (u64)prod, hi = (u64)(prod >> 64);
            u64 t = r[i] - lo;
            u64 b0 = (r[i] < lo);
            r[i] = t;
            // r[i+1] -= hi + b0, split to avoid the hi+b0 wrap
            t = r[i + 1] - b0;
            u64 b1 = (r[i + 1] < b0);
            r[i + 1] = t;
            t = r[i + 1] - hi;
            b1 |= (r[i + 1] < hi);
            r[i + 1] = t;
            for (int j = i + 2; j < 5 && b1; j++) {
                t = r[j] - 1;
                b1 = (r[j] < 1);
                r[j] = t;
            }
            if (b1) neg = true;
        }
        // negative test: top bit of the 5-word value (two's
        // complement). The add-back carry must propagate into r[4].
        while ((r[4] >> 63) != 0) {
            u64 carry = 0;
            for (int i = 0; i < 4; i++) {
                u64 s = r[i] + d[i];
                u64 c1 = (s < r[i]);
                u64 s2 = s + carry;
                u64 c2 = (s2 < s);
                r[i] = s2;
                carry = c1 | c2;
            }
            r[4] += carry;
        }
        for (int rep = 0; rep < 4; rep++) {
            int ge = 1;
            for (int i = 3; i >= 0; i--) {
                if (r[i] > d[i]) break;
                if (r[i] < d[i]) { ge = 0; break; }
            }
            if (!ge) break;
            u64 borrow = 0;
            for (int i = 0; i < 4; i++) {
                u64 t = r[i] - d[i];
                u64 b1 = (r[i] < d[i]);
                u64 t2 = t - borrow;
                u64 b2 = (t < borrow);
                r[i] = t2;
                borrow = b1 | b2;
            }
        }
        u[0] = r[0]; u[1] = r[1]; u[2] = r[2]; u[3] = r[3];
        if (verbose)
            printf("  knuthbar step %llu: u=%016llx %016llx %016llx %016llx\n",
                   (unsigned long long)i,
                   (unsigned long long)u[0], (unsigned long long)u[1],
                   (unsigned long long)u[2], (unsigned long long)u[3]);
    }
    return u[0] ^ u[1] ^ u[2] ^ u[3];
}

// ---------------- c. quadratic map x = x*x + c mod 2^128 ----------------
static u64 run_quad128(u64 c, u64 words) {
    u128 x = c;
    for (u64 i = 0; i < words; i++) {
        x = x * x + c;
    }
    return (u64)(x >> 64) ^ (u64)x;
}

// ---------------- d. SHA256 chain (scalar reference) ----------------
static inline u32 rotr32(u32 x, int n) { return (x >> n) | (x << (32 - n)); }

static void sha256_compress(u32* h, const u8* block) {
    static const u32 K[64] = {
        0x428a2f98,0x71374491,0xb5c0fbcf,0xe9b5dba5,0x3956c25b,0x59f111f1,0x923f82a4,0xab1c5ed5,
        0xd807aa98,0x12835b01,0x243185be,0x550c7dc3,0x72be5d74,0x80deb1fe,0x9bdc06a7,0xc19bf174,
        0xe49b69c1,0xefbe4786,0x0fc19dc6,0x240ca1cc,0x2de92c6f,0x4a7484aa,0x5cb0a9dc,0x76f988da,
        0x983e5152,0xa831c66d,0xb00327c8,0xbf597fc7,0xc6e00bf3,0xd5a79147,0x06ca6351,0x14292967,
        0x27b70a85,0x2e1b2138,0x4d2c6dfc,0x53380d13,0x650a7354,0x766a0abb,0x81c2c92e,0x92722c85,
        0xa2bfe8a1,0xa81a664b,0xc24b8b70,0xc76c51a3,0xd192e819,0xd6990624,0xf40e3585,0x106aa070,
        0x19a4c116,0x1e376c08,0x2748774c,0x34b0bcb5,0x391c0cb3,0x4ed8aa4a,0x5b9cca4f,0x682e6ff3,
        0x748f82ee,0x78a5636f,0x84c87814,0x8cc70208,0x90befffa,0xa4506ceb,0xbef9a3f7,0xc67178f2,
    };
    u32 w[64];
    for (int i = 0; i < 16; i++)
        w[i] = (block[i * 4] << 24) | (block[i * 4 + 1] << 16) | (block[i * 4 + 2] << 8) | block[i * 4 + 3];
    for (int i = 16; i < 64; i++) {
        u32 s0 = rotr32(w[i - 15], 7) ^ rotr32(w[i - 15], 18) ^ (w[i - 15] >> 3);
        u32 s1 = rotr32(w[i - 2], 17) ^ rotr32(w[i - 2], 19) ^ (w[i - 2] >> 10);
        w[i] = w[i - 16] + s0 + w[i - 7] + s1;
    }
    u32 a = h[0], b = h[1], c2 = h[2], d2 = h[3], e = h[4], f = h[5], g = h[6], hh = h[7];
    for (int i = 0; i < 64; i++) {
        u32 S1 = rotr32(e, 6) ^ rotr32(e, 11) ^ rotr32(e, 25);
        u32 ch = (e & f) ^ (~e & g);
        u32 t1 = hh + S1 + ch + K[i] + w[i];
        u32 S0 = rotr32(a, 2) ^ rotr32(a, 13) ^ rotr32(a, 22);
        u32 maj = (a & b) ^ (a & c2) ^ (b & c2);
        u32 t2 = S0 + maj;
        hh = g; g = f; f = e; e = d2 + t1;
        d2 = c2; c2 = b; b = a; a = t1 + t2;
    }
    h[0] += a; h[1] += b; h[2] += c2; h[3] += d2;
    h[4] += e; h[5] += f; h[6] += g; h[7] += hh;
}

static u64 run_sha256chain(u64 words) {
    u32 h[8] = {0x6a09e667,0xbb67ae85,0x3c6ef372,0xa54ff53a,0x510e527f,0x9b05688c,0x1f83d9ab,0x5be0cd19};
    u8 block[64] = {0};
    for (u64 i = 0; i < words; i++) {
        u64 w0 = prng_next();
        memcpy(block, &w0, 8);
        sha256_compress(h, block);
    }
    return h[0] ^ h[7];
}

int main(int argc, char** argv) {
    u64 words = argc > 1 ? strtoull(argv[1], nullptr, 10) : (1ull << 24);
    bool verbose = argc > 2 && strcmp(argv[2], "--verbose") == 0;
    prng_state = 0x9e3779b97f4a7c15ull;

    u64 d1 = 0x9e3779b97f4a7c15ull | 1;
    u64 d4[4] = {0x243f6a8885a308d3 | 1, 0x13198a2e03707344, 0xa4093822299f31d0, 0x082efa98ec4e6c89 | (1ull << 63)};

    double secs[6];
    u64 res[6];
    const char* names[6] = {
        "M=1 div hw    ", "M=1 barrett   ", "M=4 knuth div ",
        "M=4 knuth barr", "x^2+c mod 2^128", "sha256 chain  ",
    };
    const u64 seed = 0x9e3779b97f4a7c15ull;
    auto run = [&](int idx, u64 r) { res[idx] = r; };

    prng_state = seed;
    auto t0 = std::chrono::steady_clock::now();
    run(0, run_div1(d1, words, verbose));
    auto t1 = std::chrono::steady_clock::now();
    prng_state = seed;
    run(1, run_div1_barrett(d1, words, verbose));
    auto t2 = std::chrono::steady_clock::now();
    prng_state = seed;
    run(2, run_knuth4_div(d4, words, verbose));
    auto t3 = std::chrono::steady_clock::now();
    prng_state = seed;
    run(3, run_knuth4_barrett(d4, words, verbose));
    auto t4 = std::chrono::steady_clock::now();
    prng_state = seed;
    run(4, run_quad128(d1, words));
    auto t5 = std::chrono::steady_clock::now();
    prng_state = seed;
    run(5, run_sha256chain(words));
    auto t6 = std::chrono::steady_clock::now();

    secs[0] = std::chrono::duration<double>(t1 - t0).count();
    secs[1] = std::chrono::duration<double>(t2 - t1).count();
    secs[2] = std::chrono::duration<double>(t3 - t2).count();
    secs[3] = std::chrono::duration<double>(t4 - t3).count();
    secs[4] = std::chrono::duration<double>(t5 - t4).count();
    secs[5] = std::chrono::duration<double>(t6 - t5).count();

    printf("words %llu\n", (unsigned long long)words);
    for (int i = 0; i < 6; i++)
        printf("%s: %8.3f s  %6.2f ns/step  (%llx)\n", names[i], secs[i],
               secs[i] / words * 1e9, (unsigned long long)res[i]);
    return 0;
}