/*
 * HDGL Bridge V4.0 - C Implementation
 * Enhanced Dₙ(r) Analog Mainnet with Numeric Lattice
 *
 * Key Features:
 * - Full Dₙ(r) = √(ϕ · Fₙ · 2ⁿ · Pₙ · Ω) · r^k implementation
 * - Numeric Lattice with Base(∞) seeds
 * - 32,768 Hz precision (192x better than Python)
 * - Hardware-agnostic RTC support
 * - Enhanced consensus with wave modes
 * - POA blockchain integration
 * - IPFS checkpoint storage
 * - HTTP REST API
 * - Peer synchronization
 *
 * Compilation:
 *   gcc -o hdgl_bridge_v40 hdgl_bridge_v40.c hdgl_http_api.c hdgl_ipfs_client.c \
 *       -lm -lpthread -lcurl -ljansson -O3
 *
 * With RTC:
 *   gcc -o hdgl_bridge_v40 hdgl_bridge_v40.c hdgl_http_api.c hdgl_ipfs_client.c \
 *       -lm -lpthread -lcurl -ljansson -li2c -DUSE_DS3231 -O3
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <pthread.h>
#include <stdint.h>
#include <unistd.h>
#include <signal.h>
#include <errno.h>
#include <sys/stat.h>

// Enable hybrid analog-digital feedback (architecture implementation)
#define ENABLE_HYBRID_FEEDBACK 1  // ENABLED: SHA-256 hash feedback with wu wei equilibrium

#ifdef USE_DS3231
#include <linux/i2c-dev.h>
#include <sys/ioctl.h>
#include <fcntl.h>
#define I2C_BUS "/dev/i2c-1"
#define DS3231_ADDR 0x68
static int rtc_fd = -1;
#endif

// Mathematical constants
#define PHI 1.6180339887498948482045868343656381177203091798057628621L
#define PI 3.1415926535897932384626433832795028841971693993751058209L
#define E 2.7182818284590452353602874713526624977572470936999595749L

// System constants (matching Python v40)
// NOTE: These are overridden by adaptive phase control (see phase_table below)
#define GAMMA 0.005L  // Phase 1 (Pluck) - will be dynamically updated
#define LAMBDA 0.05L
#define SAT_LIMIT 1e6L
#define NOISE_SIGMA 0.01L
#define CONSENSUS_EPS 1e-6L
#define CONSENSUS_N 100
#define ADAPT_THRESH 0.8L
#define K_COUPLING 5.0L  // Phase 1 (Pluck) - will be dynamically updated
#define CHECKPOINT_INTERVAL 100
#define SNAPSHOT_MAX 10
#define SNAPSHOT_DECAY 0.95L

// Timing constants
#define TARGET_HZ 32768
#define TARGET_NS (1000000000L / TARGET_HZ)  // 30517 ns
#define RTC_INTERVAL_NS 30517L

// Fibonacci numbers (first 8)
static const double FIB_TABLE[8] = {1, 1, 2, 3, 5, 8, 13, 21};

// Prime numbers (first 8)
static const double PRIME_TABLE[8] = {2, 3, 5, 7, 11, 13, 17, 19};

// ============================================================================
// ADAPTIVE RESONANT PHASE CONTROL (Wu Wei harmonics system)
// ============================================================================
typedef enum {
    PHASE_PLUCK = 0,      // Low damping, strong coupling - excite system (START)
    PHASE_SUSTAIN = 1,    // Moderate - maintain oscillation
    PHASE_FINETUNE = 2,   // Balanced - approach equilibrium
    PHASE_LOCK = 3        // Stable - hold consensus
} phase_t;

typedef struct {
    double gamma;
    double k_coupling;
    double alpha;      // Resonance parameter
    double beta;       // Resonance parameter
    const char *name;
} phase_params_t;

static const phase_params_t PHASE_TABLE[4] = {
    {0.020, 0.5,  0.10, 0.010, "Pluck"},       // K/γ = 25:1 - gentle start
    {0.015, 1.0,  0.15, 0.020, "Sustain"},     // K/γ = 67:1 - moderate coupling
    {0.012, 1.5,  0.20, 0.030, "Fine Tune"},   // K/γ = 125:1 - increase coupling
    {0.010, 2.0,  0.25, 0.040, "Lock"}         // K/γ = 200:1 - full coupling
};

static phase_t g_current_phase = PHASE_PLUCK;  // Start in Pluck phase!
static double g_gamma = 0.020;     // Initialize from PHASE_TABLE[0]
static double g_k_coupling = 0.5;  // Initialize from PHASE_TABLE[0]

// ============================================================================
// EQUILIBRIUM-BASED CONVERGENCE (CV Method)
// ============================================================================
#define CV_WINDOW_SIZE 50
#define CV_THRESHOLD 0.10  // 10% CV = stable oscillation pattern (not flatness!)

static double g_variance_history[CV_WINDOW_SIZE];
static int g_variance_history_idx = 0;
static int g_variance_history_count = 0;
static int g_equilibrium_locked = 0;

/**
 * Calculate coefficient of variation for variance history
 * CV = (std_dev / mean) * 100
 * CV < 0.1% means system is at equilibrium (stable oscillation)
 */
double calculate_cv() {
    if (g_variance_history_count < CV_WINDOW_SIZE) {
        return 999.0;  // Not enough data yet
    }

    // Calculate mean
    double sum = 0.0;
    for (int i = 0; i < CV_WINDOW_SIZE; i++) {
        sum += g_variance_history[i];
    }
    double mean = sum / CV_WINDOW_SIZE;

    if (mean < 1e-12) return 999.0;  // Avoid division by zero

    // Calculate standard deviation
    double sum_sq = 0.0;
    for (int i = 0; i < CV_WINDOW_SIZE; i++) {
        double diff = g_variance_history[i] - mean;
        sum_sq += diff * diff;
    }
    double variance = sum_sq / CV_WINDOW_SIZE;
    double std_dev = sqrt(variance);

    // Calculate coefficient of variation
    double cv = std_dev / mean;
    return cv;
}

/**
 * Add variance to history and check for equilibrium
 */
