/* Copyright 2026 Rhett Creighton - Apache License 2.0
 *
 * BLS12-381 pairing — pure C23 implementation.
 * 381-bit prime, 6 x 64-bit limbs, Montgomery multiplication.
 * Full tower: Fp -> Fp2 -> Fp6 -> Fp12, G1, G2, optimal Ate pairing. */

#include "sapling/bls12_381.h"
#include <string.h>
#include <stdlib.h>
#include "util/safe_alloc.h"
#include "util/log_macros.h"

/* Runtime-dispatched Montgomery multiply (CPUID → BMI2+ADX or portable) */
extern void fp_mont_mul_accel(uint64_t r[6], const uint64_t a[6], const uint64_t b[6]);
#define fp_mont_mul fp_mont_mul_accel

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wpedantic"

/* ===== Fp constants ===== */

/* q = 0x1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffaaab */
static const uint64_t FP_Q[6] = {
    0xb9feffffffffaaabULL, 0x1eabfffeb153ffffULL,
    0x6730d2a0f6b0f624ULL, 0x64774b84f38512bfULL,
    0x4b1ba7b6434bacd7ULL, 0x1a0111ea397fe69aULL
};

/* FP_INV moved to fr_avx512.c (runtime dispatch) */

/* R = 2^384 mod q */
static const uint64_t FP_R[6] = {
    0x760900000002fffdULL, 0xebf4000bc40c0002ULL,
    0x5f48985753c758baULL, 0x77ce585370525745ULL,
    0x5c071a97a256ec6dULL, 0x15f65ec3fa80e493ULL
};

/* R^2 mod q */
static const uint64_t FP_R2[6] = {
    0xf4df1f341c341746ULL, 0x0a76e6a609d104f1ULL,
    0x8de5476c4c95b6d5ULL, 0x67eb88a9939d83c0ULL,
    0x9a793e85b519952dULL, 0x11988fe592cae3aaULL
};

/* ===== Fp arithmetic ===== */

static bool fp_gte(const uint64_t a[6], const uint64_t b[6])
{
    for (int i = 5; i >= 0; i--) {
        if (a[i] > b[i]) return true;
        if (a[i] < b[i]) return false;
    }
    return true;
}

static void fp_sub_noborrow(uint64_t r[6], const uint64_t a[6], const uint64_t b[6])
{
    unsigned __int128 borrow = 0;
    for (int i = 0; i < 6; i++) {
        unsigned __int128 tmp = (unsigned __int128)a[i] - b[i] - borrow;
        r[i] = (uint64_t)tmp;
        borrow = (tmp >> 127) & 1;
    }
}

/* Montgomery multiply implementation is in fr_avx512.c (runtime dispatch) */

void fp_zero(struct fp *r) { memset(r->d, 0, 48); }

void fp_one(struct fp *r) { memcpy(r->d, FP_R, 48); }

bool fp_is_zero(const struct fp *a)
{
    for (int i = 0; i < 6; i++)
        if (a->d[i] != 0) return false;
    return true;
}

bool fp_eq(const struct fp *a, const struct fp *b)
{
    for (int i = 0; i < 6; i++)
        if (a->d[i] != b->d[i]) return false;
    return true;
}

void fp_add(struct fp *r, const struct fp *a, const struct fp *b)
{
    unsigned __int128 carry = 0;
    uint64_t tmp[6];
    for (int i = 0; i < 6; i++) {
        unsigned __int128 sum = (unsigned __int128)a->d[i] + b->d[i] + carry;
        tmp[i] = (uint64_t)sum;
        carry = sum >> 64;
    }
    if (carry || fp_gte(tmp, FP_Q))
        fp_sub_noborrow(r->d, tmp, FP_Q);
    else
        memcpy(r->d, tmp, 48);
}

void fp_sub(struct fp *r, const struct fp *a, const struct fp *b)
{
    if (fp_gte(a->d, b->d)) {
        fp_sub_noborrow(r->d, a->d, b->d);
    } else {
        uint64_t tmp[6];
        unsigned __int128 carry = 0;
        for (int i = 0; i < 6; i++) {
            unsigned __int128 sum = (unsigned __int128)a->d[i] + FP_Q[i] + carry;
            tmp[i] = (uint64_t)sum;
            carry = sum >> 64;
        }
        fp_sub_noborrow(r->d, tmp, b->d);
    }
}

void fp_neg(struct fp *r, const struct fp *a)
{
    if (fp_is_zero(a))
        fp_zero(r);
    else
        fp_sub_noborrow(r->d, FP_Q, a->d);
}

void fp_mul(struct fp *r, const struct fp *a, const struct fp *b)
{
    fp_mont_mul(r->d, a->d, b->d);
}

void fp_sq(struct fp *r, const struct fp *a)
{
    fp_mont_mul(r->d, a->d, a->d);
}

/* fp_pow: a^exp where exp is an array of 6 uint64_t (little-endian limbs) */
static void fp_pow(struct fp *r, const struct fp *a, const uint64_t exp[6])
{
    struct fp result;
    fp_one(&result);
    struct fp base = *a;

    for (int i = 0; i < 6; i++) {
        for (int bit = 0; bit < 64; bit++) {
            if ((exp[i] >> bit) & 1)
                fp_mul(&result, &result, &base);
            fp_sq(&base, &base);
        }
    }
    *r = result;
}

void fp_inv(struct fp *r, const struct fp *a)
{
    /* q - 2 */
    uint64_t exp[6];
    memcpy(exp, FP_Q, 48);
    /* Subtract 2 from little-endian */
    exp[0] -= 2; /* FP_Q[0] ends in ...aaab, so no borrow */
    fp_pow(r, a, exp);
}

bool fp_sqrt(struct fp *r, const struct fp *a)
{
    /* q ≡ 3 (mod 4), so sqrt(a) = a^((q+1)/4) */
    /* (q+1)/4 */
    uint64_t exp[6];
    memcpy(exp, FP_Q, 48);
    /* q+1: add 1 to ...aaab -> ...aaac */
    unsigned __int128 carry = 1;
    for (int i = 0; i < 6; i++) {
        unsigned __int128 s = (unsigned __int128)exp[i] + carry;
        exp[i] = (uint64_t)s;
        carry = s >> 64;
    }
    /* Divide by 4: right shift by 2 */
    for (int i = 0; i < 5; i++)
        exp[i] = (exp[i] >> 2) | (exp[i + 1] << 62);
    exp[5] >>= 2;

    fp_pow(r, a, exp);

    /* Verify: r^2 == a */
    struct fp check;
    fp_sq(&check, r);
    return fp_eq(&check, a);
}

bool fp_lexicographically_largest(const struct fp *a)
{
    /* Convert from Montgomery to get raw value, then check > (q-1)/2 */
    uint64_t one[6] = {1, 0, 0, 0, 0, 0};
    uint64_t raw[6];
    fp_mont_mul(raw, a->d, one);

    /* (q-1)/2 */
    uint64_t half[6];
    for (int i = 0; i < 6; i++) half[i] = FP_Q[i];
    /* Right shift by 1 */
    for (int i = 0; i < 5; i++)
        half[i] = (half[i] >> 1) | (half[i + 1] << 63);
    half[5] >>= 1;

    /* raw > half means lexicographically largest */
    for (int i = 5; i >= 0; i--) {
        if (raw[i] > half[i]) return true;
        if (raw[i] < half[i]) return false;
    }
    return false;
}

bool fp_from_bytes(struct fp *r, const uint8_t s[48])
{
    uint64_t raw[6];
    for (int i = 0; i < 6; i++) {
        raw[5 - i] = 0;
        for (int j = 0; j < 8; j++)
            raw[5 - i] |= (uint64_t)s[i * 8 + j] << (8 * (7 - j));
    }

    bool ok = true;
    for (int i = 5; i >= 0; i--) {
        if (raw[i] > FP_Q[i]) { ok = false; break; }
        if (raw[i] < FP_Q[i]) break;
    }

    fp_mont_mul(r->d, raw, FP_R2);
    return ok;
}

void fp_to_bytes(uint8_t s[48], const struct fp *a)
{
    uint64_t one[6] = {1, 0, 0, 0, 0, 0};
    uint64_t raw[6];
    fp_mont_mul(raw, a->d, one);

    for (int i = 0; i < 6; i++)
        for (int j = 0; j < 8; j++)
            s[i * 8 + j] = (uint8_t)(raw[5 - i] >> (8 * (7 - j)));
}

/* ===== Fp2 = Fp[u] / (u^2 + 1) ===== */

void fp2_zero(struct fp2 *r) { fp_zero(&r->c0); fp_zero(&r->c1); }

void fp2_one(struct fp2 *r) { fp_one(&r->c0); fp_zero(&r->c1); }

bool fp2_is_zero(const struct fp2 *a) { return fp_is_zero(&a->c0) && fp_is_zero(&a->c1); }

bool fp2_eq(const struct fp2 *a, const struct fp2 *b)
{
    return fp_eq(&a->c0, &b->c0) && fp_eq(&a->c1, &b->c1);
}

void fp2_add(struct fp2 *r, const struct fp2 *a, const struct fp2 *b)
{
    fp_add(&r->c0, &a->c0, &b->c0);
    fp_add(&r->c1, &a->c1, &b->c1);
}

void fp2_sub(struct fp2 *r, const struct fp2 *a, const struct fp2 *b)
{
    fp_sub(&r->c0, &a->c0, &b->c0);
    fp_sub(&r->c1, &a->c1, &b->c1);
}

void fp2_neg(struct fp2 *r, const struct fp2 *a)
{
    fp_neg(&r->c0, &a->c0);
    fp_neg(&r->c1, &a->c1);
}

void fp2_mul(struct fp2 *r, const struct fp2 *a, const struct fp2 *b)
{
    struct fp t0, t1, t2, t3, rc0, rc1;
    fp_mul(&t0, &a->c0, &b->c0);
    fp_mul(&t1, &a->c1, &b->c1);
    fp_sub(&rc0, &t0, &t1);
    fp_add(&t2, &a->c0, &a->c1);
    fp_add(&t3, &b->c0, &b->c1);
    fp_mul(&rc1, &t2, &t3);
    fp_sub(&rc1, &rc1, &t0);
    fp_sub(&rc1, &rc1, &t1);
    r->c0 = rc0;
    r->c1 = rc1;
}

