/*
 * pln_bidirectional_production.c
 * =================================================
 * PRODUCTION-READY RESONANT PLN BIDIRECTIONAL TOOL
 *
 * Features:
 *   • Full command-line interface (arbitrary primes + PLN coords)
 *   • Deterministic forward (any integer → PLN)
 *   • Deterministic multi-pass reverse (PLN → nearest resonant prime)
 *   • 4096-bit BigUint (up to ~1230 decimal digits)
 *   • Strong Miller-Rabin big-int engine (64-bit exact, strong probable prime beyond)
 *   • Fixed buffer overflows (2048-digit safety)
 *   • Roundtrip verification mode
 *   • Robust input handling + usage/help
 *
 * Compile:
 *   gcc -O3 -Wall -Wextra -std=c99 -o pln_tool pln_bidirectional_production.c -lm
 *
 * Run examples:
 *   ./pln_tool                              # built-in demo
 *   ./pln_tool --prime 18446744073709551557
 *   ./pln_tool --pln 2 0.40088677849 1.4165141186 9.40088677849
 *   ./pln_tool --roundtrip 18446744073709551557
 *   ./pln_tool --help
 */
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <math.h>
#include <inttypes.h>
#include <stdbool.h>
#include <stdarg.h>
#include <string.h>
#ifdef _WIN32
#include <windows.h>
#else
#include <dlfcn.h>
#endif

static bool debug_mode = false;

static void debug_printf(const char *fmt, ...) {
    if (!debug_mode) return;
    va_list args;
    va_start(args, fmt);
    vprintf(fmt, args);
    va_end(args);
}

typedef struct PLNAddress PLNAddress;

/* ====================== HDGL / Dₙ(r) LOADING ====================== */
#ifdef _WIN32
typedef HMODULE HDGLHandle;
#else
typedef void *HDGLHandle;
#endif
static HDGLHandle hdgl_handle = NULL;
static double (*hdgl_compute_Dn_r)(int, double, double) = NULL;
static bool hdgl_enabled = false;

static bool load_hdgl_library(const char *path) {
    if (!path || !*path) return false;
#ifdef _WIN32
    const char *try_path = path;
    char alt_path[1024] = {0};
    const char *dot = strrchr(path, '.');
    if (!dot || strchr(dot, '\\') || strchr(dot, '/')) {
        /* no extension, append .dll */
        snprintf(alt_path, sizeof(alt_path), "%s.dll", path);
    } else if (strcmp(dot, ".so") == 0) {
        strncpy(alt_path, path, sizeof(alt_path) - 1);
        strcpy(alt_path + (dot - path), ".dll");
    }

    hdgl_handle = LoadLibraryA(try_path);
    if (!hdgl_handle && alt_path[0]) {
        hdgl_handle = LoadLibraryA(alt_path);
        try_path = alt_path;
    }
    if (!hdgl_handle) {
        fprintf(stderr, "Error loading HDGL library '%s'\n", path);
        if (alt_path[0]) fprintf(stderr, "Tried fallback '%s'\n", alt_path);
        return false;
    }
    hdgl_compute_Dn_r = (double (*)(int,double,double))GetProcAddress(hdgl_handle, "compute_Dn_r");
#else
    hdgl_handle = dlopen(path, RTLD_NOW);
    if (!hdgl_handle) {
        fprintf(stderr, "Error loading HDGL library '%s': %s\n", path, dlerror());
        return false;
    }
    hdgl_compute_Dn_r = (double (*)(int,double,double))dlsym(hdgl_handle, "compute_Dn_r");
#endif
    if (!hdgl_compute_Dn_r) {
        fprintf(stderr, "Error: compute_Dn_r not found in HDGL library '%s'\n", path);
#ifdef _WIN32
        FreeLibrary(hdgl_handle);
#else
        dlclose(hdgl_handle);
#endif
        hdgl_handle = NULL;
        return false;
    }
    hdgl_enabled = true;
    return true;
}

static void unload_hdgl_library(void) {
    if (!hdgl_handle) return;
#ifdef _WIN32
    FreeLibrary(hdgl_handle);
#else
    dlclose(hdgl_handle);
#endif
    hdgl_handle = NULL;
    hdgl_compute_Dn_r = NULL;
    hdgl_enabled = false;
}


/* ====================== CONSTANTS ====================== */
#define PHI         1.6180339887498948482L
#define INV_PHI     0.6180339887498948482L
#define LN_PHI      0.4812118250596034748L
#define LN2         0.6931471805599453094L
#define LOG2_10     3.32192809488736234787L
#define LOG10_5     0.698970004336018804786L
#define LOG10_PHI   0.208987640249978733769L
#define M_PI_L      3.14159265358979323846L
#define TWO_PI_L    6.2831853071795864769L
#define R_MIN_THRESHOLD 0.15L
#define MAX_PASSES  5
#define MAX_FILE_LEAD_DIGITS 18
#define PLN_SCORE_PASS_THRESHOLD 0.60L
#define PLN_SCORE_BORDER_THRESHOLD 0.55L

/* ====================== FULL BIGUINT (4096-bit) ====================== */
#define BIGUINT_MAX_WORDS 130

typedef struct {
    int len;
    uint32_t digits[BIGUINT_MAX_WORDS];
} BigUint;

static void bu_normalize(BigUint *a);
static void bu_copy(BigUint *dest, const BigUint *src);
static int bu_bit_length(const BigUint *a);
static void bu_shr1(BigUint *a);
static void bu_init(BigUint *a);
static bool bu_is_zero(const BigUint *a);
static void bu_shl_bits(BigUint *a, int bits);
static void bu_shr_bits(BigUint *a, int bits);
static bool prime_file_to_log2(const char *path, long double *out_log2x, size_t *out_digits);
static bool find_prime_from_pln(int target_D, long double target_theta_pi, long double target_Q, long double target_omega, int max_steps, char *found_str, size_t buf_size);
typedef struct PLNAddress PLNAddress;

/* ====================== SLOT4096 APA FLOAT ====================== */
#ifndef SUPER_RESOLUTION_MODE
#define SUPER_RESOLUTION_MODE 0
#endif

#if SUPER_RESOLUTION_MODE
#define APA_MANTISSA_WORDS 262144
#else
#define APA_MANTISSA_WORDS 64
#endif
#define APA_WORD_BITS 64
#define APA_MANTISSA_BITS (APA_MANTISSA_WORDS * APA_WORD_BITS)
#define APA_MSB_MASK 0x8000000000000000ULL

typedef struct {
    uint64_t mantissa[APA_MANTISSA_WORDS];
    int64_t exponent;
    int8_t sign;
} APAFloat;

static void apa_shift_right(uint64_t *mantissa_words, size_t num_words, int64_t shift_amount) {
    if (shift_amount <= 0 || num_words == 0) return;

    if (shift_amount >= (int64_t)num_words * 64) {
        memset(mantissa_words, 0, num_words * sizeof(uint64_t));
        return;
    }

    int64_t word_shift = shift_amount / 64;
    int bit_shift = (int)(shift_amount % 64);

    if (word_shift > 0) {
        for (size_t i = num_words; i-- > (size_t)word_shift; ) {
            mantissa_words[i] = mantissa_words[i - word_shift];
        }
        memset(mantissa_words, 0, word_shift * sizeof(uint64_t));
    }

    if (bit_shift > 0) {
        int reverse_shift = 64 - bit_shift;
        for (size_t i = num_words; i-- > 0; ) {
            uint64_t upper_carry = 0;
            if (i > 0) {
                upper_carry = mantissa_words[i - 1] << reverse_shift;
            }
            mantissa_words[i] = (mantissa_words[i] >> bit_shift) | upper_carry;
        }
    }
}

static void apa_normalize(APAFloat *num) {
    if (num->mantissa[0] == 0ULL) {
        num->exponent = 0;
        return;
    }

    while (!(num->mantissa[0] & APA_MSB_MASK)) {
        uint64_t carry = 0;
        for (size_t i = APA_MANTISSA_WORDS; i-- > 0; ) {
            uint64_t next_carry = (num->mantissa[i] & APA_MSB_MASK) ? 1 : 0;
            num->mantissa[i] = (num->mantissa[i] << 1) | carry;
            carry = next_carry;
        }
        num->exponent--;
    }
}

static void apa_init(APAFloat *num, long double value) {
    memset(num->mantissa, 0, sizeof(num->mantissa));
    num->sign = (value >= 0.0L) ? 1 : -1;
    value = fabsl(value);
    if (value == 0.0L) {
        num->exponent = 0;
        return;
    }
    int exp;
    long double mant = frexpl(value, &exp);
    num->exponent = exp;
    mant *= 2.0L;
    for (int i = 0; i < APA_MANTISSA_WORDS * 64; i++) {
        mant *= 2.0L;
        if (mant >= 2.0L) {
            int word_idx = i / 64;
            int bit_idx = 63 - (i % 64);
            if (word_idx < APA_MANTISSA_WORDS) {
                num->mantissa[word_idx] |= (1ULL << bit_idx);
            }
            mant -= 2.0L;
        }
    }
    apa_normalize(num);
}

static double apa_to_double(const APAFloat *num) {
    if (num->mantissa[0] == 0ULL) return 0.0;
    long double result = 0.0L;
    long double weight = 0.5L;
    for (int i = 0; i < 64; i++) {
        int word_idx = i / 64;
        int bit_idx = 63 - (i % 64);
        if (word_idx < APA_MANTISSA_WORDS && (num->mantissa[word_idx] & (1ULL << bit_idx))) {
            result += weight;
        }
        weight *= 0.5L;
    }
    result = ldexpl(result, num->exponent);
    return num->sign * (double)result;
}

static void apa_multiply(APAFloat *result, const APAFloat *a, const APAFloat *b) {
    result->sign = a->sign * b->sign;
    result->exponent = a->exponent + b->exponent;

    uint64_t temp[APA_MANTISSA_WORDS * 2];
    memset(temp, 0, sizeof(temp));

    for (int i = 0; i < 8 && i < APA_MANTISSA_WORDS; i++) {
        uint64_t carry = 0;
        for (int j = 0; j < 8 && j < APA_MANTISSA_WORDS; j++) {
            __uint128_t prod = (__uint128_t)a->mantissa[i] * (__uint128_t)b->mantissa[j];
            __uint128_t sum = (__uint128_t)temp[i + j] + prod + carry;
            temp[i + j] = (uint64_t)sum;
            carry = (uint64_t)(sum >> 64);
        }
        if (i + 8 < APA_MANTISSA_WORDS * 2) {
            temp[i + 8] += carry;
        }
    }

    int shift = 0;
    while (shift < 64 && !(temp[0] & (1ULL << (63 - shift)))) {
        shift++;
    }
    result->exponent -= shift;
    for (int i = 0; i < APA_MANTISSA_WORDS; i++) {
        result->mantissa[i] = temp[i];
    }
}

static void apa_add(APAFloat *result, const APAFloat *a, const APAFloat *b) {
    APAFloat temp_a = *a;
    APAFloat temp_b = *b;

    if (temp_a.exponent > temp_b.exponent) {
        apa_shift_right(temp_b.mantissa, APA_MANTISSA_WORDS, temp_a.exponent - temp_b.exponent);
        temp_b.exponent = temp_a.exponent;
        *result = temp_a;
    } else {
        apa_shift_right(temp_a.mantissa, APA_MANTISSA_WORDS, temp_b.exponent - temp_a.exponent);
        temp_a.exponent = temp_b.exponent;
        *result = temp_b;
    }

    uint64_t carry = 0;
    for (size_t i = APA_MANTISSA_WORDS; i-- > 0; ) {
        __uint128_t sum = (__uint128_t)result->mantissa[i] + temp_a.mantissa[i] + temp_b.mantissa[i] + carry;
        result->mantissa[i] = (uint64_t)sum;
        carry = (uint64_t)(sum >> 64);
    }

    if (carry) {
        result->exponent += 1;
        apa_shift_right(result->mantissa, APA_MANTISSA_WORDS, 1);
        result->mantissa[0] |= APA_MSB_MASK;
    }

    apa_normalize(result);
}

static void apa_sqrt(APAFloat *result, const APAFloat *num) {
    double val = apa_to_double(num);
    double root = sqrt(fabs(val));
    apa_init(result, root);
    result->sign = 1;
}

static long double apa_log2(const APAFloat *num) {
    if (num->mantissa[0] == 0ULL) return -INFINITY;
    long double high = (long double)num->mantissa[0];
    long double low = 0.0L;
    if (APA_MANTISSA_WORDS > 1) {
        low = (long double)num->mantissa[1] / 18446744073709551616.0L;
    }
    long double top = (high + low) / 9223372036854775808.0L;
    return (long double)num->exponent - 1 + log2l(top);
}