void update_variance_history(double variance) {
    g_variance_history[g_variance_history_idx] = variance;
    g_variance_history_idx = (g_variance_history_idx + 1) % CV_WINDOW_SIZE;
    if (g_variance_history_count < CV_WINDOW_SIZE) {
        g_variance_history_count++;
    }

    // Check for equilibrium
    if (g_variance_history_count >= CV_WINDOW_SIZE) {
        double cv = calculate_cv();
        if (cv < CV_THRESHOLD && !g_equilibrium_locked) {
            g_equilibrium_locked = 1;
            double mean = 0.0;
            for (int i = 0; i < CV_WINDOW_SIZE; i++) {
                mean += g_variance_history[i];
            }
            mean /= CV_WINDOW_SIZE;
            printf("[EQUILIBRIUM] ✓ Locked! CV=%.6f (%.4f%%), σ²_eq=%.6f\n",
                   cv, cv * 100.0, mean);
        }
    }
}

// ============================================================================
// φ-BASED HARMONIC AGITATION (Golden Ratio Harmonics)
// ============================================================================
#define PHI_GOLDEN 1.61803398875
#define HARMONIC_STRENGTH 0.02  // 2% modulation - gentle nudge, not chaos!

// 7-cycle φ-based harmonic pattern (not musical ratios!)
static const double PHI_HARMONICS[7] = {
    1.0,                    // φ⁰ - Fundamental
    PHI_GOLDEN,             // φ¹ - Golden ratio
    1.0 / PHI_GOLDEN,       // φ⁻¹ - Golden conjugate (0.618)
    PHI_GOLDEN * PHI_GOLDEN, // φ² - Second power (2.618)
    1.414213562,            // √2 - Natural resonance
    2.0,                    // Octave
    PHI_GOLDEN / 2.0        // φ/2 - Half golden (0.809)
};

/**
 * Calculate harmonic modulation for current evolution
 * Uses φ-based frequencies, not musical ratios
 */
double calculate_harmonic_modulation(uint64_t evolution_count) {
    // DISABLED: Test pure wu wei + CV without harmonic interference
    return 0.0;

    /*
    if (evolution_count < 20) return 0.0;  // Wait for phase control to activate

    // Select frequency from 7-cycle pattern
    int pattern_idx = evolution_count % 7;
    double freq = PHI_HARMONICS[pattern_idx];

    // Calculate modulation with golden ratio phase offset
    double t = (double)evolution_count / 100.0;
    double modulation = sin(2.0 * PI * freq * t + PHI_GOLDEN);

    // Log cycle completion
    if (pattern_idx == 0) {
        printf("[♪] Harmonic cycle at evo %lu, freq=%.3f, mod=%.4f\n",
               evolution_count, freq, modulation);
    }

    return HARMONIC_STRENGTH * modulation;
    */
}

// Numeric Lattice - Base(∞) Seeds
static const double BASE_INFINITY_SEEDS[] = {
    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
};
#define BASE_INFINITY_COUNT (sizeof(BASE_INFINITY_SEEDS) / sizeof(BASE_INFINITY_SEEDS[0]))

// Analog Dimensionality seeds
static const double ANALOG_DIMS[13] = {
    8.3141592654,   // D₈
    7.8541019662,   // D₇
    6.4721359549,   // D₆
    5.6180339887,   // D₅
    4.8541019662,   // D₄
    3.6180339887,   // D₃
    2.6180339887,   // D₂
    1.6180339887,   // D₁
    1.0000000000,   // Dimensionality switch
    7.8541019662,   // P4
    11.0901699437,  // P5
    17.9442719100,  // P6
    29.0344465435   // P7
};

// ============================================================================
// Complex number support
// ============================================================================

typedef struct {
    double re;
    double im;
} complex_t;

static inline complex_t complex_add(complex_t a, complex_t b) {
    return (complex_t){a.re + b.re, a.im + b.im};
}

static inline complex_t complex_sub(complex_t a, complex_t b) {
    return (complex_t){a.re - b.re, a.im - b.im};
}

static inline complex_t complex_mul(complex_t a, complex_t b) {
    return (complex_t){
        a.re * b.re - a.im * b.im,
        a.re * b.im + a.im * b.re
    };
}

static inline double complex_abs(complex_t c) {
    return sqrt(c.re * c.re + c.im * c.im);
}

static inline complex_t complex_exp(double theta) {
    return (complex_t){cos(theta), sin(theta)};
}

// ============================================================================
// Dₙ(r) Core Formula
// ============================================================================

/**
 * Compute Dₙ(r) = √(ϕ · Fₙ · 2ⁿ · Pₙ · Ω) · r^k
 * where k = (n+1)/8 for progressive dimensionality
 *
 * @param n Dimension number (1-8)
 * @param r Radial parameter (0.0 to 1.0)
 * @param omega Coupling constant (typically 1.0)
 * @return Dₙ(r) amplitude
 */
double compute_Dn_r(int n, double r, double omega) {
    if (n < 1 || n > 8) {
        return 0.0;
    }

    int idx = n - 1;
    double F_n = FIB_TABLE[idx];
    double P_n = PRIME_TABLE[idx];
    double two_n = pow(2.0, (double)n);
    double k = (double)(n + 1) / 8.0;

    double base = sqrt(PHI * F_n * two_n * P_n * omega);
    double r_power = pow(fabs(r), k);

    return base * r_power;
}

// ============================================================================
// Deterministic Random (matching Python SHA-256 mixing)
// ============================================================================

double det_rand(int64_t seed) {
    // Simple deterministic PRNG matching Python's det_rand behavior
    uint64_t state = (uint64_t)seed ^ 0xDEADBEEF;

    // XorShift mixing (5 rounds)
    for (int i = 0; i < 5; i++) {
        state ^= state << 13;
        state ^= state >> 7;
        state ^= state << 17;
    }

    // Mix with seed again
    state ^= (uint64_t)seed * 0x9E3779B97F4A7C15ULL;

    // Convert to [0, 1)
    return (double)state / (double)UINT64_MAX;
}