void fp2_sq(struct fp2 *r, const struct fp2 *a)
{
    struct fp t0, t1, rc0, rc1;
    fp_add(&t0, &a->c0, &a->c1);
    fp_sub(&t1, &a->c0, &a->c1);
    fp_mul(&rc0, &t0, &t1);
    fp_mul(&t0, &a->c0, &a->c1);
    fp_add(&rc1, &t0, &t0);
    r->c0 = rc0;
    r->c1 = rc1;
}

void fp2_inv(struct fp2 *r, const struct fp2 *a)
{
    struct fp t0, t1, inv;
    fp_sq(&t0, &a->c0);
    fp_sq(&t1, &a->c1);
    fp_add(&t0, &t0, &t1);
    fp_inv(&inv, &t0);
    fp_mul(&r->c0, &a->c0, &inv);
    fp_neg(&r->c1, &a->c1);
    fp_mul(&r->c1, &r->c1, &inv);
}

/* (a0 + a1*u) * (1 + u) = (a0 - a1) + (a0 + a1)*u */
void fp2_mul_by_nonresidue(struct fp2 *r, const struct fp2 *a)
{
    struct fp t0;
    fp_sub(&t0, &a->c0, &a->c1);
    fp_add(&r->c1, &a->c0, &a->c1);
    r->c0 = t0;
}

/* FROBENIUS_COEFF_FQ2_C1[1] = -1 in Montgomery form (q - R mod q) */
static const struct fp FP_NEGATIVE_ONE = {{
    0x43f5fffffffcaaaeULL, 0x32b7fff2ed47fffdULL,
    0x07e83a49a2e99d69ULL, 0xeca8f3318332bb7aULL,
    0xef148d1ea0f4c069ULL, 0x040ab3263eff0206ULL
}};

void fp2_frobenius_map(struct fp2 *a, int power)
{
    if (power % 2 != 0)
        fp_mul(&a->c1, &a->c1, &FP_NEGATIVE_ONE);
}

static void fp2_pow(struct fp2 *r, const struct fp2 *a, const uint64_t exp[6])
{
    struct fp2 result;
    fp2_one(&result);
    struct fp2 base = *a;

    for (int i = 0; i < 6; i++) {
        for (int bit = 0; bit < 64; bit++) {
            if ((exp[i] >> bit) & 1)
                fp2_mul(&result, &result, &base);
            fp2_sq(&base, &base);
        }
    }
    *r = result;
}

bool fp2_sqrt(struct fp2 *r, const struct fp2 *a)
{
    if (fp2_is_zero(a)) {
        fp2_zero(r);
        return true;
    }

    /* Algorithm 9, https://eprint.iacr.org/2012/685.pdf */
    /* a1 = a^((q - 3) / 4) */
    static const uint64_t QM3D4[6] = {
        0xee7fbfffffffeaaaULL, 0x07aaffffac54ffffULL,
        0xd9cc34a83dac3d89ULL, 0xd91dd2e13ce144afULL,
        0x92c6e9ed90d2eb35ULL, 0x0680447a8e5ff9a6ULL
    };
    struct fp2 a1;
    fp2_pow(&a1, a, QM3D4);

    struct fp2 alpha;
    fp2_sq(&alpha, &a1);
    fp2_mul(&alpha, &alpha, a);

    struct fp2 a0 = alpha;
    fp2_frobenius_map(&a0, 1);
    fp2_mul(&a0, &a0, &alpha);

    struct fp2 neg1;
    neg1.c0 = FP_NEGATIVE_ONE;
    fp_zero(&neg1.c1);

    if (fp2_eq(&a0, &neg1))
        return false;

    fp2_mul(&a1, &a1, a);

    if (fp2_eq(&alpha, &neg1)) {
        struct fp2 i_val;
        fp_zero(&i_val.c0);
        fp_one(&i_val.c1);
        fp2_mul(&a1, &a1, &i_val);
    } else {
        struct fp2 one_val;
        fp2_one(&one_val);
        fp2_add(&alpha, &alpha, &one_val);

        /* alpha = alpha^((q - 1) / 2) */
        static const uint64_t QM1D2[6] = {
            0xdcff7fffffffd555ULL, 0x0f55ffff58a9ffffULL,
            0xb39869507b587b12ULL, 0xb23ba5c279c2895fULL,
            0x258dd3db21a5d66bULL, 0x0d0088f51cbff34dULL
        };
        fp2_pow(&alpha, &alpha, QM1D2);
        fp2_mul(&a1, &a1, &alpha);
    }

    *r = a1;

    /* Verify */
    struct fp2 check;
    fp2_sq(&check, r);
    return fp2_eq(&check, a);
}

bool fp2_lexicographically_largest(const struct fp2 *a)
{
    /* Lexicographic order: compare c1 first, then c0 */
    if (!fp_is_zero(&a->c1))
        return fp_lexicographically_largest(&a->c1);
    return fp_lexicographically_largest(&a->c0);
}

/* ===== Frobenius coefficients for Fp6 and Fp12 ===== */

/* FROBENIUS_COEFF_FQ6_C1[i] — all in Montgomery form */
static const struct fp2 FROB_FQ6_C1[6] = {
    /* [0] = 1 */
    {{{0x760900000002fffdULL, 0xebf4000bc40c0002ULL,
       0x5f48985753c758baULL, 0x77ce585370525745ULL,
       0x5c071a97a256ec6dULL, 0x15f65ec3fa80e493ULL}},
     {{0,0,0,0,0,0}}},
    /* [1] */
    {{{0,0,0,0,0,0}},
     {{0xcd03c9e48671f071ULL, 0x5dab22461fcda5d2ULL,
       0x587042afd3851b95ULL, 0x8eb60ebe01bacb9eULL,
       0x03f97d6e83d050d2ULL, 0x18f0206554638741ULL}}},
    /* [2] */
    {{{0x30f1361b798a64e8ULL, 0xf3b8ddab7ece5a2aULL,
       0x16a8ca3ac61577f7ULL, 0xc26a2ff874fd029bULL,
       0x3636b76660701c6eULL, 0x051ba4ab241b6160ULL}},
     {{0,0,0,0,0,0}}},
    /* [3] */
    {{{0,0,0,0,0,0}},
     {{0x760900000002fffdULL, 0xebf4000bc40c0002ULL,
       0x5f48985753c758baULL, 0x77ce585370525745ULL,
       0x5c071a97a256ec6dULL, 0x15f65ec3fa80e493ULL}}},
    /* [4] */
    {{{0xcd03c9e48671f071ULL, 0x5dab22461fcda5d2ULL,
       0x587042afd3851b95ULL, 0x8eb60ebe01bacb9eULL,
       0x03f97d6e83d050d2ULL, 0x18f0206554638741ULL}},
     {{0,0,0,0,0,0}}},
    /* [5] */
    {{{0,0,0,0,0,0}},
     {{0x30f1361b798a64e8ULL, 0xf3b8ddab7ece5a2aULL,
       0x16a8ca3ac61577f7ULL, 0xc26a2ff874fd029bULL,
       0x3636b76660701c6eULL, 0x051ba4ab241b6160ULL}}}
};

static const struct fp2 FROB_FQ6_C2[6] = {
    /* [0] = 1 */
    {{{0x760900000002fffdULL, 0xebf4000bc40c0002ULL,
       0x5f48985753c758baULL, 0x77ce585370525745ULL,
       0x5c071a97a256ec6dULL, 0x15f65ec3fa80e493ULL}},
     {{0,0,0,0,0,0}}},
    /* [1] */
    {{{0x890dc9e4867545c3ULL, 0x2af322533285a5d5ULL,
       0x50880866309b7e2cULL, 0xa20d1b8c7e881024ULL,
       0x14e4f04fe2db9068ULL, 0x14e56d3f1564853aULL}},
     {{0,0,0,0,0,0}}},
    /* [2] */
    {{{0xcd03c9e48671f071ULL, 0x5dab22461fcda5d2ULL,
       0x587042afd3851b95ULL, 0x8eb60ebe01bacb9eULL,
       0x03f97d6e83d050d2ULL, 0x18f0206554638741ULL}},
     {{0,0,0,0,0,0}}},
    /* [3] */
    {{{0x43f5fffffffcaaaeULL, 0x32b7fff2ed47fffdULL,
       0x07e83a49a2e99d69ULL, 0xeca8f3318332bb7aULL,
       0xef148d1ea0f4c069ULL, 0x040ab3263eff0206ULL}},
     {{0,0,0,0,0,0}}},
    /* [4] */
    {{{0x30f1361b798a64e8ULL, 0xf3b8ddab7ece5a2aULL,
       0x16a8ca3ac61577f7ULL, 0xc26a2ff874fd029bULL,
       0x3636b76660701c6eULL, 0x051ba4ab241b6160ULL}},
     {{0,0,0,0,0,0}}},
    /* [5] */
    {{{0xecfb361b798dba3aULL, 0xc100ddb891865a2cULL,
       0x0ec08ff1232bda8eULL, 0xd5c13cc6f1ca4721ULL,
       0x47222a47bf7b5c04ULL, 0x0110f184e51c5f59ULL}},
     {{0,0,0,0,0,0}}}
};