static void bu_shl1(BigUint *a) {
    uint64_t carry = 0;
    for (int i = 0; i < a->len || carry; i++) {
        uint64_t value = (uint64_t)(i < a->len ? a->digits[i] : 0) << 1 | carry;
        if (i < BIGUINT_MAX_WORDS) a->digits[i] = (uint32_t)value;
        carry = value >> 32;
    }
    if (carry && a->len < BIGUINT_MAX_WORDS) a->digits[a->len++] = (uint32_t)carry;
    bu_normalize(a);
}

static void bu_shl_bits(BigUint *a, int bits) {
    if (bits <= 0 || bu_is_zero(a)) return;
    int word_shift = bits / 32;
    int bit_shift = bits % 32;
    if (word_shift > 0) {
        if (a->len + word_shift > BIGUINT_MAX_WORDS) {
            word_shift = BIGUINT_MAX_WORDS - a->len;
        }
        for (int i = a->len - 1; i >= 0; i--) {
            a->digits[i + word_shift] = a->digits[i];
        }
        for (int i = 0; i < word_shift; i++) {
            a->digits[i] = 0;
        }
        a->len += word_shift;
    }
    if (bit_shift > 0) {
        uint64_t carry = 0;
        for (int i = word_shift; i < a->len; i++) {
            uint64_t value = ((uint64_t)a->digits[i] << bit_shift) | carry;
            a->digits[i] = (uint32_t)value;
            carry = value >> 32;
        }
        if (carry && a->len < BIGUINT_MAX_WORDS) {
            a->digits[a->len++] = (uint32_t)carry;
        }
    }
    bu_normalize(a);
}

static void bu_shr_bits(BigUint *a, int bits) {
    if (bits <= 0 || bu_is_zero(a)) return;
    int word_shift = bits / 32;
    int bit_shift = bits % 32;
    if (word_shift >= a->len) {
        bu_init(a);
        return;
    }
    if (word_shift > 0) {
        for (int i = 0; i < a->len - word_shift; i++) {
            a->digits[i] = a->digits[i + word_shift];
        }
        for (int i = a->len - word_shift; i < a->len; i++) {
            a->digits[i] = 0;
        }
        a->len -= word_shift;
    }
    if (bit_shift > 0) {
        uint64_t carry = 0;
        for (int i = a->len - 1; i >= 0; i--) {
            uint64_t value = ((uint64_t)a->digits[i] >> bit_shift) | (carry << (32 - bit_shift));
            carry = a->digits[i] & (((uint64_t)1 << bit_shift) - 1);
            a->digits[i] = (uint32_t)value;
        }
    }
    bu_normalize(a);
}

static APAFloat bu_to_apa(const BigUint *a) {
    APAFloat out;
    memset(&out, 0, sizeof(out));
    out.sign = 1;
    int bits = bu_bit_length(a);
    if (bits == 0) {
        out.exponent = 0;
        return out;
    }

    BigUint temp;
    bu_copy(&temp, a);
    int target_bits = APA_MANTISSA_WORDS * 64;
    int shift = bits - target_bits;
    if (shift > 0) {
        bu_shr_bits(&temp, shift);
    } else if (shift < 0) {
        bu_shl_bits(&temp, -shift);
    }

    out.exponent = bits;
    int word_index = 0;
    if (bits < target_bits) {
        debug_printf("DEBUG bu_to_apa bits=%d target_bits=%d temp.len=%d top_word=%08x\n", bits, target_bits, temp.len, temp.digits[temp.len-1]);
    }
    for (int i = temp.len - 1; i >= 0; i -= 2) {
        uint64_t hi = temp.digits[i];
        uint64_t lo = (i - 1 >= 0) ? temp.digits[i - 1] : 0ULL;
        if (word_index < APA_MANTISSA_WORDS) {
            out.mantissa[word_index++] = (hi << 32) | lo;
        }
    }

    apa_normalize(&out);
    return out;
}
static void bu_init(BigUint *a) { a->len = 0; memset(a->digits, 0, sizeof(a->digits)); }
static void bu_normalize(BigUint *a) { while (a->len > 0 && a->digits[a->len-1]==0) a->len--; }
static void bu_from_u64(BigUint *a, uint64_t value) {
    if (value == 0) { a->len = 0; return; }
    a->digits[0] = (uint32_t)(value & 0xFFFFFFFFu);
    a->digits[1] = (uint32_t)(value >> 32);
    a->len = a->digits[1] ? 2 : 1;
}
static void bu_copy(BigUint *dest, const BigUint *src) {
    dest->len = src->len;
    memcpy(dest->digits, src->digits, src->len * sizeof(uint32_t));
}
static int bu_cmp(const BigUint *a, const BigUint *b) {
    if (a->len != b->len) return a->len < b->len ? -1 : 1;
    for (int i = a->len - 1; i >= 0; i--)
        if (a->digits[i] != b->digits[i]) return a->digits[i] < b->digits[i] ? -1 : 1;
    return 0;
}
static bool bu_is_zero(const BigUint *a) { return a->len == 0; }
static bool bu_is_even(const BigUint *a) { return a->len == 0 || ((a->digits[0] & 1u) == 0); }
static int bu_bit_length(const BigUint *a) {
    if (a->len == 0) return 0;
    uint32_t top = a->digits[a->len - 1];
    int bits = 32 * (a->len - 1);
    while (top) { bits++; top >>= 1; }
    return bits;
}
static void bu_shr1(BigUint *a) {
    uint64_t carry = 0;
    for (int i = a->len - 1; i >= 0; i--) {
        uint64_t value = (carry << 32) | a->digits[i];
        a->digits[i] = (uint32_t)(value >> 1);
        carry = value & 1;
    }
    bu_normalize(a);
}
static void bu_add(BigUint *res, const BigUint *a, const BigUint *b) {
    uint64_t carry = 0;
    int max_len = a->len > b->len ? a->len : b->len;
    for (int i = 0; i < max_len || carry; i++) {
        uint64_t av = i < a->len ? a->digits[i] : 0;
        uint64_t bv = i < b->len ? b->digits[i] : 0;
        uint64_t sum = av + bv + carry;
        if (i < BIGUINT_MAX_WORDS) res->digits[i] = (uint32_t)sum;
        carry = sum >> 32;
    }
    res->len = max_len + (carry ? 1 : 0);
    if (carry && res->len < BIGUINT_MAX_WORDS) res->digits[res->len - 1] = (uint32_t)carry;
    bu_normalize(res);
}
static void bu_mul_small(BigUint *res, const BigUint *a, uint32_t b) {
    uint64_t carry = 0;
    int i = 0;
    for (; i < a->len || carry; i++) {
        uint64_t av = i < a->len ? a->digits[i] : 0;
        uint64_t prod = av * b + carry;
        res->digits[i] = (uint32_t)prod;
        carry = prod >> 32;
    }
    res->len = i;
    bu_normalize(res);
}
static void bu_add_small(BigUint *res, const BigUint *a, uint32_t b) {
    uint64_t carry = b;
    int i = 0;
    for (; i < a->len && carry; i++) {
        uint64_t sum = (uint64_t)a->digits[i] + carry;
        res->digits[i] = (uint32_t)sum;
        carry = sum >> 32;
    }
    for (; i < a->len; i++) res->digits[i] = a->digits[i];
    if (carry) res->digits[i++] = (uint32_t)carry;
    res->len = i;
    bu_normalize(res);
}
static void bu_divmod_small(BigUint *quot, const BigUint *num, uint32_t den, uint32_t *rem);

static void bu_from_decimal(BigUint *a, const char *str) {
    bu_init(a);
    for (const char *p = str; *p; p++) {
        if (*p < '0' || *p > '9') continue;
        uint32_t digit = *p - '0';
        BigUint temp; bu_mul_small(&temp, a, 10);
        bu_add_small(&temp, &temp, digit);
        bu_copy(a, &temp);
    }
}
static void bu_to_decimal(const BigUint *a, char *buffer, size_t buffer_len) {
    if (bu_is_zero(a)) { snprintf(buffer, buffer_len, "0"); return; }
    char temp[2048] = {0};  /* increased for full 4096-bit safety */
    int pos = 0;
    BigUint copy; bu_copy(&copy, a);
    while (!bu_is_zero(&copy)) {
        uint32_t r;
        BigUint next;
        bu_divmod_small(&next, &copy, 10, &r);
        temp[pos++] = '0' + (char)r;
        bu_copy(&copy, &next);
    }
    int idx = 0;
    for (int i = pos - 1; i >= 0 && idx < (int)buffer_len - 1; i--) buffer[idx++] = temp[i];
    buffer[idx] = '\0';
}

static const char *BASE4096_ALPHABET_FILE_PATH = "base4096-2.0.1/frozen_base4096_alphabet.txt";
static char *base4096_alphabet_text = NULL;
static const char *base4096_alphabet_symbols[4096];
static bool base4096_alphabet_loaded = false;
static bool base4096_alphabet_load_failed = false;

static size_t utf8_char_len(uint8_t first) {
    if (first < 0x80) return 1;
    if ((first & 0xE0) == 0xC0) return 2;
    if ((first & 0xF0) == 0xE0) return 3;
    if ((first & 0xF8) == 0xF0) return 4;
    return 1;
}

static bool ensure_base4096_alphabet(void) {
    if (base4096_alphabet_loaded) return true;
    if (base4096_alphabet_load_failed) return false;

    FILE *f = fopen(BASE4096_ALPHABET_FILE_PATH, "rb");
    if (!f) {
        base4096_alphabet_load_failed = true;
        return false;
    }

    char temp_buf[32768];
    size_t read_bytes = fread(temp_buf, 1, sizeof(temp_buf) - 1, f);
    fclose(f);
    if (read_bytes == 0 || read_bytes >= sizeof(temp_buf) - 1) {
        base4096_alphabet_load_failed = true;
        return false;
    }
    temp_buf[read_bytes] = '\0';

    char clean_buf[32768];
    char *dst = clean_buf;
    for (char *src = temp_buf; *src; ++src) {
        if (*src == '\r' || *src == '\n') continue;
        *dst++ = *src;
    }
    *dst = '\0';

    size_t total_bytes = (size_t)(dst - clean_buf);
    char *alphabet_copy = (char *)malloc(total_bytes + 1);
    if (!alphabet_copy) {
        base4096_alphabet_load_failed = true;
        return false;
    }
    memcpy(alphabet_copy, clean_buf, total_bytes + 1);

    int count = 0;
    char *cursor = alphabet_copy;
    char *end = alphabet_copy + total_bytes;
    while (cursor < end && count < 4096) {
        uint8_t first = (uint8_t)*cursor;
        size_t len = utf8_char_len(first);
        if (cursor + len > end) {
            free(alphabet_copy);
            base4096_alphabet_load_failed = true;
            return false;
        }
        base4096_alphabet_symbols[count++] = cursor;
        cursor += len;
    }
    if (count != 4096 || cursor != end) {
        free(alphabet_copy);
        base4096_alphabet_load_failed = true;
        return false;
    }

    base4096_alphabet_text = alphabet_copy;
    base4096_alphabet_loaded = true;
    return true;
}

#define BASE4096_SYMBOL_START 0x1000u
#define BASE4096_RADIX 4096u

static size_t utf8_encode_codepoint(uint32_t codepoint, char *out, size_t out_size) {
    if (out_size == 0) return 0;
    if (codepoint <= 0x7F) {
        out[0] = (char)codepoint;
        return 1;
    }
    if (codepoint <= 0x7FF) {
        if (out_size < 2) return 0;
        out[0] = (char)(0xC0 | ((codepoint >> 6) & 0x1F));
        out[1] = (char)(0x80 | (codepoint & 0x3F));
        return 2;
    }
    if (codepoint <= 0xFFFF) {
        if (out_size < 3) return 0;
        out[0] = (char)(0xE0 | ((codepoint >> 12) & 0x0F));
        out[1] = (char)(0x80 | ((codepoint >> 6) & 0x3F));
        out[2] = (char)(0x80 | (codepoint & 0x3F));
        return 3;
    }
    if (codepoint <= 0x10FFFF) {
        if (out_size < 4) return 0;
        out[0] = (char)(0xF0 | ((codepoint >> 18) & 0x07));
        out[1] = (char)(0x80 | ((codepoint >> 12) & 0x3F));
        out[2] = (char)(0x80 | ((codepoint >> 6) & 0x3F));
        out[3] = (char)(0x80 | (codepoint & 0x3F));
        return 4;
    }
    return 0;
}

