/*
 * pln_analog_hdgl_final_v12.c
 * PLN + HDGL v3.0 — Final Mersenne Resonance Explorer
 * Clean, stable, and ready for continued hunting
 */

#ifndef _CRT_SECURE_NO_WARNINGS
#define _CRT_SECURE_NO_WARNINGS
#endif

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <math.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)
#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)
#endif

#define PHI_L        1.6180339887498948482L
#define LN_PHI_L     0.4812118250596034748L
#define LN2_L        0.6931471805599453094L
#define LOG2_10_L    3.32192809488736234787L
#define M_PI_L       3.14159265358979323846L
#define TWO_PI_L     6.2831853071795864769L

typedef struct { int D; long double theta; long double Q; long double omega; } PLNAddress;

static double (*hdgl_Dn)(int, double, double) = NULL;
static DLHandle hdgl_h = NULL;

static 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_L*PHI_L*(1.0L - powl(PHI_L, -D_cont)) + 1.0L/D_cont;
    lat = fmaxl(1.0L, fminl(lat, PHI_L*PHI_L));
    long double Rre = cosl(theta) - 1.0L + lat;
    long double Rim = sinl(theta);
    return (double)sqrtl(Rre*Rre + Rim*Rim);
}

static PLNAddress pln_from_omega(long double 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;
    double r = orbit_r(omega);
    long double Q = 1.0L / (r + 1e-20L);
    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;
    long double dw = fabsl(fmodl(a->omega, 1.0L) - fmodl(b->omega, 1.0L));
    if (dw > 0.5L) dw = 1.0L - dw;
    long double dD_penalty = (a->D == b->D) ? 0.0L : 0.25L;   // penalize D-band jumps
    return 0.5L/(1.0L + dt/M_PI_L) + 0.3L/(1.0L + fabsl(a->Q - b->Q))
         + 0.2L/(1.0L + dw) - dD_penalty;
}

static uint64_t mersenne_p_from_omega(long double omega) {
    if (omega <= 0.0L) return 0;
    long double log2x = (LN_PHI_L * expl(omega * LN_PHI_L)) / LN2_L;
    return (uint64_t)llroundl(log2x);
}

static long double omega_from_log2(long double log2p) {
    long double arg = log2p * LN2_L / LN_PHI_L;
    if (arg <= 1.0L) return 0.0L;
    return logl(arg) / LN_PHI_L;
}

static bool get_prime_log_info(const char *path, size_t *digits, long double *log10x, char *lead_buf, size_t lead_len) {
    FILE *f = fopen(path, "r");
    if (!f) return false;
    *digits = 0;
    char lead[64] = {0};
    size_t pos = 0;
    int c;
    while ((c = fgetc(f)) != EOF) {
        if (c >= '0' && c <= '9') {
            (*digits)++;
            if (pos < lead_len - 1) lead[pos++] = (char)c;
        }
    }
    fclose(f);
    if (*digits == 0 || pos == 0) return false;
    lead[pos] = '\0';
    strncpy(lead_buf, lead, lead_len);
    long double lead_val = strtold(lead, NULL);
    *log10x = (long double)(*digits - pos) + log10l(lead_val);
    return true;
}

static double default_hdgl_dn(int D, double r, double omega) {
    return sqrt(PHI_L * (D + 1.0) * pow(2.0, D) * (D + 2.0) * omega) * pow(r, (D + 1.0) / 8.0);
}