static const struct fp2 FROB_FQ12_C1[12] = {
    /* [0] = 1 */
    {{{0x760900000002fffdULL, 0xebf4000bc40c0002ULL,
       0x5f48985753c758baULL, 0x77ce585370525745ULL,
       0x5c071a97a256ec6dULL, 0x15f65ec3fa80e493ULL}},
     {{0,0,0,0,0,0}}},
    /* [1] */
    {{{0x07089552b319d465ULL, 0xc6695f92b50a8313ULL,
       0x97e83cccd117228fULL, 0xa35baecab2dc29eeULL,
       0x1ce393ea5daace4dULL, 0x08f2220fb0fb66ebULL}},
     {{0xb2f66aad4ce5d646ULL, 0x5842a06bfc497cecULL,
       0xcf4895d42599d394ULL, 0xc11b9cba40a8e8d0ULL,
       0x2e3813cbe5a0de89ULL, 0x110eefda88847fafULL}}},
    /* [2] */
    {{{0xecfb361b798dba3aULL, 0xc100ddb891865a2cULL,
       0x0ec08ff1232bda8eULL, 0xd5c13cc6f1ca4721ULL,
       0x47222a47bf7b5c04ULL, 0x0110f184e51c5f59ULL}},
     {{0,0,0,0,0,0}}},
    /* [3] */
    {{{0x3e2f585da55c9ad1ULL, 0x4294213d86c18183ULL,
       0x382844c88b623732ULL, 0x92ad2afd19103e18ULL,
       0x1d794e4fac7cf0b9ULL, 0x0bd592fc7d825ec8ULL}},
     {{0x7bcfa7a25aa30fdaULL, 0xdc17dec12a927e7cULL,
       0x2f088dd86b4ebef1ULL, 0xd1ca2087da74d4a7ULL,
       0x2da2596696cebc1dULL, 0x0e2b7eedbbfd87d2ULL}}},
    /* [4] */
    {{{0x30f1361b798a64e8ULL, 0xf3b8ddab7ece5a2aULL,
       0x16a8ca3ac61577f7ULL, 0xc26a2ff874fd029bULL,
       0x3636b76660701c6eULL, 0x051ba4ab241b6160ULL}},
     {{0,0,0,0,0,0}}},
    /* [5] */
    {{{0x3726c30af242c66cULL, 0x7c2ac1aad1b6fe70ULL,
       0xa04007fbba4b14a2ULL, 0xef517c3266341429ULL,
       0x0095ba654ed2226bULL, 0x02e370eccc86f7ddULL}},
     {{0x82d83cf50dbce43fULL, 0xa2813e53df9d018fULL,
       0xc6f0caa53c65e181ULL, 0x7525cf528d50fe95ULL,
       0x4a85ed50f4798a6bULL, 0x171da0fd6cf8eebdULL}}},
    /* [6] = -1 */
    {{{0x43f5fffffffcaaaeULL, 0x32b7fff2ed47fffdULL,
       0x07e83a49a2e99d69ULL, 0xeca8f3318332bb7aULL,
       0xef148d1ea0f4c069ULL, 0x040ab3263eff0206ULL}},
     {{0,0,0,0,0,0}}},
    /* [7] */
    {{{0xb2f66aad4ce5d646ULL, 0x5842a06bfc497cecULL,
       0xcf4895d42599d394ULL, 0xc11b9cba40a8e8d0ULL,
       0x2e3813cbe5a0de89ULL, 0x110eefda88847fafULL}},
     {{0x07089552b319d465ULL, 0xc6695f92b50a8313ULL,
       0x97e83cccd117228fULL, 0xa35baecab2dc29eeULL,
       0x1ce393ea5daace4dULL, 0x08f2220fb0fb66ebULL}}},
    /* [8] */
    {{{0xcd03c9e48671f071ULL, 0x5dab22461fcda5d2ULL,
       0x587042afd3851b95ULL, 0x8eb60ebe01bacb9eULL,
       0x03f97d6e83d050d2ULL, 0x18f0206554638741ULL}},
     {{0,0,0,0,0,0}}},
    /* [9] */
    {{{0x7bcfa7a25aa30fdaULL, 0xdc17dec12a927e7cULL,
       0x2f088dd86b4ebef1ULL, 0xd1ca2087da74d4a7ULL,
       0x2da2596696cebc1dULL, 0x0e2b7eedbbfd87d2ULL}},
     {{0x3e2f585da55c9ad1ULL, 0x4294213d86c18183ULL,
       0x382844c88b623732ULL, 0x92ad2afd19103e18ULL,
       0x1d794e4fac7cf0b9ULL, 0x0bd592fc7d825ec8ULL}}},
    /* [10] */
    {{{0x890dc9e4867545c3ULL, 0x2af322533285a5d5ULL,
       0x50880866309b7e2cULL, 0xa20d1b8c7e881024ULL,
       0x14e4f04fe2db9068ULL, 0x14e56d3f1564853aULL}},
     {{0,0,0,0,0,0}}},
    /* [11] */
    {{{0x82d83cf50dbce43fULL, 0xa2813e53df9d018fULL,
       0xc6f0caa53c65e181ULL, 0x7525cf528d50fe95ULL,
       0x4a85ed50f4798a6bULL, 0x171da0fd6cf8eebdULL}},
     {{0x3726c30af242c66cULL, 0x7c2ac1aad1b6fe70ULL,
       0xa04007fbba4b14a2ULL, 0xef517c3266341429ULL,
       0x0095ba654ed2226bULL, 0x02e370eccc86f7ddULL}}}
};

/* ===== Fp6 = Fp2[v] / (v^3 - (u+1)) ===== */

void fp6_zero(struct fp6 *r) { fp2_zero(&r->c0); fp2_zero(&r->c1); fp2_zero(&r->c2); }
void fp6_one(struct fp6 *r) { fp2_one(&r->c0); fp2_zero(&r->c1); fp2_zero(&r->c2); }
bool fp6_is_zero(const struct fp6 *a) { return fp2_is_zero(&a->c0) && fp2_is_zero(&a->c1) && fp2_is_zero(&a->c2); }

void fp6_add(struct fp6 *r, const struct fp6 *a, const struct fp6 *b)
{
    fp2_add(&r->c0, &a->c0, &b->c0);
    fp2_add(&r->c1, &a->c1, &b->c1);
    fp2_add(&r->c2, &a->c2, &b->c2);
}

void fp6_sub(struct fp6 *r, const struct fp6 *a, const struct fp6 *b)
{
    fp2_sub(&r->c0, &a->c0, &b->c0);
    fp2_sub(&r->c1, &a->c1, &b->c1);
    fp2_sub(&r->c2, &a->c2, &b->c2);
}

void fp6_neg(struct fp6 *r, const struct fp6 *a)
{
    fp2_neg(&r->c0, &a->c0);
    fp2_neg(&r->c1, &a->c1);
    fp2_neg(&r->c2, &a->c2);
}

/* Multiply by v: (c0, c1, c2) -> (c2*beta, c0, c1) where beta = (1+u) */
void fp6_mul_by_nonresidue(struct fp6 *r, const struct fp6 *a)
{
    struct fp2 t = a->c2;
    fp2_mul_by_nonresidue(&t, &a->c2);
    r->c2 = a->c1;
    r->c1 = a->c0;
    r->c0 = t;
}

void fp6_mul(struct fp6 *r, const struct fp6 *a, const struct fp6 *b)
{
    struct fp2 aa, bb, cc;
    fp2_mul(&aa, &a->c0, &b->c0);
    fp2_mul(&bb, &a->c1, &b->c1);
    fp2_mul(&cc, &a->c2, &b->c2);

    /* t1 = ((b1+b2)*(a1+a2) - bb - cc) * beta + aa */
    struct fp2 t1, tmp;
    fp2_add(&t1, &b->c1, &b->c2);
    fp2_add(&tmp, &a->c1, &a->c2);
    fp2_mul(&t1, &t1, &tmp);
    fp2_sub(&t1, &t1, &bb);
    fp2_sub(&t1, &t1, &cc);
    fp2_mul_by_nonresidue(&t1, &t1);
    fp2_add(&t1, &t1, &aa);

    /* t3 = (b0+b2)*(a0+a2) - aa + bb - cc */
    struct fp2 t3;
    fp2_add(&t3, &b->c0, &b->c2);
    fp2_add(&tmp, &a->c0, &a->c2);
    fp2_mul(&t3, &t3, &tmp);
    fp2_sub(&t3, &t3, &aa);
    fp2_add(&t3, &t3, &bb);
    fp2_sub(&t3, &t3, &cc);

    /* t2 = (b0+b1)*(a0+a1) - aa - bb + cc*beta */
    struct fp2 t2;
    fp2_add(&t2, &b->c0, &b->c1);
    fp2_add(&tmp, &a->c0, &a->c1);
    fp2_mul(&t2, &t2, &tmp);
    fp2_sub(&t2, &t2, &aa);
    fp2_sub(&t2, &t2, &bb);
    struct fp2 cc_nr;
    fp2_mul_by_nonresidue(&cc_nr, &cc);
    fp2_add(&t2, &t2, &cc_nr);

    r->c0 = t1;
    r->c1 = t2;
    r->c2 = t3;
}

void fp6_sq(struct fp6 *r, const struct fp6 *a)
{
    struct fp2 s0, ab, s1, s2, bc, s3, s4;

    fp2_sq(&s0, &a->c0);
    fp2_mul(&ab, &a->c0, &a->c1);
    fp2_add(&s1, &ab, &ab);

    struct fp2 t;
    fp2_sub(&t, &a->c0, &a->c1);
    fp2_add(&s2, &t, &a->c2);
    fp2_sq(&s2, &s2);

    fp2_mul(&bc, &a->c1, &a->c2);
    fp2_add(&s3, &bc, &bc);

    fp2_sq(&s4, &a->c2);

    /* c0 = s3 * beta + s0 */
    fp2_mul_by_nonresidue(&r->c0, &s3);
    fp2_add(&r->c0, &r->c0, &s0);

    /* c1 = s4 * beta + s1 */
    fp2_mul_by_nonresidue(&r->c1, &s4);
    fp2_add(&r->c1, &r->c1, &s1);

    /* c2 = s1 + s2 + s3 - s0 - s4 */
    fp2_add(&r->c2, &s1, &s2);
    fp2_add(&r->c2, &r->c2, &s3);
    fp2_sub(&r->c2, &r->c2, &s0);
    fp2_sub(&r->c2, &r->c2, &s4);
}

