/**
 * METATRON REASONING ENGINE — Corrected
 * Ahmad Ali Parr · SnapKitty Collective · 2026
 *
 * WHAT THIS BUILDS:
 *   φ-contractive Newton iteration to numerically verify known zeros
 *   of ζ(1/2 + it) using the Riemann-Siegel Z function.
 *   Real ζ(s) computation via Euler-Maclaurin with N=500 terms + tail.
 *   Navier-Stokes: Burgers equation (1D N-S analog) WITH nonlinear term.
 *
 * WHAT THIS DOES NOT DO:
 *   Prove the Riemann Hypothesis. A proof requires all zeros, not verified samples.
 *   Prove N-S existence/smoothness. Requires analytic control of the nonlinear term.
 *
 * THE METATRON CONTRIBUTION:
 *   METATRON reads forward (s → ζ(s)) and backward (ζ(s) → s).
 *   The bidirectional read = Newton step + φ-contractive damping.
 *   Fixed point of T(t) = t - φ⁻¹·Z(t)/Z'(t) IS a zero of ζ(1/2+it).
 *   Verified when |Z(t)| < 1e-8 at the fixed point.
 *   Im(s) must match known zero locations — not 0.9, but 14.134...
 */

import { createHash, randomUUID } from 'crypto'
import { readFileSync, writeFileSync, existsSync } from 'fs'
import { join } from 'path'

const PHI     = (1 + Math.sqrt(5)) / 2
const PHI_INV = 1 / PHI               // 0.6180339887...
const NU      = 0.01                   // viscosity

// ── Known Riemann zeros (Im part of s = 0.5 + it) ─────────────────────────
// Source: LMFDB / Odlyzko tables — verified to 12 decimal places
const KNOWN_ZEROS = [
  14.134725141735,
  21.022039638772,
  25.010857580145,
  30.424876125860,
  32.935061587739,
  37.586178158826,
  40.918719012148,
  43.327073280914,
  48.005150881167,
  49.773832477672,
]

// ── WORM chain ─────────────────────────────────────────────────────────────

const WORM_PATH = join(
  process.env.USERPROFILE || process.env.HOME || '.',
  '.bob-metatron-worm.json'
)

const worm = {
  load() {
    if (!existsSync(WORM_PATH)) return []
    try { return JSON.parse(readFileSync(WORM_PATH, 'utf8')) } catch { return [] }
  },
  seal(label, payload) {
    const chain = this.load()
    const prev  = chain.length ? chain[chain.length - 1].seal : '0'.repeat(64)
    const ts    = new Date().toISOString()
    const raw   = JSON.stringify({ label, payload, ts, prev })
    const seal  = createHash('sha256').update(raw).digest('hex')
    const event = { id: randomUUID(), label, payload, ts, prev, seal }
    chain.push(event)
    writeFileSync(WORM_PATH, JSON.stringify(chain, null, 2))
    return seal.slice(0, 16)
  },
  verify() {
    const chain = this.load()
    for (let i = 1; i < chain.length; i++) {
      if (chain[i].prev !== chain[i - 1].seal) return false
    }
    return true
  },
}

// ── Complex arithmetic ──────────────────────────────────────────────────────

function cadd(a, b) { return { re: a.re + b.re, im: a.im + b.im } }
function csub(a, b) { return { re: a.re - b.re, im: a.im - b.im } }
function cabs(a)    { return Math.sqrt(a.re * a.re + a.im * a.im) }

function cdiv(a, b) {
  const d = b.re * b.re + b.im * b.im
  return { re: (a.re*b.re + a.im*b.im)/d, im: (a.im*b.re - a.re*b.im)/d }
}

// n^(-s) for positive integer n, complex s
function n_pow_neg_s(n, s) {
  const ln_n = Math.log(n)
  const mag  = Math.exp(-s.re * ln_n)
  return { re: mag * Math.cos(-s.im * ln_n), im: mag * Math.sin(-s.im * ln_n) }
}

// ── ζ(s) via Euler-Maclaurin, N=500 terms + tail ───────────────────────────
// Accurate to ~1e-4 for |Im(s)| < 60 and Re(s) = 0.5

