/*
 * bigg_primality.c — Adapted Mersenne verification using the BigG/HDGL framework
 * ============================================================================
 *
 * This program reuses the fingerprint and PLN machinery from the BigG overfit
 * framework, then performs an exact Mersenne Lucas-Lehmer verification using GMP.
 *
 * Usage:
 *   gcc -O3 -std=c99 -DUSE_GMP -o bigg_primality.exe bigg_primality.c -lgmp -lm
 *   bigg_primality --verify 22514117
 *   bigg_primality --score 22514117
 *
 * NOTE: Lucas-Lehmer for p=22,514,117 is extremely large and may take a long time.
 */

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

#ifdef USE_GMP
#include <gmp.h>
#endif

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#define PHI      1.6180339887498948482
#define LN_PHI   0.4812118250596034748
#define LN2      0.6931471805599453094
#define LOG10_2  0.30102999566398119521
#define SQRT5    2.2360679774997896964
#define GOLDEN_ANG 3.8832220774530017295

#define NUM_DN 16
static const double FIB_TABLE[NUM_DN] = {
    1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987
};
static const double PRIME_TABLE[NUM_DN] = {
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53
};

static const double UPPER_FIELD[7] = {
    170.6180339887, 150.9442719100, 12.6180339887,
    8.8541019662, 4.2360679775, 3.6180339887, 1.6180339887
};
static const double ANALOG_DIMS[13] = {
    8.3141592654, 7.8541019662, 6.4721359549, 5.6180339887,
    4.8541019662, 3.6180339887, 2.6180339887, 1.6180339887,
    1.0000000000, 7.8541019662, 11.0901699437, 17.9442719100,
    29.0344465435
};
static const double LOWER_FIELD[8] = {
    0.0000000001, 0.0344465435, 0.0557280900, 0.0901699437,
    0.1458980338, 0.2360679775, 0.3819660113, 0.6180339887
};
static const double SIBLING_HARMONICS[8] = {
    0.0901699437, 0.1458980338, 0.2360679775, 0.3090169944,
    0.3819660113, 0.4721359549, 0.6545084972, 0.8729833462
};
static const double BASE_INF_SEEDS[64] = {
    0.6180339887, 1.6180339887, 2.6180339887, 3.6180339887, 4.8541019662,
    5.6180339887, 6.4721359549, 7.8541019662, 8.3141592654, 0.0901699437,
    0.1458980338, 0.2360679775, 0.3090169944, 0.3819660113, 0.4721359549,
    0.6545084972, 0.8729833462, 1.0000000000, 1.2360679775, 1.6180339887,
    2.2360679775, 2.6180339887, 3.1415926535, 3.6180339887, 4.2360679775,
    4.8541019662, 5.6180339887, 6.4721359549, 7.2360679775, 7.8541019662,
    8.6180339887, 9.2360679775, 9.8541019662, 10.6180339887, 11.0901699437,
    11.9442719100, 12.6180339887, 13.6180339887, 14.2360679775, 14.8541019662,
    15.6180339887, 16.4721359549, 17.2360679775, 17.9442719100, 18.6180339887,
    19.2360679775, 19.8541019662, 20.6180339887, 21.0901699437, 21.9442719100,
    22.6180339887, 23.6180339887, 24.2360679775, 24.8541019662, 25.6180339887,
    26.4721359549, 27.2360679775, 27.9442719100, 28.6180339887, 29.0344465435,
    29.6180339887, 30.2360679775, 30.8541019662, 31.6180339887
};

typedef struct {
    uint64_t p;
    double omega;
    int    D;
    double theta;
    double r;
} PLNCoord;

static double gen_fib(double x) {
    return (pow(PHI, x) - cos(M_PI * x) * pow(PHI, -x)) / SQRT5;
}

static double prismatic_omega_dim(double omega_p, int n) {
    return 0.5 + 0.5 * sin(omega_p + (double)n * GOLDEN_ANG);
}

static double compute_Dn_r(int n, double r, double omega_p) {
    if (n < 1 || n > NUM_DN) return 0.0;
    int idx = n - 1;
    double beta = omega_p - floor(omega_p);
    double x = (double)n + beta;
    double k = (double)(n + 1) / 8.0;
    double F_x = fabs(gen_fib(x));
    double Omega_n = prismatic_omega_dim(omega_p, n);
    double P_n = PRIME_TABLE[idx];
    double two_n = pow(2.0, (double)n);
    double phi_kx = pow(PHI, k * x);
    double coeff = sqrt(PHI * F_x * two_n * P_n * Omega_n * phi_kx);
    double delta = fabs(cos(M_PI * beta * PHI)) * log(P_n) / pow(PHI, x);
    double one_eff = 1.0 + delta;
    double scale = coeff * pow(fabs(r), k);
    double phase_re = one_eff * cos(M_PI * beta);
    double phase_im = one_eff * sin(M_PI * beta);
    double re = scale + phase_re;
    return sqrt(re * re + phase_im * phase_im);
}

