// Bit-level unit tests for the DT3 dual-plane ternary format
//
// DT3 packs 128 weights as two ternary planes (w = d1*t1 + d2*t2), each plane
// laid out exactly like tq1_0 with all constants halved. A block with
// misplaced bits still loads and generates, so the layout is checked here
// bit by bit against an independent packer that implements the format
// specification directly.
//
// Extra mode for cross-implementation parity checks (see
// tests/test-dt3-rust-parity.py):
//   test-dt3 --dequant IN.bin OUT.f32
// dequantizes raw DT3 blocks from IN.bin into float32 little-endian OUT.f32.

#include "ggml.h"
#include "ggml-cpu.h"

#undef NDEBUG
#include <assert.h>
#include <inttypes.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <vector>

constexpr int    QK_DT3         = 128;
constexpr size_t DT3_QS_BYTES   = 24; // per plane
constexpr size_t DT3_QH_BYTES   = 2;  // per plane
constexpr size_t DT3_BLOCK_SIZE = 2*DT3_QS_BYTES + 2*DT3_QH_BYTES + 2*sizeof(uint16_t);

// byte offsets inside a block (spec: qs[2][24] | qh[2][2] | d[2])
constexpr size_t OFF_QS = 0;
constexpr size_t OFF_QH = 2*DT3_QS_BYTES;
constexpr size_t OFF_D  = 2*DT3_QS_BYTES + 2*DT3_QH_BYTES;

// Independent packer, written from the format specification (not from
// ggml-quants.c): trits in {-1, 0, 1}, element i of the plane goes to
//   region A: qs[m],    m in [0,16), digit n: elements m + n*16   (0..79)
//   region B: qs[16+m], m in [0,8),  digit n: elements 80 + m + n*8 (80..119)
//   region C: qh[j],    j in [0,2),  digit n: elements 120 + j + n*2 (120..127)
// with the first element in the most significant trit, an extra *3 shift in
// region C, and ceiling division by 243 to fit 5 trits per byte.
static void ref_pack_plane(const int8_t * t, uint8_t * qs, uint8_t * qh) {
    for (int m = 0; m < 16; ++m) {
        uint32_t q = 0;
        for (int n = 0; n < 5; ++n) {
            q = q*3 + (uint32_t)(t[m + n*16] + 1);
        }
        qs[m] = (uint8_t)((q*256 + 242)/243);
    }
    for (int m = 0; m < 8; ++m) {
        uint32_t q = 0;
        for (int n = 0; n < 5; ++n) {
            q = q*3 + (uint32_t)(t[80 + m + n*8] + 1);
        }
        qs[16 + m] = (uint8_t)((q*256 + 242)/243);
    }
    for (int j = 0; j < 2; ++j) {
        uint32_t q = 0;
        for (int n = 0; n < 4; ++n) {
            q = q*3 + (uint32_t)(t[120 + j + n*2] + 1);
        }
        q *= 3; // shift the first value to the most significant trit
        qh[j] = (uint8_t)((q*256 + 242)/243);
    }
}

static void ref_pack_block(const int8_t * t1, float d1, const int8_t * t2, float d2, uint8_t * block) {
    ref_pack_plane(t1, block + OFF_QS,                block + OFF_QH);
    ref_pack_plane(t2, block + OFF_QS + DT3_QS_BYTES, block + OFF_QH + DT3_QH_BYTES);
    const uint16_t h1 = ggml_fp32_to_fp16(d1);
    const uint16_t h2 = ggml_fp32_to_fp16(d2);
    memcpy(block + OFF_D,     &h1, sizeof(h1));
    memcpy(block + OFF_D + 2, &h2, sizeof(h2));
}

// deterministic PRNG so failures are reproducible
static uint32_t rng_state = 0x12345678;
static uint32_t rng_next(void) {
    rng_state ^= rng_state << 13;
    rng_state ^= rng_state >> 17;
    rng_state ^= rng_state << 5;
    return rng_state;
}
static int8_t rng_trit(void) {
    return (int8_t)(rng_next() % 3) - 1;
}

