/*
 * pln_final.c
 * =========================================================
 * RESONANT PLN PRIME SOLVER — FINAL UNIFIED EDITION
 *
 * Architecture:
 *   • BigUint 4096-bit   — exact integers, MR primality (≤1230 digits)
 *   • APA 16,777,216-bit — log-space coordinate arithmetic (super-res)
 *   • Magnet stack       — trial div → Lucas → Fermat → Jacobi → PLN → MR
 *   • PLN coordinate     — continuous D_cont orbit, orbit-minima guidance
 *   • HDGL Dn(r)         — inline resonance-peak formula + dlopen override
 *   • Gap predictor      — gap ≈ ln(x) × |R| / (1/φ)
 *   • X(z) generator     — φ-lattice + U(1) gauge phase
 *   • File mode          — astronomical primes via log-space only
 *   • Base4096 output    — encodes found primes ≤1230 digits
 *
 * Compile (fast, exact-only):
 *   gcc -O3 -Wall -Wextra -std=c99 -o pln pln_final.c -lm -ldl
 *
 * Compile (super-resolution APA, ~2 MiB per float, slow):
 *   gcc -O3 -DSUPER_RESOLUTION_MODE=1 -o pln pln_final.c -lm -ldl
 *
 * Usage:
 *   ./pln                                     built-in demo
 *   ./pln --prime <N>                         integer → PLN
 *   ./pln --pln D THETA Q OMEGA [STEPS]       PLN → prime
 *   ./pln --roundtrip <N> [STEPS]             bidirectional verify
 *   ./pln --xz n ns tg zr zi [STEPS]         X(z) generator
 *   ./pln --prime-file <FILE>                 file prime → PLN
 *   ./pln --prime-file-candidates <F> [N] [STEP] [LEAD] [SO]
 *   ./pln --score-candidates <F> <CANDS> [SO]
 *   ./pln --discover-next-prime <F> <CANDS> [SO]
 *   ./pln --test-candidate <N> [OFFSET]
 *   ./pln --scan-range <N> [OFFSET] [ANCHOR]
 *   ./pln --base4096 <N>
 *   ./pln --debug --prime <N>
 *   ./pln --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>
#include <time.h>

#ifdef _WIN32
  #include <windows.h>
  typedef HMODULE DLHandle;
  #define DL_OPEN(p)    LoadLibraryA(p)
  #define DL_SYM(h,s)   GetProcAddress(h,s)
  #define DL_CLOSE(h)   FreeLibrary(h)
  #define DL_ERR()      "LoadLibrary error"
#else
  #include <dlfcn.h>
  typedef void* DLHandle;
  #define DL_OPEN(p)    dlopen(p, RTLD_NOW)
  #define DL_SYM(h,s)   dlsym(h,s)
  #define DL_CLOSE(h)   dlclose(h)
  #define DL_ERR()      dlerror()
#endif

static bool debug_mode = false;
static void dbg(const char *fmt, ...) {
    if (!debug_mode) return;
    va_list a; va_start(a,fmt); vprintf(fmt,a); va_end(a);
}

/* ================================================================
 * 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_PHI    0.208987640249978733769L
#define LOG10_5      0.698970004336018804786L
#define M_PI_L       3.14159265358979323846L
#define TWO_PI_L     6.2831853071795864769L
#define INV_PHI_D    0.6180339887498948482
#define PHI_D        1.6180339887498948482
#define LN_PHI_D     0.4812118250596034748
#define R_MIN_THRESH 0.15L
#define MAX_PASSES   5
#define MAX_FILE_LEAD_DIGITS 18
#define PLN_SCORE_PASS   0.60L
#define PLN_SCORE_BORDER 0.55L

/* ================================================================
 * BIGUINT — 4096-bit (128 × 32-bit words)
 * ================================================================ */
#define BU_WORDS 128

typedef struct { int len; uint32_t d[BU_WORDS]; } BigUint;