void fp6_inv(struct fp6 *r, const struct fp6 *a)
{
    /* c0 = a2*a1*beta - a0^2, negated later */
    struct fp2 c0, c1, c2;

    /* c0 = a2*beta*a1, negated, + a0^2 */
    struct fp2 tmp;
    fp2_mul_by_nonresidue(&tmp, &a->c2);
    fp2_mul(&c0, &tmp, &a->c1);
    fp2_neg(&c0, &c0);
    struct fp2 a0s;
    fp2_sq(&a0s, &a->c0);
    fp2_add(&c0, &c0, &a0s);

    /* c1 = a2^2 * beta - a0*a1 */
    fp2_sq(&c1, &a->c2);
    fp2_mul_by_nonresidue(&c1, &c1);
    fp2_mul(&tmp, &a->c0, &a->c1);
    fp2_sub(&c1, &c1, &tmp);

    /* c2 = a1^2 - a0*a2 */
    fp2_sq(&c2, &a->c1);
    fp2_mul(&tmp, &a->c0, &a->c2);
    fp2_sub(&c2, &c2, &tmp);

    /* tmp1 = a2*c1 + a1*c2, tmp1 *= beta, tmp1 += a0*c0 */
    struct fp2 tmp1, tmp2;
    fp2_mul(&tmp1, &a->c2, &c1);
    fp2_mul(&tmp2, &a->c1, &c2);
    fp2_add(&tmp1, &tmp1, &tmp2);
    fp2_mul_by_nonresidue(&tmp1, &tmp1);
    fp2_mul(&tmp2, &a->c0, &c0);
    fp2_add(&tmp1, &tmp1, &tmp2);

    struct fp2 inv;
    fp2_inv(&inv, &tmp1);

    fp2_mul(&r->c0, &c0, &inv);
    fp2_mul(&r->c1, &c1, &inv);
    fp2_mul(&r->c2, &c2, &inv);
}

/* Multiply Fp6 element by (0, c1, 0) */
void fp6_mul_by_1(struct fp6 *r, const struct fp6 *a, const struct fp2 *c1)
{
    struct fp2 bb;
    fp2_mul(&bb, &a->c1, c1);

    struct fp2 t1, tmp;
    fp2_add(&tmp, &a->c1, &a->c2);
    fp2_mul(&t1, c1, &tmp);
    fp2_sub(&t1, &t1, &bb);
    fp2_mul_by_nonresidue(&t1, &t1);

    struct fp2 t2;
    fp2_add(&tmp, &a->c0, &a->c1);
    fp2_mul(&t2, c1, &tmp);
    fp2_sub(&t2, &t2, &bb);

    r->c0 = t1;
    r->c1 = t2;
    r->c2 = bb;
}

/* Multiply Fp6 element by (c0, c1, 0) */
void fp6_mul_by_01(struct fp6 *r, const struct fp6 *a, const struct fp2 *c0, const struct fp2 *c1)
{
    struct fp2 aa, bb;
    fp2_mul(&aa, &a->c0, c0);
    fp2_mul(&bb, &a->c1, c1);

    struct fp2 t1, tmp;
    fp2_add(&tmp, &a->c1, &a->c2);
    fp2_mul(&t1, c1, &tmp);
    fp2_sub(&t1, &t1, &bb);
    fp2_mul_by_nonresidue(&t1, &t1);
    fp2_add(&t1, &t1, &aa);

    struct fp2 t3;
    fp2_add(&tmp, &a->c0, &a->c2);
    fp2_mul(&t3, c0, &tmp);
    fp2_sub(&t3, &t3, &aa);
    fp2_add(&t3, &t3, &bb);

    struct fp2 t2, c01;
    fp2_add(&c01, c0, c1);
    fp2_add(&tmp, &a->c0, &a->c1);
    fp2_mul(&t2, &c01, &tmp);
    fp2_sub(&t2, &t2, &aa);
    fp2_sub(&t2, &t2, &bb);

    r->c0 = t1;
    r->c1 = t2;
    r->c2 = t3;
}

void fp6_frobenius_map(struct fp6 *a, int power)
{
    fp2_frobenius_map(&a->c0, power);
    fp2_frobenius_map(&a->c1, power);
    fp2_frobenius_map(&a->c2, power);

    fp2_mul(&a->c1, &a->c1, &FROB_FQ6_C1[power % 6]);
    fp2_mul(&a->c2, &a->c2, &FROB_FQ6_C2[power % 6]);
}

/* ===== Fp12 = Fp6[w] / (w^2 - v) ===== */

void fp12_zero(struct fp12 *r) { fp6_zero(&r->c0); fp6_zero(&r->c1); }
void fp12_one(struct fp12 *r) { fp6_one(&r->c0); fp6_zero(&r->c1); }
bool fp12_is_zero(const struct fp12 *a) { return fp6_is_zero(&a->c0) && fp6_is_zero(&a->c1); }

void fp12_add(struct fp12 *r, const struct fp12 *a, const struct fp12 *b)
{
    fp6_add(&r->c0, &a->c0, &b->c0);
    fp6_add(&r->c1, &a->c1, &b->c1);
}

void fp12_sub(struct fp12 *r, const struct fp12 *a, const struct fp12 *b)
{
    fp6_sub(&r->c0, &a->c0, &b->c0);
    fp6_sub(&r->c1, &a->c1, &b->c1);
}

void fp12_neg(struct fp12 *r, const struct fp12 *a)
{
    fp6_neg(&r->c0, &a->c0);
    fp6_neg(&r->c1, &a->c1);
}

void fp12_mul(struct fp12 *r, const struct fp12 *a, const struct fp12 *b)
{
    struct fp6 aa, bb;
    fp6_mul(&aa, &a->c0, &b->c0);
    fp6_mul(&bb, &a->c1, &b->c1);

    struct fp6 o;
    fp6_add(&o, &b->c0, &b->c1);

    struct fp6 c1_tmp;
    fp6_add(&c1_tmp, &a->c1, &a->c0);
    fp6_mul(&r->c1, &c1_tmp, &o);
    fp6_sub(&r->c1, &r->c1, &aa);
    fp6_sub(&r->c1, &r->c1, &bb);

    struct fp6 bb_nr;
    fp6_mul_by_nonresidue(&bb_nr, &bb);
    fp6_add(&r->c0, &bb_nr, &aa);
}

void fp12_sq(struct fp12 *r, const struct fp12 *a)
{
    struct fp6 ab;
    fp6_mul(&ab, &a->c0, &a->c1);

    struct fp6 c0c1;
    fp6_add(&c0c1, &a->c0, &a->c1);

    struct fp6 c0_tmp;
    fp6_mul_by_nonresidue(&c0_tmp, &a->c1);
    fp6_add(&c0_tmp, &c0_tmp, &a->c0);
    fp6_mul(&c0_tmp, &c0_tmp, &c0c1);
    fp6_sub(&c0_tmp, &c0_tmp, &ab);

    r->c1 = ab;
    fp6_add(&r->c1, &r->c1, &ab);

    struct fp6 ab_nr;
    fp6_mul_by_nonresidue(&ab_nr, &ab);
    fp6_sub(&r->c0, &c0_tmp, &ab_nr);
}

void fp12_inv(struct fp12 *r, const struct fp12 *a)
{
    struct fp6 c0s, c1s;
    fp6_sq(&c0s, &a->c0);
    fp6_sq(&c1s, &a->c1);
    fp6_mul_by_nonresidue(&c1s, &c1s);
    fp6_sub(&c0s, &c0s, &c1s);

    struct fp6 t;
    fp6_inv(&t, &c0s);

    fp6_mul(&r->c0, &a->c0, &t);
    fp6_mul(&r->c1, &a->c1, &t);
    fp6_neg(&r->c1, &r->c1);
}

void fp12_conjugate(struct fp12 *a)
{
    fp6_neg(&a->c1, &a->c1);
}

void fp12_frobenius_map(struct fp12 *a, int power)
{
    fp6_frobenius_map(&a->c0, power);
    fp6_frobenius_map(&a->c1, power);

    const struct fp2 *coeff = &FROB_FQ12_C1[power % 12];
    fp2_mul(&a->c1.c0, &a->c1.c0, coeff);
    fp2_mul(&a->c1.c1, &a->c1.c1, coeff);
    fp2_mul(&a->c1.c2, &a->c1.c2, coeff);
}

void fp12_mul_by_014(struct fp12 *f, const struct fp2 *c0, const struct fp2 *c1, const struct fp2 *c4)
{
    struct fp6 aa, bb;
    fp6_mul_by_01(&aa, &f->c0, c0, c1);
    fp6_mul_by_1(&bb, &f->c1, c4);

    struct fp2 o;
    fp2_add(&o, c1, c4);

    struct fp6 t;
    fp6_add(&t, &f->c1, &f->c0);
    fp6_mul_by_01(&t, &t, c0, &o);
    fp6_sub(&t, &t, &aa);
    fp6_sub(&t, &t, &bb);
    f->c1 = t;

    struct fp6 bb_nr;
    fp6_mul_by_nonresidue(&bb_nr, &bb);
    fp6_add(&f->c0, &bb_nr, &aa);
}

static void fp12_pow_u64(struct fp12 *r, const struct fp12 *a, uint64_t exp)
{
    struct fp12 result;
    fp12_one(&result);
    struct fp12 base = *a;

    while (exp > 0) {
        if (exp & 1)
            fp12_mul(&result, &result, &base);
        fp12_sq(&base, &base);
        exp >>= 1;
    }
    *r = result;
}

/* ===== G1 point operations (Jacobian, y^2 = x^3 + 4) ===== */

/* B coefficient = 4 in Montgomery form */
static const struct fp G1_B = {{
    0xaa270000000cfff3ULL, 0x53cc0032fc34000aULL,
    0x478fe97a6b0a807fULL, 0xb1d37ebee6ba24d7ULL,
    0x8ec9733bbf78ab2fULL, 0x09d645513d83de7eULL
}};

/* Generator of G1 in Montgomery form (used by g1_generator) */
const struct fp G1_GEN_X = {{
    0x5cb38790fd530c16ULL, 0x7817fc679976fff5ULL,
    0x154f95c7143ba1c1ULL, 0xf0ae6acdf3d0e747ULL,
    0xedce6ecc21dbf440ULL, 0x120177419e0bfb75ULL
}};
const struct fp G1_GEN_Y = {{
    0xbaac93d50ce72271ULL, 0x8c22631a7918fd8eULL,
    0xdd595f13570725ceULL, 0x51ac582950405194ULL,
    0x0e1c8c3fad0059c0ULL, 0x0bbc3efc5008a26aULL
}};

void g1_identity(struct g1_point *p)
{
    fp_zero(&p->x);
    fp_one(&p->y);
    fp_zero(&p->z);
}

bool g1_is_identity(const struct g1_point *p)
{
    return fp_is_zero(&p->z);
}

void g1_neg(struct g1_point *r, const struct g1_point *p)
{
    r->x = p->x;
    fp_neg(&r->y, &p->y);
    r->z = p->z;
}