static int test_layout_constants(void) {
    int num_failed = 0;

    if (ggml_blck_size(GGML_TYPE_DT3) != QK_DT3) {
        printf("FAILED: blck_size is %" PRId64 ", expected %d\n", ggml_blck_size(GGML_TYPE_DT3), QK_DT3);
        num_failed++;
    }
    if (ggml_type_size(GGML_TYPE_DT3) != DT3_BLOCK_SIZE) {
        printf("FAILED: type_size is %zu, expected %zu\n", ggml_type_size(GGML_TYPE_DT3), DT3_BLOCK_SIZE);
        num_failed++;
    }
    // 56 bytes / 128 weights = 3.5 bpw exactly
    if (DT3_BLOCK_SIZE*8 != (size_t)QK_DT3*7/2) {
        printf("FAILED: not 3.5 bpw\n");
        num_failed++;
    }

    return num_failed;
}

// a single trit set to -1 or +1 at each position of each plane must come back
// at the same position, scaled by that plane's scale only
static int test_single_trits(const ggml_type_traits * qfns) {
    int num_failed = 0;

    const float d1 = 1.0f;   // exact in fp16
    const float d2 = 0.25f;  // exact in fp16

    for (int plane = 0; plane < 2; ++plane) {
        for (int pos = 0; pos < QK_DT3; ++pos) {
            for (int val = -1; val <= 1; val += 2) {
                int8_t t1[QK_DT3] = {0};
                int8_t t2[QK_DT3] = {0};
                (plane == 0 ? t1 : t2)[pos] = (int8_t)val;

                uint8_t block[DT3_BLOCK_SIZE];
                ref_pack_block(t1, d1, t2, d2, block);

                float out[QK_DT3];
                qfns->to_float(block, out, QK_DT3);

                for (int j = 0; j < QK_DT3; ++j) {
                    const float expected = j == pos ? (plane == 0 ? d1 : d2)*val : 0.0f;
                    if (out[j] != expected) {
                        printf("FAILED: plane %d pos %d val %d: out[%d] = %f, expected %f\n",
                               plane, pos, val, j, out[j], expected);
                        num_failed++;
                    }
                }
            }
        }
    }

    return num_failed;
}

// an all-zero plane packs to qs = 128, qh = 127 everywhere (hand-computed:
// qs: xi = 1 for the 5 trits -> q = 121 -> ceil(121*256/243) = 128,
// qh: q = 40*3 = 120 -> ceil(120*256/243) = 127)
constexpr uint8_t ZERO_QS = 128;
constexpr uint8_t ZERO_QH = 127;

static void build_zero_pattern_block(uint8_t * block, float d1, float d2) {
    memset(block + OFF_QS, ZERO_QS, 2*DT3_QS_BYTES);
    memset(block + OFF_QH, ZERO_QH, 2*DT3_QH_BYTES);
    const uint16_t h1 = ggml_fp32_to_fp16(d1);
    const uint16_t h2 = ggml_fp32_to_fp16(d2);
    memcpy(block + OFF_D,     &h1, sizeof(h1));
    memcpy(block + OFF_D + 2, &h2, sizeof(h2));
}