static double orbit_r(double omega) {
    double frac = omega - floor(omega);
    double theta = frac * 6.2831853071795864769;
    double D_cont = fmod(omega, 8.0) + 1.0;
    double lat = PHI*PHI*(1.0 - pow(PHI, -D_cont)) + 1.0/D_cont;
    if (lat < 1.0) lat = 1.0;
    if (lat > PHI*PHI) lat = PHI*PHI;
    double Rre = cos(theta) - 1.0 + lat;
    double Rim = sin(theta);
    return sqrt(Rre*Rre + Rim*Rim);
}

static PLNCoord compute_pln(uint64_t p) {
    PLNCoord m;
    m.p = p;
    double log2x = (double)p;
    double lnx = log2x * LN2;
    double arg = lnx / LN_PHI;
    m.omega = (arg > 1.0) ? log(arg) / LN_PHI : 0.0;
    m.D = (int)fmod(floor(m.omega), 8.0) + 1;
    m.theta = fmod(m.omega, 1.0);
    m.r = orbit_r(m.omega);
    return m;
}

#define FP_LEN 48

typedef struct {
    uint64_t p;
    double features[FP_LEN];
    double score;
} BigGFingerprint;

static BigGFingerprint compute_fingerprint(uint64_t p) {
    BigGFingerprint fp;
    memset(&fp, 0, sizeof(fp));
    fp.p = p;
    PLNCoord pln = compute_pln(p);
    int fi = 0;
    for (int n = 1; n <= NUM_DN && fi < FP_LEN; n++) {
        fp.features[fi++] = compute_Dn_r(n, pln.r, pln.omega);
    }
    fp.features[fi++] = pln.omega;
    fp.features[fi++] = pln.theta;
    fp.features[fi++] = pln.r;
    fp.features[fi++] = (double)pln.D;
    fp.features[fi++] = fmod(pln.omega * UPPER_FIELD[0], 1.0);
    fp.features[fi++] = fmod(pln.omega * UPPER_FIELD[6], 1.0);
    fp.features[fi++] = fmod(pln.r * UPPER_FIELD[3], 1.0);
    fp.features[fi++] = fmod(pln.theta + LOWER_FIELD[7], 1.0);
    fp.features[fi++] = fmod(pln.omega * LOWER_FIELD[5], 1.0);
    fp.features[fi++] = fabs(sin(pln.omega * LOWER_FIELD[6] * M_PI));
    for (int i = 0; i < 4 && fi < FP_LEN; i++) {
        fp.features[fi++] = fabs(cos(pln.omega * SIBLING_HARMONICS[i*2] * M_PI));
    }
    for (int i = 0; i < 4 && fi < FP_LEN; i++) {
        double phase = fmod(pln.omega, ANALOG_DIMS[i]);
        fp.features[fi++] = phase / ANALOG_DIMS[i];
    }
    double min4[4] = {1e30, 1e30, 1e30, 1e30};
    for (int s = 0; s < 64 && fi < FP_LEN; s++) {
        double dist = fabs(fmod(pln.omega, 32.0) - BASE_INF_SEEDS[s]);
        for (int j = 0; j < 4; j++) {
            if (dist < min4[j]) {
                for (int k = 3; k > j; k--) min4[k] = min4[k-1];
                min4[j] = dist;
                break;
            }
        }
    }
    for (int i = 0; i < 4 && fi < FP_LEN; i++) {
        fp.features[fi++] = min4[i];
    }
    double d1 = compute_Dn_r(1, pln.r, pln.omega);
    double d8 = compute_Dn_r(8, pln.r, pln.omega);
    double d4 = compute_Dn_r(4, pln.r, pln.omega);
    double d5 = compute_Dn_r(5, pln.r, pln.omega);
    fp.features[fi++] = (d4*d5 > 1e-15) ? (d1*d8)/(d4*d5) : 0.0;
    double golden_dist = fabs(pln.r - 1.0/PHI);
    fp.features[fi++] = exp(-golden_dist * golden_dist * 10.0);
    fp.features[fi++] = fmod((double)p * PHI, 1.0);
    fp.features[fi++] = fmod((double)p * PHI * PHI, 1.0);
    fp.features[fi++] = fmod((double)p * 1.41421356237, 1.0);
    fp.features[fi++] = fmod((double)p * M_PI, 1.0);
    fp.features[fi++] = (double)(p % 6) / 6.0;
    fp.features[fi++] = (double)(p % 8) / 8.0;
    fp.features[fi++] = sqrt((double)p) - floor(sqrt((double)p));
    fp.features[fi++] = log2((double)p) - floor(log2((double)p));
    return fp;
}