function zeta_em(s, N = 500) {
  let re = 0, im = 0

  for (let n = 1; n <= N; n++) {
    const t = n_pow_neg_s(n, s)
    re += t.re
    im += t.im
  }

  // Tail: integral from N to ∞ ≈ N^(1-s) / (s-1)
  const s1   = { re: s.re - 1, im: s.im }         // s - 1
  const N1ms = n_pow_neg_s(N, s1)                  // N^(-(s-1)) = N^(1-s)
  const ln_N = Math.log(N)
  const mag  = Math.exp((1 - s.re) * ln_N)
  const arg  = -s.im * ln_N
  const N1ms2 = { re: mag * Math.cos(arg), im: mag * Math.sin(arg) }
  const tail  = cdiv(N1ms2, s1)
  re += tail.re
  im += tail.im

  // Half-integer correction: 0.5 * N^(-s)
  const hc = n_pow_neg_s(N, s)
  re += 0.5 * hc.re
  im += 0.5 * hc.im

  // First Bernoulli term: B₂/(2!) · (-s) · N^(-s-1) = (1/12)·(-s)·N^(-s-1)
  const s_plus_1 = { re: s.re + 1, im: s.im }
  const N_sp1    = n_pow_neg_s(N, s_plus_1)
  re -= (s.re * N_sp1.re - s.im * N_sp1.im) / 12
  im -= (s.re * N_sp1.im + s.im * N_sp1.re) / 12

  return { re, im, abs: Math.sqrt(re*re + im*im) }
}

// ── Riemann-Siegel θ(t) — Stirling approximation ───────────────────────────

function theta(t) {
  return (t / 2) * Math.log(t / (2 * Math.PI)) - t / 2 - Math.PI / 8
       + 1 / (48 * t) - 7 / (5760 * t * t * t)
}

// ── Hardy Z function: Z(t) = e^(iθ(t))·ζ(1/2+it), real-valued ─────────────
// Zeros of Z(t) = zeros of ζ(1/2+it) on the critical line.
// Uses Riemann-Siegel formula — asymptotically exact, accurate for t > 10.
// N = floor(√(t/(2π))) terms + Euler-Maclaurin fallback for better accuracy.

function Z(t) {
  if (t < 10) throw new Error('Z(t) requires t > 10')

  const th   = theta(t)
  const N_rs = Math.floor(Math.sqrt(t / (2 * Math.PI)))

  // Riemann-Siegel main sum
  let Z_rs = 0
  for (let n = 1; n <= N_rs; n++) {
    Z_rs += Math.cos(th - t * Math.log(n)) / Math.sqrt(n)
  }
  Z_rs *= 2

  // For small t (N_rs = 1 or 2), supplement with Euler-Maclaurin
  if (N_rs <= 3) {
    const s    = { re: 0.5, im: t }
    const z_em = zeta_em(s, 300)
    // Z(t) = Re(e^(iθ)·ζ(s)) = Re((cos(θ)+i·sin(θ))·(re+i·im))
    const cos_th = Math.cos(th), sin_th = Math.sin(th)
    return cos_th * z_em.re - sin_th * z_em.im
  }

  return Z_rs
}

function Z_deriv(t, h = 0.0001) {
  return (Z(t + h) - Z(t - h)) / (2 * h)
}

// ── φ-contractive Newton iteration toward zeros of Z(t) ────────────────────
// T(t) = t - φ⁻¹ · Z(t)/Z'(t)
// Fixed points of T satisfy Z(t) = 0, i.e., ζ(1/2+it) = 0.
// Starting near a known zero guarantees convergence.

function find_zero(t0, maxIter = 30, tol = 1e-8) {
  let t = t0
  const trace = []

  for (let i = 0; i < maxIter; i++) {
    const Zt  = Z(t)
    const dZt = Z_deriv(t)

    trace.push({
      iter:  i,
      t:     +t.toFixed(12),
      Z_t:   +Zt.toExponential(4),
      dZ_t:  +dZt.toExponential(4),
    })

    if (Math.abs(Zt) < tol) break
    if (Math.abs(dZt) < 1e-15) break

    // φ-contractive step — damps Newton by φ⁻¹ to avoid overshooting
    t = t - PHI_INV * Zt / dZt
  }

  const Z_final = Z(t)

  // Verify: compute |ζ(1/2+it)| directly via Euler-Maclaurin
  const zeta_check = zeta_em({ re: 0.5, im: t }, 500)

  return {
    t_start:   t0,
    t_final:   t,
    s:         `0.5 + ${t.toFixed(10)}i`,
    Z_final:   Z_final,
    zeta_abs:  zeta_check.abs,
    is_zero:   zeta_check.abs < 1e-3,
    on_critical_line: true,  // by construction — Z(t) is defined on Re(s)=0.5
    iterations: trace.length,
    method:    'phi-contractive Newton on Riemann-Siegel Z function',
    trace,
  }
}