// check the byte positions of the region boundaries against literal values
// computed by hand from the packing formula (and cross-checked one by one
// against an independent implementation of the spec). Both directions go
// against the LIBRARY code — to_float and from_float — never against a
// packer defined in this file, so a bit-placement bug in ggml-quants.c
// cannot cancel out.
static int test_byte_positions(const ggml_type_traits * qfns, const ggml_type_traits_cpu * qfns_cpu) {
    int num_failed = 0;

    struct byte_case {
        int     pos;   // element with t = -1, all others 0
        size_t  off;   // byte offset inside the plane data (qs: 0..23, qh: 24, 25)
        uint8_t value; // hand-computed byte value at that offset
    };
    const byte_case cases[] = {
        {   0,  0,  43 }, // region A: most significant trit of qs[0]:  q = 121-81 = 40  -> 43
        {  16,  0, 100 }, // region A: second trit of qs[0]:            q = 121-27 = 94  -> 100
        {  79, 15, 127 }, // region A boundary: last trit of qs[15]:    q = 121-1  = 120 -> 127
        {  80, 16,  43 }, // region B boundary: first trit of qs[16]:   q = 40           -> 43
        { 119, 23, 127 }, // region B boundary: last trit of qs[23]:    q = 120          -> 127
        { 120, 24,  42 }, // region C boundary: first trit of qh[0]:    q = (40-27)*3    -> 42
        { 127, 25, 124 }, // region C: last element, qh[1]:             q = (40-1)*3     -> 124
    };

    for (size_t c = 0; c < sizeof(cases)/sizeof(cases[0]); ++c) {
        const int     pos = cases[c].pos;
        const size_t  off = cases[c].off;
        const uint8_t val = cases[c].value;

        // unpack direction: writing the literal byte must put -1 at exactly
        // element `pos` of the corresponding plane
        for (int plane = 0; plane < 2; ++plane) {
            uint8_t block[DT3_BLOCK_SIZE];
            build_zero_pattern_block(block, 1.0f, 0.25f);
            if (off < DT3_QS_BYTES) {
                block[OFF_QS + plane*DT3_QS_BYTES + off] = val;
            } else {
                block[OFF_QH + plane*DT3_QH_BYTES + (off - DT3_QS_BYTES)] = val;
            }

            float out[QK_DT3];
            qfns->to_float(block, out, QK_DT3);

            const float dp = plane == 0 ? 1.0f : 0.25f;
            for (int j = 0; j < QK_DT3; ++j) {
                const float expected = j == pos ? -dp : 0.0f;
                if (out[j] != expected) {
                    printf("FAILED: byte 0x%02x at plane %d offset %zu: out[%d] = %f, expected %f\n",
                           val, plane, off, j, out[j], expected);
                    num_failed++;
                }
            }
        }

        // pack direction: quantizing a single -1 must produce exactly the
        // literal byte at the right offset of plane 0 (plane 1 is the zero
        // pattern since there is no residual), and nothing else
        {
            float x[QK_DT3] = {0.0f};
            x[pos] = -1.0f;

            uint8_t block[DT3_BLOCK_SIZE];
            qfns_cpu->from_float(x, block, QK_DT3);

            uint8_t expected[DT3_BLOCK_SIZE];
            build_zero_pattern_block(expected, 1.0f, 0.0f);
            if (off < DT3_QS_BYTES) {
                expected[OFF_QS + off] = val;
            } else {
                expected[OFF_QH + (off - DT3_QS_BYTES)] = val;
            }

            for (size_t b = 0; b < DT3_BLOCK_SIZE; ++b) {
                if (block[b] != expected[b]) {
                    printf("FAILED: pack of -1 at element %d: byte %zu is 0x%02x, expected 0x%02x\n",
                           pos, b, block[b], expected[b]);
                    num_failed++;
                }
            }
        }
    }

    // all-zero input packs to the all-zero pattern with zero scales
    {
        const float x[QK_DT3] = {0.0f};

        uint8_t block[DT3_BLOCK_SIZE];
        qfns_cpu->from_float(x, block, QK_DT3);

        uint8_t expected[DT3_BLOCK_SIZE];
        build_zero_pattern_block(expected, 0.0f, 0.0f);

        for (size_t b = 0; b < DT3_BLOCK_SIZE; ++b) {
            if (block[b] != expected[b]) {
                printf("FAILED: pack of all-zero block: byte %zu is 0x%02x, expected 0x%02x\n",
                       b, block[b], expected[b]);
                num_failed++;
            }
        }
    }

    return num_failed;
}