void g1_double(struct g1_point *r, const struct g1_point *p)
{
    if (g1_is_identity(p)) { *r = *p; return; }

    struct fp a, b, c, d, e, f;
    fp_sq(&a, &p->x);
    fp_sq(&b, &p->y);
    fp_sq(&c, &b);

    fp_add(&d, &p->x, &b);
    fp_sq(&d, &d);
    fp_sub(&d, &d, &a);
    fp_sub(&d, &d, &c);
    fp_add(&d, &d, &d);

    fp_add(&e, &a, &a);
    fp_add(&e, &e, &a);

    fp_sq(&f, &e);

    fp_mul(&r->z, &p->y, &p->z);
    fp_add(&r->z, &r->z, &r->z);

    fp_sub(&r->x, &f, &d);
    fp_sub(&r->x, &r->x, &d);

    fp_sub(&r->y, &d, &r->x);
    fp_mul(&r->y, &r->y, &e);
    fp_add(&c, &c, &c);
    fp_add(&c, &c, &c);
    fp_add(&c, &c, &c);
    fp_sub(&r->y, &r->y, &c);
}

void g1_add(struct g1_point *r, const struct g1_point *a, const struct g1_point *b)
{
    if (g1_is_identity(a)) { *r = *b; return; }
    if (g1_is_identity(b)) { *r = *a; return; }

    struct fp z1z1, z2z2, u1, u2, s1, s2;
    fp_sq(&z1z1, &a->z);
    fp_sq(&z2z2, &b->z);
    fp_mul(&u1, &a->x, &z2z2);
    fp_mul(&u2, &b->x, &z1z1);
    fp_mul(&s1, &a->y, &b->z);
    fp_mul(&s1, &s1, &z2z2);
    fp_mul(&s2, &b->y, &a->z);
    fp_mul(&s2, &s2, &z1z1);

    if (fp_eq(&u1, &u2) && fp_eq(&s1, &s2)) {
        g1_double(r, a);
        return;
    }

    struct fp h, i, j, rr, v;
    fp_sub(&h, &u2, &u1);
    fp_add(&i, &h, &h);
    fp_sq(&i, &i);
    fp_mul(&j, &h, &i);
    fp_sub(&rr, &s2, &s1);
    fp_add(&rr, &rr, &rr);
    fp_mul(&v, &u1, &i);

    fp_sq(&r->x, &rr);
    fp_sub(&r->x, &r->x, &j);
    fp_sub(&r->x, &r->x, &v);
    fp_sub(&r->x, &r->x, &v);

    fp_sub(&r->y, &v, &r->x);
    fp_mul(&r->y, &r->y, &rr);
    fp_mul(&s1, &s1, &j);
    fp_add(&s1, &s1, &s1);
    fp_sub(&r->y, &r->y, &s1);

    fp_add(&r->z, &a->z, &b->z);
    fp_sq(&r->z, &r->z);
    fp_sub(&r->z, &r->z, &z1z1);
    fp_sub(&r->z, &r->z, &z2z2);
    fp_mul(&r->z, &r->z, &h);
}

void g1_to_affine(struct fp *ax, struct fp *ay, const struct g1_point *p)
{
    struct fp zinv;
    fp_inv(&zinv, &p->z);
    struct fp zinv2;
    fp_sq(&zinv2, &zinv);
    fp_mul(ax, &p->x, &zinv2);
    struct fp zinv3;
    fp_mul(&zinv3, &zinv2, &zinv);
    fp_mul(ay, &p->y, &zinv3);
}

bool g1_from_compressed(struct g1_point *p, const uint8_t in[48])
{
    uint8_t copy[48];
    memcpy(copy, in, 48);

    bool compressed = (copy[0] >> 7) & 1;
    bool infinity = (copy[0] >> 6) & 1;
    bool greatest = (copy[0] >> 5) & 1;
    copy[0] &= 0x1f;

    if (!compressed) return false;

    if (infinity) {
        g1_identity(p);
        return true;
    }

    struct fp x;
    if (!fp_from_bytes(&x, copy)) return false;

    /* y^2 = x^3 + 4 */
    struct fp x3b;
    fp_sq(&x3b, &x);
    fp_mul(&x3b, &x3b, &x);
    fp_add(&x3b, &x3b, &G1_B);

    struct fp y;
    if (!fp_sqrt(&y, &x3b)) return false;

    struct fp negy;
    fp_neg(&negy, &y);

    bool y_largest = fp_lexicographically_largest(&y);
    if (y_largest != greatest)
        y = negy;

    p->x = x;
    p->y = y;
    fp_one(&p->z);
    return true;
}

bool g1_from_uncompressed(struct g1_point *p, const uint8_t in[96])
{
    uint8_t copy[96];
    memcpy(copy, in, 96);

    bool compressed = (copy[0] >> 7) & 1;
    bool infinity = (copy[0] >> 6) & 1;
    copy[0] &= 0x1f;

    if (compressed) return false;

    if (infinity) {
        g1_identity(p);
        return true;
    }

    struct fp x, y;
    if (!fp_from_bytes(&x, copy)) return false;
    if (!fp_from_bytes(&y, copy + 48)) return false;

    /* Verify on curve: y^2 == x^3 + 4 */
    struct fp lhs, rhs;
    fp_sq(&lhs, &y);
    fp_sq(&rhs, &x);
    fp_mul(&rhs, &rhs, &x);
    fp_add(&rhs, &rhs, &G1_B);
    if (!fp_eq(&lhs, &rhs)) return false;

    p->x = x;
    p->y = y;
    fp_one(&p->z);
    return true;
}

/* ===== G2 point operations (Jacobian, y^2 = x^3 + 4(u+1)) ===== */

/* B coefficient for G2 = 4*(1+u) in Montgomery form */
static const struct fp2 G2_B = {
    {{0xaa270000000cfff3ULL, 0x53cc0032fc34000aULL,
      0x478fe97a6b0a807fULL, 0xb1d37ebee6ba24d7ULL,
      0x8ec9733bbf78ab2fULL, 0x09d645513d83de7eULL}},
    {{0xaa270000000cfff3ULL, 0x53cc0032fc34000aULL,
      0x478fe97a6b0a807fULL, 0xb1d37ebee6ba24d7ULL,
      0x8ec9733bbf78ab2fULL, 0x09d645513d83de7eULL}}
};

void g2_identity(struct g2_point *p)
{
    fp2_zero(&p->x);
    fp2_one(&p->y);
    fp2_zero(&p->z);
}

bool g2_is_identity(const struct g2_point *p)
{
    return fp2_is_zero(&p->z);
}

void g2_neg(struct g2_point *r, const struct g2_point *p)
{
    r->x = p->x;
    fp2_neg(&r->y, &p->y);
    r->z = p->z;
}

void g2_double(struct g2_point *r, const struct g2_point *p)
{
    if (g2_is_identity(p)) { *r = *p; return; }

    struct fp2 a, b, c, d, e, f;
    fp2_sq(&a, &p->x);
    fp2_sq(&b, &p->y);
    fp2_sq(&c, &b);

    fp2_add(&d, &p->x, &b);
    fp2_sq(&d, &d);
    fp2_sub(&d, &d, &a);
    fp2_sub(&d, &d, &c);
    fp2_add(&d, &d, &d);

    fp2_add(&e, &a, &a);
    fp2_add(&e, &e, &a);

    fp2_sq(&f, &e);

    fp2_mul(&r->z, &p->y, &p->z);
    fp2_add(&r->z, &r->z, &r->z);

    fp2_sub(&r->x, &f, &d);
    fp2_sub(&r->x, &r->x, &d);

    fp2_sub(&r->y, &d, &r->x);
    fp2_mul(&r->y, &r->y, &e);
    fp2_add(&c, &c, &c);
    fp2_add(&c, &c, &c);
    fp2_add(&c, &c, &c);
    fp2_sub(&r->y, &r->y, &c);
}

void g2_add(struct g2_point *r, const struct g2_point *a, const struct g2_point *b)
{
    if (g2_is_identity(a)) { *r = *b; return; }
    if (g2_is_identity(b)) { *r = *a; return; }

    struct fp2 z1z1, z2z2, u1, u2, s1, s2;
    fp2_sq(&z1z1, &a->z);
    fp2_sq(&z2z2, &b->z);
    fp2_mul(&u1, &a->x, &z2z2);
    fp2_mul(&u2, &b->x, &z1z1);
    fp2_mul(&s1, &a->y, &b->z);
    fp2_mul(&s1, &s1, &z2z2);
    fp2_mul(&s2, &b->y, &a->z);
    fp2_mul(&s2, &s2, &z1z1);

    if (fp2_eq(&u1, &u2) && fp2_eq(&s1, &s2)) {
        g2_double(r, a);
        return;
    }

    struct fp2 h, i, j, rr, v;
    fp2_sub(&h, &u2, &u1);
    fp2_add(&i, &h, &h);
    fp2_sq(&i, &i);
    fp2_mul(&j, &h, &i);
    fp2_sub(&rr, &s2, &s1);
    fp2_add(&rr, &rr, &rr);
    fp2_mul(&v, &u1, &i);

    fp2_sq(&r->x, &rr);
    fp2_sub(&r->x, &r->x, &j);
    fp2_sub(&r->x, &r->x, &v);
    fp2_sub(&r->x, &r->x, &v);

    fp2_sub(&r->y, &v, &r->x);
    fp2_mul(&r->y, &r->y, &rr);
    fp2_mul(&s1, &s1, &j);
    fp2_add(&s1, &s1, &s1);
    fp2_sub(&r->y, &r->y, &s1);

    fp2_add(&r->z, &a->z, &b->z);
    fp2_sq(&r->z, &r->z);
    fp2_sub(&r->z, &r->z, &z1z1);
    fp2_sub(&r->z, &r->z, &z2z2);
    fp2_mul(&r->z, &r->z, &h);
}

void g2_to_affine(struct fp2 *ax, struct fp2 *ay, const struct g2_point *p)
{
    struct fp2 zinv;
    fp2_inv(&zinv, &p->z);
    struct fp2 zinv2;
    fp2_sq(&zinv2, &zinv);
    fp2_mul(ax, &p->x, &zinv2);
    struct fp2 zinv3;
    fp2_mul(&zinv3, &zinv2, &zinv);
    fp2_mul(ay, &p->y, &zinv3);
}