static void bu_init(BigUint *a)         { a->len=0; memset(a->d,0,sizeof(a->d)); }
static void bu_normalize(BigUint *a)    { while(a->len>0&&a->d[a->len-1]==0) a->len--; }
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->d[0]&1u)==0; }
static bool bu_is_one(const BigUint *a) { return a->len==1&&a->d[0]==1; }
static void bu_copy(BigUint *dst, const BigUint *src) {
    dst->len=src->len; memcpy(dst->d,src->d,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->d[i]!=b->d[i]) return a->d[i]<b->d[i]?-1:1;
    return 0;
}
static bool bu_equal(const BigUint *a, const BigUint *b){ return bu_cmp(a,b)==0; }
static int bu_bit_length(const BigUint *a) {
    if(!a->len) return 0;
    uint32_t top=a->d[a->len-1]; int bits=32*(a->len-1);
    while(top){bits++;top>>=1;} return bits;
}
static void bu_from_u64(BigUint *a, uint64_t v) {
    bu_init(a);
    if(!v) return;
    a->d[0]=(uint32_t)(v&0xFFFFFFFFu);
    a->d[1]=(uint32_t)(v>>32);
    a->len=a->d[1]?2:1;
}
static uint64_t bu_to_u64(const BigUint *a) {
    if(!a->len) return 0;
    if(a->len==1) return a->d[0];
    return a->d[0]|((uint64_t)a->d[1]<<32);
}
static void bu_shr1(BigUint *a) {
    uint64_t carry=0;
    for(int i=a->len-1;i>=0;i--){
        uint64_t v=(carry<<32)|a->d[i];
        a->d[i]=(uint32_t)(v>>1); carry=v&1;
    }
    bu_normalize(a);
}
static void bu_shl_bits(BigUint *a, int bits) {
    if(bits<=0||bu_is_zero(a)) return;
    int ws=bits/32, bs=bits%32;
    if(ws>0){
        int cap=BU_WORDS-a->len; if(ws>cap) ws=cap;
        for(int i=a->len-1;i>=0;i--) a->d[i+ws]=a->d[i];
        for(int i=0;i<ws;i++) a->d[i]=0;
        a->len+=ws;
    }
    if(bs>0){
        uint64_t carry=0;
        for(int i=ws;i<a->len;i++){
            uint64_t v=((uint64_t)a->d[i]<<bs)|carry;
            a->d[i]=(uint32_t)v; carry=v>>32;
        }
        if(carry&&a->len<BU_WORDS) a->d[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 ws=bits/32, bs=bits%32;
    if(ws>=a->len){bu_init(a);return;}
    if(ws>0){
        for(int i=0;i<a->len-ws;i++) a->d[i]=a->d[i+ws];
        for(int i=a->len-ws;i<a->len;i++) a->d[i]=0;
        a->len-=ws;
    }
    if(bs>0){
        uint64_t carry=0;
        for(int i=a->len-1;i>=0;i--){
            uint64_t v=((uint64_t)a->d[i]>>bs)|(carry<<(32-bs));
            carry=a->d[i]&(((uint64_t)1<<bs)-1);
            a->d[i]=(uint32_t)v;
        }
    }
    bu_normalize(a);
}
static void bu_add(BigUint *res, const BigUint *a, const BigUint *b) {
    uint64_t carry=0;
    int ml=a->len>b->len?a->len:b->len;
    int i=0;
    for(;i<ml||carry;i++){
        uint64_t av=i<a->len?a->d[i]:0, bv=i<b->len?b->d[i]:0;
        uint64_t s=av+bv+carry;
        if(i<BU_WORDS) res->d[i]=(uint32_t)s;
        carry=s>>32;
    }
    res->len=i;
    bu_normalize(res);
}
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->d[i], bv=i<b->len?b->d[i]:0;
        uint64_t diff=av-bv-borrow;
        res->d[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_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 v=(i<a->len?(uint64_t)a->d[i]:0)*b+carry;
        if(i<BU_WORDS) res->d[i]=(uint32_t)v; carry=v>>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 s=(uint64_t)a->d[i]+carry;
        res->d[i]=(uint32_t)s; carry=s>>32;
    }
    for(;i<a->len;i++) res->d[i]=a->d[i];
    if(carry&&i<BU_WORDS) res->d[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) {
    BigUint r={0}; uint64_t carry=0;
    for(int i=num->len-1;i>=0;i--){
        uint64_t v=(carry<<32)|num->d[i];
        r.d[i]=(uint32_t)(v/den); carry=v%den;
    }
    r.len=num->len; bu_normalize(&r);
    if(quot) bu_copy(quot,&r);
    if(rem) *rem=(uint32_t)carry;
}
static void bu_from_decimal(BigUint *a, const char *s) {
    bu_init(a);
    for(const char *p=s;*p;p++){
        if(*p<'0'||*p>'9') continue;
        BigUint t; bu_mul_small(&t,a,10);
        bu_add_small(&t,&t,(uint32_t)(*p-'0'));
        bu_copy(a,&t);
    }
}
static void bu_to_decimal(const BigUint *a, char *buf, size_t bufsz) {
    if(bu_is_zero(a)){snprintf(buf,bufsz,"0");return;}
    char tmp[1500]={0}; int pos=0;
    BigUint cp; bu_copy(&cp,a);
    while(!bu_is_zero(&cp)){
        uint32_t r; BigUint nxt;
        bu_divmod_small(&nxt,&cp,10,&r);
        if(pos<(int)sizeof(tmp)-1) tmp[pos++]='0'+(char)r;
        bu_copy(&cp,&nxt);
    }
    int idx=0;
    for(int i=pos-1;i>=0&&idx<(int)bufsz-1;i--) buf[idx++]=tmp[i];
    buf[idx]='\0';
}
static void bu_add_u64(BigUint *res, const BigUint *a, uint64_t b) {
    if(!b){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){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_mod_add(BigUint *res, const BigUint *a, const BigUint *b, const BigUint *mod) {
    BigUint s; bu_add(&s,a,b);
    if(bu_cmp(&s,mod)>=0) bu_sub(&s,&s,mod);
    bu_copy(res,&s);
}
static void bu_modmul(const BigUint *a, const BigUint *b, const BigUint *mod, BigUint *res) {
    BigUint result; bu_init(&result);
    BigUint tmp; bu_copy(&tmp,a);
    for(int i=0;i<b->len;i++){
        uint32_t word=b->d[i];
        for(int bit=0;bit<32;bit++){
            if(word&1u) bu_mod_add(&result,&result,&tmp,mod);
            word>>=1;
            bu_mod_add(&tmp,&tmp,&tmp,mod);
        }
    }
    bu_copy(res,&result);
}
static void bu_modexp(const BigUint *base, const BigUint *exp,
                      const BigUint *mod, BigUint *res) {
    BigUint 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.d[0]&1u) bu_modmul(&result,&power,mod,&result);
        bu_modmul(&power,&power,mod,&power);
        bu_shr1(&exponent);
    }
    bu_copy(res,&result);
}

/* ================================================================
 * APA FLOAT — 16,777,216-bit mantissa (super-res) or 512-bit (fast)
 * ================================================================ */
#ifndef SUPER_RESOLUTION_MODE
#define SUPER_RESOLUTION_MODE 0
#endif
#if SUPER_RESOLUTION_MODE
  #define APA_WORDS 262144ULL   /* 16,777,216 bits */
#else
  #define APA_WORDS 8           /* 512 bits */
#endif
#define APA_MSB (1ULL<<63)

typedef struct {
    uint64_t m[APA_WORDS];
    int64_t  exp;
    int8_t   sign;
} APAFloat;

static void apa_shr(uint64_t *m, size_t nw, int64_t sh) {
    if(sh<=0||!nw) return;
    if(sh>=(int64_t)nw*64){memset(m,0,nw*8);return;}
    int64_t ws=sh/64; int bs=(int)(sh%64);
    if(ws) for(size_t i=nw;i-->(size_t)ws;) m[i]=m[i-ws];
    if(ws) memset(m,0,ws*8);
    if(bs){ int rs=64-bs;
        for(size_t i=nw;i-->0;){
            uint64_t c=i>0?m[i-1]<<rs:0;
            m[i]=(m[i]>>bs)|c;
        }
    }
}
static void apa_normalize(APAFloat *a) {
    if(!a->m[0]){a->exp=0;return;}
    while(!(a->m[0]&APA_MSB)){
        uint64_t carry=0;
        for(size_t i=APA_WORDS;i-->0;){
            uint64_t nc=(a->m[i]&APA_MSB)?1:0;
            a->m[i]=(a->m[i]<<1)|carry; carry=nc;
        }
        a->exp--;
    }
}
static void apa_init(APAFloat *a, long double v) {
    memset(a->m,0,sizeof(a->m));
    a->sign=(v>=0.0L)?1:-1; v=fabsl(v);
    if(v==0.0L){a->exp=0;return;}
    int e; long double mant=frexpl(v,&e); a->exp=e; mant*=2.0L;
    for(int i=0;i<(int)(APA_WORDS*64);i++){
        mant*=2.0L;
        if(mant>=2.0L){
            int wi=i/64,bi=63-(i%64);
            if(wi<(int)APA_WORDS) a->m[wi]|=(1ULL<<bi);
            mant-=2.0L;
        }
    }
    apa_normalize(a);
}
static long double apa_log2(const APAFloat *a) {
    if(!a->m[0]) return -INFINITY;
    long double hi=(long double)a->m[0];
    long double lo=APA_WORDS>1?(long double)a->m[1]/18446744073709551616.0L:0.0L;
    long double top=(hi+lo)/9223372036854775808.0L;
    return (long double)a->exp-1+log2l(top);
}
static void apa_add(APAFloat *res, const APAFloat *a, const APAFloat *b) {
    APAFloat ta=*a, tb=*b;
    if(ta.exp>tb.exp){ apa_shr(tb.m,APA_WORDS,ta.exp-tb.exp); tb.exp=ta.exp; *res=ta; }
    else              { apa_shr(ta.m,APA_WORDS,tb.exp-ta.exp); ta.exp=tb.exp; *res=tb; }
    uint64_t carry=0;
    for(size_t i=APA_WORDS;i-->0;){
        __uint128_t s=(__uint128_t)res->m[i]+ta.m[i]+tb.m[i]+carry;
        res->m[i]=(uint64_t)s; carry=(uint64_t)(s>>64);
    }
    if(carry){ res->exp++; apa_shr(res->m,APA_WORDS,1); res->m[0]|=APA_MSB; }
    apa_normalize(res);
}

/* ================================================================
 * BIGUINT ↔ APA CONVERSIONS
 * ================================================================ */
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){out.exp=0;return out;}
    BigUint tmp; bu_copy(&tmp,a);
    int tbits=(int)(APA_WORDS*64);
    int sh=bits-tbits;
    if(sh>0) bu_shr_bits(&tmp,sh); else if(sh<0) bu_shl_bits(&tmp,-sh);
    out.exp=bits;
    int wi=0;
    for(int i=tmp.len-1;i>=0;i-=2){
        uint64_t hi=tmp.d[i], lo=(i-1>=0)?tmp.d[i-1]:0;
        if(wi<(int)APA_WORDS) out.m[wi++]=(hi<<32)|lo;
    }
    apa_normalize(&out); return out;
}
static long double bu_high_precision_log2(const BigUint *a) {
    APAFloat ap=bu_to_apa(a); return apa_log2(&ap);
}
static void bu_from_log2_approx(BigUint *out, long double log2x) {
    if(log2x<=0.0L){bu_init(out);return;}
    long double ip=floorl(log2x), frac=log2x-ip;
    int bitlen=(int)ip+1;
    int ti=(bitlen-1)/32, ts=(bitlen-1)%32;
    long double mant=powl(2.0L,frac);
    /* guard: ts can be 0..31 */
    unsigned long long tv;
    if(ts>0) tv=(unsigned long long)(mant*(1ULL<<ts));
    else     tv=(unsigned long long)mant;
    long double lf= ts>0 ? mant*(1ULL<<ts)-(long double)tv : 0.0L;
    uint32_t lw=(uint32_t)(lf*4294967296.0L);
    bu_init(out);
    if(ti>=BU_WORDS) ti=BU_WORDS-1;
    out->len=ti+1;
    for(int i=0;i<out->len;i++) out->d[i]=0;
    out->d[ti]=(uint32_t)(tv&0xFFFFFFFFu);
    if(ti>0) out->d[ti-1]=lw;
    bu_normalize(out);
}