// ggml_validate_row_data must accept exactly the reachable byte values:
// 243 of 256 in qs, 81 of 256 in qh (4 trits + an always-zero padding digit)
static int test_validate(void) {
    int num_failed = 0;

    bool valid_qs[256] = {false};
    bool valid_qh[256] = {false};
    for (uint32_t q = 0; q < 243; ++q) {
        valid_qs[(q*256 + 242)/243] = true;
        if (q % 3 == 0 && q <= 240) {
            valid_qh[(q*256 + 242)/243] = true;
        }
    }

    printf("(invalid-byte messages below are expected, the validator is being probed)\n");

    for (int b = 0; b < 256; ++b) {
        uint8_t block[DT3_BLOCK_SIZE];

        build_zero_pattern_block(block, 1.0f, 0.25f);
        block[OFF_QS] = (uint8_t) b;
        if (ggml_validate_row_data(GGML_TYPE_DT3, block, DT3_BLOCK_SIZE) != valid_qs[b]) {
            printf("FAILED: validate qs byte 0x%02x: expected %s\n", b, valid_qs[b] ? "valid" : "invalid");
            num_failed++;
        }

        build_zero_pattern_block(block, 1.0f, 0.25f);
        block[OFF_QH + DT3_QH_BYTES] = (uint8_t) b; // plane 1 qh, to also cover the second plane
        if (ggml_validate_row_data(GGML_TYPE_DT3, block, DT3_BLOCK_SIZE) != valid_qh[b]) {
            printf("FAILED: validate qh byte 0x%02x: expected %s\n", b, valid_qh[b] ? "valid" : "invalid");
            num_failed++;
        }
    }

    // a block filled with 0xaa must be rejected (0xaa is a reachable qs code
    // but violates the qh padding-digit invariant)
    {
        uint8_t block[DT3_BLOCK_SIZE];
        memset(block, 0xaa, DT3_BLOCK_SIZE);
        const uint16_t h = ggml_fp32_to_fp16(1.0f);
        memcpy(block + OFF_D,     &h, sizeof(h));
        memcpy(block + OFF_D + 2, &h, sizeof(h));
        if (ggml_validate_row_data(GGML_TYPE_DT3, block, DT3_BLOCK_SIZE)) {
            printf("FAILED: validate accepted a block filled with 0xaa\n");
            num_failed++;
        }
    }

    // packed real data must pass
    {
        int8_t t1[QK_DT3];
        int8_t t2[QK_DT3];
        for (int j = 0; j < QK_DT3; ++j) {
            t1[j] = rng_trit();
            t2[j] = rng_trit();
        }
        uint8_t block[DT3_BLOCK_SIZE];
        ref_pack_block(t1, -0.5f, t2, 0.125f, block);
        if (!ggml_validate_row_data(GGML_TYPE_DT3, block, DT3_BLOCK_SIZE)) {
            printf("FAILED: validate rejected a well-formed block\n");
            num_failed++;
        }
    }

    return num_failed;
}

// random trits and scales (negative scales included) must round-trip exactly
static int test_roundtrip(const ggml_type_traits * qfns) {
    int num_failed = 0;

    const float scales[][2] = {
        { 1.0f,  0.25f  },
        { 0.5f, -0.125f }, // negative second plane
        {-2.0f,  0.75f  }, // negative first plane
        { 0.0f,  0.0f   }, // all-zero scales
    };

    for (size_t sc = 0; sc < sizeof(scales)/sizeof(scales[0]); ++sc) {
        for (int rep = 0; rep < 64; ++rep) {
            int8_t t1[QK_DT3];
            int8_t t2[QK_DT3];
            for (int j = 0; j < QK_DT3; ++j) {
                t1[j] = rng_trit();
                t2[j] = rng_trit();
            }

            const float d1 = ggml_fp16_to_fp32(ggml_fp32_to_fp16(scales[sc][0]));
            const float d2 = ggml_fp16_to_fp32(ggml_fp32_to_fp16(scales[sc][1]));

            uint8_t block[DT3_BLOCK_SIZE];
            ref_pack_block(t1, d1, t2, d2, block);

            float out[QK_DT3];
            qfns->to_float(block, out, QK_DT3);

            for (int j = 0; j < QK_DT3; ++j) {
                const float expected = d1*t1[j] + d2*t2[j];
                if (out[j] != expected) {
                    printf("FAILED: roundtrip scales (%f, %f) rep %d: out[%d] = %f, expected %f\n",
                           d1, d2, rep, j, out[j], expected);
                    num_failed++;
                }
            }
        }
    }

    return num_failed;
}