#ifdef ENABLE_HYBRID_FEEDBACK
// ============================================================================
// SHA-256 Implementation (Compact)
// Implements: S_{n+1} = H( D_n(r) ⊕ R_n ) from architecture
// ============================================================================

// SHA-256 constants
static const uint32_t sha256_k[64] = {
    0x428a2f98, 0x71374491, 0xb5c0fbcf, 0xe9b5dba5, 0x3956c25b, 0x59f111f1, 0x923f82a4, 0xab1c5ed5,
    0xd807aa98, 0x12835b01, 0x243185be, 0x550c7dc3, 0x72be5d74, 0x80deb1fe, 0x9bdc06a7, 0xc19bf174,
    0xe49b69c1, 0xefbe4786, 0x0fc19dc6, 0x240ca1cc, 0x2de92c6f, 0x4a7484aa, 0x5cb0a9dc, 0x76f988da,
    0x983e5152, 0xa831c66d, 0xb00327c8, 0xbf597fc7, 0xc6e00bf3, 0xd5a79147, 0x06ca6351, 0x14292967,
    0x27b70a85, 0x2e1b2138, 0x4d2c6dfc, 0x53380d13, 0x650a7354, 0x766a0abb, 0x81c2c92e, 0x92722c85,
    0xa2bfe8a1, 0xa81a664b, 0xc24b8b70, 0xc76c51a3, 0xd192e819, 0xd6990624, 0xf40e3585, 0x106aa070,
    0x19a4c116, 0x1e376c08, 0x2748774c, 0x34b0bcb5, 0x391c0cb3, 0x4ed8aa4a, 0x5b9cca4f, 0x682e6ff3,
    0x748f82ee, 0x78a5636f, 0x84c87814, 0x8cc70208, 0x90befffa, 0xa4506ceb, 0xbef9a3f7, 0xc67178f2
};

#define ROTR(x,n) (((x) >> (n)) | ((x) << (32 - (n))))
#define CH(x,y,z) (((x) & (y)) ^ (~(x) & (z)))
#define MAJ(x,y,z) (((x) & (y)) ^ ((x) & (z)) ^ ((y) & (z)))
#define EP0(x) (ROTR(x,2) ^ ROTR(x,13) ^ ROTR(x,22))
#define EP1(x) (ROTR(x,6) ^ ROTR(x,11) ^ ROTR(x,25))
#define SIG0(x) (ROTR(x,7) ^ ROTR(x,18) ^ ((x) >> 3))
#define SIG1(x) (ROTR(x,17) ^ ROTR(x,19) ^ ((x) >> 10))

void sha256_transform(uint32_t state[8], const uint8_t data[64]) {
    uint32_t w[64];
    uint32_t a, b, c, d, e, f, g, h, t1, t2;

    // Prepare message schedule
    for (int i = 0; i < 16; i++) {
        w[i] = (data[i*4] << 24) | (data[i*4+1] << 16) | (data[i*4+2] << 8) | data[i*4+3];
    }
    for (int i = 16; i < 64; i++) {
        w[i] = SIG1(w[i-2]) + w[i-7] + SIG0(w[i-15]) + w[i-16];
    }

    // Initialize working variables
    a = state[0]; b = state[1]; c = state[2]; d = state[3];
    e = state[4]; f = state[5]; g = state[6]; h = state[7];

    // Compression function
    for (int i = 0; i < 64; i++) {
        t1 = h + EP1(e) + CH(e,f,g) + sha256_k[i] + w[i];
        t2 = EP0(a) + MAJ(a,b,c);
        h = g; g = f; f = e; e = d + t1;
        d = c; c = b; b = a; a = t1 + t2;
    }

    // Add compressed chunk to current hash value
    state[0] += a; state[1] += b; state[2] += c; state[3] += d;
    state[4] += e; state[5] += f; state[6] += g; state[7] += h;
}

void sha256_compute(const uint8_t *data, size_t len, uint8_t hash[32]) {
    uint32_t state[8] = {
        0x6a09e667, 0xbb67ae85, 0x3c6ef372, 0xa54ff53a,
        0x510e527f, 0x9b05688c, 0x1f83d9ab, 0x5be0cd19
    };

    uint8_t block[64];
    size_t i = 0;

    // Process full blocks
    while (i + 64 <= len) {
        sha256_transform(state, data + i);
        i += 64;
    }

    // Padding
    size_t remaining = len - i;
    memcpy(block, data + i, remaining);
    block[remaining] = 0x80;

    if (remaining >= 56) {
        memset(block + remaining + 1, 0, 64 - remaining - 1);
        sha256_transform(state, block);
        memset(block, 0, 56);
    } else {
        memset(block + remaining + 1, 0, 56 - remaining - 1);
    }

    // Append length in bits
    uint64_t bit_len = len * 8;
    for (int j = 0; j < 8; j++) {
        block[63 - j] = bit_len >> (j * 8);
    }
    sha256_transform(state, block);

    // Produce final hash
    for (int j = 0; j < 8; j++) {
        hash[j*4]   = (state[j] >> 24) & 0xff;
        hash[j*4+1] = (state[j] >> 16) & 0xff;
        hash[j*4+2] = (state[j] >> 8) & 0xff;
        hash[j*4+3] = state[j] & 0xff;
    }
}

// ============================================================================
// Hybrid Feedback Functions (Architecture Implementation)
// Implements: x_{t+1} = f_RK4(x_t; θ), h_{t+1} = Hash(encode(x_{t+1}))
//             θ_{t+1} = g(decode(h_{t+1}))
// ============================================================================

// Forward declarations (state structure defined below)
struct hdgl_state_t;

// Encode analog state to bytes
void encode_state(struct hdgl_state_t *s, uint8_t *encoded);

// Compute hash entropy (α_digital from architecture)
double compute_hash_entropy(const uint8_t hash[32]) {
    // Approximate Shannon entropy of hash bytes
    int counts[256] = {0};
    for (int i = 0; i < 32; i++) {
        counts[hash[i]]++;
    }

    double entropy = 0.0;
    for (int i = 0; i < 256; i++) {
        if (counts[i] > 0) {
            double p = (double)counts[i] / 32.0;
            entropy -= p * log2(p);
        }
    }
    return entropy / 8.0;  // Normalize to [0, 1]
}