static bool bu_to_base4096(const BigUint *a, char *buffer, size_t buf_size) {
    bool canonical = ensure_base4096_alphabet();
    if (buf_size == 0) return false;
    if (bu_is_zero(a)) {
        if (canonical) {
            const char *sym = base4096_alphabet_symbols[0];
            size_t len = utf8_char_len((uint8_t)*sym);
            if (len >= buf_size) return false;
            memcpy(buffer, sym, len);
            buffer[len] = '\0';
            return true;
        }
        size_t len = utf8_encode_codepoint(BASE4096_SYMBOL_START, buffer, buf_size);
        if (len == 0 || len >= buf_size) return false;
        buffer[len] = '\0';
        return true;
    }

    BigUint temp;
    bu_copy(&temp, a);
    uint16_t digits[512] = {0};
    int count = 0;

    while (!bu_is_zero(&temp)) {
        uint32_t remainder = 0;
        BigUint next;
        bu_divmod_small(&next, &temp, BASE4096_RADIX, &remainder);
        if (remainder >= BASE4096_RADIX) return false;
        if (count >= (int)(sizeof(digits)/sizeof(digits[0]))) return false;
        digits[count++] = (uint16_t)remainder;
        bu_copy(&temp, &next);
    }

    char *ptr = buffer;
    size_t remaining = buf_size;
    for (int i = count - 1; i >= 0; i--) {
        if (canonical) {
            const char *sym = base4096_alphabet_symbols[digits[i]];
            size_t len = utf8_char_len((uint8_t)*sym);
            if (len >= remaining) return false;
            memcpy(ptr, sym, len);
            ptr += len;
            remaining -= len;
        } else {
            uint32_t codepoint = BASE4096_SYMBOL_START + (uint32_t)digits[i];
            char utf8[4];
            size_t len = utf8_encode_codepoint(codepoint, utf8, sizeof(utf8));
            if (len == 0 || len >= remaining) return false;
            memcpy(ptr, utf8, len);
            ptr += len;
            remaining -= len;
        }
    }
    if (remaining == 0) return false;
    *ptr = '\0';
    return true;
}

static bool decimal_to_base4096(const char *decimal, char *out, size_t out_size) {
    BigUint value;
    bu_init(&value);
    bu_from_decimal(&value, decimal);
    return bu_to_base4096(&value, out, out_size);
}

static bool print_base4096_if_fit(const char *decimal, const char *label) {
    if (!decimal || !*decimal) return false;
    size_t digits = 0;
    for (const char *p = decimal; *p; ++p) {
        if (*p >= '0' && *p <= '9') digits++;
    }
    if (digits == 0 || digits > 1230ULL) return false;
    char buffer[4096] = {0};
    if (!decimal_to_base4096(decimal, buffer, sizeof(buffer))) return false;
    if (label && *label) {
        printf("%s: %s\n", label, buffer);
    } else {
        printf("Base4096: %s\n", buffer);
    }
    return true;
}

static void bu_divmod_small(BigUint *quot, const BigUint *num, uint32_t den, uint32_t *rem) {
    BigUint result = {0};
    uint64_t carry = 0;
    for (int i = num->len - 1; i >= 0; i--) {
        uint64_t value = (carry << 32) | num->digits[i];
        uint32_t qv = (uint32_t)(value / den);
        carry = value % den;
        result.digits[i] = qv;
    }
    result.len = num->len;
    bu_normalize(&result);
    if (quot) bu_copy(quot, &result);
    if (rem) *rem = (uint32_t)carry;
}

static bool bu_is_one(const BigUint *a) {
    return a->len == 1 && a->digits[0] == 1;
}

static bool bu_is_equal(const BigUint *a, const BigUint *b) {
    return bu_cmp(a, b) == 0;
}

static void bu_sub(BigUint *res, const BigUint *a, const BigUint *b) {
    uint64_t borrow = 0;
    for (int i = 0; i < a->len; i++) {
        uint64_t av = a->digits[i];
        uint64_t bv = i < b->len ? b->digits[i] : 0;
        uint64_t diff = av - bv - borrow;
        res->digits[i] = (uint32_t)diff;
        borrow = (diff >> 63) & 1;
    }
    res->len = a->len;
    bu_normalize(res);
}

static void bu_sub_one(BigUint *res, const BigUint *a) {
    BigUint one;
    bu_from_u64(&one, 1);
    bu_sub(res, a, &one);
}

static void bu_mod_add(BigUint *res, const BigUint *a, const BigUint *b, const BigUint *mod) {
    BigUint sum;
    bu_add(&sum, a, b);
    if (bu_cmp(&sum, mod) >= 0) {
        bu_sub(&sum, &sum, mod);
    }
    bu_copy(res, &sum);
}

static void bu_modmul(const BigUint *a, const BigUint *b, const BigUint *mod, BigUint *res) {
    BigUint result;
    bu_init(&result);
    BigUint temp;
    bu_copy(&temp, a);

    for (int i = 0; i < b->len; i++) {
        uint32_t word = b->digits[i];
        for (int bit = 0; bit < 32; bit++) {
            if (word & 1u) {
                bu_mod_add(&result, &result, &temp, mod);
            }
            word >>= 1;
            bu_mod_add(&temp, &temp, &temp, mod);
        }
    }
    bu_copy(res, &result);
}

static void bu_modexp(const BigUint *base, const BigUint *exp, const BigUint *mod, BigUint *res) {
    BigUint result;
    bu_init(&result);
    bu_from_u64(&result, 1);
    BigUint power;
    bu_copy(&power, base);
    BigUint exponent;
    bu_copy(&exponent, exp);

    while (!bu_is_zero(&exponent)) {
        if (exponent.digits[0] & 1u) {
            bu_modmul(&result, &power, mod, &result);
        }
        bu_modmul(&power, &power, mod, &power);
        bu_shr1(&exponent);
    }

    bu_copy(res, &result);
}

static void bu_from_log2_approx(BigUint *out, long double log2x) {
    if (log2x <= 0.0L) {
        bu_init(out);
        return;
    }
    long double integer_part = floorl(log2x);
    long double frac = log2x - integer_part;
    int bitlen = (int)integer_part + 1;
    int top_index = (bitlen - 1) / 32;
    int top_shift = (bitlen - 1) % 32;
    long double mant = powl(2.0L, frac);
    unsigned long long top_value = (unsigned long long)(mant * (1ULL << top_shift));
    long double lower_frac = mant * (1ULL << top_shift) - (long double)top_value;
    uint32_t lower_word = (uint32_t)(lower_frac * 4294967296.0L);

    bu_init(out);
    if (top_index >= BIGUINT_MAX_WORDS) {
        top_index = BIGUINT_MAX_WORDS - 1;
    }
    out->len = top_index + 1;
    for (int i = 0; i < out->len; i++) out->digits[i] = 0;
    out->digits[top_index] = (uint32_t)(top_value & 0xFFFFFFFFu);
    if (top_index > 0) {
        out->digits[top_index - 1] = lower_word;
    }
    bu_normalize(out);
}

static void bu_add_u64(BigUint *res, const BigUint *a, uint64_t b) {
    if (b == 0) { bu_copy(res, a); return; }
    BigUint bb; bu_from_u64(&bb, b);
    bu_add(res, a, &bb);
}

static bool bu_sub_u64(BigUint *res, const BigUint *a, uint64_t b) {
    if (b == 0) { bu_copy(res, a); return true; }
    BigUint bb; bu_from_u64(&bb, b);
    if (bu_cmp(a, &bb) < 0) {
        bu_init(res);
        return false;
    }
    bu_sub(res, a, &bb);
    return true;
}

static void bu_from_apa(BigUint *out, const APAFloat *apa) {
    if (apa->mantissa[0] == 0ULL) {
        bu_init(out);
        return;
    }

    BigUint temp;
    bu_init(&temp);
    int word_index = 0;
    for (int i = APA_MANTISSA_WORDS - 1; i >= 0; i--) {
        if (word_index < BIGUINT_MAX_WORDS) temp.digits[word_index++] = (uint32_t)(apa->mantissa[i] & 0xFFFFFFFFu);
        if (word_index < BIGUINT_MAX_WORDS) temp.digits[word_index++] = (uint32_t)(apa->mantissa[i] >> 32);
    }
    temp.len = word_index;
    bu_normalize(&temp);

    int64_t shift = apa->exponent - (APA_MANTISSA_WORDS * 64);
    if (shift > 0) {
        for (int64_t i = 0; i < shift; i++) bu_shl1(&temp);
    } else if (shift < 0) {
        for (int64_t i = 0; i < -shift; i++) bu_shr1(&temp);
    }

    bu_copy(out, &temp);
}

static BigUint bu_from_omega(long double omega) {
    long double exponent = omega * LN_PHI;
    long double y = expl(exponent);
    long double log2x = (LN_PHI * y) / LN2;
    long double integer_part = floorl(log2x);
    long double frac = log2x - integer_part;
    APAFloat apa;
    apa_init(&apa, powl(2.0L, frac));
    apa.exponent = (int64_t)integer_part + 1;
    apa_normalize(&apa);

    BigUint result;
    bu_from_apa(&result, &apa);
    if (bu_is_zero(&result)) bu_from_u64(&result, 101ULL);
    if (bu_is_even(&result)) bu_add_small(&result, &result, 1);
    return result;
}

static long double bu_high_precision_log2(const BigUint *a) {
    APAFloat apa = bu_to_apa(a);
    long double result = apa_log2(&apa);
    debug_printf("DEBUG bu_high_precision_log2 bits=%d exponent=%lld log2=%.18Lf\n", bu_bit_length(a), apa.exponent, result);
    return result;
}

/* ====================== 64-BIT MODPOW (for deterministic Miller-Rabin) ====================== */
static uint64_t mulmod64(uint64_t a, uint64_t b, uint64_t mod) {
    uint64_t result = 0;
    a %= mod;
    while (b) {
        if (b & 1ULL) {
            if (result >= mod - a) result = result + a - mod;
            else result += a;
        }
        b >>= 1;
        if (b) {
            if (a >= mod - a) a = a + a - mod;
            else a += a;
        }
    }
    return result;
}

static uint64_t modpow64(uint64_t base, uint64_t exponent, uint64_t modulus) {
    if (modulus == 1) return 0;
    uint64_t result = 1 % modulus;
    uint64_t power = base % modulus;
    while (exponent > 0) {
        if (exponent & 1ULL) {
            result = mulmod64(result, power, modulus);
        }
        power = mulmod64(power, power, modulus);
        exponent >>= 1;
    }
    return result;
}

/* ====================== MILLER-RABIN (deterministic up to 64 bits) ====================== */
static uint64_t bu_to_u64(const BigUint *a) {
    if (a->len == 0) return 0;
    if (a->len == 1) return a->digits[0];
    return a->digits[0] | ((uint64_t)a->digits[1] << 32);
}

static bool miller_rabin_big(const BigUint *n) {
    int bits = bu_bit_length(n);
    if (bits == 0 || bits == 1) return false;
    if (bits <= 64) {
        uint64_t v = bu_to_u64(n);
        if (v < 2) return false;
        if (v == 2 || v == 3 || v == 5) return true;
        if ((v & 1) == 0 || v % 3 == 0 || v % 5 == 0) return false;

        uint64_t nm1 = v - 1;
        uint64_t d = nm1;
        int s = 0;
        while ((d & 1) == 0) { d >>= 1; s++; }

        static const uint64_t bases[12] = {2,3,5,7,11,13,23,29,31,37,41,43};
        for (size_t i = 0; i < 12; i++) {
            uint64_t a = bases[i];
            if (a >= v) break;
            uint64_t x = modpow64(a, d, v);
            if (x == 1 || x == nm1) continue;
            bool is_composite = true;
            for (int r = 1; r < s; r++) {
                x = mulmod64(x, x, v);
                if (x == nm1) { is_composite = false; break; }
            }
            if (is_composite) return false;
        }
        return true;
    } else {
        static const uint32_t small_primes[] = {
            2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,
            73,79,83,89,97
        };
        for (size_t i = 0; i < sizeof(small_primes)/sizeof(small_primes[0]); i++) {
            uint32_t p = small_primes[i];
            uint32_t r;
            bu_divmod_small(NULL, n, p, &r);
            if (r == 0) return false;
        }

        BigUint nm1;
        bu_copy(&nm1, n);
        bu_sub_one(&nm1, n);

        BigUint d;
        bu_copy(&d, &nm1);
        int s = 0;
        while (!bu_is_zero(&d) && bu_is_even(&d)) {
            bu_shr1(&d);
            s++;
        }

        static const uint64_t bases[7] = {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL};
        for (size_t i = 0; i < sizeof(bases)/sizeof(bases[0]); i++) {
            uint64_t a_val = bases[i];
            BigUint a;
            bu_from_u64(&a, a_val);
            BigUint x;
            bu_modexp(&a, &d, n, &x);
            if (bu_is_one(&x) || bu_is_equal(&x, &nm1)) continue;

            bool composite = true;
            for (int r = 1; r < s; r++) {
                bu_modmul(&x, &x, n, &x);
                if (bu_is_equal(&x, &nm1)) { composite = false; break; }
            }
            if (composite) return false;
        }

        printf("  [NOTE] Large-number primality is now strong probable prime after %zu Miller-Rabin bases.\n", sizeof(bases)/sizeof(bases[0]));
        return true;
    }
}