static const uint64_t ALL_MERSENNE_P[] = {
    2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127,
    521, 607, 1279, 2203, 2281, 3217, 4253, 4423,
    9689, 9941, 11213, 19937, 21701, 23209, 44497,
    86243, 110503, 132049, 216091, 756839, 859433,
    1257787, 1398269, 2976221, 3021377, 6972593,
    13466917, 20996011, 24036583, 25964951, 30402457,
    32582657, 37156667, 42643801, 43112609,
    57885161, 74207281, 77232917, 82589933, 136279841
};
#define N_KNOWN (int)(sizeof(ALL_MERSENNE_P)/sizeof(ALL_MERSENNE_P[0]))

typedef struct {
    BigGFingerprint fps[N_KNOWN];
    double feat_std[FP_LEN];
    double feat_weight[FP_LEN];
} OverfitModel;

static OverfitModel fit_model(void) {
    OverfitModel model;
    for (int i = 0; i < N_KNOWN; i++) {
        model.fps[i] = compute_fingerprint(ALL_MERSENNE_P[i]);
    }
    for (int f = 0; f < FP_LEN; f++) {
        double sum = 0, sum2 = 0;
        for (int i = 0; i < N_KNOWN; i++) {
            sum += model.fps[i].features[f];
            sum2 += model.fps[i].features[f] * model.fps[i].features[f];
        }
        double mean = sum / N_KNOWN;
        double var = sum2 / N_KNOWN - mean * mean;
        model.feat_std[f] = (var > 1e-20) ? sqrt(var) : 1.0;
        if (var < 1e-10) {
            model.feat_weight[f] = 0.0;
        } else {
            model.feat_weight[f] = 1.0 / model.feat_std[f];
        }
    }
    double wsum = 0;
    for (int f = 0; f < FP_LEN; f++) wsum += model.feat_weight[f];
    if (wsum > 0) {
        for (int f = 0; f < FP_LEN; f++) model.feat_weight[f] /= wsum;
    }
    return model;
}

static double score_candidate(const OverfitModel *model, const BigGFingerprint *fp) {
    double max_sim = 0.0;
    for (int i = 0; i < N_KNOWN; i++) {
        double wdist2 = 0.0;
        for (int f = 0; f < FP_LEN; f++) {
            if (model->feat_weight[f] < 1e-15) continue;
            double diff = fp->features[f] - model->fps[i].features[f];
            double z = diff / model->feat_std[f];
            wdist2 += model->feat_weight[f] * z * z;
        }
        double sim = exp(-wdist2 / 2.0);
        if (sim > max_sim) max_sim = sim;
    }
    return max_sim;
}

static inline uint64_t mulmod64(uint64_t a, uint64_t b, uint64_t m) {
    return (uint64_t)((__uint128_t)a * b % m);
}

static uint64_t powmod64(uint64_t base, uint64_t exp, uint64_t mod) {
    if (mod == 1) return 0;
    uint64_t result = 1;
    base %= mod;
    while (exp) {
        if (exp & 1) result = mulmod64(result, base, mod);
        base = mulmod64(base, base, mod);
        exp >>= 1;
    }
    return result;
}

static bool is_prime_u64(uint64_t n) {
    if (n < 2) return false;
    if (n == 2 || n == 3 || n == 5 || n == 7) return true;
    if (n % 2 == 0 || n % 3 == 0 || n % 5 == 0) return false;
    for (uint64_t d = 7; d <= 1000 && d * d <= n; d += 2)
        if (n % d == 0) return false;
    uint64_t d = n - 1; int s = 0;
    while (!(d & 1)) { d >>= 1; s++; }
    static const uint64_t bases[] = {2,325,9375,28178,450775,9780504,1795265022};
    for (int i = 0; i < 7; i++) {
        uint64_t a = bases[i] % n;
        if (a == 0) return true;
        uint64_t x = powmod64(a, d, n);
        if (x == 1 || x == n - 1) continue;
        bool composite = true;
        for (int r = 1; r < s; r++) {
            x = mulmod64(x, x, n);
            if (x == n - 1) { composite = false; break; }
        }
        if (composite) return false;
    }
    return true;
}

static uint64_t trial_factor(uint64_t p, uint64_t trial_limit) {
    for (uint64_t k = 1; ; k++) {
        if (k > (trial_limit - 1) / (2 * p)) break;
        uint64_t q = 2 * k * p + 1;
        if (q > trial_limit) break;
        uint64_t qmod8 = q % 8;
        if (qmod8 != 1 && qmod8 != 7) continue;
        if (q % 3 == 0 || q % 5 == 0 || q % 7 == 0 || q % 11 == 0 || q % 13 == 0) {
            if (q > 13) continue;
        }
        if (powmod64(2, p, q) == 1) return q;
    }
    return 0;
}