bool g2_from_compressed(struct g2_point *p, const uint8_t in[96])
{
    uint8_t copy[96];
    memcpy(copy, in, 96);

    bool compressed = (copy[0] >> 7) & 1;
    bool infinity = (copy[0] >> 6) & 1;
    bool greatest = (copy[0] >> 5) & 1;
    copy[0] &= 0x1f;

    if (!compressed) return false;

    if (infinity) {
        g2_identity(p);
        return true;
    }

    /* Fp2 x: c1 is in first 48 bytes, c0 in next 48 bytes */
    struct fp2 x;
    if (!fp_from_bytes(&x.c1, copy)) return false;
    if (!fp_from_bytes(&x.c0, copy + 48)) return false;

    /* y^2 = x^3 + 4(1+u) */
    struct fp2 x3b;
    fp2_sq(&x3b, &x);
    fp2_mul(&x3b, &x3b, &x);
    fp2_add(&x3b, &x3b, &G2_B);

    struct fp2 y;
    if (!fp2_sqrt(&y, &x3b)) return false;

    struct fp2 negy;
    fp2_neg(&negy, &y);

    bool y_largest = fp2_lexicographically_largest(&y);
    if (y_largest != greatest)
        y = negy;

    p->x = x;
    p->y = y;
    fp2_one(&p->z);
    return true;
}

bool g2_from_uncompressed(struct g2_point *p, const uint8_t in[192])
{
    uint8_t copy[192];
    memcpy(copy, in, 192);

    bool compressed = (copy[0] >> 7) & 1;
    bool infinity = (copy[0] >> 6) & 1;
    copy[0] &= 0x1f;

    if (compressed) return false;

    if (infinity) {
        g2_identity(p);
        return true;
    }

    /* Fp2 x: c1 first 48 bytes, c0 next 48 */
    struct fp2 x, y;
    if (!fp_from_bytes(&x.c1, copy)) return false;
    if (!fp_from_bytes(&x.c0, copy + 48)) return false;
    if (!fp_from_bytes(&y.c1, copy + 96)) return false;
    if (!fp_from_bytes(&y.c0, copy + 144)) return false;

    /* Verify: y^2 == x^3 + 4(1+u) */
    struct fp2 lhs, rhs;
    fp2_sq(&lhs, &y);
    fp2_sq(&rhs, &x);
    fp2_mul(&rhs, &rhs, &x);
    fp2_add(&rhs, &rhs, &G2_B);
    if (!fp2_eq(&lhs, &rhs)) return false;

    p->x = x;
    p->y = y;
    fp2_one(&p->z);
    return true;
}

/* ===== Optimal Ate pairing ===== */

/* BLS parameter x = -0xd201000000010000 */
#define BLS_X 0xd201000000010000ULL

/* Fr modulus is 255-bit, so a scalar safely holds 254 bits of payload
 * (BLS12-381). Was previously 253 (BN-254), which broke Sapling Groth16
 * verification — bit 253 of the nullifier landed in the wrong packed
 * scalar. Used by both multipack_bytes_to_fr variants. */
#define BLS12_381_FR_CAPACITY 254

/* Doubling step for Miller loop (Algorithm 26 from https://eprint.iacr.org/2010/354.pdf) */
static void doubling_step(struct g2_point *r, struct fp2 *c0, struct fp2 *c1, struct fp2 *c2)
{
    struct fp2 tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, zsquared;

    fp2_sq(&tmp0, &r->x);
    fp2_sq(&tmp1, &r->y);
    fp2_sq(&tmp2, &tmp1);

    fp2_add(&tmp3, &tmp1, &r->x);
    fp2_sq(&tmp3, &tmp3);
    fp2_sub(&tmp3, &tmp3, &tmp0);
    fp2_sub(&tmp3, &tmp3, &tmp2);
    fp2_add(&tmp3, &tmp3, &tmp3);

    fp2_add(&tmp4, &tmp0, &tmp0);
    fp2_add(&tmp4, &tmp4, &tmp0);

    fp2_add(&tmp6, &r->x, &tmp4);

    fp2_sq(&tmp5, &tmp4);

    fp2_sq(&zsquared, &r->z);

    r->x = tmp5;
    fp2_sub(&r->x, &r->x, &tmp3);
    fp2_sub(&r->x, &r->x, &tmp3);

    fp2_add(&r->z, &r->y, &r->z);
    fp2_sq(&r->z, &r->z);
    fp2_sub(&r->z, &r->z, &tmp1);
    fp2_sub(&r->z, &r->z, &zsquared);

    r->y = tmp3;
    fp2_sub(&r->y, &r->y, &r->x);
    fp2_mul(&r->y, &r->y, &tmp4);

    fp2_add(&tmp2, &tmp2, &tmp2);
    fp2_add(&tmp2, &tmp2, &tmp2);
    fp2_add(&tmp2, &tmp2, &tmp2);
    fp2_sub(&r->y, &r->y, &tmp2);

    fp2_mul(&tmp3, &tmp4, &zsquared);
    fp2_add(&tmp3, &tmp3, &tmp3);
    fp2_neg(&tmp3, &tmp3);

    fp2_sq(&tmp6, &tmp6);
    fp2_sub(&tmp6, &tmp6, &tmp0);
    fp2_sub(&tmp6, &tmp6, &tmp5);
    fp2_add(&tmp1, &tmp1, &tmp1);
    fp2_add(&tmp1, &tmp1, &tmp1);
    fp2_sub(&tmp6, &tmp6, &tmp1);

    fp2_mul(&tmp0, &r->z, &zsquared);
    fp2_add(&tmp0, &tmp0, &tmp0);

    *c0 = tmp0;
    *c1 = tmp3;
    *c2 = tmp6;
}

/* Addition step for Miller loop (Algorithm 27) */
static void addition_step(struct g2_point *r, const struct fp2 *qx, const struct fp2 *qy,
                           struct fp2 *c0, struct fp2 *c1, struct fp2 *c2)
{
    struct fp2 zsquared, ysquared, t0, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;

    fp2_sq(&zsquared, &r->z);
    fp2_sq(&ysquared, qy);

    fp2_mul(&t0, &zsquared, qx);

    fp2_add(&t1, qy, &r->z);
    fp2_sq(&t1, &t1);
    fp2_sub(&t1, &t1, &ysquared);
    fp2_sub(&t1, &t1, &zsquared);
    fp2_mul(&t1, &t1, &zsquared);

    fp2_sub(&t2, &t0, &r->x);

    fp2_sq(&t3, &t2);

    fp2_add(&t4, &t3, &t3);
    fp2_add(&t4, &t4, &t4);

    fp2_mul(&t5, &t4, &t2);

    fp2_sub(&t6, &t1, &r->y);
    fp2_sub(&t6, &t6, &r->y);

    fp2_mul(&t9, &t6, qx);

    fp2_mul(&t7, &t4, &r->x);

    fp2_sq(&r->x, &t6);
    fp2_sub(&r->x, &r->x, &t5);
    fp2_sub(&r->x, &r->x, &t7);
    fp2_sub(&r->x, &r->x, &t7);

    fp2_add(&r->z, &r->z, &t2);
    fp2_sq(&r->z, &r->z);
    fp2_sub(&r->z, &r->z, &zsquared);
    fp2_sub(&r->z, &r->z, &t3);

    fp2_add(&t10, qy, &r->z);

    fp2_sub(&t8, &t7, &r->x);
    fp2_mul(&t8, &t8, &t6);

    fp2_mul(&t0, &r->y, &t5);
    fp2_add(&t0, &t0, &t0);

    r->y = t8;
    fp2_sub(&r->y, &r->y, &t0);

    fp2_sq(&t10, &t10);
    fp2_sub(&t10, &t10, &ysquared);
    struct fp2 ztsquared;
    fp2_sq(&ztsquared, &r->z);
    fp2_sub(&t10, &t10, &ztsquared);

    fp2_add(&t9, &t9, &t9);
    fp2_sub(&t9, &t9, &t10);

    struct fp2 t10_out;
    t10_out = r->z;
    fp2_add(&t10_out, &t10_out, &t10_out);

    fp2_neg(&t6, &t6);
    struct fp2 t1_out;
    fp2_add(&t1_out, &t6, &t6);

    *c0 = t10_out;
    *c1 = t1_out;
    *c2 = t9;
}

/* Twist: apply line function evaluation to G1 affine point */
static void ell(struct fp12 *f, const struct fp2 *c0_coeff, const struct fp2 *c1_coeff,
                const struct fp2 *c2_coeff, const struct fp *px, const struct fp *py)
{
    struct fp2 c0_eval = *c0_coeff;
    struct fp2 c1_eval = *c1_coeff;

    /* c0 *= py (both components) */
    fp_mul(&c0_eval.c0, &c0_eval.c0, py);
    fp_mul(&c0_eval.c1, &c0_eval.c1, py);

    /* c1 *= px (both components) */
    fp_mul(&c1_eval.c0, &c1_eval.c0, px);
    fp_mul(&c1_eval.c1, &c1_eval.c1, px);

    fp12_mul_by_014(f, c2_coeff, &c1_eval, &c0_eval);
}

/* Miller loop only (no final exponentiation).
 * Returns the Miller loop result f for pairing(P, Q). */
static void bls12_381_miller_loop(struct fp12 *result, const struct g1_point *p, const struct g2_point *q)
{
    if (g1_is_identity(p) || g2_is_identity(q)) {
        fp12_one(result);
        return;
    }

    /* Convert P to affine */
    struct fp px, py;
    g1_to_affine(&px, &py, p);

    /* Convert Q to affine for addition steps */
    struct fp2 qx, qy;
    g2_to_affine(&qx, &qy, q);

    /* Miller loop R starts at Q */
    struct g2_point r_point;
    r_point.x = qx;
    r_point.y = qy;
    fp2_one(&r_point.z);

    struct fp12 f;
    fp12_one(&f);

    /* Iterate over bits of BLS_X >> 1 = 0x6900800000008000 */
    uint64_t x_half = BLS_X >> 1;
    bool found_one = false;

    for (int i = 63; i >= 0; i--) {
        bool bit = (x_half >> i) & 1;
        if (!found_one) {
            found_one = bit;
            continue;
        }

        struct fp2 c0, c1, c2;
        doubling_step(&r_point, &c0, &c1, &c2);
        ell(&f, &c0, &c1, &c2, &px, &py);

        if (bit) {
            addition_step(&r_point, &qx, &qy, &c0, &c1, &c2);
            ell(&f, &c0, &c1, &c2, &px, &py);
        }

        fp12_sq(&f, &f);
    }

    /* Final doubling step */
    {
        struct fp2 c0, c1, c2;
        doubling_step(&r_point, &c0, &c1, &c2);
        ell(&f, &c0, &c1, &c2, &px, &py);
    }

    /* BLS_X is negative, so conjugate */
    fp12_conjugate(&f);
    *result = f;
}