static void mersenne_tuner(const char *target, long double radius, int steps, const char *so_path) {
    if (so_path && *so_path) {
        hdgl_h = DL_OPEN(so_path);
        if (hdgl_h) hdgl_Dn = (double(*)(int,double,double))DL_SYM(hdgl_h, "compute_Dn_r");
    }

    size_t digits = 0;
    long double log10_anchor = 0.0L;
    char lead[32] = {0};
    long double anchor_omega = 0.0L;
    uint64_t known_p = 0;

    bool is_file = get_prime_log_info(target, &digits, &log10_anchor, lead, sizeof(lead));

    if (is_file) {
        long double log2_approx = log10_anchor * LOG2_10_L;
        known_p = (uint64_t)llroundl(log2_approx);
        anchor_omega = omega_from_log2(log2_approx);
        printf("Anchor: %s (%zu digits)  p=%llu\n", target, digits, (unsigned long long)known_p);
    } else {
        anchor_omega = strtold(target, NULL);
        known_p = mersenne_p_from_omega(anchor_omega);
    }

    PLNAddress anchor = pln_from_omega(anchor_omega);
    printf("Anchor PLN = [D_%d  θ=%.12Lfπ  ω=%.12Lf]\n", anchor.D, anchor.theta/M_PI_L, anchor.omega);

    if (radius <= 0.0L) radius = 5.0L;
    if (steps < 1000) steps = 80000;

    printf("Scanning ±%.6Lf ω ...\n\n", radius);

    struct Cand { long double omega; long double score; double dn; uint64_t p; } top[15] = {0};
    int top_cnt = 0;
    time_t start = time(NULL);

    // Force the exact anchor into the list
    {
        double r = orbit_r(anchor_omega);
        double dn = hdgl_Dn ? hdgl_Dn(anchor.D, r, (double)anchor_omega) : default_hdgl_dn(anchor.D, r, (double)anchor_omega);
        long double score = pln_resonance_score(&anchor, &anchor) + 0.45L * ((long double)dn / (1.0L + dn));
        top[top_cnt++] = (struct Cand){anchor_omega, score, dn, known_p};
    }

    for (int i = 0; i < steps; i++) {
        if (i % 2000 == 0 && i > 0) {
            double pct = 100.0 * i / steps;
            double elapsed = difftime(time(NULL), start);
            double eta = (i > 0) ? elapsed * (steps - i) / (double)i : 0;
            printf("\rProgress: %5.1f%%  ETA ~ %.0fs", pct, eta);
            fflush(stdout);
        }

        long double t = (long double)i / (steps - 1);
        long double omega = anchor_omega - radius + 2.0L * radius * t;
        if (omega <= 0.0L) continue;

        PLNAddress cand = pln_from_omega(omega);
        if (cand.D != anchor.D) continue;   // Strong D-band filter

        double r = orbit_r(omega);
        double dn = hdgl_Dn ? hdgl_Dn(cand.D, r, (double)omega) : default_hdgl_dn(cand.D, r, (double)omega);
        long double score = pln_resonance_score(&cand, &anchor) + 0.45L * ((long double)dn / (1.0L + dn));

        uint64_t p = mersenne_p_from_omega(omega);

        if (top_cnt < 15 || score > top[14].score) {
            if (top_cnt == 15) top_cnt--;
            top[top_cnt++] = (struct Cand){omega, score, dn, p};
            for (int a = 0; a < top_cnt - 1; a++) {
                for (int b = a + 1; b < top_cnt; b++) {
                    if (top[a].score < top[b].score) {
                        struct Cand tmp = top[a]; top[a] = top[b]; top[b] = tmp;
                    }
                }
            }
        }
    }

    printf("\n\n=== TOP EXACT MERSENNE CANDIDATES ===\n");

    FILE *f = fopen("mersenne_candidates.txt", "w");
    if (f) {
        fprintf(f, "# PLN+HDGL v12 Candidates\n");
        fprintf(f, "# Anchor p = %llu\n\n", (unsigned long long)known_p);

        for (int i = 0; i < (top_cnt < 12 ? top_cnt : 12); i++) {
            uint64_t p = top[i].p;
            printf("%2d  ω=%.8Lf  p=%llu  Dn=%.5f  score=%.6Lf\n",
                   i+1, top[i].omega, (unsigned long long)p, top[i].dn, top[i].score);
            fprintf(f, "p = %llu    # ω=%.8Lf  score=%.6Lf  Dn=%.5f\n",
                    (unsigned long long)p, top[i].omega, top[i].score, top[i].dn);
            fprintf(f, "   gpuowl -p %llu\n\n", (unsigned long long)p);
        }
        fclose(f);
        printf("\nSaved to mersenne_candidates.txt\n");
    }

    if (hdgl_h) DL_CLOSE(hdgl_h);
}

int main(int argc, char *argv[]) {
    printf("=== PLN + HDGL v3.0 — Exact Mersenne Explorer (v12) ===\n\n");

    if (argc >= 3 && strcmp(argv[1], "--mersenne-tune") == 0) {
        long double r = (argc > 3) ? strtold(argv[3], NULL) : 5.0L;
        int s = (argc > 4) ? atoi(argv[4]) : 80000;
        const char *so = (argc > 5) ? argv[5] : "hdgl_analog_v30.so";
        mersenne_tuner(argv[2], r, s, so);
        return 0;
    }

    mersenne_tuner("M136279841.txt", 5.0L, 80000, "hdgl_analog_v30.so");
    return 0;
}