/*
 * mersenne_gen.c
 * ============================================================
 * MERSENNE PRIME FULL INTEGER GENERATOR
 *
 * Generates 2^p - 1 as:
 *   • Binary file  (little-endian, exact, always fast)
 *   • Decimal text (requires GMP, ~1 min for p=136279879)
 *   • Hex text     (no GMP, instant)
 *   • PLN info     (omega, D, theta, Q, |R|)
 *
 * Compile (with GMP — recommended):
 *   gcc -O3 -std=c99 -o mersenne_gen mersenne_gen.c -lgmp -lm
 *
 * Compile (without GMP — binary+hex only):
 *   gcc -O3 -std=c99 -DNO_GMP -o mersenne_gen mersenne_gen.c -lm
 *
 * Usage:
 *   ./mersenne_gen --info   136279841
 *   ./mersenne_gen --binary 136279841 [output.bin]
 *   ./mersenne_gen --hex    136279841 [output.hex]
 *   ./mersenne_gen --decimal 136279841 [output.txt]   (needs GMP)
 *   ./mersenne_gen --next-candidates 136279841 [radius] [steps]
 *   ./mersenne_gen --verify-pln 136279879 39.681394074241442733
 *
 * The binary format is big-endian (most significant byte first),
 * matching standard cryptographic conventions.
 * For p=136279841 the binary is exactly 17,034,981 bytes.
 *
 * Notes:
 *   • --decimal for p~136M takes ~60s with GMP on modern hardware
 *   • --binary is instant (just writes 0xFF bytes + top byte)
 *   • --hex is instant (same logic, hexadecimal output)
 *   • Leading digits shown in --info match file-based PLN anchor
 */

#ifndef _GNU_SOURCE
#define _GNU_SOURCE
#endif

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <errno.h>

#ifndef NO_GMP
#include <gmp.h>
#endif

/* ================================================================
 * CONSTANTS
 * ================================================================ */
#define PHI         1.6180339887498948482L
#define LN_PHI      0.4812118250596034748L
#define LN2         0.6931471805599453094L
#define LOG10_2     0.30102999566398119521L
#define M_PI_L      3.14159265358979323846L
#define TWO_PI_L    6.2831853071795864769L
#define INV_PHI     0.6180339887498948482L

/* ================================================================
 * PLN COORDINATE — log-space, no big integer needed
 * ================================================================ */
typedef struct {
    uint64_t    p;          /* Mersenne exponent */
    long double omega;      /* PLN omega (phi-lattice coordinate) */
    int         D;          /* Dimension class 1..8 */
    long double theta_pi;   /* theta / pi */
    long double Q;          /* resonance quality */
    long double r;          /* orbit radius |R| */
    long double log10x;     /* log10(2^p - 1) ≈ p * log10(2) */
    uint64_t    digits;     /* decimal digit count */
} MersennePLN;