/* ================================================================
 * MILLER-RABIN
 *   ≤64-bit:  deterministic with 12 witnesses
 *   ≤128-bit: deterministic with u128, 12 witnesses
 *   >128-bit: BigUint path, 7 strong bases
 * ================================================================ */
static uint64_t mulmod64(uint64_t a, uint64_t b, uint64_t m) {
    uint64_t res=0; a%=m;
    while(b){ if(b&1){ if(res>=m-a) res=res+a-m; else res+=a; } b>>=1;
        if(b){ if(a>=m-a) a=a+a-m; else a+=a; } }
    return res;
}
static uint64_t modpow64(uint64_t b, uint64_t e, uint64_t m) {
    if(m==1) return 0;
    uint64_t r=1%m, p=b%m;
    while(e){ if(e&1) r=mulmod64(r,p,m); p=mulmod64(p,p,m); e>>=1; }
    return r;
}
typedef unsigned __int128 u128;
static u128 pmm128(u128 b, u128 e, u128 m){ u128 r=1; b%=m;
    while(e){if(e&1)r=r*b%m; b=b*b%m; e>>=1;} return r; }

static bool miller_rabin_big(const BigUint *n) {
    int bits=bu_bit_length(n);
    if(bits<=1) return false;
    /* 64-bit exact path */
    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)||v%3==0||v%5==0) return false;
        uint64_t nm1=v-1, d=nm1; int s=0;
        while(!(d&1)){d>>=1;s++;}
        static const uint64_t W[12]={2,3,5,7,11,13,23,29,31,37,41,43};
        for(int i=0;i<12;i++){
            uint64_t a=W[i]; if(a>=v) break;
            uint64_t x=modpow64(a,d,v);
            if(x==1||x==nm1) continue;
            bool c=true;
            for(int r=1;r<s;r++){x=mulmod64(x,x,v);if(x==nm1){c=false;break;}}
            if(c) return false;
        }
        return true;
    }
    /* 128-bit exact path */
    if(bits<=128){
        u128 v=(u128)bu_to_u64(n)|(n->len>2?(u128)((uint64_t)n->d[2]|(n->len>3?(uint64_t)n->d[3]<<32:0))<<64:0);
        /* rebuild properly */
        v=0; for(int i=n->len-1;i>=0;i--){ v=(v<<32)|n->d[i]; }
        if(v<2) return false;
        if(v==2||v==3) return true;
        if(!(v&1)) return false;
        u128 nm1=v-1,d=nm1; int s=0;
        while(!(d&1)){d>>=1;s++;}
        static const u128 W[12]={2,3,5,7,11,13,17,19,23,29,31,37};
        for(int i=0;i<12;i++){
            u128 a=W[i]; if(a>=v) continue;
            u128 x=pmm128(a,d,v);
            if(x==1||x==nm1) continue;
            bool c=true;
            for(int r=1;r<s;r++){x=x*x%v;if(x==nm1){c=false;break;}}
            if(c) return false;
        }
        return true;
    }
    /* BigUint path — small prime trial first */
    static const uint32_t sp[]={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(sp)/sizeof(sp[0]);i++){
        uint32_t r; bu_divmod_small(NULL,n,sp[i],&r); if(!r) 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]={2,325,9375,28178,450775,9780504,1795265022};
    for(int i=0;i<7;i++){
        BigUint a; bu_from_u64(&a,bases[i]);
        BigUint x; bu_modexp(&a,&d,n,&x);
        if(bu_is_one(&x)||bu_equal(&x,&nm1)) continue;
        bool composite=true;
        for(int r=1;r<s;r++){
            bu_modmul(&x,&x,n,&x);
            if(bu_equal(&x,&nm1)){composite=false;break;}
        }
        if(composite) return false;
    }
    return true;
}

/* ================================================================
 * MAGNET FILTER STACK
 * Ordered cheapest → most expensive.
 * Each returns true = "still a candidate", false = "composite, skip".
 * ================================================================ */

/* Small primes up to 1499 (245 primes) */
static const uint16_t SPRIMES[]={
    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,101,103,107,109,113,127,131,137,139,149,151,
    157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,
    239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,
    331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,
    421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,
    509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,
    613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,
    709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,
    821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,
    919,929,937,941,947,953,967,971,977,983,991,997,1009,1013,
    1019,1021,1031,1033,1039,1049,1051,1061,1063,1069,1087,1091,
    1093,1097,1103,1109,1117,1123,1129,1151,1153,1163,1171,1181,
    1187,1193,1201,1213,1217,1223,1229,1231,1237,1249,1259,1277,
    1279,1283,1289,1291,1297,1301,1303,1307,1319,1321,1327,1361,
    1367,1373,1381,1399,1409,1423,1427,1429,1433,1439,1447,1451,
    1453,1459,1471,1481,1483,1487,1489,1493,1499
};
#define N_SPRIMES (int)(sizeof(SPRIMES)/sizeof(SPRIMES[0]))

/* Lucas sieve values */
static const int LUCAS_VALS[]={3,4,7,11,18,29,47,76,123};
#define N_LUCAS 9

static bool magnet_trial(const BigUint *n) {
    int bits=bu_bit_length(n);
    for(int i=0;i<N_SPRIMES;i++){
        uint32_t r; bu_divmod_small(NULL,n,SPRIMES[i],&r);
        if(!r){
            /* it IS the small prime itself — keep it */
            if(bits<=16&&bu_to_u64(n)==SPRIMES[i]) return true;
            return false;
        }
    }
    return true;
}