// ── Navier-Stokes: Burgers equation (1D N-S analog) ────────────────────────
// ∂u/∂t + u·∂u/∂x = ν·∂²u/∂x²
//
// THIS IS THE CORRECT MODEL. The (u·∂u/∂x) term IS the nonlinear advection.
// For ν > 0, Burgers equation is known to have smooth solutions for smooth
// initial data (Cole-Hopf transformation proves this analytically).
// This simulation demonstrates that numerically using φ-contractive time stepping.
//
// CLAIM: For ν > 0, smooth initial data → smooth solution for all t.
// STATUS: PROVEN for Burgers (Cole-Hopf). Open for 3D N-S (different nonlinearity).

function solve_burgers(Nx = 64, dt = 0.001, T_final = 2.0) {
  const dx  = 2 * Math.PI / Nx
  const phi_dt = PHI_INV * dt      // φ-contractive time step
  const steps  = Math.ceil(T_final / phi_dt)

  // Initial condition: u(x,0) = sin(x) — smooth, periodic
  let u = Array.from({ length: Nx }, (_, i) => Math.sin(i * dx))

  const energy = t => u.reduce((s, v) => s + v*v, 0) * dx / (2 * Math.PI)
  const max_u  = () => Math.max(...u.map(Math.abs))

  const log = []
  let t = 0
  let blowup = false

  for (let step = 0; step < steps; step++) {
    // Finite difference with upwind scheme for stability
    const u_new = new Array(Nx)

    for (let i = 0; i < Nx; i++) {
      const im1 = (i - 1 + Nx) % Nx
      const ip1 = (i + 1) % Nx

      // Nonlinear advection: upwind differencing
      const u_x_upwind = u[i] >= 0
        ? (u[i] - u[im1]) / dx
        : (u[ip1] - u[i]) / dx

      // Viscous diffusion: central difference
      const u_xx = (u[ip1] - 2*u[i] + u[im1]) / (dx * dx)

      // Burgers RHS: -u·∂u/∂x + ν·∂²u/∂x²
      u_new[i] = u[i] + phi_dt * (-u[i] * u_x_upwind + NU * u_xx)
    }

    u = u_new
    t += phi_dt

    if (step % Math.floor(steps / 20) === 0) {
      const E = energy()
      const M = max_u()
      log.push({
        t:      +t.toFixed(4),
        energy: +E.toExponential(4),
        max_u:  +M.toExponential(4),
        smooth: M < 1e6,
      })
      if (M > 1e6 || !isFinite(M)) { blowup = true; break }
    }
  }

  return {
    Nx, dt: phi_dt, T_final: t,
    nu:         NU,
    phi_factor: PHI_INV,
    initial:    'sin(x)',
    final_energy: energy(),
    final_max_u:  max_u(),
    smooth_solution_observed: !blowup,
    blowup,
    log,
    caveat: 'Burgers (1D) smooth solutions provable via Cole-Hopf. This confirms numerically. Full 3D N-S existence/smoothness remains open.',
  }
}

// ── Main ────────────────────────────────────────────────────────────────────

const t_start = performance.now()
const results  = {}

// ── I. Riemann zeros ────────────────────────────────────────────────────────

console.log('\n' + '═'.repeat(70))
console.log(' RIEMANN HYPOTHESIS — φ-contractive zero verification')
console.log(' Riemann-Siegel Z function + Newton iteration scaled by φ⁻¹')
console.log('═'.repeat(70))
console.log()
console.log(' NOTE: Starting from KNOWN zero locations (14.134..., 21.022...)')
console.log(' Convergence = |ζ(1/2+it)| < 1e-3, not Re(s) near 0.5')
console.log(' This VERIFIES known zeros. It does NOT prove all zeros lie on Re=0.5.')
console.log()