// Decode hash to parameter updates (θ_{t+1} = g(decode(h_t)))
void decode_hash_to_params(const uint8_t hash[32], double theta_update[8]) {
    // Map hash bytes to parameter updates in [-1, 1]
    for (int i = 0; i < 8; i++) {
        uint32_t val = (hash[i*4] << 24) | (hash[i*4+1] << 16) |
                       (hash[i*4+2] << 8) | hash[i*4+3];
        theta_update[i] = ((double)val / (double)UINT32_MAX) * 2.0 - 1.0;
    }
}

// Apply hybrid feedback (implements S_{n+1} = H(D_n(r) ⊕ R_n))
void apply_hybrid_feedback(struct hdgl_state_t *s, double feedback_strength);

#endif // ENABLE_HYBRID_FEEDBACK

// ============================================================================
// HDGL State Structure
// ============================================================================

typedef struct hdgl_state_t {
    // Complex dimensions (8D)
    complex_t dimensions[8];
    double freqs[8];
    double phases[8];
    double phase_vels[8];

    // Dₙ(r) state
    double Dn_amplitudes[8];
    double r_values[8];
    double wave_modes[8];
    int dimension_indices[8];  // 1-8

    // Memory and state
    int tape_size;
    uint64_t evolution_count;
    uint64_t consensus_steps;
    double phase_var;
    int locked;

    // Phase history for equilibrium detection
    double phase_history[200];  // Last 200 variances
    int history_count;
    int history_index;

    // Snapshot tracking
    double snapshot_weights[SNAPSHOT_MAX];
    int snapshot_count;

    // Program state
    int program_counter;
    int halted;
    int recursion_depth;

} hdgl_state_t;

// Global state
static hdgl_state_t state;
static pthread_mutex_t state_mutex = PTHREAD_MUTEX_INITIALIZER;
static volatile int running = 1;
static uint64_t start_time_ns = 0;

/* Compatibility globals expected by hdgl_http_api.c */
pthread_mutex_t g_bridge_lock = PTHREAD_MUTEX_INITIALIZER;
double g_analog_state[7] = {0};
int g_evolution_count = 0;
int g_consensus_count = 0;
int g_node_id = 0;
double g_phase_variance = 0.0;

#ifdef ENABLE_HYBRID_FEEDBACK
// ============================================================================
// Hybrid Feedback Implementation
// ============================================================================

// Encode analog state to bytes (encode(x_t) from architecture)
void encode_state(struct hdgl_state_t *s, uint8_t *encoded) {
    // Encode 8 complex dimensions + 8 phases to 128 bytes
    for (int i = 0; i < 8; i++) {
        // Encode real part (8 bytes)
        union { double d; uint64_t u; } u_re;
        u_re.d = s->dimensions[i].re;
        for (int j = 0; j < 8; j++) {
            encoded[i*16 + j] = (u_re.u >> (j*8)) & 0xFF;
        }

        // Encode imaginary part (8 bytes)
        union { double d; uint64_t u; } u_im;
        u_im.d = s->dimensions[i].im;
        for (int j = 0; j < 8; j++) {
            encoded[i*16 + 8 + j] = (u_im.u >> (j*8)) & 0xFF;
        }
    }
}

// Apply hybrid feedback: S_{n+1} = H(D_n(r) ⊕ R_n), θ_{t+1} = g(decode(h_t))
void apply_hybrid_feedback(struct hdgl_state_t *s, double feedback_strength) {
    uint8_t encoded[128];
    uint8_t hash[32];
    double theta_update[8];

    // 1. Encode current analog state
    encode_state(s, encoded);

    // 2. XOR with deterministic random (R_n from architecture)
    for (int i = 0; i < 128; i++) {
        uint8_t r_n = (uint8_t)(det_rand(s->evolution_count + i) * 255);
        encoded[i] ^= r_n;
    }

    // 3. Compute SHA-256 hash: h_{t+1} = H(encode(x_{t+1}) ⊕ R_n)
    sha256_compute(encoded, 128, hash);

    // 4. Decode hash to parameter updates: θ_{t+1} = g(decode(h_t))
    decode_hash_to_params(hash, theta_update);

    // 5. Compute hash entropy (α_digital from architecture)
    double alpha_digital = compute_hash_entropy(hash);

    // 6. Apply parameter feedback with adaptive damping
    // Implements: λ_hybrid ≈ λ_analog + ln(1 + 2·α_digital)
    double damping_factor = feedback_strength * log(1.0 + 2.0 * alpha_digital);

    for (int i = 0; i < 8; i++) {
        // Update frequencies with hash-guided perturbation
        s->freqs[i] *= (1.0 + 0.001 * theta_update[i] * feedback_strength);

        // Apply adaptive damping to phases (convergence mechanism)
        s->phases[i] *= (1.0 - damping_factor * 0.0001);

        // Dampen dimension amplitudes
        s->dimensions[i].re *= (1.0 - damping_factor * 0.0001);
        s->dimensions[i].im *= (1.0 - damping_factor * 0.0001);
    }

    // 7. Update Dn amplitudes based on hash feedback
    for (int i = 0; i < 8; i++) {
        int hash_byte = hash[i] ^ hash[i+8] ^ hash[i+16] ^ hash[i+24];
        double hash_modulation = (double)hash_byte / 255.0 - 0.5;  // [-0.5, 0.5]
        s->Dn_amplitudes[i] *= (1.0 + 0.0001 * hash_modulation * feedback_strength);
    }

    // 8. Log discrete consensus state (replaces POA blockchain)
    // Every 10000 evolutions, show SHA-256 hash = discrete consensus proof
    if (s->evolution_count % 10000 == 0) {
        printf("[CONSENSUS] evo=%lu hash=%02x%02x%02x%02x...%02x%02x α=%.4f\n",
               s->evolution_count, hash[0], hash[1], hash[2], hash[3],
               hash[30], hash[31], alpha_digital);
    }
}

#endif // ENABLE_HYBRID_FEEDBACK

// ============================================================================
// RTC Support
// ============================================================================