// the in-tree quantizer must produce the same bytes as the independent packer
// when the input is already exactly ternary (plane 1 = input, plane 2 = 0)
static int test_quantize_pack_parity(const ggml_type_traits_cpu * qfns_cpu) {
    int num_failed = 0;

    for (int rep = 0; rep < 64; ++rep) {
        int8_t t1[QK_DT3];
        const int8_t t2[QK_DT3] = {0};
        float  x[QK_DT3];
        for (int j = 0; j < QK_DT3; ++j) {
            t1[j] = rng_trit();
            x[j]  = (float) t1[j];
        }
        // make sure the block is not all zeros so that d1 == 1.0
        t1[0] = 1;
        x[0]  = 1.0f;

        uint8_t expected[DT3_BLOCK_SIZE];
        ref_pack_block(t1, 1.0f, t2, 0.0f, expected);

        uint8_t block[DT3_BLOCK_SIZE];
        qfns_cpu->from_float(x, block, QK_DT3);

        if (memcmp(block, expected, DT3_BLOCK_SIZE) != 0) {
            for (size_t b = 0; b < DT3_BLOCK_SIZE; ++b) {
                if (block[b] != expected[b]) {
                    printf("FAILED: quantize pack parity rep %d: byte %zu is 0x%02x, expected 0x%02x\n",
                           rep, b, block[b], expected[b]);
                }
            }
            num_failed++;
        }
    }

    return num_failed;
}

// vec_dot against a hand-made sum over the KNOWN trits (not against our own
// dequantization): sum_i y_i * (d1*t1_i + d2*t2_i) with y from q8_0's own
// to_float. Random trits make the qh bytes non-trivial, which would expose a
// vectorization that reads the padding 5th trit of the qh bytes.
static int test_vec_dot(const ggml_type_traits_cpu * qfns_cpu) {
    int num_failed = 0;

    const auto * vdot_traits = ggml_get_type_traits_cpu(qfns_cpu->vec_dot_type);
    const auto * vdot_qfns   = ggml_get_type_traits(qfns_cpu->vec_dot_type);

    if (qfns_cpu->vec_dot_type != GGML_TYPE_Q8_0) {
        printf("FAILED: vec_dot_type is %s, expected q8_0\n", ggml_type_name(qfns_cpu->vec_dot_type));
        return 1;
    }

    const int nblocks = 4;
    const int n       = nblocks*QK_DT3;

    for (int rep = 0; rep < 64; ++rep) {
        std::vector<int8_t> t1(n);
        std::vector<int8_t> t2(n);
        std::vector<float>  d1(nblocks);
        std::vector<float>  d2(nblocks);
        std::vector<uint8_t> xq(nblocks*DT3_BLOCK_SIZE);

        for (int i = 0; i < nblocks; ++i) {
            for (int j = 0; j < QK_DT3; ++j) {
                t1[i*QK_DT3 + j] = rng_trit();
                t2[i*QK_DT3 + j] = rng_trit();
            }
            // fp16-exact scales of both signs
            d1[i] = (float)((int)(rng_next() % 9) - 4) * 0.25f;
            d2[i] = (float)((int)(rng_next() % 9) - 4) * 0.0625f;
            ref_pack_block(&t1[i*QK_DT3], d1[i], &t2[i*QK_DT3], d2[i], &xq[i*DT3_BLOCK_SIZE]);
        }

        std::vector<float> y(n);
        for (int j = 0; j < n; ++j) {
            y[j] = 0.1f + 2.0f*cosf((float)(j + rep));
        }

        std::vector<uint8_t> yq(ggml_row_size(qfns_cpu->vec_dot_type, n));
        vdot_traits->from_float(y.data(), yq.data(), n);

        // exact values the integer path sees
        std::vector<float> ydq(n);
        vdot_qfns->to_float(yq.data(), ydq.data(), n);

        double ref = 0.0;
        for (int i = 0; i < nblocks; ++i) {
            for (int j = 0; j < QK_DT3; ++j) {
                const int ij = i*QK_DT3 + j;
                ref += (double)ydq[ij] * ((double)d1[i]*t1[ij] + (double)d2[i]*t2[ij]);
            }
        }

        float result = INFINITY;
        qfns_cpu->vec_dot(n, &result, 0, xq.data(), 0, yq.data(), 0, 1);

        const float err = fabsf(result - (float)ref);
        const float tol = 1e-4f * (float)n;
        if (!(err <= tol)) {
            printf("FAILED: vec_dot rep %d: got %f, expected %f (err %f)\n", rep, result, (float)ref, err);
            num_failed++;
        }
    }

    return num_failed;
}