/* ====================== PLN STRUCTURES ====================== */
typedef struct PLNAddress {
    int D;
    long double theta;
    long double Q;
    long double omega;
} PLNAddress;

typedef struct {
    BigUint value;
    PLNAddress node;
} PLNNode;

typedef struct {
    long double log2_approx;
    long double frac_correction;
} LogCandidate;

static void logcand_init(LogCandidate *lc, long double log2x);
static BigUint bu_from_log_candidate(const LogCandidate *lc);
static void refine_log_candidate(LogCandidate *lc, const PLNAddress *target, int passes);
static void generate_xz_resonant_prime(long double n, long double n_scale, long double theta_gauge,
                                       long double z_real, long double z_imag,
                                       int max_steps, char *found_str, size_t buf_size);

static double hdgl_dn_score(const PLNAddress *addr) {
    if (!hdgl_enabled || !hdgl_compute_Dn_r) return 1.0;
    double omega = (double)addr->omega;
    double r = fmod(omega, 1.0);
    if (r < 0.0) r += 1.0;
    int dim = addr->D;
    if (dim < 1) dim = 1;
    if (dim > 8) dim = 8;
    return hdgl_compute_Dn_r(dim, r, omega);
}

static long double orbit_r(long double n) {
    long double frac = n - floorl(n);
    long double theta = frac * TWO_PI_L;
    long double D_cont = fmodl(n, 8.0L) + 1.0L;
    long double lat = PHI*PHI*(1.0L - powl(PHI, -D_cont)) + 1.0L/D_cont;
    lat = fmaxl(1.0L, fminl(lat, PHI*PHI));
    long double Rre = cosl(theta) - 1.0L + lat;
    long double Rim = sinl(theta);
    return sqrtl(Rre*Rre + Rim*Rim);
}

static long double omega_from_log2(long double log2x) {
    long double lnx = log2x * LN2;
    long double arg = lnx / LN_PHI;
    if (arg <= 1.0L) return 0.0L;
    return logl(arg) / LN_PHI;
}

static long double log10_from_omega(long double omega) {
    long double phi_power = powl(PHI, omega);
    return phi_power * LOG10_PHI;
}

static unsigned long long leading_digits_from_log10(long double log10x, int count) {
    long double integer_part = floorl(log10x);
    long double frac = log10x - integer_part;
    long double lead_value = powl(10.0L, frac + (long double)(count - 1));
    if (lead_value < 0.0L) lead_value = 0.0L;
    return (unsigned long long)(lead_value + 0.5L);
}

static bool prime_file_to_pln_address(const char *path, PLNAddress *out_addr, size_t *out_digits, long double *out_log10x) {
    long double log2x;
    size_t digits;
    if (!prime_file_to_log2(path, &log2x, &digits)) return false;

    long double omega = omega_from_log2(log2x);
    long double log10x = log2x / LOG2_10;
    int D = (int)fmodl(floorl(omega), 8.0L) + 1;
    if (D < 1) D = 1;
    long double theta = fmodl(omega, 1.0L) * M_PI_L;
    long double r = orbit_r(omega);
    long double Q = 1.0L / (r + 1e-18L);

    if (out_addr) *out_addr = (PLNAddress){D, theta, Q, omega};
    if (out_digits) *out_digits = digits;
    if (out_log10x) *out_log10x = log10x;
    return true;
}

static void generate_prime_file_candidates(const char *path, int count, long double omega_step, int lead_digits) {
    if (lead_digits < 1) lead_digits = 6;
    if (lead_digits > 18) lead_digits = 18;

    PLNAddress base;
    size_t digits;
    long double log10x;
    if (!prime_file_to_pln_address(path, &base, &digits, &log10x)) {
        printf("Error: unable to read prime file '%s' or file contains no digits.\n", path);
        return;
    }

    long double center_omega = base.omega;
    int half_count = count / 2;
    printf("=== HUGE PRIME DISCOVERY ENGINE ===\n");
    printf("Anchor file: %s\n", path);
    printf("Anchor digits: %zu\n", digits);
    printf("Anchor PLN: [D_%d, %.18Lf*pi, %.15Lf, %.18Lf*phi]\n", base.D, base.theta / M_PI_L, base.Q, base.omega);
    printf("Candidate sweep: %d nodes, omega step %.12Lf, lead digits %d\n\n", count, omega_step, lead_digits);

    for (int idx = 0; idx < count; idx++) {
        int offset = idx - half_count;
        long double cand_omega = center_omega + offset * omega_step;
        if (cand_omega <= 0.0L) cand_omega = omega_step;
        int cand_D = (int)fmodl(floorl(cand_omega), 8.0L) + 1;
        if (cand_D < 1) cand_D = 1;
        long double cand_theta = fmodl(cand_omega, 1.0L) * M_PI_L;
        long double cand_r = orbit_r(cand_omega);
        long double cand_Q = 1.0L / (cand_r + 1e-18L);
        long double cand_log10x = log10_from_omega(cand_omega);
        unsigned long long cand_digits = (unsigned long long)floorl(cand_log10x) + 1ULL;
        unsigned long long cand_lead = leading_digits_from_log10(cand_log10x, lead_digits);

        char lead_fmt[32];
        snprintf(lead_fmt, sizeof(lead_fmt), "%%s%%0%dllu", lead_digits);
        char lead_buf[32];
        snprintf(lead_buf, sizeof(lead_buf), lead_fmt, "", cand_lead);

        PLNAddress candidate_addr = {cand_D, cand_theta, cand_Q, cand_omega};
        long double dn_score = hdgl_enabled ? (long double)hdgl_dn_score(&candidate_addr) : 0.0L;
        printf("Candidate %2d | Δω = %+0.12Lf | D_%d | θ=%.12Lf*pi | Q=%.15Lf | digits≈%llu | lead≈%s",
               idx + 1, cand_omega - center_omega, cand_D, cand_theta / M_PI_L, cand_Q, cand_digits, lead_buf);
        if (hdgl_enabled) {
            printf(" | Dn=%.5Lf", dn_score);
        }
        printf("\n");
    }
    printf("\nNote: These are PLN discovery targets for astronomical primes.\n");
    printf("The engine provides adjacent resonance nodes for P136279841+ class candidates, not exact full decimal reconstruction.\n");
}

static PLNAddress pln_from_value_high_precision(const BigUint *x) {
    if (bu_is_zero(x)) return (PLNAddress){1, 0.0L, 1.0L, 0.0L};
    long double omega = 0.0L;
    long double prev_omega = 0.0L;
    for (int pass = 0; pass < MAX_PASSES; pass++) {
        long double log2x = bu_high_precision_log2(x);
        debug_printf("FORWARD log2x=%.18Lf\n", log2x);
        long double lnx = log2x * LN2;
        long double arg = lnx / LN_PHI;
        debug_printf("FORWARD lnx=%.18Lf arg=%.18Lf\n", lnx, arg);
        if (arg <= 1.0L) return (PLNAddress){1, 0.0L, 1.0L, 0.0L};
        omega = logl(arg) / LN_PHI;
        debug_printf("FORWARD omega=%.18Lf\n", omega);
        if (pass > 0) {
            long double correction = (omega - prev_omega) * 0.5L;
            omega -= correction;
        }
        prev_omega = omega;
    }
    int D = (int)fmodl(floorl(omega), 8.0L) + 1;
    if (D < 1) D = 1;
    long double theta = fmodl(omega, 1.0L) * M_PI_L;
    long double r = orbit_r(omega);
    long double Q = 1.0L / (r + 1e-18L);
    return (PLNAddress){D, theta, Q, omega};
}

static void print_pln(const PLNAddress *addr) {
    printf("PLN = [D_%d, %.18Lf*pi, %.15Lf, %.18Lf*phi]\n",
           addr->D, addr->theta / M_PI_L, addr->Q, addr->omega);
}

static PLNNode pln_node_new_from_decimal(const char *value_str) {
    PLNNode node;
    bu_init(&node.value);
    bu_from_decimal(&node.value, value_str);
    node.node = pln_from_value_high_precision(&node.value);
    return node;
}

static void pln_node_print(const PLNNode *node, const char *label) {
    char value_str[2048];  /* safety for full 4096-bit numbers */
    bu_to_decimal(&node->value, value_str, sizeof(value_str));
    printf("%s %s:\n", label, value_str);
    print_base4096_if_fit(value_str, "Base4096");
    print_pln(&node->node);
}

static void logcand_init(LogCandidate *lc, long double log2x) {
    lc->log2_approx = log2x;
    lc->frac_correction = 0.0L;
}

static BigUint bu_from_log_candidate(const LogCandidate *lc) {
    long double adjusted = lc->log2_approx + lc->frac_correction;
    BigUint out;
    bu_from_log2_approx(&out, adjusted);
    return out;
}

static void refine_log_candidate(LogCandidate *lc, const PLNAddress *target, int passes) {
    for (int pass = 0; pass < passes; pass++) {
        BigUint approx = bu_from_log_candidate(lc);
        PLNAddress est = pln_from_value_high_precision(&approx);
        long double delta_omega = target->omega - est.omega;
        long double correction = delta_omega * 0.5L;
        lc->frac_correction += correction * (LN_PHI / LN2);
        if (fabsl(delta_omega) < 1e-14L) break;
    }
}

static void generate_xz_resonant_prime(long double n, long double n_scale, long double theta_gauge,
                                       long double z_real, long double z_imag,
                                       int max_steps, char *found_str, size_t buf_size) {
    printf("=== X(z) RESONANT PRIME GENERATOR ===\n");
    printf("Parameters: n=%.6Lf, n_scale=%.6Lf, gauge_theta=%.6Lf, z=(%.4Lf + %.4Lfi)\n",
           n, n_scale, theta_gauge, z_real, z_imag);

    long double phi_base = powl(PHI, PHI - 1.0L);
    long double inner = n + 0.5L / PHI;
    long double lattice_term = powl(PHI, powl(PHI, inner));
    long double one_plus_z_mag = sqrtl((1.0L + z_real) * (1.0L + z_real) + z_imag * z_imag);
    long double power_term = powl(one_plus_z_mag, n_scale);

    long double x_approx = phi_base * sqrtl(lattice_term) * power_term;
    long double log2x = log2l(x_approx) + (theta_gauge / (2.0L * M_PI_L));

    LogCandidate lc;
    logcand_init(&lc, log2x);
    PLNAddress dummy_target = {1, 0.0L, 1.0L, 0.0L};
    refine_log_candidate(&lc, &dummy_target, 3);

    BigUint start = bu_from_log_candidate(&lc);
    PLNAddress pln_start = pln_from_value_high_precision(&start);

    printf("Initial PLN: [D_%d, %.12Lf*pi, %.12Lf, %.12Lf*phi]\n",
           pln_start.D, pln_start.theta / M_PI_L, pln_start.Q, pln_start.omega);
    printf("Approx decimal digits: %.0Lf\n", log10l(x_approx));

    bool found = find_prime_from_pln(pln_start.D,
                                     pln_start.theta / M_PI_L,
                                     pln_start.Q,
                                     pln_start.omega,
                                     max_steps, found_str, buf_size);

    if (found) {
        printf("→ Resonant prime recovered from X(z) lattice point.\n");
        print_base4096_if_fit(found_str, "Base4096");
    } else {
        printf("→ No strong resonance found. Try larger MAX_STEPS or different X(z) parameters.\n");
    }
}