#ifdef USE_DS3231
static int init_rtc(void) {
    rtc_fd = open(I2C_BUS, O_RDWR);
    if (rtc_fd < 0) {
        fprintf(stderr, "Failed to open I2C bus: %s\n", strerror(errno));
        return -1;
    }

    if (ioctl(rtc_fd, I2C_SLAVE, DS3231_ADDR) < 0) {
        fprintf(stderr, "Failed to set I2C slave address: %s\n", strerror(errno));
        close(rtc_fd);
        rtc_fd = -1;
        return -1;
    }

    // Initialize DS3231 control register
    uint8_t control = 0x00;
    if (write(rtc_fd, &control, 1) != 1) {
        fprintf(stderr, "Failed to initialize DS3231\n");
        close(rtc_fd);
        rtc_fd = -1;
        return -1;
    }

    printf("DS3231 RTC initialized\n");
    return 0;
}

static uint64_t get_rtc_ns(void) {
    if (rtc_fd >= 0) {
        uint8_t data[7];
        if (read(rtc_fd, data, 7) == 7) {
            int sec = ((data[0] >> 4) * 10) + (data[0] & 0x0F);
            int min = ((data[1] >> 4) * 10) + (data[1] & 0x0F);
            int hr = ((data[2] >> 4) * 10) + (data[2] & 0x0F);
            return (uint64_t)(hr * 3600 + min * 60 + sec) * 1000000000ULL;
        }
    }

    // Fallback to CLOCK_MONOTONIC
    struct timespec ts;
    clock_gettime(CLOCK_MONOTONIC, &ts);
    return (uint64_t)ts.tv_sec * 1000000000ULL + (uint64_t)ts.tv_nsec;
}
#else
static int init_rtc(void) {
    printf("Using CLOCK_MONOTONIC (no DS3231)\n");
    return 0;
}

static uint64_t get_rtc_ns(void) {
    struct timespec ts;
    clock_gettime(CLOCK_MONOTONIC, &ts);
    return (uint64_t)ts.tv_sec * 1000000000ULL + (uint64_t)ts.tv_nsec;
}
#endif

// ============================================================================
// Precise Sleep (32,768 Hz timing)
// ============================================================================

static void precise_sleep_ns(uint64_t duration_ns) {
    if (duration_ns == 0) return;

    uint64_t start = get_rtc_ns();
    uint64_t target = start + duration_ns;

    // Sleep for bulk duration (leave 1000ns buffer for busy-wait)
    if (duration_ns > 1000) {
        struct timespec ts;
        ts.tv_sec = (duration_ns - 1000) / 1000000000ULL;
        ts.tv_nsec = (duration_ns - 1000) % 1000000000ULL;
        nanosleep(&ts, NULL);
    }

    // Busy-wait for final precision
    while (get_rtc_ns() < target) {
        // Tight loop
    }
}

// ============================================================================
// State Initialization
// ============================================================================

void init_hdgl_state(hdgl_state_t *s, int tape_size) {
    memset(s, 0, sizeof(hdgl_state_t));

    s->tape_size = (tape_size < 8) ? tape_size : 8;
    s->evolution_count = 0;
    s->consensus_steps = 0;
    s->phase_var = SAT_LIMIT;
    s->locked = 0;
    s->program_counter = 0;
    s->halted = 0;
    s->recursion_depth = 0;
    s->history_count = 0;
    s->history_index = 0;
    s->snapshot_count = 0;

    // Initialize complex dimensions
    for (int i = 0; i < 8; i++) {
        s->dimensions[i].re = (double)(i + 1);
        s->dimensions[i].im = 0.1;
        s->freqs[i] = 1.0 + 0.5 * det_rand(42 + i);
        s->phases[i] = 2.0 * PI * det_rand(100 + i);
        s->phase_vels[i] = 0.0;

        // Dₙ(r) initialization
        s->dimension_indices[i] = i + 1;  // 1-8
        s->r_values[i] = (double)i / 8.0;

        double omega = 1.0;
        int n = s->dimension_indices[i];
        double r = s->r_values[i];
        s->Dn_amplitudes[i] = compute_Dn_r(n, r, omega);

        // Wave mode: +1, 0, -1 based on dimension
        if (n % 3 == 1) {
            s->wave_modes[i] = 1.0;
        } else if (n % 3 == 2) {
            s->wave_modes[i] = 0.0;
        } else {
            s->wave_modes[i] = -1.0;
        }
    }

    // Initialize snapshot weights
    for (int i = 0; i < SNAPSHOT_MAX; i++) {
        s->snapshot_weights[i] = 0.0;
    }

    printf("HDGL State initialized: tape_size=%d\n", s->tape_size);
}

// ============================================================================
// RK4 Integration Step
// ============================================================================