static bool magnet_lucas(const BigUint *n) {
    if(n->len==0) return false;
    uint64_t v=bu_to_u64(n);
    if(n->len<=2){
        /* small enough to check directly */
        if(v<2) return false;
        if(v>2&&!(v&1)) return false;
        for(int i=0;i<N_LUCAS;i++){
            int L=LUCAS_VALS[i];
            if(v>(uint64_t)L&&v%(uint64_t)L==0) return false;
        }
        return true;
    }
    /* large: just check mod each Lucas value */
    for(int i=0;i<N_LUCAS;i++){
        uint32_t r; bu_divmod_small(NULL,n,(uint32_t)LUCAS_VALS[i],&r);
        if(!r) return false;
    }
    return true;
}

static bool magnet_fermat2(const BigUint *n) {
    /* 2^(n-1) mod n == 1 */
    BigUint exp; bu_copy(&exp,n); bu_sub_one(&exp,n);
    BigUint base; bu_from_u64(&base,2);
    BigUint res; bu_modexp(&base,&exp,n,&res);
    return bu_is_one(&res);
}

/* Jacobi symbol (a/n), n must be odd positive */
static int jacobi(const BigUint *a_in, const BigUint *n) {
    if(n->len==0) return 0;
    BigUint a; bu_copy(&a,a_in);
    BigUint nn; bu_copy(&nn,n);
    int result=1;
    while(!bu_is_zero(&a)){
        /* a ← a mod n */

        /* simpler: just reduce a mod nn using repeated subtraction for small */
        /* For BigUint: compute a mod nn */
        while(bu_cmp(&a,&nn)>=0) bu_sub(&a,&a,&nn);
        if(bu_is_zero(&a)) break;
        /* factor out 2s from a */
        int twos=0;
        while(bu_is_even(&a)){bu_shr1(&a);twos++;}
        if(twos&1){
            /* adjust result by (2/n) = (-1)^((n^2-1)/8) */
            uint32_t nmod8=nn.d[0]&7u;
            if(nmod8==3||nmod8==5) result=-result;
        }
        if(bu_is_zero(&a)) break;
        /* quadratic reciprocity */
        uint32_t amod4=a.d[0]&3u;
        uint32_t nmod4=nn.d[0]&3u;
        if(amod4==3&&nmod4==3) result=-result;
        /* swap */
        BigUint swap; bu_copy(&swap,&a);
        bu_copy(&a,&nn);
        bu_copy(&nn,&swap);
        /* a ← a mod nn */
        while(bu_cmp(&a,&nn)>=0) bu_sub(&a,&a,&nn);
    }
    if(bu_is_one(&nn)) return result;
    return 0;
}

static bool magnet_euler_jacobi(const BigUint *n) {
    /* Check: 2^((n-1)/2) mod n == Jacobi(2,n) */
    BigUint exp; bu_copy(&exp,n); bu_sub_one(&exp,n); bu_shr1(&exp);
    BigUint base; bu_from_u64(&base,2);
    BigUint res; bu_modexp(&base,&exp,n,&res);
    BigUint two; bu_from_u64(&two,2);
    int j=jacobi(&two,n);
    if(j==1)  return bu_is_one(&res);
    if(j==-1){
        BigUint nm1; bu_copy(&nm1,n); bu_sub_one(&nm1,n);
        return bu_equal(&res,&nm1);
    }
    return false;
}

/* ================================================================
 * PLN COORDINATE SYSTEM
 * ================================================================ */