#ifdef USE_GMP
static void reduce_mod_mersenne(mpz_t x, uint64_t p, const mpz_t M) {
    mpz_t low, high, tmp;
    mpz_inits(low, high, tmp, NULL);
    while ((uint64_t)mpz_sizeinbase(x, 2) > p) {
        mpz_tdiv_r_2exp(low, x, p);
        mpz_tdiv_q_2exp(high, x, p);
        mpz_add(tmp, low, high);
        mpz_set(x, tmp);
    }
    if (mpz_cmp(x, M) >= 0) mpz_sub(x, x, M);
    mpz_clears(low, high, tmp, NULL);
}

static bool lucas_lehmer_mersenne(uint64_t p) {
    if (p < 2) return false;
    mpz_t M, s, t;
    mpz_inits(M, s, t, NULL);
    mpz_set_ui(M, 1);
    mpz_mul_2exp(M, M, p);
    mpz_sub_ui(M, M, 1);
    mpz_set_ui(s, 4);
    for (uint64_t i = 0; i < p - 2; i++) {
        mpz_mul(t, s, s);
        mpz_sub_ui(t, t, 2);
        reduce_mod_mersenne(t, p, M);
        mpz_set(s, t);
    }
    bool is_prime = (mpz_cmp_ui(s, 0) == 0);
    mpz_clears(M, s, t, NULL);
    return is_prime;
}
#endif

static void print_usage(void) {
    printf("BigG Primality Adapter\n");
    printf("Usage:\n");
    printf("  bigg_primality --score <p1> [p2 ...]      Score candidates with HDGL fingerprint\n");
    printf("  bigg_primality --verify <p>               Verify candidate with trial factor + GMP LL\n");
    printf("  bigg_primality --help\n");
    printf("\n");
    printf("This program adapts the BigG fingerprint framework into a candidate verifier.\n");
    printf("Exact Lucas-Lehmer requires GMP and may be very slow for p >> 10^6.\n");
}

static void cmd_score(uint64_t *exponents, int n) {
    OverfitModel model = fit_model();
    for (int i = 0; i < n; i++) {
        uint64_t p = exponents[i];
        BigGFingerprint fp = compute_fingerprint(p);
        fp.score = score_candidate(&model, &fp);
        PLNCoord pln = compute_pln(p);
        printf("p=%-12llu  score=%.6f  omega=%.10f  D=%d\n",
               (unsigned long long)p, fp.score, pln.omega, pln.D);
    }
}

static void cmd_verify(uint64_t p) {
    OverfitModel model = fit_model();
    printf("=== HDGL candidate verification for p=%llu ===\n", (unsigned long long)p);
    printf("1) exponent primality: %s\n", is_prime_u64(p) ? "prime" : "composite");
    uint64_t factor = trial_factor(p, (uint64_t)1 << 24);
    if (factor) {
        printf("2) small factor found: q=%llu, so M_%llu is composite\n",
               (unsigned long long)factor, (unsigned long long)p);
    } else {
        printf("2) no factor found below 2^24; candidate survives trial factoring\n");
    }
    BigGFingerprint fp = compute_fingerprint(p);
    fp.score = score_candidate(&model, &fp);
    PLNCoord pln = compute_pln(p);
    printf("3) BigG fingerprint score: %.6f\n", fp.score);
    printf("   omega=%.10f  D=%d  r=%.10f\n", pln.omega, pln.D, pln.r);
#ifdef USE_GMP
    printf("4) performing Lucas-Lehmer...\n");
    bool ll = lucas_lehmer_mersenne(p);
    printf("   Lucas-Lehmer result: %s\n", ll ? "prime" : "composite");
#else
    printf("4) GMP not enabled: exact Lucas-Lehmer unavailable. Recompile with -DUSE_GMP -lgmp.\n");
#endif
}

int main(int argc, char *argv[]) {
    if (argc < 2) { print_usage(); return 1; }
    if (strcmp(argv[1], "--help") == 0 || strcmp(argv[1], "-h") == 0) {
        print_usage();
        return 0;
    }
    if (strcmp(argv[1], "--score") == 0 && argc >= 3) {
        uint64_t *exponents = malloc((argc - 2) * sizeof(uint64_t));
        int n = 0;
        for (int i = 2; i < argc; i++) {
            uint64_t p = strtoull(argv[i], NULL, 10);
            if (p > 1) exponents[n++] = p;
        }
        cmd_score(exponents, n);
        free(exponents);
        return 0;
    }
    if (strcmp(argv[1], "--verify") == 0 && argc == 3) {
        uint64_t p = strtoull(argv[2], NULL, 10);
        if (p < 2) {
            fprintf(stderr, "Invalid exponent: %s\n", argv[2]);
            return 1;
        }
        cmd_verify(p);
        return 0;
    }
    print_usage();
    return 1;
}