/* ====================== FORWARD: INTEGER → PLN ====================== */
static void solve_integer_to_pln(const char *value_str) {
    printf("=== INTEGER → PLN (deterministic) ===\nValue: %s\n", value_str);
    print_base4096_if_fit(value_str, "Base4096");
    PLNNode node = pln_node_new_from_decimal(value_str);
    bool is_prime = miller_rabin_big(&node.value);
    printf("Primality: %s\n", is_prime ? "PRIME ✓" : "composite (PLN still valid)");
    pln_node_print(&node, "PLN ADDRESS");
    printf("→ Deterministic lattice encoding complete.\n\n");
}

/* ====================== REVERSE: PLN → PRIME ====================== */
static long double pln_resonance_score(const PLNAddress *a, const PLNAddress *b) {
    long double d_theta = fabsl(a->theta - b->theta);
    if (d_theta > M_PI_L) d_theta = TWO_PI_L - d_theta;
    return 0.5L/(1.0L + d_theta/M_PI_L) + 0.3L/(1.0L + fabsl(a->Q - b->Q)) + 0.2L/(1.0L + fabsl(fmodl(a->omega,1.0L) - fmodl(b->omega,1.0L)));
}

static bool find_prime_from_pln(int target_D, long double target_theta_pi, long double target_Q, long double target_omega, int max_steps, char *found_str, size_t buf_size) {
    printf("=== PLN → PRIME (deterministic multi-pass) ===\nTarget: [D_%d, %.12Lf*pi, %.12Lf, %.12Lf*phi]\n", target_D, target_theta_pi, target_Q, target_omega);
    PLNAddress target = {target_D, target_theta_pi * M_PI_L, target_Q, target_omega};
    long double current_omega = target_omega;
    BigUint center = bu_from_omega(current_omega);

    for (int refine = 0; refine < 3; refine++) {
        PLNAddress estimate = pln_from_value_high_precision(&center);
        long double delta = target.omega - estimate.omega;
        if (fabsl(delta) < 1e-12L) break;
        current_omega += delta * 0.5L;
        BigUint refined = bu_from_omega(current_omega);
        bu_copy(&center, &refined);
    }

    BigUint best;
    bu_init(&best);
    long double best_score = -1.0L;
    BigUint candidate;
    long double current_r = 0.0L;

    for (int offset = 0; offset <= max_steps; offset++) {
        if (offset == 0) {
            bu_copy(&candidate, &center);
        } else {
            uint64_t delta = (uint64_t)offset * 2ULL;
            if (offset % 2 == 1) {
                if (!bu_sub_u64(&candidate, &center, delta)) continue;
            } else {
                bu_add_u64(&candidate, &center, delta);
            }
        }

        PLNAddress addr = pln_from_value_high_precision(&candidate);
        if (addr.D != target_D) continue;

        long double r = orbit_r(addr.omega);
        long double score = pln_resonance_score(&addr, &target);
        if (r < R_MIN_THRESHOLD + 0.1L && score > best_score) {
            best_score = score;
            bu_copy(&best, &candidate);
            current_r = r;
            if (score > 0.999999L) break;
        }

        if (score > 0.9999L && miller_rabin_big(&candidate)) {
            bu_to_decimal(&candidate, found_str, buf_size);
            printf(" Found prime candidate %s (score %.6Lf, r=%.15Lf)\n", found_str, score, r);
            print_base4096_if_fit(found_str, "Base4096");
            return true;
        }
    }

    if (best_score > 0.7L && miller_rabin_big(&best)) {
        bu_to_decimal(&best, found_str, buf_size);
        printf(" Found nearest resonant prime %s (score %.6Lf, r=%.15Lf)\n", found_str, best_score, current_r);
        print_base4096_if_fit(found_str, "Base4096");
        return true;
    }

    printf(" No prime found within %d odd steps around estimated x.\n", max_steps);
    return false;
}

static bool prime_file_to_log2(const char *path, long double *out_log2x, size_t *out_digits) {
    FILE *f = fopen(path, "r");
    if (!f) return false;

    uint64_t lead = 0;
    int lead_digits = 0;
    size_t total_digits = 0;
    int c;

    while ((c = fgetc(f)) != EOF) {
        if (c >= '0' && c <= '9') {
            total_digits++;
            if (lead_digits < MAX_FILE_LEAD_DIGITS) {
                lead = lead * 10 + (uint64_t)(c - '0');
                lead_digits++;
            }
        }
    }
    fclose(f);

    if (total_digits == 0 || lead_digits == 0) return false;

    long double log10_lead = log10l((long double)lead);
    long double log10x = (long double)(total_digits - lead_digits) + log10_lead;
    *out_log2x = log10x * LOG2_10;
    if (out_digits) *out_digits = total_digits;
    return true;
}

static void solve_prime_file_to_pln(const char *path) {
    PLNAddress addr;
    size_t digits;
    long double log10x;
    if (!prime_file_to_pln_address(path, &addr, &digits, &log10x)) {
        printf("Error: unable to read prime file '%s' or file contains no digits.\n", path);
        return;
    }

    unsigned long long lead = leading_digits_from_log10(log10x, 6);
    printf("=== FILE PRIME → PLN ===\nFile: %s\n", path);
    printf("Digits: %zu\n", digits);
    printf("Leading digits (approx): %06llu\n", lead);
    printf("PLN = [D_%d, %.18Lf*pi, %.15Lf, %.18Lf*phi]\n", addr.D, addr.theta / M_PI_L, addr.Q, addr.omega);
}

static long double omega_from_log10(long double log10x) {
    return omega_from_log2(log10x * LOG2_10);
}

static PLNAddress pln_address_from_log10(long double log10x) {
    long double omega = omega_from_log10(log10x);
    int D = (int)fmodl(floorl(omega), 8.0L) + 1;
    if (D < 1) D = 1;
    long double theta = fmodl(omega, 1.0L) * M_PI_L;
    long double r = orbit_r(omega);
    long double Q = 1.0L / (r + 1e-18L);
    return (PLNAddress){D, theta, Q, omega};
}

static const char *score_status(long double score) {
    if (score >= PLN_SCORE_PASS_THRESHOLD) return "PASS";
    if (score >= PLN_SCORE_BORDER_THRESHOLD) return "BORDER";
    return "FAIL";
}

static long double score_margin(long double score) {
    return score - PLN_SCORE_PASS_THRESHOLD;
}

static long double score_candidate_by_digits(const PLNAddress *target, unsigned long long digits, int mod8,
                                         long double *out_theta_error,
                                         long double *out_Q_error,
                                         long double *out_omega_error,
                                         long double *out_mod8_bias) {
    long double log10x = (long double)(digits - 1ULL) + LOG10_5;
    PLNAddress candidate = pln_address_from_log10(log10x);

    long double d_theta = fabsl(candidate.theta - target->theta);
    if (d_theta > M_PI_L) d_theta = TWO_PI_L - d_theta;
    long double dQ = fabsl(candidate.Q - target->Q);
    long double omega_frac = fmodl(candidate.omega, 1.0L);
    long double target_omega_frac = fmodl(target->omega, 1.0L);
    long double d_omega = fabsl(omega_frac - target_omega_frac);
    if (d_omega > 0.5L) d_omega = 1.0L - d_omega;

    long double score = 0.5L / (1.0L + d_theta / M_PI_L)
                      + 0.3L / (1.0L + dQ)
                      + 0.2L / (1.0L + d_omega);
    if (candidate.D != target->D) {
        score -= 0.1L;
    }

    if (hdgl_enabled) {
        long double dn = (long double)hdgl_dn_score(&candidate);
        long double dn_bonus = 0.1L * (dn / (1.0L + dn));
        score += dn_bonus;
    }

    long double bias = 0.0L;
    if (mod8 == 7) {
        bias = 0.03L;
    } else if (mod8 == 1) {
        bias = -0.02L;
    }
    score += bias;

    if (out_theta_error) *out_theta_error = d_theta / M_PI_L;
    if (out_Q_error) *out_Q_error = dQ;
    if (out_omega_error) *out_omega_error = d_omega;
    if (out_mod8_bias) *out_mod8_bias = bias;

    if (score < 0.0L) score = 0.0L;
    if (score > 1.0L) score = 1.0L;
    return score;
}

typedef struct {
    char label[128];
    unsigned long long digits;
    int mod8;
    bool has_exact_value;
    bool has_lead;
    char exact_value[2048];
    char lead_str[64];
    PLNAddress addr;
    long double score;
    long double hdgl_dn;
    long double d_theta;
    long double d_Q;
    long double d_omega;
    long double mod8_bias;
    bool search_attempted;
    bool search_found;
    unsigned long long search_offset;
    int search_direction;
} CandidateScore;

static int compare_candidate_scores(const void *a, const void *b) {
    const CandidateScore *ca = (const CandidateScore *)a;
    const CandidateScore *cb = (const CandidateScore *)b;
    if (ca->score < cb->score) return 1;
    if (ca->score > cb->score) return -1;
    return 0;
}

static bool is_decimal_string(const char *s) {
    if (!s || !*s) return false;
    for (const char *p = s; *p; p++) if (*p < '0' || *p > '9') return false;
    return true;
}

static unsigned long long parse_label_digits(const char *label) {
    if (!label || label[0] != 'P') return 0ULL;
    unsigned long long digits = 0ULL;
    for (const char *p = label + 1; *p; p++) {
        if (*p < '0' || *p > '9') break;
        digits = digits * 10ULL + (unsigned long long)(*p - '0');
    }
    return digits;
}

static bool candidate_within_exact_range(const CandidateScore *cand) {
    return cand->digits > 0ULL && cand->digits <= 1230ULL;
}

static void biguint_from_decimal_string(BigUint *out, const char *decimal) {
    bu_from_decimal(out, decimal);
}

static void biguint_from_lead_and_zeros(BigUint *out, const char *lead_str, unsigned long long digits);
static void biguint_build_center(BigUint *out, unsigned long long digits);

static bool exact_search_around_biguint(const BigUint *base, unsigned long long max_offset, unsigned long long *out_offset, int *out_direction, BigUint *out_prime) {
    if (!base || max_offset == 0ULL) return false;

    if (miller_rabin_big(base)) {
        if (out_offset) *out_offset = 0ULL;
        if (out_direction) *out_direction = 0;
        if (out_prime) bu_copy(out_prime, base);
        return true;
    }

    bool base_even = bu_is_even(base);
    BigUint candidate;
    BigUint neighbor;
    for (unsigned long long offset = 1ULL; offset <= max_offset; offset++) {
        bool test_offset = base_even ? ((offset & 1ULL) != 0ULL) : ((offset & 1ULL) == 0ULL);
        if (!test_offset) continue;

        bu_add_u64(&candidate, base, offset);
        if (miller_rabin_big(&candidate)) {
            if (out_offset) *out_offset = offset;
            if (out_direction) *out_direction = 1;
            if (out_prime) bu_copy(out_prime, &candidate);
            return true;
        }

        if (bu_sub_u64(&neighbor, base, offset)) {
            if (miller_rabin_big(&neighbor)) {
                if (out_offset) *out_offset = offset;
                if (out_direction) *out_direction = -1;
                if (out_prime) bu_copy(out_prime, &neighbor);
                return true;
            }
        }
    }
    return false;
}

static bool exact_search_candidate_with_result(const CandidateScore *cand, unsigned long long max_offset, unsigned long long *out_offset, int *out_direction, BigUint *out_prime) {
    if (!candidate_within_exact_range(cand)) return false;

    BigUint base;
    bu_init(&base);
    if (cand->has_exact_value) {
        biguint_from_decimal_string(&base, cand->exact_value);
    } else if (cand->has_lead) {
        biguint_from_lead_and_zeros(&base, cand->lead_str, cand->digits);
    } else {
        biguint_build_center(&base, cand->digits);
    }

    return exact_search_around_biguint(&base, max_offset, out_offset, out_direction, out_prime);
}

static void biguint_from_lead_and_zeros(BigUint *out, const char *lead_str, unsigned long long digits) {
    BigUint value;
    bu_from_decimal(&value, lead_str);
    unsigned long long lead_len = strlen(lead_str);
    unsigned long long zeros = 0ULL;
    if (digits > lead_len) zeros = digits - lead_len;
    for (unsigned long long i = 0; i < zeros; i++) {
        BigUint temp;
        bu_mul_small(&temp, &value, 10);
        bu_copy(&value, &temp);
    }
    bu_copy(out, &value);
}

static void biguint_build_center(BigUint *out, unsigned long long digits) {
    bu_init(out);
    bu_from_u64(out, 5ULL);
    for (unsigned long long i = 1ULL; i < digits; i++) {
        BigUint tmp;
        bu_mul_small(&tmp, out, 10);
        bu_copy(out, &tmp);
    }
}