/* Final exponentiation: f^((q^12 - 1) / r) */
static void bls12_381_final_exp(struct fp12 *result, const struct fp12 *f_in)
{
    struct fp12 f = *f_in;

    /* Final exponentiation */
    /* Easy part: f = f^(q^6-1) * f^(q^2+1) */
    struct fp12 f1 = f;
    fp12_conjugate(&f1);

    struct fp12 f2;
    fp12_inv(&f2, &f);
    fp12_mul(&f, &f1, &f2);
    f2 = f;
    fp12_frobenius_map(&f, 2);
    fp12_mul(&f, &f, &f2);

    /* Hard part */
    struct fp12 y0, y1, y2, y3;

    struct fp12 r_save = f;

    fp12_sq(&y0, &f);

    fp12_pow_u64(&y1, &y0, BLS_X);
    fp12_conjugate(&y1); /* exp_by_x conjugates since BLS_X_IS_NEGATIVE */

    fp12_pow_u64(&y2, &y1, BLS_X >> 1);
    fp12_conjugate(&y2);

    y3 = r_save;
    fp12_conjugate(&y3);
    fp12_mul(&y1, &y1, &y3);
    fp12_conjugate(&y1);
    fp12_mul(&y1, &y1, &y2);

    fp12_pow_u64(&y2, &y1, BLS_X);
    fp12_conjugate(&y2);

    fp12_pow_u64(&y3, &y2, BLS_X);
    fp12_conjugate(&y3);

    fp12_conjugate(&y1);
    fp12_mul(&y3, &y3, &y1);

    fp12_conjugate(&y1);
    fp12_frobenius_map(&y1, 3);

    fp12_frobenius_map(&y2, 2);
    fp12_mul(&y1, &y1, &y2);

    fp12_pow_u64(&y2, &y3, BLS_X);
    fp12_conjugate(&y2);
    fp12_mul(&y2, &y2, &y0);
    fp12_mul(&y2, &y2, &r_save);
    fp12_mul(&y1, &y1, &y2);

    y2 = y3;
    fp12_frobenius_map(&y2, 1);
    fp12_mul(&y1, &y1, &y2);

    *result = y1;
}

void bls12_381_pairing(struct fp12 *result, const struct g1_point *p, const struct g2_point *q)
{
    struct fp12 ml;
    bls12_381_miller_loop(&ml, p, q);
    bls12_381_final_exp(result, &ml);
}

/* Multi-pairing: product of Miller loops, single final exponentiation.
 * Checks if product of pairings equals 1 in GT.
 * This is both faster and more numerically robust than computing
 * individual pairings and multiplying in GT. */
bool bls12_381_multi_pairing_check(const struct g1_point *a_pts, const struct g2_point *b_pts,
                                    size_t n_pairs)
{
    struct fp12 acc;
    fp12_one(&acc);

    /* Accumulate Miller loop products (before final exponentiation) */
    for (size_t i = 0; i < n_pairs; i++) {
        if (g1_is_identity(&a_pts[i]) || g2_is_identity(&b_pts[i]))
            continue;

        struct fp12 ml;
        bls12_381_miller_loop(&ml, &a_pts[i], &b_pts[i]);
        fp12_mul(&acc, &acc, &ml);
    }

    /* Single final exponentiation over the accumulated product */
    struct fp12 result;
    bls12_381_final_exp(&result, &acc);

    struct fp12 one;
    fp12_one(&one);
    return memcmp(&result, &one, sizeof(struct fp12)) == 0;
}

/* Constant-time conditional move for Fp.  If mask == all-ones,
 * dst = src; if mask == 0, dst is unchanged.  No branches, no
 * secret-dependent memory access. */
static inline void fp_cmov(struct fp *dst, const struct fp *src, uint64_t mask)
{
    for (int i = 0; i < 6; i++)
        dst->d[i] = (dst->d[i] & ~mask) | (src->d[i] & mask);
}

static inline void g1_cmov(struct g1_point *dst, const struct g1_point *src,
                           uint64_t mask)
{
    fp_cmov(&dst->x, &src->x, mask);
    fp_cmov(&dst->y, &src->y, mask);
    fp_cmov(&dst->z, &src->z, mask);
}

/* G1 scalar multiplication: r = scalar * p (double-and-add, 256 bits).
 *
 * Threat model: callers include Groth16 proving, where the scalar
 * is a secret blinding factor (r_blind, s_blind — groth16_prover.c:906,
 * 934).  An attacker measuring wall-time of the prover must not learn
 * the Hamming weight of the scalar.
 *
 * CT properties:
 *   - Every iteration unconditionally runs g1_add and g1_double, so the
 *     total number of field-arithmetic ops is fixed (256 adds, 256
 *     doubles) independent of scalar content.
 *   - The bit-dependent update of `result` is done via `g1_cmov` masked
 *     on the bit, not a branch.
 *
 * Any future reviewer: do not reintroduce `if (bit)` without revisiting
 * the side-channel threat model. */
void g1_scalar_mul(struct g1_point *r, const struct g1_point *p, const uint64_t scalar[4])
{
    struct g1_point result;
    g1_identity(&result);
    struct g1_point base = *p;

    for (int i = 0; i < 4; i++) {
        for (int bit = 0; bit < 64; bit++) {
            struct g1_point sum;
            g1_add(&sum, &result, &base);
            uint64_t mask = (uint64_t)0 - ((scalar[i] >> bit) & 1ULL);
            g1_cmov(&result, &sum, mask);
            g1_double(&base, &base);
        }
    }
    *r = result;
}

/* G2 scalar multiplication: r = scalar * p (double-and-add, 256 bits).
 *
 * Unlike g1_scalar_mul, this is NOT constant-time: its only caller is the
 * subgroup membership check below, whose input is the PUBLIC proof point
 * (no secret material), so a side channel reveals nothing. Keeping it
 * branch-on-bit makes the [r]P==O check trivially auditable. */
static void g2_scalar_mul(struct g2_point *r, const struct g2_point *p,
                          const uint64_t scalar[4])
{
    struct g2_point result;
    g2_identity(&result);
    struct g2_point base = *p;

    for (int i = 0; i < 4; i++) {
        for (int bit = 0; bit < 64; bit++) {
            if ((scalar[i] >> bit) & 1ULL) {
                struct g2_point sum;
                g2_add(&sum, &result, &base);
                result = sum;
            }
            struct g2_point dbl;
            g2_double(&dbl, &base);
            base = dbl;
        }
    }
    *r = result;
}

/* BLS12-381 subgroup order r (little-endian 64-bit limbs):
 * r = 0x73eda753299d7d483339d80809a1d80553bda402fffe5bfeffffffff00000001
 * (matches FR_MODULUS in fr.c / fr_avx512.c). */
static const uint64_t BLS12_381_R[4] = {
    0xffffffff00000001ULL, 0x53bda402fffe5bfeULL,
    0x3339d80809a1d805ULL, 0x73eda753299d7d48ULL
};

/* Subgroup membership: G1/G2 have large cofactors, so an on-curve point can
 * still carry a torsion (non-prime-order) component that the pairing accepts —
 * a Groth16 soundness gap (a prover could forge an accepting "proof").  A point
 * P lies in the prime-order-r subgroup iff [r]P == O.  These checks ADD ONLY
 * rejections: they cannot accept anything the on-curve check already rejected,
 * and cannot cause a UAF/split.
 *
 * The identity (point at infinity) IS in the subgroup and MUST pass:
 * [r]O == O.  Callers must never feed the identity branch of
 * g1/g2_from_compressed through anything that would reject it — but [r]O==O
 * holds here regardless, so the function is correct for the identity too. */
bool g1_in_subgroup(const struct g1_point *p)
{
    struct g1_point rp;
    g1_scalar_mul(&rp, p, BLS12_381_R);
    return g1_is_identity(&rp);
}

bool g2_in_subgroup(const struct g2_point *p)
{
    struct g2_point rp;
    g2_scalar_mul(&rp, p, BLS12_381_R);
    return g2_is_identity(&rp);
}

/* Groth16 proof deserialization: A (G1, 48B) + B (G2, 96B) + C (G1, 48B) = 192B.
 *
 * Each point is checked for (1) on-curve (inside g{1,2}_from_compressed) AND
 * (2) prime-order subgroup membership.  Without (2), a malicious prover could
 * submit on-curve-but-torsion points that the pairing still accepts, breaking
 * Groth16 soundness on the consensus Sapling proof path. */
bool groth16_proof_read(struct groth16_proof *proof, const uint8_t data[192])
{
    if (!g1_from_compressed(&proof->a, data))
        LOG_FAIL("groth16", "proof_read: g1_from_compressed(A) failed");
    if (!g1_in_subgroup(&proof->a))
        LOG_FAIL("groth16", "proof_read: A not in G1 prime-order subgroup");
    if (!g2_from_compressed(&proof->b, data + 48))
        LOG_FAIL("groth16", "proof_read: g2_from_compressed(B) failed");
    if (!g2_in_subgroup(&proof->b))
        LOG_FAIL("groth16", "proof_read: B not in G2 prime-order subgroup");
    if (!g1_from_compressed(&proof->c, data + 144))
        LOG_FAIL("groth16", "proof_read: g1_from_compressed(C) failed");
    if (!g1_in_subgroup(&proof->c))
        LOG_FAIL("groth16", "proof_read: C not in G1 prime-order subgroup");
    return true;
}