static long double orbit_r_ld(long double omega) {
    long double frac = omega - floorl(omega);
    long double theta = frac * TWO_PI_L;
    long double D_cont = fmodl(omega, 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 MersennePLN compute_pln(uint64_t p) {
    MersennePLN m;
    m.p = p;
    /* log2(2^p - 1) ≈ p for large p */
    long double log2x = (long double)p;
    long double lnx   = log2x * LN2;
    long double arg   = lnx / LN_PHI;
    m.omega     = (arg > 1.0L) ? logl(arg) / LN_PHI : 0.0L;
    m.D         = (int)fmodl(floorl(m.omega), 8.0L) + 1;
    m.theta_pi  = fmodl(m.omega, 1.0L);
    m.r         = orbit_r_ld(m.omega);
    m.Q         = 1.0L / (m.r + 1e-30L);
    m.log10x    = (long double)p * LOG10_2;
    m.digits    = (uint64_t)floorl(m.log10x) + 1;
    return m;
}

static void print_pln(const MersennePLN *m) {
    printf("  p                 = %llu\n", (unsigned long long)m->p);
    printf("  Decimal digits    = %llu\n", (unsigned long long)m->digits);
    printf("  log10(2^p-1)      = %.15Lf\n", m->log10x);
    printf("  omega             = %.18Lf\n", m->omega);
    printf("  D                 = %d  (D_cont = %.15Lf)\n",
           m->D, fmodl(m->omega, 8.0L) + 1.0L);
    printf("  theta / pi        = %.18Lf\n", m->theta_pi);
    printf("  theta (rad)       = %.15Lf\n", m->theta_pi * M_PI_L);
    printf("  Q                 = %.15Lf\n", m->Q);
    printf("  |R|               = %.15Lf\n", m->r);

    /* Leading 18 decimal digits via log10 fractional part */
    long double frac18 = fmodl(m->log10x, 1.0L) + 17.0L;
    unsigned long long lead18 = (unsigned long long)powl(10.0L, frac18);
    printf("  Leading 18 digits = %018llu\n", lead18);
}

/* ================================================================
 * BINARY OUTPUT — 2^p-1 in big-endian bytes
 * ================================================================ */
static int write_binary(uint64_t p, const char *path) {
    FILE *f = path ? fopen(path, "wb") : stdout;
    if (!f) { perror(path); return 1; }

    /* 2^p - 1 in binary:
       - top byte has (p % 8) ones in its low bits, or 0xFF if p%8==0
       - all remaining bytes are 0xFF
       Total bytes: ceil(p/8) */
    uint64_t total_bytes = (p + 7) / 8;
    int top_bits = (int)(p % 8);
    uint8_t top_byte = top_bits ? (uint8_t)((1u << top_bits) - 1u) : 0xFFu;

    /* Write top byte first (big-endian) */
    fwrite(&top_byte, 1, 1, f);

    /* Write remaining 0xFF bytes in chunks */
    uint8_t buf[65536];
    memset(buf, 0xFF, sizeof(buf));
    uint64_t remaining = total_bytes - 1;
    while (remaining > 0) {
        size_t chunk = remaining < sizeof(buf) ? (size_t)remaining : sizeof(buf);
        fwrite(buf, 1, chunk, f);
        remaining -= chunk;
    }

    if (path) {
        fclose(f);
        printf("Binary written: %s (%llu bytes = %.2f MB)\n",
               path, (unsigned long long)total_bytes,
               (double)total_bytes / 1048576.0);
    }
    return 0;
}

/* ================================================================
 * HEX OUTPUT — instant, no GMP
 * ================================================================ */
static int write_hex(uint64_t p, const char *path) {
    FILE *f = path ? fopen(path, "w") : stdout;
    if (!f) { perror(path); return 1; }

    /* 2^p - 1: top hex digits, then all Fs */
    uint64_t total_hex_digits = (p + 3) / 4;
    int top_bits = (int)(p % 4);
    char top_nibble = top_bits ? "0123456789ABCDEF"[(1u << top_bits) - 1u] : 'F';

    fprintf(f, "%c", top_nibble);
    uint64_t remaining = total_hex_digits - 1;
    uint8_t buf[80];
    memset(buf, 'F', 80);
    while (remaining > 0) {
        size_t chunk = remaining < 80 ? (size_t)remaining : 80;
        fwrite(buf, 1, chunk, f);
        remaining -= chunk;
    }
    fprintf(f, "\n");

    if (path) {
        fclose(f);
        printf("Hex written: %s (%llu hex digits)\n",
               path, (unsigned long long)total_hex_digits);
    }
    return 0;
}

/* ================================================================
 * DECIMAL OUTPUT — requires GMP
 * ================================================================ */
#ifndef NO_GMP
static int write_decimal(uint64_t p, const char *path) {
    FILE *f = path ? fopen(path, "w") : stdout;
    if (!f) { perror(path); return 1; }

    printf("Computing 2^%llu - 1 in decimal (GMP)...\n",
           (unsigned long long)p);
    printf("Expected: ~%llu digits, ~%.0f MB output\n",
           (unsigned long long)((uint64_t)(p * 0.30103) + 1),
           p * 0.30103 / 1048576.0);
    fflush(stdout);

    struct timespec t0, t1;
    clock_gettime(CLOCK_MONOTONIC, &t0);

    mpz_t x;
    mpz_init(x);
    mpz_ui_pow_ui(x, 2, (unsigned long)p);
    mpz_sub_ui(x, x, 1);

    clock_gettime(CLOCK_MONOTONIC, &t1);
    double compute_s = (t1.tv_sec - t0.tv_sec) + (t1.tv_nsec - t0.tv_nsec) * 1e-9;
    printf("Computed in %.2f s. Writing...\n", compute_s);
    fflush(stdout);

    clock_gettime(CLOCK_MONOTONIC, &t0);
    gmp_fprintf(f, "%Zd\n", x);
    clock_gettime(CLOCK_MONOTONIC, &t1);
    double write_s = (t1.tv_sec - t0.tv_sec) + (t1.tv_nsec - t0.tv_nsec) * 1e-9;

    size_t actual_digits = mpz_sizeinbase(x, 10);
    mpz_clear(x);

    if (path) {
        fclose(f);
        printf("Decimal written: %s (%zu digits) in %.2f s\n",
               path, actual_digits, write_s);
    }
    return 0;
}
#else
static int write_decimal(uint64_t p, const char *path) {
    (void)p; (void)path;
    fprintf(stderr,
        "Error: decimal output requires GMP.\n"
        "Recompile without -DNO_GMP:\n"
        "  gcc -O3 -std=c99 -o mersenne_gen mersenne_gen.c -lgmp -lm\n");
    return 1;
}
#endif

/* ================================================================
 * NEXT CANDIDATE SWEEP — PLN-based, no GMP needed
 * ================================================================ */
static void next_candidates(uint64_t anchor_p, long double radius, int steps) {
    MersennePLN anchor = compute_pln(anchor_p);
    printf("=== NEXT MERSENNE CANDIDATE SWEEP ===\n");
    printf("Anchor: p=%llu  omega=%.12Lf  D=%d\n\n",
           (unsigned long long)anchor_p, anchor.omega, anchor.D);

    /* Score function: PLN resonance similarity to anchor */
    typedef struct { uint64_t p; long double score; long double omega; int D; } Cand;
    #define MAX_CANDS 20
    Cand top[MAX_CANDS];
    memset(top, 0, sizeof(top));
    int tcnt = 0;

    for (int i = 0; i < steps; i++) {
        long double t = (long double)i / (steps - 1);
        long double omega = anchor.omega - radius + 2.0L * radius * t;
        if (omega <= 0.0L) continue;

        /* This omega → p value */
        long double log2x = expl(omega * LN_PHI) * LN_PHI / LN2;
        uint64_t p = (uint64_t)llroundl(log2x);
        if (p <= anchor_p) continue; /* only look forward */

        /* PLN at this omega */
        MersennePLN cand = compute_pln(p);

        /* Resonance score vs anchor */
        long double dth = fabsl(cand.theta_pi - anchor.theta_pi);
        if (dth > 0.5L) dth = 1.0L - dth;
        long double dQ  = fabsl(cand.Q - anchor.Q);
        long double dw  = fabsl(fmodl(cand.omega,1.0L) - fmodl(anchor.omega,1.0L));
        if (dw > 0.5L) dw = 1.0L - dw;
        long double score = 0.5L/(1.0L + dth)
                          + 0.3L/(1.0L + dQ)
                          + 0.2L/(1.0L + dw);
        /* Bonus for same D class */
        if (cand.D == anchor.D) score += 0.05L;

        /* Deduplicate and maintain top-N */
        bool dup = false;
        for (int k = 0; k < tcnt; k++) if (top[k].p == p) { dup = true; break; }
        if (!dup && (tcnt < MAX_CANDS || score > top[MAX_CANDS-1].score)) {
            if (tcnt == MAX_CANDS) tcnt--;
            top[tcnt++] = (Cand){p, score, cand.omega, cand.D};
            /* insertion sort descending */
            for (int a = tcnt-1; a > 0 && top[a].score > top[a-1].score; a--) {
                Cand tmp = top[a]; top[a] = top[a-1]; top[a-1] = tmp;
            }
        }
    }

    printf("  Rank  p                  omega          D   score    digits\n");
    printf("  ----  -----------------  -------------  --  -------  ------\n");
    int show = tcnt < 15 ? tcnt : 15;
    for (int i = 0; i < show; i++) {
        MersennePLN m = compute_pln(top[i].p);
        printf("  %4d  %-17llu  %.12Lf  D%d  %.5Lf  %llu\n",
               i+1, (unsigned long long)top[i].p,
               top[i].omega, top[i].D, top[i].score,
               (unsigned long long)m.digits);
    }

    /* Write candidates file */
    FILE *f = fopen("next_mersenne_candidates.txt", "w");
    if (f) {
        fprintf(f, "# PLN Next-Mersenne Candidates\n");
        fprintf(f, "# Anchor: p=%llu  omega=%.12Lf\n\n",
                (unsigned long long)anchor_p, anchor.omega);
        for (int i = 0; i < show; i++) {
            MersennePLN m = compute_pln(top[i].p);
            fprintf(f, "p = %-17llu # omega=%.10Lf  score=%.5Lf  D%d  digits=%llu\n",
                    (unsigned long long)top[i].p, top[i].omega,
                    top[i].score, top[i].D, (unsigned long long)m.digits);
            fprintf(f, "   gpuowl -p %llu\n\n", (unsigned long long)top[i].p);
        }
        fclose(f);
        printf("\nSaved: next_mersenne_candidates.txt\n");
    }
    #undef MAX_CANDS
}

/* ================================================================
 * VERIFY PLN — check claimed values against our formula
 * ================================================================ */
static void verify_pln(uint64_t p, long double claimed_omega) {
    MersennePLN m = compute_pln(p);
    printf("=== PLN VERIFICATION: p=%llu ===\n\n", (unsigned long long)p);
    printf("Our computed values:\n");
    print_pln(&m);
    printf("\nClaimed omega = %.18Lf\n", claimed_omega);
    printf("Our omega     = %.18Lf\n", m.omega);
    printf("Difference    = %.6Le\n", m.omega - claimed_omega);
    long double rel = fabsl(m.omega - claimed_omega) / fabsl(claimed_omega + 1e-30L);
    printf("Relative diff = %.6Le  (%s)\n", rel,
           rel < 1e-10L ? "MATCH (within double precision)" :
           rel < 1e-6L  ? "CLOSE (formula variant)" : "MISMATCH");
    printf("\nD match: %s  theta/pi match: %s\n",
           m.D == (int)fmodl(floorl(claimed_omega), 8.0L)+1 ? "YES" : "NO",
           fabsl(m.theta_pi - fmodl(claimed_omega, 1.0L)) < 1e-10L ? "YES" : "NO");
}

/* ================================================================
 * MAIN
 * ================================================================ */
static void print_usage(void) {
    printf("Usage:\n");
    printf("  mersenne_gen --info    <p>\n");
    printf("  mersenne_gen --binary  <p> [output.bin]\n");
    printf("  mersenne_gen --hex     <p> [output.hex]\n");
    printf("  mersenne_gen --decimal <p> [output.txt]    (needs GMP)\n");
    printf("  mersenne_gen --next-candidates <p> [radius] [steps]\n");
    printf("  mersenne_gen --verify-pln <p> <claimed_omega>\n");
    printf("  mersenne_gen --help\n\n");
    printf("Examples:\n");
    printf("  mersenne_gen --info 136279841\n");
    printf("  mersenne_gen --binary 136279841 M136279841.bin\n");
    printf("  mersenne_gen --decimal 136279841 M136279841.txt\n");
    printf("  mersenne_gen --next-candidates 136279841 8.0 500000\n");
    printf("  mersenne_gen --verify-pln 136279879 39.681394074241442733\n\n");
#ifdef NO_GMP
    printf("GMP: NOT compiled in (--decimal unavailable)\n");
#else
    printf("GMP: compiled in (--decimal available)\n");
#endif
}

int main(int argc, char *argv[]) {
    printf("MERSENNE PRIME GENERATOR — PLN Edition\n");
#ifdef NO_GMP
    printf("Mode: binary+hex+PLN only (no GMP)\n\n");
#else
    printf("Mode: full (binary + hex + decimal + PLN, GMP enabled)\n\n");
#endif

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

    if (strcmp(argv[1], "--info") == 0 && argc >= 3) {
        uint64_t p = strtoull(argv[2], NULL, 10);
        MersennePLN m = compute_pln(p);
        printf("=== 2^%llu - 1 ===\n\n", (unsigned long long)p);
        print_pln(&m);
        return 0;
    }

    if (strcmp(argv[1], "--binary") == 0 && argc >= 3) {
        uint64_t p = strtoull(argv[2], NULL, 10);
        const char *path = argc > 3 ? argv[3] : NULL;
        if (!path) {
            /* auto-name */
            static char auto_path[64];
            snprintf(auto_path, sizeof(auto_path), "M%llu.bin",
                     (unsigned long long)p);
            path = auto_path;
        }
        MersennePLN m = compute_pln(p);
        printf("Generating 2^%llu - 1\n", (unsigned long long)p);
        printf("  Digits  : %llu\n", (unsigned long long)m.digits);
        printf("  Bytes   : %llu\n", (unsigned long long)((p + 7) / 8));
        return write_binary(p, path);
    }

    if (strcmp(argv[1], "--hex") == 0 && argc >= 3) {
        uint64_t p = strtoull(argv[2], NULL, 10);
        const char *path = argc > 3 ? argv[3] : NULL;
        if (!path) {
            static char auto_path[64];
            snprintf(auto_path, sizeof(auto_path), "M%llu.hex",
                     (unsigned long long)p);
            path = auto_path;
        }
        return write_hex(p, path);
    }

    if (strcmp(argv[1], "--decimal") == 0 && argc >= 3) {
        uint64_t p = strtoull(argv[2], NULL, 10);
        const char *path = argc > 3 ? argv[3] : NULL;
        if (!path) {
            static char auto_path[64];
            snprintf(auto_path, sizeof(auto_path), "M%llu.txt",
                     (unsigned long long)p);
            path = auto_path;
        }
        return write_decimal(p, path);
    }

    if (strcmp(argv[1], "--next-candidates") == 0 && argc >= 3) {
        uint64_t p = strtoull(argv[2], NULL, 10);
        long double radius = argc > 3 ? strtold(argv[3], NULL) : 8.0L;
        int steps = argc > 4 ? atoi(argv[4]) : 200000;
        next_candidates(p, radius, steps);
        return 0;
    }

    if (strcmp(argv[1], "--verify-pln") == 0 && argc >= 4) {
        uint64_t p = strtoull(argv[2], NULL, 10);
        long double claimed = strtold(argv[3], NULL);
        verify_pln(p, claimed);
        return 0;
    }

    print_usage();
    return 1;
}