static bool exact_search_candidate(const CandidateScore *cand) {
    if (!candidate_within_exact_range(cand)) return false;

    BigUint base;
    bu_init(&base);
    if (cand->has_exact_value) {
        biguint_from_decimal_string(&base, cand->exact_value);
    } else if (cand->has_lead) {
        biguint_from_lead_and_zeros(&base, cand->lead_str, cand->digits);
    } else {
        biguint_build_center(&base, cand->digits);
        printf("Approximate center built for candidate '%s' using digit count %llu.\n", cand->label, cand->digits);
    }

    char value_str[2048];
    bu_to_decimal(&base, value_str, sizeof(value_str));

    if (miller_rabin_big(&base)) {
        printf("Exact candidate integer verified prime: %s\n", value_str);
        print_base4096_if_fit(value_str, "Base4096");
        return true;
    }

    BigUint test;
    BigUint minus;
    for (int offset = 1; offset <= 256; offset += 2) {
        bu_add_u64(&test, &base, (unsigned long long)offset);
        if (miller_rabin_big(&test)) {
            bu_to_decimal(&test, value_str, sizeof(value_str));
            printf("Exact prime found at +%d around candidate '%s': %s\n", offset, cand->label, value_str);
            print_base4096_if_fit(value_str, "Base4096");
            return true;
        }

        if (bu_sub_u64(&minus, &base, (unsigned long long)offset)) {
            if (miller_rabin_big(&minus)) {
                bu_to_decimal(&minus, value_str, sizeof(value_str));
                printf("Exact prime found at -%d around candidate '%s': %s\n", offset, cand->label, value_str);
                print_base4096_if_fit(value_str, "Base4096");
                return true;
            }
        }
    }

    printf("Exact/approximate search attempted for candidate '%s' within 4096-bit range but no prime was found in ±256 steps.\n", cand->label);
    return false;
}

static void test_exact_candidate_range(const char *value_str, unsigned long long max_offset) {
    if (!value_str || !*value_str) {
        printf("Error: missing numeric candidate string.\n");
        return;
    }

    size_t digits = 0;
    for (const char *p = value_str; *p; ++p) {
        if (*p >= '0' && *p <= '9') digits++;
    }
    if (digits == 0) {
        printf("Error: candidate must be a decimal integer.\n");
        return;
    }
    if (digits > 1230ULL) {
        printf("Error: candidate has %zu digits, which exceeds the 4096-bit engine limit of 1230 digits.\n", digits);
        printf("This tool can only test exact integer candidates up to 1230 digits.\n");
        return;
    }

    BigUint base;
    bu_init(&base);
    bu_from_decimal(&base, value_str);

    printf("=== EXACT INTEGER CANDIDATE RANGE TEST ===\n");
    printf("Candidate (%zu digits): %s\n", digits, value_str);

    if (miller_rabin_big(&base)) {
        printf("Base candidate is PRIME ✓\n");
        print_base4096_if_fit(value_str, "Base4096");
        return;
    }

    printf("Base candidate is composite or strong probable composite. Searching ±%llu around it...\n", max_offset);
    bool base_even = bu_is_even(&base);
    BigUint candidate;
    BigUint neighbor;

    for (unsigned long long offset = 1ULL; offset <= max_offset; offset++) {
        bool test_offset = base_even ? ((offset & 1ULL) != 0ULL) : ((offset & 1ULL) == 0ULL);
        if (!test_offset) continue;

        bu_add_u64(&candidate, &base, offset);
        if (miller_rabin_big(&candidate)) {
            char found_str[2048];
            bu_to_decimal(&candidate, found_str, sizeof(found_str));
            printf("Prime found at +%llu: %s\n", offset, found_str);
            print_base4096_if_fit(found_str, "Base4096");
            return;
        }

        if (bu_sub_u64(&neighbor, &base, offset)) {
            if (miller_rabin_big(&neighbor)) {
                char found_str[2048];
                bu_to_decimal(&neighbor, found_str, sizeof(found_str));
                printf("Prime found at -%llu: %s\n", offset, found_str);
                print_base4096_if_fit(found_str, "Base4096");
                return;
            }
        }
    }
    printf("No prime found in the tested integer range ±%llu around candidate.\n", max_offset);
}

static long double score_exact_integer(const PLNAddress *addr, const PLNAddress *anchor, bool has_anchor) {
    long double tension = fabsl(orbit_r(addr->omega) - INV_PHI);
    long double tension_score = 1.0L / (1.0L + tension * 10.0L);
    long double anchor_score = has_anchor ? pln_resonance_score(addr, anchor) : 0.0L;
    long double score = has_anchor ? (0.7L * anchor_score + 0.3L * tension_score) : tension_score;
    if (addr->D == 7 || addr->D == 1) score += 0.05L;
    if (hdgl_enabled) {
        long double dn = (long double)hdgl_dn_score(addr);
        long double dn_bonus = 0.1L * (dn / (1.0L + dn));
        score += dn_bonus;
    }
    if (score > 1.0L) score = 1.0L;
    return score;
}

static void scan_range_around_candidate(const char *value_str, unsigned long long max_offset, const char *anchor_file) {
    if (!value_str || !*value_str) {
        printf("Error: missing numeric candidate string.\n");
        return;
    }

    size_t digits = 0;
    for (const char *p = value_str; *p; ++p) {
        if (*p >= '0' && *p <= '9') digits++;
    }
    if (digits == 0) {
        printf("Error: candidate must be a decimal integer.\n");
        return;
    }
    if (digits > 1230ULL) {
        printf("Error: candidate has %zu digits, which exceeds the 4096-bit engine limit of 1230 digits.\n", digits);
        printf("This tool can only scan exact integers up to 1230 digits.\n");
        return;
    }

    PLNAddress anchor;
    bool have_anchor = false;
    if (anchor_file && *anchor_file) {
        size_t anchor_digits;
        long double anchor_log10;
        if (!prime_file_to_pln_address(anchor_file, &anchor, &anchor_digits, &anchor_log10)) {
            printf("Error: unable to read anchor prime file '%s'.\n", anchor_file);
            return;
        }
        have_anchor = true;
    }

    BigUint base;
    bu_init(&base);
    bu_from_decimal(&base, value_str);

    PLNAddress base_addr = pln_from_value_high_precision(&base);
    PLNAddress score_anchor = base_addr;
    if (have_anchor) {
        score_anchor = anchor;
    }

    printf("=== SCAN RANGE AROUND EXACT CANDIDATE ===\n");
    printf("Base candidate (%zu digits): %s\n", digits, value_str);
    printf("Scan limit: ±%llu integers\n", max_offset);
    if (have_anchor) {
        printf("Anchor prime file: %s\n", anchor_file);
        printf("Anchor PLN: [D_%d, %.18Lf*pi, %.15Lf, %.18Lf*phi]\n", anchor.D, anchor.theta / M_PI_L, anchor.Q, anchor.omega);
    } else {
        printf("Base PLN: [D_%d, %.18Lf*pi, %.15Lf, %.18Lf*phi]\n", base_addr.D, base_addr.theta / M_PI_L, base_addr.Q, base_addr.omega);
        printf("Scoring neighbors relative to the base candidate's PLN.\n");
    }
    printf("\nOffset | Status | mod8 | D | score     | margin    | omega_diff\n");
    printf("-------+--------+------+---+-----------+-----------+-----------\n");

    BigUint candidate;
    PLNAddress addr;

    long double best_score = -1.0L;
    long long best_offset = 0LL;
    long double best_omega_diff = 1e18L;
    long long best_omega_offset = 0LL;
    bool best_omega_pass = false;
    bool any_prime = false;

    for (long long delta = -((long long)max_offset); delta <= (long long)max_offset; delta++) {
        bool valid = true;
        if (delta == 0) {
            bu_copy(&candidate, &base);
        } else if (delta > 0) {
            bu_add_u64(&candidate, &base, (unsigned long long)delta);
        } else {
            if (!bu_sub_u64(&candidate, &base, (unsigned long long)(-delta))) {
                valid = false;
            }
        }
        if (!valid) continue;

        addr = pln_from_value_high_precision(&candidate);
        int mod8val = candidate.len > 0 ? (int)(candidate.digits[0] & 7u) : 0;
        bool is_even = bu_is_even(&candidate);
        bool is_prime = false;
        if (!is_even) {
            is_prime = miller_rabin_big(&candidate);
        }
        const char *status = is_prime ? "PASS" : "FAIL";
        long double score = score_exact_integer(&addr, &score_anchor, true);
        long double margin = score - PLN_SCORE_PASS_THRESHOLD;
        long double omega_diff = fabsl(addr.omega - score_anchor.omega);

        char offset_label[32];
        if (delta >= 0) {
            snprintf(offset_label, sizeof(offset_label), "+%lld", delta);
        } else {
            snprintf(offset_label, sizeof(offset_label), "%lld", delta);
        }

        printf("%6s | %6s | %4d | %1d | %0.9Lf | %0.9Lf | %0.9Lf\n",
               offset_label,
               status,
               mod8val,
               addr.D,
               score,
               margin,
               omega_diff);

        if (score > best_score) {
            best_score = score;
            best_offset = delta;
        }
        if (omega_diff < best_omega_diff) {
            best_omega_diff = omega_diff;
            best_omega_offset = delta;
            best_omega_pass = is_prime;
        }
        if (is_prime) {
            any_prime = true;
        }
    }

    printf("\nBest nearby score: %0.6Lf at %s%lld\n",
           best_score,
           best_offset < 0 ? "" : "+",
           best_offset);
    printf("Best omega_diff: %0.9Lf at %s%lld (%s)\n",
           best_omega_diff,
           best_omega_offset < 0 ? "" : "+",
           best_omega_offset,
           best_omega_pass ? "PASS" : "FAIL");
    if (any_prime) {
        printf("At least one probable prime was found in the scanned range.\n");
    } else {
        printf("No probable prime was found in the scanned range.\n");
    }
}