/* Multipack: pack bytes into Fr scalars (253 bits per scalar, LE bit order) */
void multipack_bytes_to_fr(uint64_t (*out)[4], size_t *n_out,
                           const uint8_t *bytes, size_t n_bytes)
{
    const size_t capacity = BLS12_381_FR_CAPACITY;
    size_t n_bits = n_bytes * 8;
    size_t n_scalars = (n_bits + capacity - 1) / capacity;

    for (size_t s = 0; s < n_scalars; s++) {
        memset(out[s], 0, 32);
        uint64_t cur[4] = {0, 0, 0, 0};
        uint64_t coeff[4] = {1, 0, 0, 0};

        size_t bit_start = s * capacity;
        size_t bit_end = bit_start + capacity;
        if (bit_end > n_bits) bit_end = n_bits;

        for (size_t b = bit_start; b < bit_end; b++) {
            size_t byte_idx = b / 8;
            int bit_idx = b % 8;
            bool bit_val = (bytes[byte_idx] >> bit_idx) & 1;

            if (bit_val) {
                /* cur += coeff (simple 256-bit add, no modular reduction needed
                 * since we only accumulate 253 bits) */
                unsigned __int128 carry = 0;
                for (int i = 0; i < 4; i++) {
                    unsigned __int128 sum = (unsigned __int128)cur[i] + coeff[i] + carry;
                    cur[i] = (uint64_t)sum;
                    carry = sum >> 64;
                }
            }
            /* coeff *= 2 */
            uint64_t carry2 = 0;
            for (int i = 0; i < 4; i++) {
                uint64_t new_carry = coeff[i] >> 63;
                coeff[i] = (coeff[i] << 1) | carry2;
                carry2 = new_carry;
            }
        }
        memcpy(out[s], cur, 32);
    }
    *n_out = n_scalars;
}

/* Multipack: pack bytes into Fr scalars (253 bits per scalar, BE bit order).
 * Sprout uses MSB-first bit ordering (bytes_to_bits in Rust). */
void multipack_bytes_to_fr_be(uint64_t (*out)[4], size_t *n_out,
                               const uint8_t *bytes, size_t n_bytes)
{
    /* Sprout Groth16 uses BLS12-381, BE bit order. */
    const size_t capacity = BLS12_381_FR_CAPACITY;
    size_t n_bits = n_bytes * 8;
    size_t n_scalars = (n_bits + capacity - 1) / capacity;

    for (size_t s = 0; s < n_scalars; s++) {
        uint64_t cur[4] = {0, 0, 0, 0};
        uint64_t coeff[4] = {1, 0, 0, 0};

        size_t bit_start = s * capacity;
        size_t bit_end = bit_start + capacity;
        if (bit_end > n_bits) bit_end = n_bits;

        for (size_t b = bit_start; b < bit_end; b++) {
            /* BE bit order: bit b maps to byte b/8, bit (7 - b%8) */
            size_t byte_idx = b / 8;
            int bit_idx = 7 - (int)(b % 8);
            bool bit_val = (bytes[byte_idx] >> bit_idx) & 1;

            if (bit_val) {
                unsigned __int128 carry = 0;
                for (int i = 0; i < 4; i++) {
                    unsigned __int128 sum = (unsigned __int128)cur[i] + coeff[i] + carry;
                    cur[i] = (uint64_t)sum;
                    carry = sum >> 64;
                }
            }
            uint64_t carry2 = 0;
            for (int i = 0; i < 4; i++) {
                uint64_t new_carry = coeff[i] >> 63;
                coeff[i] = (coeff[i] << 1) | carry2;
                carry2 = new_carry;
            }
        }
        memcpy(out[s], cur, 32);
    }
    *n_out = n_scalars;
}

/* Groth16 verification:
 * e(A, B) == e(alpha, beta) * e(vk_x, gamma) * e(C, delta)
 * Rearranged: e(A, B) * e(vk_x, -gamma) * e(C, -delta) == e(alpha, beta) */
bool groth16_verify(const struct groth16_vk *vk,
                    const struct groth16_proof *proof,
                    const uint64_t (*public_inputs)[4],
                    size_t n_inputs)
{
    if (n_inputs + 1 != vk->ic_len)
        LOG_FAIL("groth16",
                 "verify: public input count mismatch: n_inputs=%zu ic_len=%zu (want %zu)",
                 n_inputs, vk->ic_len, vk->ic_len - 1);

    /* Compute vk_x = IC[0] + sum(IC[i+1] * input[i]) */
    struct g1_point vk_x = vk->ic[0];

    for (size_t i = 0; i < n_inputs; i++) {
        /* Skip zero inputs (optimization + avoids identity issues) */
        if (public_inputs[i][0] == 0 && public_inputs[i][1] == 0 &&
            public_inputs[i][2] == 0 && public_inputs[i][3] == 0)
            continue;
        struct g1_point term;
        g1_scalar_mul(&term, &vk->ic[i + 1], public_inputs[i]);
        g1_add(&vk_x, &vk_x, &term);
    }

    /* Negate gamma and delta for the pairing check */
    struct g2_point neg_gamma, neg_delta;
    g2_neg(&neg_gamma, &vk->gamma_g2);
    g2_neg(&neg_delta, &vk->delta_g2);

    /* Check: e(A, B) * e(vk_x, -gamma) * e(C, -delta) * e(-alpha, beta) == 1
     * Equivalent to: e(A, B) == e(alpha, beta) * e(vk_x, gamma) * e(C, delta) */
    struct g1_point neg_alpha;
    g1_neg(&neg_alpha, &vk->alpha_g1);

    struct g1_point pts[4] = { proof->a, vk_x, proof->c, neg_alpha };
    struct g2_point qts[4] = { proof->b, neg_gamma, neg_delta, vk->beta_g2 };

    return bls12_381_multi_pairing_check(pts, qts, 4);
}

/* ===== VK reader (bellman format) ===== */

static bool read_g1_uncompressed(struct g1_point *p, const uint8_t **data, size_t *remaining)
{
    if (*remaining < 96)
        LOG_FAIL("groth16_vk",
                 "read_g1_uncompressed: short buffer (%zu < 96)", *remaining);
    if (!g1_from_uncompressed(p, *data))
        LOG_FAIL("groth16_vk", "read_g1_uncompressed: g1_from_uncompressed failed");
    *data += 96;
    *remaining -= 96;
    return true;
}

static bool read_g2_uncompressed(struct g2_point *p, const uint8_t **data, size_t *remaining)
{
    if (*remaining < 192)
        LOG_FAIL("groth16_vk",
                 "read_g2_uncompressed: short buffer (%zu < 192)", *remaining);
    if (!g2_from_uncompressed(p, *data))
        LOG_FAIL("groth16_vk", "read_g2_uncompressed: g2_from_uncompressed failed");
    *data += 192;
    *remaining -= 192;
    return true;
}

static uint32_t read_be32(const uint8_t **data, size_t *remaining)
{
    uint32_t v = ((uint32_t)(*data)[0] << 24) | ((uint32_t)(*data)[1] << 16) |
                 ((uint32_t)(*data)[2] << 8) | (*data)[3];
    *data += 4;
    *remaining -= 4;
    return v;
}

bool groth16_vk_read_raw(struct groth16_vk *vk, const uint8_t *data, size_t len)
{
    /* VK format: alpha_g1(96) beta_g1(96) beta_g2(192) gamma_g2(192)
     *            delta_g1(96) delta_g2(192) ic_len(4) ic[](96 each)
     * Minimum: 96+96+192+192+96+192+4 = 868 bytes */
    if (len < 868)
        LOG_FAIL("groth16_vk",
                 "vk_read_raw: buffer too small: len=%zu < 868 minimum", len);

    struct g1_point beta_g1, delta_g1;
    if (!read_g1_uncompressed(&vk->alpha_g1, &data, &len))
        LOG_FAIL("groth16_vk", "vk_read_raw: alpha_g1 parse failed");
    if (!read_g1_uncompressed(&beta_g1, &data, &len))  /* beta_g1 unused */
        LOG_FAIL("groth16_vk", "vk_read_raw: beta_g1 parse failed");
    if (!read_g2_uncompressed(&vk->beta_g2, &data, &len))
        LOG_FAIL("groth16_vk", "vk_read_raw: beta_g2 parse failed");
    if (!read_g2_uncompressed(&vk->gamma_g2, &data, &len))
        LOG_FAIL("groth16_vk", "vk_read_raw: gamma_g2 parse failed");
    if (!read_g1_uncompressed(&delta_g1, &data, &len))  /* delta_g1 unused */
        LOG_FAIL("groth16_vk", "vk_read_raw: delta_g1 parse failed");
    if (!read_g2_uncompressed(&vk->delta_g2, &data, &len))
        LOG_FAIL("groth16_vk", "vk_read_raw: delta_g2 parse failed");

    if (len < 4)
        LOG_FAIL("groth16_vk", "vk_read_raw: no room for ic_len u32 (remaining=%zu)", len);
    uint32_t ic_len = read_be32(&data, &len);
    if (ic_len > 1000000)
        LOG_FAIL("groth16_vk",
                 "vk_read_raw: ic_len=%u exceeds 1M sanity cap", ic_len);
    if (len < (size_t)ic_len * 96)
        LOG_FAIL("groth16_vk",
                 "vk_read_raw: short for ic array: remaining=%zu need=%u",
                 len, ic_len * 96);

    vk->ic = zcl_malloc(ic_len * sizeof(struct g1_point), "groth16_vk_ic");
    if (!vk->ic)
        LOG_FAIL("groth16_vk",
                 "vk_read_raw: zcl_malloc failed for ic[] (ic_len=%u)", ic_len);
    vk->ic_len = ic_len;

    for (uint32_t i = 0; i < ic_len; i++) {
        if (!read_g1_uncompressed(&vk->ic[i], &data, &len)) {
            free(vk->ic);
            vk->ic = NULL;
            LOG_FAIL("groth16_vk",
                     "vk_read_raw: ic[%u] parse failed (of %u)", i, ic_len);
        }
    }

    return true;
}

bool groth16_vk_read(struct groth16_vk *vk, const uint8_t *data, size_t len)
{
    /* Parameters file starts with VK, then proving key arrays.
     * We only need the VK portion. */
    return groth16_vk_read_raw(vk, data, len);
}

#pragma GCC diagnostic pop