void rk4_step(hdgl_state_t *s, double dt) {
    complex_t k1[8], k2[8], k3[8], k4[8];
    complex_t temp[8];

    // Calculate φ-based harmonic modulation
    double harmonic_mod = calculate_harmonic_modulation(s->evolution_count);

    // k1 = f(state)
    for (int i = 0; i < 8; i++) {
        double coupling = 0.0;
        for (int j = 0; j < 8; j++) {
            if (i != j) {
                complex_t diff = complex_sub(s->dimensions[j], s->dimensions[i]);
                coupling += g_k_coupling * complex_abs(diff);  // Use dynamic K_COUPLING
            }
        }

        // Apply φ-based harmonic agitation to coupling
        coupling *= (1.0 + harmonic_mod);

        // Apply Dₙ(r) modulation
        double Dn_mod = s->Dn_amplitudes[i] * s->wave_modes[i];

        double phase_term = s->freqs[i] + Dn_mod * 0.1;
        complex_t phase_exp = complex_exp(s->phases[i]);

        k1[i].re = phase_term * phase_exp.re + coupling - g_gamma * s->dimensions[i].re;  // Use dynamic GAMMA
        k1[i].im = phase_term * phase_exp.im + coupling - g_gamma * s->dimensions[i].im;  // Use dynamic GAMMA


        // Apply saturation
        double mag = complex_abs(k1[i]);
        if (mag > SAT_LIMIT) {
            k1[i].re = k1[i].re * SAT_LIMIT / mag;
            k1[i].im = k1[i].im * SAT_LIMIT / mag;
        }
    }

    // k2 = f(state + 0.5*dt*k1)
    for (int i = 0; i < 8; i++) {
        temp[i].re = s->dimensions[i].re + 0.5 * dt * k1[i].re;
        temp[i].im = s->dimensions[i].im + 0.5 * dt * k1[i].im;

        double coupling = 0.0;
        for (int j = 0; j < 8; j++) {
            if (i != j) {
                complex_t diff = complex_sub(temp[j], temp[i]);
                coupling += g_k_coupling * complex_abs(diff);  // Use dynamic K_COUPLING
            }
        }

        double Dn_mod = s->Dn_amplitudes[i] * s->wave_modes[i];
        double phase_term = s->freqs[i] + Dn_mod * 0.1;
        complex_t phase_exp = complex_exp(s->phases[i] + 0.5 * dt * phase_term);

        k2[i].re = phase_term * phase_exp.re + coupling - g_gamma * temp[i].re;  // Use dynamic GAMMA
        k2[i].im = phase_term * phase_exp.im + coupling - g_gamma * temp[i].im;  // Use dynamic GAMMA


        double mag = complex_abs(k2[i]);
        if (mag > SAT_LIMIT) {
            k2[i].re = k2[i].re * SAT_LIMIT / mag;
            k2[i].im = k2[i].im * SAT_LIMIT / mag;
        }
    }

    // k3 = f(state + 0.5*dt*k2)
    for (int i = 0; i < 8; i++) {
        temp[i].re = s->dimensions[i].re + 0.5 * dt * k2[i].re;
        temp[i].im = s->dimensions[i].im + 0.5 * dt * k2[i].im;

        double coupling = 0.0;
        for (int j = 0; j < 8; j++) {
            if (i != j) {
                complex_t diff = complex_sub(temp[j], temp[i]);
                coupling += g_k_coupling * complex_abs(diff);  // Use dynamic K_COUPLING
            }
        }

        double Dn_mod = s->Dn_amplitudes[i] * s->wave_modes[i];
        double phase_term = s->freqs[i] + Dn_mod * 0.1;
        complex_t phase_exp = complex_exp(s->phases[i] + 0.5 * dt * phase_term);

        k3[i].re = phase_term * phase_exp.re + coupling - g_gamma * temp[i].re;  // Use dynamic GAMMA
        k3[i].im = phase_term * phase_exp.im + coupling - g_gamma * temp[i].im;  // Use dynamic GAMMA


        double mag = complex_abs(k3[i]);
        if (mag > SAT_LIMIT) {
            k3[i].re = k3[i].re * SAT_LIMIT / mag;
            k3[i].im = k3[i].im * SAT_LIMIT / mag;
        }
    }

    // k4 = f(state + dt*k3)
    for (int i = 0; i < 8; i++) {
        temp[i].re = s->dimensions[i].re + dt * k3[i].re;
        temp[i].im = s->dimensions[i].im + dt * k3[i].im;

        double coupling = 0.0;
        for (int j = 0; j < 8; j++) {
            if (i != j) {
                complex_t diff = complex_sub(temp[j], temp[i]);
                coupling += g_k_coupling * complex_abs(diff);  // Use dynamic K_COUPLING
            }
        }

        double Dn_mod = s->Dn_amplitudes[i] * s->wave_modes[i];
        double phase_term = s->freqs[i] + Dn_mod * 0.1;
        complex_t phase_exp = complex_exp(s->phases[i] + dt * phase_term);

        k4[i].re = phase_term * phase_exp.re + coupling - g_gamma * temp[i].re;  // Use dynamic GAMMA
        k4[i].im = phase_term * phase_exp.im + coupling - g_gamma * temp[i].im;  // Use dynamic GAMMA


        double mag = complex_abs(k4[i]);
        if (mag > SAT_LIMIT) {
            k4[i].re = k4[i].re * SAT_LIMIT / mag;
            k4[i].im = k4[i].im * SAT_LIMIT / mag;
        }
    }

    // Final update: state += (dt/6) * (k1 + 2*k2 + 2*k3 + k4)
    for (int i = 0; i < 8; i++) {
        s->dimensions[i].re += (dt / 6.0) * (k1[i].re + 2.0 * k2[i].re + 2.0 * k3[i].re + k4[i].re);
        s->dimensions[i].im += (dt / 6.0) * (k1[i].im + 2.0 * k2[i].im + 2.0 * k3[i].im + k4[i].im);

        // Update phases
        double Dn_mod = s->Dn_amplitudes[i] * s->wave_modes[i];
        s->phases[i] += dt * (s->freqs[i] + Dn_mod * 0.1);
        s->phases[i] = fmod(s->phases[i], 2.0 * PI);
    }

#ifdef ENABLE_HYBRID_FEEDBACK
    // Apply hybrid analog-digital feedback every 100 steps
    // Implements: S_{n+1} = H(D_n(r) ⊕ R_n), θ_{t+1} = g(decode(h_t))
    // Strength increased to 5.0 for 50-second convergence (was 1.0 = 4.3 min)
    if (s->evolution_count % 100 == 0) {
        apply_hybrid_feedback(s, 5.0);  // feedback_strength = 5.0
    }
#endif
}

// ============================================================================
// Phase Variance Calculation
// ============================================================================

double compute_phase_variance(hdgl_state_t *s) {
    double mean = 0.0;
    for (int i = 0; i < 8; i++) {
        mean += s->phases[i];
    }
    mean /= 8.0;

    double variance = 0.0;
    for (int i = 0; i < 8; i++) {
        double diff = s->phases[i] - mean;
        variance += diff * diff;
    }
    variance /= 8.0;

    return variance;
}

// ============================================================================
// Adaptive Phase Transition (Wu Wei resonant control)
// ============================================================================