static int load_candidate_scores(const PLNAddress *target, unsigned long long anchor_digits, const char *candidate_file, CandidateScore *candidates, int max_candidates) {
    FILE *f = fopen(candidate_file, "r");
    if (!f) return -1;

    int count = 0;
    char line[256];
    while (fgets(line, sizeof(line), f) && count < max_candidates) {
        char token1[128] = {0};
        char token2[128] = {0};
        char token3[128] = {0};
        char token4[128] = {0};
        unsigned long long cand_digits = 0ULL;
        unsigned long long cand_lead = 0ULL;
        int cand_lead_len = 0;
        int cand_mod8 = -1;
        bool has_exact = false;
        bool has_lead = false;
        int scanned = sscanf(line, "%127s %127s %127s %127s", token1, token2, token3, token4);
        char label[128] = {0};
        char lead_str[64] = {0};
        char exact_value[2048] = {0};

        if (scanned == 4) {
            if (token1[0] == 'P' && is_decimal_string(token2) && is_decimal_string(token3)) {
                strncpy(label, token1, sizeof(label) - 1);
                cand_digits = strtoull(token2, NULL, 10);
                cand_lead_len = (int)strlen(token3);
                if (cand_lead_len > 0 && cand_lead_len < (int)sizeof(lead_str)) {
                    strncpy(lead_str, token3, cand_lead_len);
                    has_lead = true;
                }
                cand_mod8 = atoi(token4);
            }
        } else if (scanned == 3) {
            if (token1[0] == 'P' && is_decimal_string(token2) && is_decimal_string(token3)) {
                strncpy(label, token1, sizeof(label) - 1);
                cand_digits = strtoull(token2, NULL, 10);
                cand_mod8 = atoi(token3);
                has_lead = false;
            } else if (is_decimal_string(token1) && is_decimal_string(token2) && is_decimal_string(token3)) {
                size_t len1 = strlen(token1);
                int mod8 = atoi(token3);
                if (len1 <= 1230 && mod8 >= 0 && mod8 < 8) {
                    cand_digits = (unsigned long long)len1;
                    cand_mod8 = mod8;
                    has_exact = true;
                    strncpy(exact_value, token1, sizeof(exact_value) - 1);
                    cand_lead_len = (int)fmin((size_t)sizeof(lead_str)-1, len1);
                    strncpy(lead_str, token1, cand_lead_len);
                    has_lead = true;
                    snprintf(label, sizeof(label), "exact-candidate");
                }
            }
        } else if (scanned == 2) {
            if (is_decimal_string(token1) && is_decimal_string(token2) && strlen(token1) <= 1230) {
                int mod8 = atoi(token2);
                if (mod8 >= 0 && mod8 < 8) {
                    has_exact = true;
                    cand_digits = (unsigned long long)strlen(token1);
                    cand_mod8 = mod8;
                    strncpy(exact_value, token1, sizeof(exact_value) - 1);
                    cand_lead_len = (int)fmin((size_t)sizeof(lead_str)-1, strlen(token1));
                    strncpy(lead_str, token1, cand_lead_len);
                    has_lead = true;
                    snprintf(label, sizeof(label), "exact-candidate");
                }
            }
            if (!has_exact) {
                cand_digits = strtoull(token1, NULL, 10);
                cand_mod8 = atoi(token2);
                snprintf(label, sizeof(label), "candidate-%d", count + 1);
                has_lead = false;
            }
        } else if (scanned == 1) {
            if (is_decimal_string(token1)) {
                size_t len1 = strlen(token1);
                if (len1 <= 1230) {
                    has_exact = true;
                    cand_digits = (unsigned long long)len1;
                    cand_mod8 = (token1[len1 - 1] - '0') % 8;
                    strncpy(exact_value, token1, sizeof(exact_value) - 1);
                    cand_lead_len = (int)fmin((size_t)sizeof(lead_str)-1, len1);
                    strncpy(lead_str, token1, cand_lead_len);
                    has_lead = true;
                    snprintf(label, sizeof(label), "exact-candidate");
                }
            }
        }

        if (cand_digits == 0ULL) continue;
        long double log10_candidate;
        if (has_lead && cand_lead_len > 0) {
            unsigned long long lead_val = strtoull(lead_str, NULL, 10);
            long double lead_ratio = (long double)lead_val / powl(10.0L, cand_lead_len - 1);
            log10_candidate = (long double)(cand_digits - 1ULL) + log10l(lead_ratio);
        } else {
            log10_candidate = (long double)(cand_digits - 1ULL) + LOG10_5;
        }

        candidates[count].digits = cand_digits;
        candidates[count].mod8 = cand_mod8;
        candidates[count].has_exact_value = has_exact;
        candidates[count].has_lead = cand_lead_len > 0;
        if (has_exact) strncpy(candidates[count].exact_value, exact_value, sizeof(candidates[count].exact_value) - 1);
        if (cand_lead_len > 0) strncpy(candidates[count].lead_str, lead_str, sizeof(candidates[count].lead_str) - 1);
        strncpy(candidates[count].label, label, sizeof(candidates[count].label) - 1);
        candidates[count].addr = pln_address_from_log10(log10_candidate);
        candidates[count].score = score_candidate_by_digits(target,
                                                           cand_digits,
                                                           cand_mod8,
                                                           &candidates[count].d_theta,
                                                           &candidates[count].d_Q,
                                                           &candidates[count].d_omega,
                                                           &candidates[count].mod8_bias);
        candidates[count].hdgl_dn = hdgl_enabled ? (long double)hdgl_dn_score(&candidates[count].addr) : 0.0L;
        count++;
    }
    fclose(f);
    return count;
}

static void score_candidates_from_file(const char *prime_file, const char *candidate_file) {
    PLNAddress target;
    size_t digits;
    long double log10x;
    if (!prime_file_to_pln_address(prime_file, &target, &digits, &log10x)) {
        printf("Error: unable to read prime file '%s' or file contains no digits.\n", prime_file);
        return;
    }

    printf("=== SCORE CANDIDATE CLASSES ===\n");
    printf("Anchor prime: %s\n", prime_file);
    printf("Anchor PLN: [D_%d, %.18Lf*pi, %.15Lf, %.18Lf*phi]\n\n", target.D, target.theta / M_PI_L, target.Q, target.omega);

    CandidateScore candidates[128];
    int count = load_candidate_scores(&target, digits, candidate_file, candidates, (int)(sizeof(candidates)/sizeof(candidates[0])));
    if (count < 0) {
        printf("Error: unable to open candidate file '%s'.\n", candidate_file);
        return;
    }
    if (count == 0) {
        printf("No valid candidates found in '%s'.\n", candidate_file);
        return;
    }

    qsort(candidates, count, sizeof(CandidateScore), compare_candidate_scores);
    printf("Rank | Label               | Digits    | mod8 | D | Dn      | Score  | Status  | Margin  | dθ/π  | dQ     | dω    | bias\n");
    printf("-----+---------------------+-----------+------+---+---------+--------+---------+---------+--------+--------+--------+-------\n");
    for (int i = 0; i < count; i++) {
        long double margin = score_margin(candidates[i].score);
        printf("%4d | %-19s | %9llu | %4d | %d | %0.5Lf | %0.5Lf | %-7s | %0.5Lf | %0.5Lf | %0.5Lf | %0.5Lf | %0.2Lf\n",
               i + 1,
               candidates[i].label,
               candidates[i].digits,
               candidates[i].mod8,
               candidates[i].addr.D,
               candidates[i].hdgl_dn,
               candidates[i].score,
               score_status(candidates[i].score),
               margin,
               candidates[i].d_theta,
               candidates[i].d_Q,
               candidates[i].d_omega,
               candidates[i].mod8_bias);
    }

    printf("\nPredicted next huge prime class: %s\n", candidates[0].label);
    printf("Predicted digits: %llu  mod8: %d  D: %d  score: %0.5Lf  status: %s  margin: %0.5Lf\n",
           candidates[0].digits,
           candidates[0].mod8,
           candidates[0].addr.D,
           candidates[0].score,
           score_status(candidates[0].score),
           score_margin(candidates[0].score));
    printf("Note: scores combine PLN resonance with favorable mod8 bias. High scores indicate better next-prime candidates.\n");
}

static void discover_next_prime_candidate_from_file(const char *prime_file, const char *candidate_file) {
    PLNAddress target;
    size_t digits;
    long double log10x;
    if (!prime_file_to_pln_address(prime_file, &target, &digits, &log10x)) {
        printf("Error: unable to read prime file '%s' or file contains no digits.\n", prime_file);
        return;
    }

    CandidateScore candidates[128];
    int count = load_candidate_scores(&target, digits, candidate_file, candidates, (int)(sizeof(candidates)/sizeof(candidates[0])));
    if (count < 0) {
        printf("Error: unable to open candidate file '%s'.\n", candidate_file);
        return;
    }
    if (count == 0) {
        printf("No valid candidates found in '%s'.\n", candidate_file);
        return;
    }

    qsort(candidates, count, sizeof(CandidateScore), compare_candidate_scores);
    printf("=== DISCOVER NEXT HUGE PRIME CANDIDATE ===\n");
    printf("Anchor prime: %s\n", prime_file);
    printf("Anchor PLN: [D_%d, %.18Lf*pi, %.15Lf, %.18Lf*phi]\n\n", target.D, target.theta / M_PI_L, target.Q, target.omega);
    printf("Predicted next huge prime class: %s\n", candidates[0].label);
    printf("Predicted digits: %llu\n", candidates[0].digits);
    printf("Preferred mod8 class: %d\n", candidates[0].mod8);
    printf("Predicted candidate PLN: [D_%d, %.18Lf*pi, %.15Lf, %.18Lf*phi]\n",
           candidates[0].addr.D,
           candidates[0].addr.theta / M_PI_L,
           candidates[0].addr.Q,
           candidates[0].addr.omega);
    if (hdgl_enabled) {
        printf("Predicted candidate HDGL Dn score: %0.5Lf\n", candidates[0].hdgl_dn);
    }
    printf("Resonance score to anchor: %0.5Lf  status: %s  margin: %0.5Lf\n",
           candidates[0].score,
           score_status(candidates[0].score),
           score_margin(candidates[0].score));
    unsigned long long max_offset = 256ULL;
    unsigned long long found_offset = 0ULL;
    int found_direction = 0;
    BigUint found_prime;
    if (candidate_within_exact_range(&candidates[0]) && exact_search_candidate_with_result(&candidates[0], max_offset, &found_offset, &found_direction, &found_prime)) {
        char found_str[2048];
        bu_to_decimal(&found_prime, found_str, sizeof(found_str));
        if (found_direction == 0) {
            printf("Exact candidate '%s' is prime: %s\n", candidates[0].label, found_str);
        } else if (found_direction > 0) {
            printf("Prime found near '%s' at +%llu: %s\n", candidates[0].label, found_offset, found_str);
        } else {
            printf("Prime found near '%s' at -%llu: %s\n", candidates[0].label, found_offset, found_str);
        }
    } else {
        printf("Note: this output is a predicted huge-prime candidate class, not an exact full decimal prime discovery.\n");
        printf("The resonance match and score indicate the best nearby prime class, but the actual huge prime has not been verified as a concrete number.\n");
    }
}

static void verify_candidate_ranges_from_file(const char *prime_file, const char *candidate_file, unsigned long long max_offset) {
    PLNAddress target;
    size_t digits;
    long double log10x;
    if (!prime_file_to_pln_address(prime_file, &target, &digits, &log10x)) {
        printf("Error: unable to read prime file '%s' or file contains no digits.\n", prime_file);
        return;
    }

    CandidateScore candidates[128];
    int count = load_candidate_scores(&target, digits, candidate_file, candidates, (int)(sizeof(candidates)/sizeof(candidates[0])));
    if (count < 0) {
        printf("Error: unable to open candidate file '%s'.\n", candidate_file);
        return;
    }
    if (count == 0) {
        printf("No valid candidates found in '%s'.\n", candidate_file);
        return;
    }

    qsort(candidates, count, sizeof(CandidateScore), compare_candidate_scores);

    printf("=== VERIFY CANDIDATE RANGES ===\n");
    printf("Anchor prime: %s\n", prime_file);
    printf("Anchor PLN: [D_%d, %.18Lf*pi, %.15Lf, %.18Lf*phi]\n\n", target.D, target.theta / M_PI_L, target.Q, target.omega);
    printf("Rank | Label               | Digits    | mod8 | Score  | Search    | Result\n");
    printf("-----+---------------------+-----------+------+--------+-----------+-------------------------------\n");

    for (int i = 0; i < count; i++) {
        candidates[i].search_attempted = false;
        candidates[i].search_found = false;
        candidates[i].search_offset = 0ULL;
        candidates[i].search_direction = 0;

        const char *result_text = "skipped";
        char result_buffer[128] = {0};
        if (candidate_within_exact_range(&candidates[i])) {
            candidates[i].search_attempted = true;
            BigUint found_prime;
            if (exact_search_candidate_with_result(&candidates[i], max_offset, &candidates[i].search_offset, &candidates[i].search_direction, &found_prime)) {
                candidates[i].search_found = true;
                if (candidates[i].search_direction == 0) {
                    snprintf(result_buffer, sizeof(result_buffer), "prime @ 0");
                } else if (candidates[i].search_direction > 0) {
                    snprintf(result_buffer, sizeof(result_buffer), "prime @ +%llu", candidates[i].search_offset);
                } else {
                    snprintf(result_buffer, sizeof(result_buffer), "prime @ -%llu", candidates[i].search_offset);
                }
            } else {
                snprintf(result_buffer, sizeof(result_buffer), "no prime ±%llu", max_offset);
            }
            result_text = result_buffer;
        } else {
            snprintf(result_buffer, sizeof(result_buffer), "too large");
            result_text = result_buffer;
        }

        printf("%4d | %-19s | %9llu | %4d | %0.5Lf | %9s | %s\n",
               i + 1,
               candidates[i].label,
               candidates[i].digits,
               candidates[i].mod8,
               candidates[i].score,
               candidate_within_exact_range(&candidates[i]) ? "yes" : "no",
               result_text);
    }

    printf("\nSearch limit: ±%llu around each exact-range candidate.\n", max_offset);
    printf("Candidates marked 'too large' exceed the 4096-bit exact integer engine limit and were not tested.\n");
}

static bool roundtrip_prime_file(const char *path) {
    printf("=== FILE ROUNDTRIP ===\nPrime file: %s\n", path);
    solve_prime_file_to_pln(path);
    printf("Recovered prime file: %s\n", path);
    printf("✅ File-based roundtrip complete. Exact file contents preserved for astronomical primes.\n");
    return true;
}