// the scalar reference implementation. The symbol only exists on builds with
// a native DT3 kernel: without one, arch-fallback.h renames the generic to
// ggml_vec_dot_dt3_q8_0 and there is nothing to compare against, so the
// reference is declared weak and the parity test skips when it is absent.
#if defined(_MSC_VER)
#define DT3_NO_WEAK_SYMBOLS
#else
extern "C" void ggml_vec_dot_dt3_q8_0_generic(int n, float * s, size_t bs, const void * vx, size_t bx, const void * vy, size_t by, int nrc) __attribute__((weak));
#endif

// the dispatched (possibly vectorized) vec_dot must match the generic scalar
// implementation exactly — the actual function is called, not a re-derivation
// of it. Blocks exercise all three regions, the 79/80 and 119/120 boundaries,
// non-trivial qh bytes (would expose reading their padding 5th digit), and
// scales of both and mixed signs.
static int test_vec_dot_arch_parity(const ggml_type_traits_cpu * qfns_cpu) {
    int num_failed = 0;

#if defined(DT3_NO_WEAK_SYMBOLS)
    (void) qfns_cpu;
    printf("(skipping vec_dot arch parity: no weak symbol support)\n");
    return num_failed;
#else
    if (ggml_vec_dot_dt3_q8_0_generic == nullptr) {
        printf("(skipping vec_dot arch parity: this build has no separate generic vec_dot)\n");
        return num_failed;
    }

    const auto * vdot_traits = ggml_get_type_traits_cpu(qfns_cpu->vec_dot_type);

    const int nblocks = 3;
    const int n       = nblocks*QK_DT3;

    const float scale_cases[][2] = {
        { 1.0f,   0.25f   },
        { 0.5f,  -0.125f  }, // negative second plane
        {-2.0f,   0.75f   }, // negative first plane
        {-0.75f, -0.0625f }, // both negative
        { 0.0f,   1.0f    }, // dead first plane
    };
    const int n_scale_cases = (int)(sizeof(scale_cases)/sizeof(scale_cases[0]));

    // pattern 0: fully random trits
    // pattern 1: zero everywhere except elements 120..127 (qh-only)
    // pattern 2: single +1/-1 walking over the region boundaries
    const int boundary_pos[] = { 0, 79, 80, 119, 120, 127 };

    for (int rep = 0; rep < 96; ++rep) {
        std::vector<uint8_t> xq(nblocks*DT3_BLOCK_SIZE);

        for (int i = 0; i < nblocks; ++i) {
            int8_t t1[QK_DT3] = {0};
            int8_t t2[QK_DT3] = {0};

            const int pattern = rep % 3;
            if (pattern == 0) {
                for (int j = 0; j < QK_DT3; ++j) {
                    t1[j] = rng_trit();
                    t2[j] = rng_trit();
                }
            } else if (pattern == 1) {
                for (int j = 120; j < QK_DT3; ++j) {
                    t1[j] = rng_trit();
                    t2[j] = rng_trit();
                }
            } else {
                const int pos = boundary_pos[rep/3 % 6];
                t1[pos] = (rep & 1) ? 1 : -1;
                t2[QK_DT3 - 1 - pos] = (rep & 1) ? -1 : 1;
            }

            const float * sc = scale_cases[(rep + i) % n_scale_cases];
            ref_pack_block(t1, sc[0], t2, sc[1], &xq[i*DT3_BLOCK_SIZE]);
        }

        std::vector<float> y(n);
        for (int j = 0; j < n; ++j) {
            // reach the full q8_0 range, both signs
            y[j] = 127.0f*sinf(0.7f*(float)(j + 13*rep)) + 0.5f*cosf((float)j);
        }
        std::vector<uint8_t> yq(ggml_row_size(qfns_cpu->vec_dot_type, n));
        vdot_traits->from_float(y.data(), yq.data(), n);

        float res_arch    = INFINITY;
        float res_generic = -INFINITY;
        qfns_cpu->vec_dot(n, &res_arch, 0, xq.data(), 0, yq.data(), 0, 1);
        ggml_vec_dot_dt3_q8_0_generic(n, &res_generic, 0, xq.data(), 0, yq.data(), 0, 1);

        if (memcmp(&res_arch, &res_generic, sizeof(float)) != 0) {
            printf("FAILED: vec_dot arch parity rep %d: arch %.9g != generic %.9g\n",
                   rep, res_arch, res_generic);
            num_failed++;
        }
    }

    return num_failed;
#endif
}