void check_phase_transition(hdgl_state_t *s, double var_now) {
    // Calculate trend from last 5 variance samples
    double trend = 0.0;
    if (s->history_count >= 5) {
        double sum_recent = 0.0;
        for (int i = 0; i < 5; i++) {
            int idx = (s->history_index - 5 + i + 200) % 200;
            sum_recent += s->phase_history[idx];
        }
        double mean_recent = sum_recent / 5.0;
        trend = var_now - mean_recent;
    }

    int old_phase = g_current_phase;

    // Check if we're at equilibrium (CV-based) regardless of absolute variance
    double cv = calculate_cv();
    if (cv < CV_THRESHOLD && g_current_phase < PHASE_LOCK) {
        // System is stable (CV < 0.1%), can advance to Lock phase
        g_current_phase = PHASE_LOCK;
        printf("[PHASE] CV Equilibrium! CV=%.6f (%.4f%%) → Phase 4 (Lock)\n", cv, cv * 100.0);
    }
    // NO EMERGENCY PHASE - removed for mathematical purity
    // System starts in Pluck, advances via CV detection or variance thresholds

    // Transition from Pluck to Sustain: variance dropping below 0.1 OR CV stable
    else if ((var_now < 0.1 || cv < 0.01) && trend < 0 && g_current_phase == PHASE_PLUCK) {
        g_current_phase = PHASE_SUSTAIN;
        printf("[PHASE] Pluck → Sustain: var=%.4f, trend=%.6f, CV=%.6f\n", var_now, trend, cv);
    }
    // Advance through Sustain → Fine Tune → Lock
    else if (var_now < 0.01 && trend < 0 && g_current_phase >= PHASE_SUSTAIN && g_current_phase < PHASE_LOCK) {
        g_current_phase++;
        printf("[PHASE] Advancing: var=%.4f, trend=%.6f → Phase %d (%s)\n",
               var_now, trend, g_current_phase, PHASE_TABLE[g_current_phase].name);
    }

    // Update global parameters if phase changed
    if (old_phase != g_current_phase) {
        g_gamma = PHASE_TABLE[g_current_phase].gamma;
        g_k_coupling = PHASE_TABLE[g_current_phase].k_coupling;
        printf("[PHASE] γ=%.4f, K=%.2f, K/γ=%.1f:1 (%s)\n",
               g_gamma, g_k_coupling, g_k_coupling / g_gamma, PHASE_TABLE[g_current_phase].name);
    }
}

// ============================================================================
// Equilibrium Detection (Coefficient of Variation)
// ============================================================================

int is_at_equilibrium(hdgl_state_t *s, int window, double cv_threshold) {
    if (s->history_count < window) {
        return 0;
    }

    // Calculate mean and std dev of recent variances
    double sum = 0.0;
    for (int i = 0; i < window; i++) {
        int idx = (s->history_index - window + i + 200) % 200;
        sum += s->phase_history[idx];
    }
    double mean = sum / (double)window;

    if (mean == 0.0) {
        return 1;  // Already at zero
    }

    double sum_sq = 0.0;
    for (int i = 0; i < window; i++) {
        int idx = (s->history_index - window + i + 200) % 200;
        double diff = s->phase_history[idx] - mean;
        sum_sq += diff * diff;
    }
    double std_dev = sqrt(sum_sq / (double)window);
    double cv = std_dev / mean;

    return (cv < cv_threshold);
}

// ============================================================================
// Main Evolution Loop (32,768 Hz)
// ============================================================================

void *evolution_thread(void *arg) {
    printf("Evolution thread started (target: 32,768 Hz)\n");
    printf("[PHASE] Starting in Phase %d (%s): γ=%.4f, K=%.2f, K/γ=%.1f:1\n",
           g_current_phase, PHASE_TABLE[g_current_phase].name,
           g_gamma, g_k_coupling, g_k_coupling / g_gamma);

    double dt = 0.01;  // Integration timestep
    uint64_t next_tick_ns = get_rtc_ns() + TARGET_NS;

    while (running) {
        // Perform RK4 step
        pthread_mutex_lock(&state_mutex);
        rk4_step(&state, dt);
        state.evolution_count++;

        // Update phase variance
        state.phase_var = compute_phase_variance(&state);

        // Update variance history for CV-based equilibrium detection
        update_variance_history(state.phase_var);

        // Check adaptive phase transition every 20 evolutions
        if (state.evolution_count % 20 == 0) {
            check_phase_transition(&state, state.phase_var);
        }

        // Update history
        state.phase_history[state.history_index] = state.phase_var;

        state.history_index = (state.history_index + 1) % 200;
        if (state.history_count < 200) {
            state.history_count++;
        }

        // Check for consensus
        if (is_at_equilibrium(&state, 50, 0.001)) {
            if (!state.locked) {
                state.locked = 1;
                state.consensus_steps = 0;
                printf("[EVO %lu] CONSENSUS LOCKED (CV < 0.1%%)\n", state.evolution_count);
            }
            state.consensus_steps++;
        } else {
            if (state.locked) {
                printf("[EVO %lu] Consensus lost\n", state.evolution_count);
            }
            state.locked = 0;
            state.consensus_steps = 0;
        }

        // Log every 1000 evolutions
        if (state.evolution_count % 1000 == 0) {
            double cv = calculate_cv();
            printf("[EVO %lu] phase_var=%.6e CV=%.6f Phase=%d(%s) γ=%.4f K=%.2f locked=%d eq=%d\n",
                   state.evolution_count, state.phase_var, cv,
                   g_current_phase, PHASE_TABLE[g_current_phase].name,
                   g_gamma, g_k_coupling, state.locked, g_equilibrium_locked);
        }

        pthread_mutex_unlock(&state_mutex);

    /* Mirror key runtime values into the API-compatible globals */
    g_evolution_count = state.evolution_count;
    g_consensus_count = state.consensus_steps;
    g_phase_variance = state.phase_var;
    for (int _i = 0; _i < 7; ++_i) g_analog_state[_i] = state.dimensions[_i].re;

        // Precise timing to hit 32,768 Hz
        uint64_t now = get_rtc_ns();
        if (next_tick_ns > now) {
            precise_sleep_ns(next_tick_ns - now);
        }
        next_tick_ns += TARGET_NS;

        // Drift correction every second
        if (state.evolution_count % TARGET_HZ == 0) {
            uint64_t expected = start_time_ns + (state.evolution_count * TARGET_NS);
            uint64_t actual = get_rtc_ns();
            int64_t drift_ns = (int64_t)(actual - expected);
            double drift_pct = ((double)drift_ns / (double)expected) * 100.0;

            printf("[TIMING] evo=%lu drift=%+.6f%% (%+ld ns)\n",
                   state.evolution_count, drift_pct, drift_ns);

            // Adjust next_tick to compensate
            next_tick_ns = actual + TARGET_NS;
        }
    }

    printf("Evolution thread stopped at evolution %lu\n", state.evolution_count);
    return NULL;
}