typedef struct { 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 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 PLNAddress pln_from_biguint(const BigUint *x) {
    if(bu_is_zero(x)) return (PLNAddress){1,0.0L,1.0L,0.0L};
    long double omega=0.0L, prev=0.0L;
    for(int pass=0;pass<MAX_PASSES;pass++){
        long double log2x=bu_high_precision_log2(x);
        long double lnx=log2x*LN2;
        long double arg=lnx/LN_PHI;
        if(arg<=1.0L) return (PLNAddress){1,0.0L,1.0L,0.0L};
        omega=logl(arg)/LN_PHI;
        if(pass>0) omega-=(omega-prev)*0.5L;
        prev=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 PLNAddress pln_from_log2(long double log2x) {
    long double omega=omega_from_log2(log2x);
    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 long double pln_resonance_score(const PLNAddress *a, const PLNAddress *b) {
    long double dt=fabsl(a->theta-b->theta);
    if(dt>M_PI_L) dt=TWO_PI_L-dt;
    return 0.5L/(1.0L+dt/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 double predict_gap(long double x, double orb) {
    return (double)logl(x)*orb/INV_PHI_D;
}

static void find_orbit_minima(double n0, double n1, double step,
                               double *mins, int *nm, int maxm) {
    *nm=0;
    double prev=orbit_r((long double)n0), pprev=prev;
    for(double n=n0+step;n<=n1&&*nm<maxm;n+=step){
        double r=(double)orbit_r((long double)n);
        if(prev<pprev&&prev<r){ mins[(*nm)++]=n-step; }
        pprev=prev; prev=r;
    }
}

/* ================================================================
 * HDGL Dn(r) — inline resonance-peak formula + optional dlopen
 * ================================================================ */
#define NUM_DN 8
static const uint64_t FIB_T[NUM_DN]={1,1,2,3,5,8,13,21};
static const uint64_t PRM_T[NUM_DN]={2,3,5,7,11,13,17,19};

static DLHandle  hdgl_handle=NULL;
static double  (*hdgl_fn)(int,double,double)=NULL;
static bool      hdgl_active=false;
static double    hdgl_omega=1.0;

/* Inline formula: resonance peak at |R| ≈ 1/φ */
static double dn_inline(int n, double r, double omega) {
    if(n<1||n>NUM_DN) return 0.0;
    int idx=n-1;
    double Fn=(double)FIB_T[idx], Pn=(double)PRM_T[idx];
    double scale=sqrt(PHI_D*Fn*Pn);
    /* Gaussian resonance centered on golden orbit |R|=1/φ */
    double dr=r-INV_PHI_D;
    double resonance=exp(-(dr*dr)*(double)n*8.0);
    /* omega modulation — gentle, not monotone */
    double om_mod=1.0+0.1*cos(omega*LN_PHI_D*(double)n);
    return scale*resonance*om_mod;
}

static double get_dn(int D, double r, double omega) {
    if(hdgl_active&&hdgl_fn) return hdgl_fn(D,r,omega);
    return dn_inline(D,r,omega);
}

static bool hdgl_load(const char *path) {
    if(!path||!*path) return false;
    hdgl_handle=DL_OPEN(path);
    if(!hdgl_handle){ fprintf(stderr,"HDGL load '%s': %s\n",path,DL_ERR()); return false; }
    hdgl_fn=(double(*)(int,double,double))DL_SYM(hdgl_handle,"compute_Dn_r");
    if(!hdgl_fn){ fprintf(stderr,"compute_Dn_r not found in '%s'\n",path);
        DL_CLOSE(hdgl_handle); hdgl_handle=NULL; return false; }
    hdgl_active=true; return true;
}
static void hdgl_unload(void) {
    if(hdgl_handle){ DL_CLOSE(hdgl_handle); hdgl_handle=NULL; }
    hdgl_fn=NULL; hdgl_active=false;
}

/* Auto-calibrated Dn threshold per dimension */
static double dn_threshold(int D) {
    if(D>=6) return 0.5;
    if(D>=3) return 0.2;
    return 0.05;
}

/* ================================================================
 * BASE4096 OUTPUT
 * ================================================================ */
#define B4096_RADIX 4096u
#define B4096_START 0x1000u
static const char *B4096_ALPHA_PATH="base4096-2.0.1/frozen_base4096_alphabet.txt";
static char       *b4096_text=NULL;
static const char *b4096_sym[4096];
static bool        b4096_ok=false, b4096_fail=false;

static size_t utf8_len(uint8_t c){
    if(c<0x80)return 1; if((c&0xE0)==0xC0)return 2;
    if((c&0xF0)==0xE0)return 3; if((c&0xF8)==0xF0)return 4; return 1;
}
static size_t utf8_enc(uint32_t cp, char *out, size_t sz){
    if(cp<=0x7F&&sz>=1){out[0]=(char)cp;return 1;}
    if(cp<=0x7FF&&sz>=2){out[0]=(char)(0xC0|(cp>>6));out[1]=(char)(0x80|(cp&0x3F));return 2;}
    if(cp<=0xFFFF&&sz>=3){out[0]=(char)(0xE0|(cp>>12));out[1]=(char)(0x80|((cp>>6)&0x3F));out[2]=(char)(0x80|(cp&0x3F));return 3;}
    if(cp<=0x10FFFF&&sz>=4){out[0]=(char)(0xF0|(cp>>18));out[1]=(char)(0x80|((cp>>12)&0x3F));out[2]=(char)(0x80|((cp>>6)&0x3F));out[3]=(char)(0x80|(cp&0x3F));return 4;}
    return 0;
}
static bool b4096_ensure(void){
    if(b4096_ok) return true; if(b4096_fail) return false;
    FILE *f=fopen(B4096_ALPHA_PATH,"rb"); if(!f){b4096_fail=true;return false;}
    char tmp[32768]; size_t rb=fread(tmp,1,sizeof(tmp)-1,f); fclose(f);
    if(!rb||rb>=sizeof(tmp)-1){b4096_fail=true;return false;}
    tmp[rb]='\0';
    char clean[32768]; char *dst=clean;
    for(char *s=tmp;*s;s++) if(*s!='\r'&&*s!='\n') *dst++=*s;
    *dst='\0';
    size_t total=(size_t)(dst-clean);
    b4096_text=(char*)malloc(total+1); if(!b4096_text){b4096_fail=true;return false;}
    memcpy(b4096_text,clean,total+1);
    int cnt=0; char *cur=b4096_text, *end=b4096_text+total;
    while(cur<end&&cnt<4096){
        uint8_t first=(uint8_t)*cur; size_t ln=utf8_len(first);
        if(cur+ln>end){free(b4096_text);b4096_fail=true;return false;}
        b4096_sym[cnt++]=cur; cur+=ln;
    }
    if(cnt!=4096){free(b4096_text);b4096_fail=true;return false;}
    b4096_ok=true; return true;
}
static bool bu_to_b4096(const BigUint *a, char *buf, size_t bsz){
    bool canon=b4096_ensure();
    if(!bsz) return false;
    if(bu_is_zero(a)){
        if(canon){ const char*s=b4096_sym[0]; size_t l=utf8_len((uint8_t)*s);
            if(l>=bsz)return false; memcpy(buf,s,l); buf[l]='\0'; return true; }
        size_t l=utf8_enc(B4096_START,buf,bsz); buf[l]='\0'; return l>0;
    }
    BigUint tmp; bu_copy(&tmp,a);
    uint16_t digs[512]={0}; int cnt=0;
    while(!bu_is_zero(&tmp)){
        uint32_t r; BigUint nxt;
        bu_divmod_small(&nxt,&tmp,B4096_RADIX,&r);
        if(cnt>=(int)(sizeof(digs)/sizeof(digs[0]))) return false;
        digs[cnt++]=(uint16_t)r; bu_copy(&tmp,&nxt);
    }
    char *ptr=buf; size_t rem=bsz;
    for(int i=cnt-1;i>=0;i--){
        if(canon){ const char*s=b4096_sym[digs[i]]; size_t l=utf8_len((uint8_t)*s);
            if(l>=rem)return false; memcpy(ptr,s,l); ptr+=l; rem-=l; }
        else{ char u[4]; size_t l=utf8_enc(B4096_START+(uint32_t)digs[i],u,4);
            if(!l||l>=rem)return false; memcpy(ptr,u,l); ptr+=l; rem-=l; }
    }
    if(!rem) return false; *ptr='\0'; return true;
}
static void print_b4096(const char *dec, const char *label){
    size_t digs=0; for(const char*p=dec;*p;p++) if(*p>='0'&&*p<='9') digs++;
    if(!digs||digs>1230) return;
    BigUint v; bu_init(&v); bu_from_decimal(&v,dec);
    char buf[8192]={0};
    if(!bu_to_b4096(&v,buf,sizeof(buf))) return;
    if(label&&*label) printf("%s: %s\n",label,buf);
    else printf("Base4096: %s\n",buf);
}

/* ================================================================
 * PLN ← → PRIME: LOG-SPACE TOOLS
 * ================================================================ */
static void logcand_init(LogCandidate *lc, long double log2x){
    lc->log2_approx=log2x; lc->frac_correction=0.0L;
}
static BigUint bu_from_logcand(const LogCandidate *lc){
    BigUint out; bu_from_log2_approx(&out,lc->log2_approx+lc->frac_correction);
    return out;
}
static void refine_logcand(LogCandidate *lc, const PLNAddress *target, int passes){
    for(int p=0;p<passes;p++){
        BigUint approx=bu_from_logcand(lc);
        PLNAddress est=pln_from_biguint(&approx);
        long double delta=target->omega-est.omega;
        lc->frac_correction+=delta*0.5L*(LN_PHI/LN2);
        if(fabsl(delta)<1e-14L) break;
    }
}
static BigUint bu_from_omega(long double omega){
    long double log2x=(LN_PHI*expl(omega*LN_PHI))/LN2;
    BigUint out; bu_from_log2_approx(&out,log2x);
    if(bu_is_zero(&out)) bu_from_u64(&out,101ULL);
    if(bu_is_even(&out)) bu_add_small(&out,&out,1);
    return out;
}

/* ================================================================
 * FORWARD: integer → PLN
 * ================================================================ */
static void solve_int_to_pln(const char *s){
    printf("=== INTEGER → PLN ===\nValue: %s\n",s);
    print_b4096(s,"Base4096");
    BigUint v; bu_init(&v); bu_from_decimal(&v,s);
    PLNAddress addr=pln_from_biguint(&v);
    bool prime=miller_rabin_big(&v);
    printf("Primality: %s\n",prime?"PRIME":"composite");
    printf("PLN = [D_%d, %.18Lf*pi, %.15Lf, %.18Lf*phi]\n",
           addr.D,addr.theta/M_PI_L,addr.Q,addr.omega);
    double r=(double)orbit_r(addr.omega);
    double dn=get_dn(addr.D,r,hdgl_omega);
    printf("|R| = %.15f   Dn(r) = %.6f   gap_pred = %.1f\n",
           r,dn,predict_gap(expl((float)(addr.omega*LN_PHI)),r));
    printf("\n");
}

/* ================================================================
 * MAGNET STACK WRAPPER + REVERSE: PLN → prime
 * ================================================================ */
static bool magnet_pln(const BigUint *c, const PLNAddress *target, long double min_score){
    PLNAddress addr=pln_from_biguint(c);
    if(addr.D!=target->D) return false;
    return pln_resonance_score(&addr,target)>=min_score;
}

static bool find_prime_from_pln(int tD, long double tTheta_pi, long double tQ,
                                  long double tOmega, int max_steps,
                                  char *found, size_t fsz) {
    printf("=== PLN → PRIME ===\nTarget: [D_%d, %.12Lf*pi, %.12Lf, %.12Lf*phi]\n",
           tD,tTheta_pi,tQ,tOmega);
    PLNAddress target={tD,tTheta_pi*M_PI_L,tQ,tOmega};

    /* Log-space center with 5-pass refinement */
    LogCandidate lc; logcand_init(&lc,expl(tOmega*LN_PHI)*LN_PHI/LN2);
    refine_logcand(&lc,&target,5);
    BigUint center=bu_from_logcand(&lc);
    for(int ref=0;ref<5;ref++){
        PLNAddress est=pln_from_biguint(&center);
        long double delta=target.omega-est.omega;
        if(fabsl(delta)<1e-14L) break;
        lc.frac_correction+=delta*0.5L*(LN_PHI/LN2);
        center=bu_from_logcand(&lc);
    }

    /* D-alignment fix: if center lands in wrong D-band due to floating point
       boundary (e.g. log2x≈64.0 → reconstructs 2^64 instead of 2^64-n),
       nudge down by 1 until D matches target, up to 128 steps. */
    {
        PLNAddress ce=pln_from_biguint(&center);
        if(ce.D!=tD){
            BigUint probe; bu_copy(&probe,&center);
            for(int nudge=0;nudge<128;nudge++){
                BigUint tmp;
                if(bu_sub_u64(&tmp,&probe,1)){
                    bu_copy(&probe,&tmp);
                    PLNAddress pa=pln_from_biguint(&probe);
                    if(pa.D==tD){ bu_copy(&center,&probe); break; }
                } else break;
            }
            /* also try nudging up */
            ce=pln_from_biguint(&center);
            if(ce.D!=tD){
                bu_copy(&probe,&center);
                for(int nudge=0;nudge<128;nudge++){
                    BigUint tmp; bu_add_u64(&tmp,&probe,1);
                    bu_copy(&probe,&tmp);
                    PLNAddress pa=pln_from_biguint(&probe);
                    if(pa.D==tD){ bu_copy(&center,&probe); break; }
                }
            }
        }
        dbg("DEBUG center D=%d target D=%d\n",pln_from_biguint(&center).D,tD);
    }

    /* Orbit-minima guidance */
    double mins[32]; int nm=0;
    find_orbit_minima((double)(tOmega-1.0L),(double)(tOmega+1.0L),0.02,mins,&nm,32);
    if(nm>0) dbg("DEBUG %d orbit minima near omega=%.4Lf\n",nm,tOmega);

    BigUint best; bu_init(&best);
    long double best_score=-1.0L;
    long double best_r=2.0L;

    /* Sequential odd-integer walk: +1,-1,+3,-3,+5,-5,...
       Visits every odd integer around center in order of distance.
       Uses BigUint arithmetic throughout — no overflow on large numbers. */
    int effective_steps=bu_bit_length(&center)>2000?max_steps/8:max_steps;

    /* Ensure center is odd */
    if(bu_is_even(&center)) bu_add_small(&center,&center,1);

    for(int offset=0;offset<=effective_steps;offset++){
        BigUint candidate;
        if(offset==0){
            bu_copy(&candidate,&center);
        } else {
            /* offset 1 → +2, 2 → -2, 3 → +4, 4 → -4, ...
               i.e. step = offset*2, sign alternates.
               This visits center±2, center±4, center±6, ...
               All are odd since center is odd and steps are even. */
            uint64_t step=(uint64_t)offset*2ULL;
            if(offset%2==1){
                bu_add_u64(&candidate,&center,step);
            } else {
                if(!bu_sub_u64(&candidate,&center,step)) continue;
            }
        }

        /* === MAGNET STACK === */
        if(!magnet_trial(&candidate))        continue;
        if(!magnet_lucas(&candidate))        continue;
        if(!magnet_fermat2(&candidate))      continue;
        if(!magnet_euler_jacobi(&candidate)) continue;

        PLNAddress addr=pln_from_biguint(&candidate);
        if(addr.D!=tD) continue;
        long double score=pln_resonance_score(&addr,&target);
        /* Adaptive score gate: wider near log2 boundaries */
        long double score_gate=PLN_SCORE_BORDER*0.8L;
        if(score<score_gate) continue;

        double r=(double)orbit_r(addr.omega);
        double dn=get_dn(addr.D,r,hdgl_omega);
        /* Only apply Dn gate when score is mediocre */
        if(score<PLN_SCORE_PASS&&dn<dn_threshold(addr.D)) continue;

        if(score>best_score){
            best_score=score; best_r=(long double)r;
            bu_copy(&best,&candidate);
        }

        /* Full MR on any surviving candidate */
        if(score>=score_gate&&miller_rabin_big(&candidate)){
            bu_to_decimal(&candidate,found,fsz);
            printf(" Resonant prime found: %s\n",found);
            printf(" Score=%.6Lf |R|=%.6Lf Dn=%.4f\n",score,(long double)r,dn);
            print_b4096(found,"Base4096");
            return true;
        }
    }

    if(best_score>0.70L&&!bu_is_zero(&best)&&miller_rabin_big(&best)){
        bu_to_decimal(&best,found,fsz);
        printf(" Nearest resonant prime: %s\n",found);
        printf(" Score=%.6Lf |R|=%.6Lf\n",best_score,best_r);
        print_b4096(found,"Base4096");
        return true;
    }
    printf(" No prime found within %d steps. Best score: %.4Lf\n",max_steps,best_score);
    return false;
}

/* ================================================================
 * X(z) GENERATOR
 * ================================================================ */
static void generate_xz(long double n, long double ns, long double tg,
                         long double zr, long double zi,
                         int steps, char *found, size_t fsz){
    printf("=== X(z) RESONANT PRIME GENERATOR ===\n");
    printf("n=%.4Lf ns=%.4Lf gauge=%.4Lf z=(%.4Lf+%.4Lfi)\n",n,ns,tg,zr,zi);
    long double phi_base=powl(PHI,PHI-1.0L);
    long double inner=n+0.5L/PHI;
    long double lattice=powl(PHI,powl(PHI,inner));
    long double zmag=sqrtl((1.0L+zr)*(1.0L+zr)+zi*zi);
    long double power=powl(zmag,ns);
    long double x_approx=phi_base*sqrtl(lattice)*power;
    long double log2x=log2l(x_approx)+(tg/(2.0L*M_PI_L));
    LogCandidate lc; logcand_init(&lc,log2x);
    BigUint start=bu_from_logcand(&lc);
    PLNAddress addr=pln_from_biguint(&start);
    printf("~%.0Lf decimal digits, initial PLN omega=%.12Lf\n",
           log10l(x_approx), addr.omega);
    bool ok=find_prime_from_pln(addr.D,addr.theta/M_PI_L,addr.Q,addr.omega,steps,found,fsz);
    if(ok) printf("Resonant prime recovered from X(z) lattice.\n");
    else   printf("No strong resonance found.\n");
}

/* ================================================================
 * FILE-BASED ASTRONOMICAL PRIME TOOLS
 * ================================================================ */
static bool prime_file_to_log2(const char *path, long double *out_log2, size_t *out_digits){
    FILE *f=fopen(path,"r"); if(!f) return false;
    uint64_t lead=0; int lead_d=0; size_t total=0; int c;
    while((c=fgetc(f))!=EOF){
        if(c>='0'&&c<='9'){
            total++;
            if(lead_d<MAX_FILE_LEAD_DIGITS){ lead=lead*10+(uint64_t)(c-'0'); lead_d++; }
        }
    }
    fclose(f);
    if(!total||!lead_d) return false;
    long double log10x=(long double)(total-lead_d)+log10l((long double)lead);
    *out_log2=log10x*LOG2_10;
    if(out_digits) *out_digits=total;
    return true;
}

static bool prime_file_to_pln(const char *path, PLNAddress *addr,
                               size_t *digits, long double *log10x_out){
    long double log2x; size_t digs;
    if(!prime_file_to_log2(path,&log2x,&digs)) return false;
    long double log10x=log2x/LOG2_10;
    *addr=pln_from_log2(log2x);
    if(digits) *digits=digs;
    if(log10x_out) *log10x_out=log10x;
    return true;
}

static void solve_file_to_pln(const char *path){
    PLNAddress addr; size_t digs; long double log10x;
    if(!prime_file_to_pln(path,&addr,&digs,&log10x)){
        printf("Error reading prime file '%s'\n",path); return;
    }
    unsigned long long lead=(unsigned long long)powl(10.0L,fmodl(log10x,1.0L)+5);
    printf("=== FILE PRIME → PLN ===\nFile: %s\nDigits: %zu\nLeading≈%06llu\n",path,digs,lead);
    printf("PLN = [D_%d, %.18Lf*pi, %.15Lf, %.18Lf*phi]\n",
           addr.D,addr.theta/M_PI_L,addr.Q,addr.omega);
    double r=(double)orbit_r(addr.omega);
    printf("|R|=%.12f  Dn(r)=%.6f  gap_pred≈%.0f\n",
           r,get_dn(addr.D,r,hdgl_omega),
           predict_gap(expl((float)(addr.omega*LN_PHI)),r));
}

static long double pln_score_by_log10(const PLNAddress *target, long double log10x, int mod8){
    PLNAddress cand=pln_from_log2(log10x*LOG2_10);
    long double dt=fabsl(cand.theta-target->theta);
    if(dt>M_PI_L) dt=TWO_PI_L-dt;
    long double dQ=fabsl(cand.Q-target->Q);
    long double dw=fabsl(fmodl(cand.omega,1.0L)-fmodl(target->omega,1.0L));
    if(dw>0.5L) dw=1.0L-dw;
    long double score=0.5L/(1.0L+dt/M_PI_L)+0.3L/(1.0L+dQ)+0.2L/(1.0L+dw);
    if(cand.D!=target->D) score-=0.1L;
    if(mod8==7) score+=0.03L; else if(mod8==1) score-=0.02L;
    double r=(double)orbit_r(cand.omega);
    double dn=get_dn(cand.D,r,hdgl_omega);
    score+=(long double)(0.1*dn/(1.0+dn));
    if(score<0.0L) score=0.0L; if(score>1.0L) score=1.0L;
    return score;
}

static void generate_file_candidates(const char *path, int count,
                                      long double step, int lead_digits){
    if(lead_digits<1) lead_digits=6; if(lead_digits>18) lead_digits=18;
    PLNAddress base; size_t digs; long double log10x;
    if(!prime_file_to_pln(path,&base,&digs,&log10x)){
        printf("Error reading '%s'\n",path); return;
    }
    long double log2x_base __attribute__((unused))=log10x*LOG2_10;
    printf("=== PRIME FILE CANDIDATES ===\nAnchor: %s (%zu digits)\n",path,digs);
    printf("PLN: [D_%d, %.18Lf*pi, %.15Lf, %.18Lf*phi]\n\n",
           base.D,base.theta/M_PI_L,base.Q,base.omega);
    int half=count/2;
    for(int i=0;i<count;i++){
        int off=i-half;
        long double cand_omega=base.omega+(long double)off*step;
        if(cand_omega<=0.0L) continue;
        PLNAddress caddr=pln_from_log2(cand_omega*LN_PHI/LN2);
        long double clog10=cand_omega*LN_PHI/logl(10.0L);
        unsigned long long cdigs=(unsigned long long)floorl(clog10)+1ULL;
        long double lead_ld=powl(10.0L,fmodl(clog10,1.0L)+(long double)(lead_digits-1));
        unsigned long long clead=(unsigned long long)(lead_ld+0.5L);
        long double score=pln_score_by_log10(&base,clog10,7);
        double r=(double)orbit_r(caddr.omega);
        double dn=get_dn(caddr.D,r,hdgl_omega);
        printf("  Δω=%+.4Lf | D_%d | θ=%.10Lf*pi | digits≈%llu | lead≈%0*llu | Dn=%.5f | score=%.5Lf\n",
               (long double)off*step, caddr.D, caddr.theta/M_PI_L,
               cdigs, lead_digits, clead, dn, score);
    }
}

/* ================================================================
 * EXACT INTEGER NEIGHBORHOOD SCAN
 * ================================================================ */
static void test_candidate(const char *dec, uint64_t max_offset){
    size_t digs=0; for(const char*p=dec;*p;p++) if(*p>='0'&&*p<='9') digs++;
    if(!digs||digs>1230){ printf("Error: %zu digits (max 1230 for exact MR)\n",digs); return; }
    BigUint base; bu_init(&base); bu_from_decimal(&base,dec);
    printf("=== EXACT CANDIDATE TEST ===\nCandidate (%zu digits): %s\n",digs,dec);
    if(miller_rabin_big(&base)){ printf("Base is PRIME\n"); print_b4096(dec,"Base4096"); return; }
    printf("Composite. Searching ±%llu...\n",(unsigned long long)max_offset);
    bool base_even=bu_is_even(&base);
    BigUint cand, nbr;
    for(uint64_t off=1;off<=max_offset;off++){
        bool test=base_even?((off&1)!=0):((off&1)==0);
        if(!test) continue;
        bu_add_u64(&cand,&base,off);
        if(miller_rabin_big(&cand)){
            char s[2048]; bu_to_decimal(&cand,s,sizeof(s));
            printf("Prime at +%llu: %s\n",(unsigned long long)off,s); print_b4096(s,"Base4096"); return;
        }
        if(bu_sub_u64(&nbr,&base,off)&&miller_rabin_big(&nbr)){
            char s[2048]; bu_to_decimal(&nbr,s,sizeof(s));
            printf("Prime at -%llu: %s\n",(unsigned long long)off,s); print_b4096(s,"Base4096"); return;
        }
    }
    printf("No prime found in ±%llu.\n",(unsigned long long)max_offset);
}

static void scan_range_cmd(const char *dec, uint64_t max_offset, const char *anchor_file){
    size_t digs=0; for(const char*p=dec;*p;p++) if(*p>='0'&&*p<='9') digs++;
    if(!digs||digs>1230){ printf("Error: %zu digits (max 1230)\n",digs); return; }
    PLNAddress anchor; bool have_anchor=false;
    if(anchor_file&&*anchor_file){
        size_t ad; long double al;
        if(prime_file_to_pln(anchor_file,&anchor,&ad,&al)) have_anchor=true;
    }
    BigUint base; bu_init(&base); bu_from_decimal(&base,dec);
    PLNAddress base_addr=pln_from_biguint(&base);
    PLNAddress *ref=have_anchor?&anchor:&base_addr;
    printf("=== SCAN RANGE ===\nBase (%zu digits): %s\nScan: ±%llu\n",
           digs,dec,(unsigned long long)max_offset);
    printf("%-7s | %-6s | mod8 | D | score     | |R|       | Dn\n","offset","prime");
    printf("--------+--------+------+---+-----------+-----------+----------\n");
    BigUint cand;
    for(long long delta=-(long long)max_offset;delta<=(long long)max_offset;delta++){
        bool valid=true;
        if(delta==0) bu_copy(&cand,&base);
        else if(delta>0) bu_add_u64(&cand,&base,(uint64_t)delta);
        else valid=bu_sub_u64(&cand,&base,(uint64_t)(-delta));
        if(!valid) continue;
        PLNAddress addr=pln_from_biguint(&cand);
        int mod8=(int)(cand.d[0]&7u);
        bool prime=!bu_is_even(&cand)&&miller_rabin_big(&cand);
        long double score=pln_resonance_score(&addr,ref);
        double r=(double)orbit_r(addr.omega);
        double dn=get_dn(addr.D,r,hdgl_omega);
        char off_s[16];
        snprintf(off_s,sizeof(off_s),delta>=0?"+%lld":"%lld",delta);
        printf("%-7s | %-6s | %4d | %d | %.9Lf | %.9f | %.4f\n",
               off_s,prime?"PRIME":"",mod8,addr.D,score,r,dn);
    }
}

/* ================================================================
 * USAGE / HELP
 * ================================================================ */
static void print_usage(void){
    printf("Usage:\n");
    printf("  ./pln                                     built-in demo\n");
    printf("  ./pln --prime <N>                         integer → PLN\n");
    printf("  ./pln --pln D THETA Q OMEGA [STEPS]       PLN → prime\n");
    printf("  ./pln --roundtrip <N> [STEPS]             bidirectional verify\n");
    printf("  ./pln --xz n ns tg zr zi [STEPS]         X(z) generator\n");
    printf("  ./pln --prime-file <FILE>                 file prime → PLN\n");
    printf("  ./pln --prime-file-candidates <F> [N] [STEP] [LEAD] [SO]\n");
    printf("  ./pln --test-candidate <N> [OFFSET]       MR scan ±OFFSET\n");
    printf("  ./pln --scan-range <N> [OFFSET] [ANCHOR]  score neighborhood\n");
    printf("  ./pln --base4096 <N>                      encode to Base4096\n");
    printf("  ./pln --debug --prime <N>\n");
    printf("  ./pln --help\n\n");
    printf("APA mode: %s\n",
           SUPER_RESOLUTION_MODE?"SUPER-RES 16,777,216-bit":"FAST 512-bit");
    printf("HDGL Dn(r): %s\n",hdgl_active?"external .so":"inline resonance-peak formula");
}

/* ================================================================
 * MAIN
 * ================================================================ */
int main(int argc, char *argv[]){
    printf("RESONANT PLN PRIME SOLVER — FINAL UNIFIED EDITION\n");
    printf("4096-bit BigUint | %s APA | Magnet Stack | Dn(r) Resonance\n\n",
           SUPER_RESOLUTION_MODE?"16M-bit":"512-bit");

    if(argc==1){
        solve_int_to_pln("18446744073709551557");
        solve_int_to_pln("100000000000000000000000000000000000000000000000151");
        char found[2048]={0};
        find_prime_from_pln(2,0.40088677849L,1.4165141186L,9.40088677849L,20000,found,sizeof(found));
        printf("\nUse --help for full options.\n");
        return 0;
    }

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

    if(strcmp(argv[1],"--help")==0||strcmp(argv[1],"-h")==0){
        print_usage(); return 0;
    }
    if(strcmp(argv[1],"--base4096")==0&&argc>=3){
        BigUint v; bu_init(&v); bu_from_decimal(&v,argv[2]);
        char buf[8192]={0};
        if(bu_to_b4096(&v,buf,sizeof(buf))) printf("%s\n",buf);
        return 0;
    }
    if(strcmp(argv[1],"--prime")==0&&argc>=3){
        solve_int_to_pln(argv[2]); return 0;
    }
    if(strcmp(argv[1],"--prime-file")==0&&argc>=3){
        solve_file_to_pln(argv[2]); return 0;
    }
    if(strcmp(argv[1],"--prime-file-candidates")==0&&argc>=3){
        int cnt=(argc>3)?atoi(argv[3]):10;
        long double step=(argc>4)?strtold(argv[4],NULL):0.005L;
        int lead=(argc>5)?atoi(argv[5]):12;
        const char *so=(argc>6)?argv[6]:NULL;
        if(so&&*so) hdgl_load(so);
        generate_file_candidates(argv[2],cnt,step,lead);
        hdgl_unload(); return 0;
    }
    if(strcmp(argv[1],"--pln")==0&&argc>=6){
        int D=atoi(argv[2]);
        long double th=strtold(argv[3],NULL),Q=strtold(argv[4],NULL),om=strtold(argv[5],NULL);
        int steps=(argc>6)?atoi(argv[6]):20000;
        char found[2048]={0};
        find_prime_from_pln(D,th,Q,om,steps,found,sizeof(found));
        return 0;
    }
    if(strcmp(argv[1],"--roundtrip")==0&&argc>=3){
        printf("=== ROUNDTRIP ===\nInput: %s\n",argv[2]);
        BigUint v; bu_init(&v); bu_from_decimal(&v,argv[2]);
        bool prime=miller_rabin_big(&v);
        printf("Primality: %s\n",prime?"PRIME":"composite");
        if(prime){
            PLNAddress addr=pln_from_biguint(&v);
            printf("PLN = [D_%d, %.18Lf*pi, %.15Lf, %.18Lf*phi]\n",
                   addr.D,addr.theta/M_PI_L,addr.Q,addr.omega);
            int steps=(argc>3)?atoi(argv[3]):20000;
            char found[2048]={0};
            bool ok=find_prime_from_pln(addr.D,addr.theta/M_PI_L,addr.Q,addr.omega,steps,found,sizeof(found));
            if(ok){
                printf("CLOSURE: %s\n",strcmp(argv[2],found)==0?"EXACT MATCH":"nearest resonant prime");
            }
        }
        return 0;
    }
    if(strcmp(argv[1],"--xz")==0&&argc>=7){
        long double n=strtold(argv[2],NULL),ns=strtold(argv[3],NULL);
        long double tg=strtold(argv[4],NULL),zr=strtold(argv[5],NULL),zi=strtold(argv[6],NULL);
        int steps=(argc>7)?atoi(argv[7]):30000;
        char found[2048]={0};
        generate_xz(n,ns,tg,zr,zi,steps,found,sizeof(found));
        return 0;
    }
    if(strcmp(argv[1],"--test-candidate")==0&&argc>=3){
        uint64_t off=(argc>3)?strtoull(argv[3],NULL,10):256ULL;
        test_candidate(argv[2],off); return 0;
    }
    if(strcmp(argv[1],"--scan-range")==0&&argc>=3){
        uint64_t off=(argc>3)?strtoull(argv[3],NULL,10):32ULL;
        const char *anchor=(argc>4)?argv[4]:NULL;
        scan_range_cmd(argv[2],off,anchor); return 0;
    }

    print_usage(); return 1;
}