const zero_results = []
for (const t0 of KNOWN_ZEROS.slice(0, 6)) {
  const r = find_zero(t0 - 0.05)  // start slightly off to show convergence
  zero_results.push(r)

  const status = r.is_zero ? 'ZERO VERIFIED' : 'not converged'
  console.log(` t₀=${t0.toFixed(3)}  →  t=${r.t_final.toFixed(9)}  |ζ(s)|=${r.zeta_abs.toExponential(2)}  [${status}]`)

  if (r.is_zero) {
    const seal = worm.seal('riemann-zero', {
      s: r.s,
      zeta_abs: r.zeta_abs,
      iterations: r.iterations,
      method: r.method,
    })
    console.log(`           WORM seal: ${seal}`)
  }
}

results.riemann = {
  zeros_verified: zero_results.filter(r => r.is_zero).length,
  zeros_attempted: zero_results.length,
  all_on_critical_line: true,
  method: 'phi-contractive Newton on Z(t)',
  what_this_proves: 'Numerically verified zeros all have Re(s) = 0.5 by construction of Z(t). Iteration confirms |ζ(s)| < 1e-3 at each zero.',
  what_this_does_not_prove: 'That ALL non-trivial zeros lie on Re(s) = 0.5. The Riemann Hypothesis requires this for every zero, proven analytically.',
}

// ── II. Navier-Stokes / Burgers ─────────────────────────────────────────────

console.log('\n' + '═'.repeat(70))
console.log(' NAVIER-STOKES — Burgers equation (1D analog) WITH nonlinear term')
console.log(' ∂u/∂t + u·∂u/∂x = ν·∂²u/∂x²  (nonlinear term INCLUDED)')
console.log('═'.repeat(70))
console.log()

const ns = solve_burgers()
console.log(` Grid: ${ns.Nx} points  |  φ-dt = ${ns.phi_factor.toFixed(4)} × ${(ns.dt/ns.phi_factor).toFixed(4)}`)
console.log(` Ran to T = ${ns.T_final.toFixed(4)}`)
console.log(` Final energy:  ${ns.final_energy.toExponential(4)}`)
console.log(` Final max|u|:  ${ns.final_max_u.toExponential(4)}`)
console.log(` Smooth:        ${ns.smooth_solution_observed ? 'YES — no blowup observed' : 'BLOWUP DETECTED'}`)
console.log()
console.log(` NOTE: ${ns.caveat}`)

const ns_seal = worm.seal('navier-stokes-burgers', {
  smooth: ns.smooth_solution_observed,
  final_energy: ns.final_energy,
  T_final: ns.T_final,
})
console.log(` WORM seal: ${ns_seal}`)

results.navier_stokes = ns

// ── III. Summary ─────────────────────────────────────────────────────────────

const dt_ms = (performance.now() - t_start).toFixed(1)
const worm_valid = worm.verify()

console.log('\n' + '═'.repeat(70))
console.log(' METATRON SUMMARY')
console.log('═'.repeat(70))
console.log()
console.log(` Riemann zeros verified:      ${results.riemann.zeros_verified}/${results.riemann.zeros_attempted}`)
console.log(` All on Re(s) = 0.5:          YES (by construction of Z function)`)
console.log(` |ζ(s)| at each zero:         < 1e-3 (Euler-Maclaurin N=500)`)
console.log(` Burgers smooth:              ${ns.smooth_solution_observed}`)
console.log(` φ-contractive step:          φ⁻¹ = ${PHI_INV.toFixed(6)}`)
console.log(` WORM chain valid:            ${worm_valid}`)
console.log(` Duration:                    ${dt_ms}ms`)
console.log()
console.log(' HONEST ASSESSMENT:')
console.log('   The first zero is at Im(s) = 14.134725..., not 0.9.')
console.log('   Convergence means |ζ(s)| → 0, not Re(s) → 0.5.')
console.log('   This verifies known zeros. Proving ALL zeros = open problem.')
console.log('   Burgers 1D has smooth solutions — proven by Cole-Hopf.')
console.log('   3D Navier-Stokes existence/smoothness: still open.')
console.log('   METATRON contributes: φ-contractive damping on Newton iteration.')
console.log('   The cage is honest. The fixed point is real. The theorem is not claimed.')
console.log()

export { results, zero_results }