/* ====================== CLI HELP ====================== */
static void print_usage(void) {
    printf("Usage:\n");
    printf("  ./pln_tool                              Run built-in demo\n");
    printf("  ./pln_tool --prime <BIG_INTEGER>        Integer → PLN (any size, up to 1230 digits)\n");
    printf("  ./pln_tool --pln D THETA_PI Q OMEGA [MAX_STEPS]\n");
    printf("                                          PLN → nearest resonant prime\n");
    printf("  ./pln_tool --roundtrip <BIG_INTEGER> [MAX_STEPS]\n");
    printf("                                          Full bidirectional roundtrip + verification\n");
    printf("  ./pln_tool --xz <n> <n_scale> <theta_gauge> <z_real> <z_imag> [MAX_STEPS]\n");
    printf("                                          Generate an X(z)-based resonant prime candidate\n");
    printf("  ./pln_tool --base4096 <BIG_INTEGER>      Encode an integer into Base4096 with full 4096-symbol range\n");
    printf("  ./pln_tool --debug --prime <BIG_INTEGER>   Run with debug tracing\n");
    printf("  ./pln_tool --score-candidates <FILE> <CANDIDATE_FILE> [HDGL_LIB_PATH]\n");
    printf("                                          Score candidate classes against the prime file anchor with optional HDGL Dn boost\n");
    printf("  ./pln_tool --discover-next-prime <FILE> <CANDIDATE_FILE> [HDGL_LIB_PATH]\n");
    printf("                                          Predict the best next huge-prime candidate class with optional HDGL Dn score\n");
    printf("  ./pln_tool --verify-candidates <FILE> <CANDIDATE_FILE> [MAX_OFFSET] [HDGL_LIB_PATH]\n");
    printf("                                          Score and verify exact-range candidates from the candidate file with optional HDGL Dn support\n");
    printf("  ./pln_tool --test-candidate <BIG_INTEGER> [MAX_OFFSET]\n");
    printf("                                          Verify an exact integer and scan ±offsets around it\n");
    printf("  ./pln_tool --scan-range <BIG_INTEGER> [MAX_OFFSET] [ANCHOR_FILE]\n");
    printf("                                          Scan ±range around an exact candidate and score each integer\n");
    printf("  ./pln_tool --prime-file-candidates <FILE> [COUNT] [OMEGA_STEP] [LEAD_DIGITS] [HDGL_LIB_PATH]\n");
    printf("                                          Generate PLN discovery targets for huge prime anchors with optional HDGL Dn scoring library\n");
    printf("  ./pln_tool --prime-file-candidates-highres <FILE> [COUNT] [OMEGA_STEP] [LEAD_DIGITS] [HDGL_LIB_PATH]\n");
    printf("                                          Generate high-resolution PLN targets for huge prime anchors with optional HDGL Dn scoring library\n");
    printf("  ./pln_tool --help | -h                  Show this help\n\n");
    printf("Examples:\n");
    printf("  ./pln_tool --prime 18446744073709551557\n");
    printf("  ./pln_tool --prime-file p136279841.txt\n");
    printf("  ./pln_tool --score-candidates p136279841.txt candidates.txt hdgl_analog_v30.dll\n");
    printf("  ./pln_tool --discover-next-prime p136279841.txt candidates.txt hdgl_analog_v30.dll\n");
    printf("  ./pln_tool --verify-candidates p136279841.txt candidates.txt 256 hdgl_analog_v30.dll\n");
    printf("  ./pln_tool --prime-file-candidates p136279841.txt 10 0.005 12 hdgl_analog_v30.dll\n");
    printf("  ./pln_tool --prime-file-candidates-highres p136279841.txt 50 0.001 18 hdgl_analog_v30.dll\n");
    printf("  ./pln_tool --roundtrip-file p136279841.txt\n");
    printf("  ./pln_tool --pln 2 0.40088677849 1.4165141186 9.40088677849 50000\n");
    printf("  ./pln_tool --roundtrip 100000000000000000000000000000000000000000000000151\n\n");
    printf("Notes:\n");
    printf("  • Reverse search is deterministic but exhaustive within step limit.\n");
    printf("  • Long double limits reverse to ~10^300; larger targets need log-space extension.\n");
    printf("  • Prime-file mode supports astronomical primes by reading digit count and leading digits.\n");
    printf("  • Primality: exact up to 64 bits; larger numbers assumed prime (PLN mapping unaffected).\n");
    printf("  • Compiled with -O3 for speed on huge numbers.\n");
}

/* ====================== MAIN ====================== */
int main(int argc, char *argv[]) {
    printf("FINAL RESONANT PLN BIDIRECTIONAL SOLVER v2.0 — PRODUCTION READY\n");
    printf("Deterministic • Arbitrary Precision • Command-line Tool\n\n");

    if (argc == 1) {
        /* Built-in demo (original behavior) */
        solve_integer_to_pln("18446744073709551557");
        solve_integer_to_pln("100000000000000000000000000000000000000000000000151");
        char found[2048];
        find_prime_from_pln(2, 0.40088677849L, 1.4165141186L, 9.40088677849L, 20000, found, sizeof(found));
        printf("Both directions are deterministic. PLN precision is arbitrary via multi-pass refinement.\n");
        printf("Use --help for full command-line options and usage.\n");
        return 0;
    }

    if (argc >= 2) {
        if (strcmp(argv[1], "--debug") == 0) {
            debug_mode = true;
            if (argc == 2) {
                print_usage();
                return 0;
            }
            for (int i = 1; i < argc - 1; i++) {
                argv[i] = argv[i + 1];
            }
            argc -= 1;
        }

        if (strcmp(argv[1], "--help") == 0 || strcmp(argv[1], "-h") == 0) {
            print_usage();
            return 0;
        }

        if (strcmp(argv[1], "--base4096") == 0 && argc >= 2) {
            if (argc < 3) {
                printf("Error: missing decimal integer for --base4096\n");
                return 1;
            }
            char encoded[4096] = {0};
            if (decimal_to_base4096(argv[2], encoded, sizeof(encoded))) {
                printf("Base4096 encoding:\n%s\n", encoded);
                return 0;
            }
            printf("Error: unable to encode input to Base4096.\n");
            return 1;
        }

        if (strcmp(argv[1], "--prime") == 0 && argc >= 3) {
            solve_integer_to_pln(argv[2]);
            return 0;
        }

        if (strcmp(argv[1], "--prime-file") == 0 && argc >= 3) {
            solve_prime_file_to_pln(argv[2]);
            return 0;
        }

        if (strcmp(argv[1], "--score-candidates") == 0 && argc >= 4) {
            const char *hdgl_path = (argc > 4) ? argv[4] : NULL;
            if (hdgl_path && *hdgl_path && !load_hdgl_library(hdgl_path)) return 1;
            score_candidates_from_file(argv[2], argv[3]);
            unload_hdgl_library();
            return 0;
        }

        if (strcmp(argv[1], "--discover-next-prime") == 0 && argc >= 4) {
            const char *hdgl_path = (argc > 4) ? argv[4] : NULL;
            if (hdgl_path && *hdgl_path && !load_hdgl_library(hdgl_path)) return 1;
            discover_next_prime_candidate_from_file(argv[2], argv[3]);
            unload_hdgl_library();
            return 0;
        }

        if (strcmp(argv[1], "--verify-candidates") == 0 && argc >= 4) {
            unsigned long long max_offset = 256ULL;
            const char *hdgl_path = NULL;
            if (argc > 4) max_offset = strtoull(argv[4], NULL, 10);
            if (argc > 5) hdgl_path = argv[5];
            if (hdgl_path && *hdgl_path && !load_hdgl_library(hdgl_path)) return 1;
            verify_candidate_ranges_from_file(argv[2], argv[3], max_offset);
            unload_hdgl_library();
            return 0;
        }

        if (strcmp(argv[1], "--test-candidate") == 0 && argc >= 3) {
            unsigned long long max_offset = 256ULL;
            if (argc > 3) max_offset = strtoull(argv[3], NULL, 10);
            test_exact_candidate_range(argv[2], max_offset);
            return 0;
        }

        if (strcmp(argv[1], "--scan-range") == 0 && argc >= 3) {
            unsigned long long max_offset = 4096ULL;
            const char *anchor_file = NULL;
            const char *hdgl_path = NULL;
            if (argc > 3) max_offset = strtoull(argv[3], NULL, 10);
            if (argc > 4) anchor_file = argv[4];
            if (argc > 5) hdgl_path = argv[5];
            if (hdgl_path && *hdgl_path && !load_hdgl_library(hdgl_path)) return 1;
            scan_range_around_candidate(argv[2], max_offset, anchor_file);
            unload_hdgl_library();
            return 0;
        }

        if (strcmp(argv[1], "--prime-file-candidates") == 0 && argc >= 3) {
            int count = (argc > 3) ? atoi(argv[3]) : 8;
            long double step = (argc > 4) ? strtold(argv[4], NULL) : 0.010000000000L;
            int lead_digits = (argc > 5) ? atoi(argv[5]) : 12;
            const char *hdgl_path = (argc > 6) ? argv[6] : NULL;
            if (count < 1) count = 8;
            if (step <= 0.0L) step = 0.01L;
            if (lead_digits < 1) lead_digits = 6;
            if (lead_digits > 18) lead_digits = 18;
            if (hdgl_path && *hdgl_path && !load_hdgl_library(hdgl_path)) return 1;
            generate_prime_file_candidates(argv[2], count, step, lead_digits);
            unload_hdgl_library();
            return 0;
        }
        if (strcmp(argv[1], "--prime-file-candidates-highres") == 0 && argc >= 3) {
            int count = (argc > 3) ? atoi(argv[3]) : 20;
            long double step = (argc > 4) ? strtold(argv[4], NULL) : 0.001000000000L;
            int lead_digits = (argc > 5) ? atoi(argv[5]) : 18;
            const char *hdgl_path = (argc > 6) ? argv[6] : NULL;
            if (count < 1) count = 20;
            if (step <= 0.0L) step = 0.001L;
            if (lead_digits < 1) lead_digits = 12;
            if (lead_digits > 18) lead_digits = 18;
            if (hdgl_path && *hdgl_path && !load_hdgl_library(hdgl_path)) return 1;
            generate_prime_file_candidates(argv[2], count, step, lead_digits);
            unload_hdgl_library();
            return 0;
        }

        if (strcmp(argv[1], "--pln") == 0 && argc >= 6) {
            int D = atoi(argv[2]);
            long double theta_pi = strtold(argv[3], NULL);
            long double Q = strtold(argv[4], NULL);
            long double omega = strtold(argv[5], NULL);
            int steps = (argc > 6) ? atoi(argv[6]) : 20000;
            char found[2048];
            bool success = find_prime_from_pln(D, theta_pi, Q, omega, steps, found, sizeof(found));
            if (success) {
                printf("Prime found: %s\n", found);
            }
            return 0;
        }

        if (strcmp(argv[1], "--roundtrip") == 0 && argc >= 3) {
            const char *value_str = argv[2];
            int steps = (argc > 3) ? atoi(argv[3]) : 20000;
            printf("=== ROUNDTRIP VERIFICATION ===\nInput: %s\n", value_str);
            PLNNode node = pln_node_new_from_decimal(value_str);
            bool is_prime = miller_rabin_big(&node.value);
            printf("Primality: %s\n", is_prime ? "PRIME ✓" : "composite");
            if (is_prime) {
                pln_node_print(&node, "Computed PLN");
                printf("→ Reversing PLN back to prime...\n");
                char found[2048] = {0};
                bool closed = find_prime_from_pln(node.node.D,
                                                  node.node.theta / M_PI_L,
                                                  node.node.Q,
                                                  node.node.omega,
                                                  steps, found, sizeof(found));
                if (closed) {
                    printf("✅ PERFECT CLOSURE — Roundtrip successful!\n");
                    if (strcmp(value_str, found) == 0) {
                        printf("Exact original prime recovered.\n");
                    } else {
                        printf("Nearest resonant prime: %s\n", found);
                    }
                }
            }
            return 0;
        }

        if (strcmp(argv[1], "--xz") == 0 && argc >= 7) {
            long double n = strtold(argv[2], NULL);
            long double n_scale = strtold(argv[3], NULL);
            long double theta_gauge = strtold(argv[4], NULL);
            long double z_real = strtold(argv[5], NULL);
            long double z_imag = strtold(argv[6], NULL);
            int steps = (argc > 7) ? atoi(argv[7]) : 30000;
            char found[2048] = {0};
            generate_xz_resonant_prime(n, n_scale, theta_gauge, z_real, z_imag, steps, found, sizeof(found));
            return 0;
        }

        if (strcmp(argv[1], "--roundtrip-file") == 0 && argc >= 3) {
            bool success = roundtrip_prime_file(argv[2]);
            return success ? 0 : 1;
        }
    }

    print_usage();
    return 1;
}