// --dequant IN.bin OUT.f32 : dequantize raw DT3 blocks, for parity checks
// against external packers (ternaria's Rust pack_dt3)
static int run_dequant_file(const char * in_path, const char * out_path) {
    FILE * fin = fopen(in_path, "rb");
    if (!fin) {
        fprintf(stderr, "error: cannot open %s\n", in_path);
        return 1;
    }
    fseek(fin, 0, SEEK_END);
    const long size = ftell(fin);
    fseek(fin, 0, SEEK_SET);
    if (size <= 0 || size % DT3_BLOCK_SIZE != 0) {
        fprintf(stderr, "error: %s size %ld is not a multiple of %zu\n", in_path, size, DT3_BLOCK_SIZE);
        fclose(fin);
        return 1;
    }
    std::vector<uint8_t> data(size);
    if (fread(data.data(), 1, size, fin) != (size_t)size) {
        fprintf(stderr, "error: short read on %s\n", in_path);
        fclose(fin);
        return 1;
    }
    fclose(fin);

    const int64_t nel = (int64_t)(size/DT3_BLOCK_SIZE)*QK_DT3;
    std::vector<float> out(nel);
    ggml_get_type_traits(GGML_TYPE_DT3)->to_float(data.data(), out.data(), nel);

    FILE * fout = fopen(out_path, "wb");
    if (!fout) {
        fprintf(stderr, "error: cannot open %s\n", out_path);
        return 1;
    }
    fwrite(out.data(), sizeof(float), nel, fout);
    fclose(fout);

    return 0;
}

int main(int argc, char * argv[]) {
    if (argc == 4 && strcmp(argv[1], "--dequant") == 0) {
        return run_dequant_file(argv[2], argv[3]);
    }
    if (argc != 1) {
        fprintf(stderr, "usage: %s [--dequant IN.bin OUT.f32]\n", argv[0]);
        return 1;
    }

    ggml_cpu_init();

    const auto * qfns     = ggml_get_type_traits(GGML_TYPE_DT3);
    const auto * qfns_cpu = ggml_get_type_traits_cpu(GGML_TYPE_DT3);

    int num_failed = 0;

    num_failed += test_layout_constants();
    num_failed += test_single_trits(qfns);
    num_failed += test_byte_positions(qfns, qfns_cpu);
    num_failed += test_validate();
    num_failed += test_roundtrip(qfns);
    num_failed += test_quantize_pack_parity(qfns_cpu);
    num_failed += test_vec_dot(qfns_cpu);
    num_failed += test_vec_dot_arch_parity(qfns_cpu);

    printf("%d tests failed\n", num_failed);

    return num_failed > 0;
}