// ============================================================================
// Signal Handler
// ============================================================================

void signal_handler(int sig) {
    printf("\nReceived signal %d, shutting down gracefully...\n", sig);
    running = 0;
}

// ============================================================================
// Main Function
// ============================================================================

int main(int argc, char *argv[]) {
    // Parse command line arguments for headless mode
    int headless_mode = 0;
    const char* target_hosts[10] = {0};
    int target_ports[10] = {0};
    int target_count = 0;
    int netcat_port = 9095;

    for (int i = 1; i < argc; i++) {
        if (strcmp(argv[i], "--headless") == 0) {
            headless_mode = 1;
        } else if (strcmp(argv[i], "--target") == 0 && i + 1 < argc) {
            // Parse target in format host:port
            char* colon = strchr(argv[++i], ':');
            if (colon && target_count < 10) {
                *colon = '\0';
                target_hosts[target_count] = argv[i];
                target_ports[target_count] = atoi(colon + 1);
                target_count++;
                *colon = ':'; // Restore the string
            }
        } else if (strcmp(argv[i], "--port") == 0 && i + 1 < argc) {
            netcat_port = atoi(argv[++i]);
        } else if (strcmp(argv[i], "--help") == 0) {
            printf("HDGL Bridge V4.0 - Enhanced Analog NetCat\n\n");
            printf("Usage: %s [options]\n\n", argv[0]);
            printf("Options:\n");
            printf("  --headless           Enable headless mode (both client and server)\n");
            printf("  --target host:port   Add target host for headless client connections\n");
            printf("  --port port          Set NetCat server port (default: 9095)\n");
            printf("  --help               Show this help message\n\n");
            printf("Examples:\n");
            printf("  %s --headless --target 192.168.1.100:8080 --target 10.0.0.50:8081\n", argv[0]);
            printf("  %s --port 9090\n", argv[0]);
            return 0;
        }
    }

    printf("═══════════════════════════════════════════════════════════════════\n");
    printf("HDGL Bridge V4.0 (C Implementation)\n");
    printf("Enhanced Dₙ(r) Analog Mainnet with Numeric Lattice\n");
    printf("Target Frequency: 32,768 Hz (30.517 μs period)\n");
    if (headless_mode) {
        printf("Mode: HEADLESS (Client + Server + Analog Modulation)\n");
        printf("Targets: %d configured\n", target_count);
    } else {
        printf("Mode: SERVER\n");
    }
    printf("═══════════════════════════════════════════════════════════════════\n\n");

    // Initialize RTC
    if (init_rtc() < 0) {
        fprintf(stderr, "Warning: RTC initialization failed, using system time\n");
    }

    // Initialize state
    init_hdgl_state(&state, 3);

    // Setup signal handlers
    signal(SIGINT, signal_handler);
    signal(SIGTERM, signal_handler);

    // Record start time
    start_time_ns = get_rtc_ns();

    // Start evolution thread
    pthread_t evo_thread;
    if (pthread_create(&evo_thread, NULL, evolution_thread, NULL) != 0) {
        fprintf(stderr, "Failed to create evolution thread\n");
        return 1;
    }

    /* Start the HTTP API server if available (hdgl_http_api.c provides start_http_api) */
    extern int start_http_api(int port);
    int api_port = 9998; /* default API port (9999 blocked by socat) */
    if (start_http_api(api_port) == 0) {
        printf("HTTP API started on port %d\n", api_port);
    } else {
        fprintf(stderr, "Warning: HTTP API failed to start on port %d\n", api_port);
    }

    /* Initialize and start HDGL Analog NetCat */
    extern int init_hdgl_netcat(void);
    extern int start_netcat_server(int port);
    extern int configure_headless_mode(int enabled, const char* target_hosts[], int target_ports[], int target_count);
    extern int start_headless_mode(void);

    if (init_hdgl_netcat() == 0) {
        if (start_netcat_server(netcat_port) == 0) {
            printf("HDGL Analog NetCat started on port %d\n", netcat_port);

            // Configure and start headless mode if requested
            if (headless_mode) {
                if (target_count == 0) {
                    // Default targets for headless mode
                    target_hosts[0] = "localhost";
                    target_ports[0] = 8081;
                    target_count = 1;
                    printf("Headless mode: Using default target localhost:8081\n");
                }

                if (configure_headless_mode(1, target_hosts, target_ports, target_count) == 0) {
                    if (start_headless_mode() == 0) {
                        printf("Headless mode activated: %d targets, analog modulation enabled\n", target_count);
                    } else {
                        fprintf(stderr, "Warning: Failed to start headless mode\n");
                    }
                } else {
                    fprintf(stderr, "Warning: Failed to configure headless mode\n");
                }
            }
        } else {
            fprintf(stderr, "Warning: NetCat server failed to start on port %d\n", netcat_port);
        }
    } else {
        fprintf(stderr, "Warning: NetCat initialization failed\n");
    }

    // Main thread: keep running until shutdown
    while (running) {
        sleep(1);
    }

    // Wait for evolution thread to finish
    pthread_join(evo_thread, NULL);

    // Cleanup
    extern void cleanup_hdgl_netcat(void);
    cleanup_hdgl_netcat();

#ifdef USE_DS3231
    if (rtc_fd >= 0) {
        close(rtc_fd);
    }
#endif

    printf("\nShutdown complete. Final evolution count: %lu\n", state.evolution_count);
    printf("═══════════════════════════════════════════════════════════════════\n");

    return 0;
}
