// Name: PhD Math Solver
// Description: Advanced symbolic and numeric mathematics toolkit with calculus, algebra, optimization, transforms, and more.
// Author: Nko00
// GitHub: Nko00

import "@johnlindquist/kit"
import { create, all, MathNode } from "mathjs"
import nerdamer from "nerdamer"
import "nerdamer/Calculus"
import "nerdamer/Algebra"
import "nerdamer/Solve"
import "nerdamer/Extra"

type Num = number
type Vec = number[]
type Mat = number[][]
type Complex = { re: number; im: number }

const math = create(all, {})

const clamp = (x: number, lo: number, hi: number) => Math.max(lo, Math.min(hi, x))
const isFiniteNumber = (x: any) => typeof x === "number" && Number.isFinite(x)
const isVec = (x: any): x is Vec => Array.isArray(x) && x.every(isFiniteNumber)
const isMat = (x: any): x is Mat => Array.isArray(x) && x.every(r => Array.isArray(r) && r.every(isFiniteNumber))
const toFixed = (x: number, d = 10) => Number.parseFloat(x.toFixed(d))

const mdHeader = (s: string) => `# ${s}`
const mdSub = (s: string) => `## ${s}`
const mdCode = (s: string, lang = "text") => `\`\`\`${lang}\n${s}\n\`\`\``
const hr = () => `\n---\n`
const mdList = (items: string[]) => items.map(s => `- ${s}`).join("\n")
const mdKV = (kv: Record<string, any>) => Object.entries(kv).map(([k, v]) => `- ${k}: ${v}`).join("\n")

// -------------------------------------------------------------------------------------
// Expression & Symbolic Helpers
// -------------------------------------------------------------------------------------
function latexOf(expr: string): string {
  try {
    const out = nerdamer(expr)
    return out.toTeX()
  } catch {
    try {
      const node = math.parse(expr)
      return node.toTex()
    } catch {
      return expr
    }
  }
}

function nerdSimplify(expr: string): { expr: string; latex: string } {
  try {
    const e = nerdamer(expr).simplify()
    return { expr: e.text(), latex: e.toTeX() }
  } catch {
    return { expr, latex: latexOf(expr) }
  }
}

function nerdExpand(expr: string): { expr: string; latex: string } {
  try {
    const e = nerdamer(`expand(${expr})`)
    return { expr: e.text(), latex: e.toTeX() }
  } catch {
    return { expr, latex: latexOf(expr) }
  }
}

function nerdFactor(expr: string): { expr: string; latex: string } {
  try {
    const e = nerdamer(`factor(${expr})`)
    return { expr: e.text(), latex: e.toTeX() }
  } catch {
    return { expr, latex: latexOf(expr) }
  }
}

function nerdDifferentiate(expr: string, v = "x"): { expr: string; latex: string } {
  try {
    const e = nerdamer(`diff(${expr}, ${v})`)
    return { expr: e.text(), latex: e.toTeX() }
  } catch {
    return { expr, latex: latexOf(expr) }
  }
}

function nerdIntegrate(expr: string, v = "x"): { expr: string; latex: string } {
  try {
    const e = nerdamer(`integrate(${expr}, ${v})`)
    return { expr: e.text(), latex: e.toTeX() }
  } catch {
    return { expr, latex: latexOf(expr) }
  }
}

function nerdSeries(expr: string, v = "x", about = 0, order = 5): { expr: string; latex: string } {
  try {
    // Nerdamer's series: series(expr, var, about, order)
    const e = nerdamer(`series(${expr}, ${v}, ${about}, ${order})`)
    return { expr: e.text(), latex: e.toTeX() }
  } catch {
    return { expr, latex: latexOf(expr) }
  }
}

function nerdLimit(expr: string, v = "x", to: number | string = "infinity"): { expr: string; latex: string } {
  try {
    const e = nerdamer(`limit(${expr}, ${v}, ${to})`)
    return { expr: e.text(), latex: e.toTeX() }
  } catch {
    return { expr, latex: latexOf(expr) }
  }
}

function nerdSolve(expr: string, v = "x"): { solutions: string[]; latex: string[] } {
  try {
    const e = nerdamer(`solve(${expr}, ${v})`)
    const arr = JSON.parse(e.toString().replace(/^\[|\]$/g, "").split(",").map(s => s.trim()).join(","))
    const sols: string[] = Array.isArray(arr) ? arr.map((x: any) => String(x)) : [String(arr)]
    const latexSols = sols.map(s => latexOf(`${v} = ${s}`))
    return { solutions: sols, latex: latexSols }
  } catch {
    return { solutions: [], latex: [] }
  }
}

function mathEval(expr: string, scope: Record<string, any> = {}): any {
  try {
    return math.evaluate(expr, scope)
  } catch (e: any) {
    throw new Error(`Evaluation error: ${e?.message || e}`)
  }
}

function parseFunction(expr: string, variable = "x"): (x: number) => number {
  const node = math.parse(expr)
  const code = node.compile()
  return (x: number) => {
    try {
      const v = code.evaluate({ [variable]: x, e: Math.E, pi: Math.PI })
      if (typeof v === "number" && Number.isFinite(v)) return v
      if (v?.entries && v.entries.length === 1) return v.entries[0]
      return Number(v)
    } catch {
      return NaN
    }
  }
}

function parseFunction2(expr: string, vars: string[] = ["x", "y"]): (...xs: number[]) => number {
  const node = math.parse(expr)
  const code = node.compile()
  return (...xs: number[]) => {
    const scope: Record<string, number> = {}
    for (let i = 0; i < vars.length; i++) scope[vars[i]] = xs[i]
    try {
      const v = code.evaluate({ ...scope, e: Math.E, pi: Math.PI })
      return Number(v)
    } catch {
      return NaN
    }
  }
}

// -------------------------------------------------------------------------------------
// Numeric Calculus
// -------------------------------------------------------------------------------------
function derivNum(f: (x: number) => number, x: number, h = 1e-6): number {
  return (f(x + h) - f(x - h)) / (2 * h)
}

function secondDerivNum(f: (x: number) => number, x: number, h = 1e-4): number {
  return (f(x + h) - 2 * f(x) + f(x - h)) / (h * h)
}

function integrateTrapezoid(f: (x: number) => number, a: number, b: number, n = 1000): number {
  const h = (b - a) / n
  let s = 0.5 * (f(a) + f(b))
  for (let i = 1; i < n; i++) s += f(a + i * h)
  return s * h
}

function integrateSimpson(f: (x: number) => number, a: number, b: number, n = 1000): number {
  if (n % 2 === 1) n += 1
  const h = (b - a) / n
  let s = f(a) + f(b)
  for (let i = 1; i < n; i++) s += f(a + i * h) * (i % 2 === 0 ? 2 : 4)
  return s * (h / 3)
}

function integrateAdaptiveSimpson(
  f: (x: number) => number,
  a: number,
  b: number,
  eps = 1e-8,
  maxDepth = 20
): number {
  const simpson = (fa: number, fm: number, fb: number, h: number) => (h / 6) * (fa + 4 * fm + fb)
  function recur(a: number, b: number, fa: number, fm: number, fb: number, I: number, depth: number): number {
    const m = 0.5 * (a + b)
    const h = 0.5 * (b - a)
    const f1 = f(0.5 * (a + m))
    const f2 = f(0.5 * (m + b))
    const I1 = simpson(fa, f1, fm, h)
    const I2 = simpson(fm, f2, fb, h)
    if (depth <= 0 || Math.abs(I1 + I2 - I) <= 15 * eps) {
      return I1 + I2 + (I1 + I2 - I) / 15
    }
    return (
      recur(a, m, fa, f1, fm, I1, depth - 1) + recur(m, b, fm, f2, fb, I2, depth - 1)
    )
  }
  const fa = f(a), fb = f(b), m = 0.5 * (a + b), fm = f(m)
  const I0 = simpson(fa, fm, fb, b - a)
  return recur(a, b, fa, fm, fb, I0, maxDepth)
}

// -------------------------------------------------------------------------------------
// Root Finding
// -------------------------------------------------------------------------------------
function bisection(f: (x: number) => number, a: number, b: number, tol = 1e-10, maxIter = 200): number {
  let fa = f(a), fb = f(b)
  if (fa * fb > 0) throw new Error("Bisection requires f(a)*f(b) < 0")
  for (let i = 0; i < maxIter; i++) {
    const m = 0.5 * (a + b)
    const fm = f(m)
    if (Math.abs(fm) < tol || (b - a) / 2 < tol) return m
    if (fa * fm < 0) {
      b = m
      fb = fm
    } else {
      a = m
      fa = fm
    }
  }
  return 0.5 * (a + b)
}

function newton(f: (x: number) => number, x0: number, tol = 1e-12, maxIter = 100, h = 1e-6): number {
  let x = x0
  for (let i = 0; i < maxIter; i++) {
    const fx = f(x)
    const dfx = derivNum(f, x, h)
    if (Math.abs(dfx) < 1e-20) break
    const x1 = x - fx / dfx
    if (Math.abs(x1 - x) < tol) return x1
    x = x1
  }
  return x
}

function secant(f: (x: number) => number, x0: number, x1: number, tol = 1e-12, maxIter = 100): number {
  let a = x0, b = x1
  let fa = f(a), fb = f(b)
  for (let i = 0; i < maxIter; i++) {
    const denom = (fb - fa)
    if (Math.abs(denom) < 1e-30) break
    const c = b - fb * (b - a) / denom
    if (Math.abs(c - b) < tol) return c
    a = b; fa = fb
    b = c; fb = f(b)
  }
  return b
}

// -------------------------------------------------------------------------------------
// Linear Algebra Utilities
// -------------------------------------------------------------------------------------
function matShape(A: Mat): [number, number] {
  return [A.length, A[0]?.length || 0]
}

function matMul(A: Mat, B: Mat): Mat {
  const [m, p] = matShape(A)
  const [p2, n] = matShape(B)
  if (p !== p2) throw new Error("Incompatible shapes in matMul")
  const C: Mat = Array.from({ length: m }, () => Array(n).fill(0))
  for (let i = 0; i < m; i++) {
    for (let k = 0; k < p; k++) {
      const aik = A[i][k]
      for (let j = 0; j < n; j++) {
        C[i][j] += aik * B[k][j]
      }
    }
  }
  return C
}

function matVec(A: Mat, x: Vec): Vec {
  return A.map(r => r.reduce((s, aij, j) => s + aij * x[j], 0))
}

function dot(a: Vec, b: Vec): number {
  let s = 0
  for (let i = 0; i < a.length; i++) s += a[i] * b[i]
  return s
}

function norm2(a: Vec): number {
  return Math.sqrt(dot(a, a))
}

function transpose(A: Mat): Mat {
  const [m, n] = matShape(A)
  const T: Mat = Array.from({ length: n }, () => Array(m).fill(0))
  for (let i = 0; i < m; i++) for (let j = 0; j < n; j++) T[j][i] = A[i][j]
  return T
}

function identity(n: number): Mat {
  const I: Mat = Array.from({ length: n }, () => Array(n).fill(0))
  for (let i = 0; i < n; i++) I[i][i] = 1
  return I
}

function gaussElim(Ain: Mat, bin: Vec): Vec {
  const A = Ain.map(r => r.slice())
  const b = bin.slice()
  const n = A.length
  for (let i = 0; i < n; i++) {
    let piv = i
    for (let r = i + 1; r < n; r++) if (Math.abs(A[r][i]) > Math.abs(A[piv][i])) piv = r
    ;[A[i], A[piv]] = [A[piv], A[i]]
    ;[b[i], b[piv]] = [b[piv], b[i]]
    const pivot = A[i][i]
    if (Math.abs(pivot) < 1e-20) throw new Error("Singular matrix")
    for (let j = i; j < n; j++) A[i][j] /= pivot
    b[i] /= pivot
    for (let r = 0; r < n; r++) {
      if (r === i) continue
      const f = A[r][i]
      for (let j = i; j < n; j++) A[r][j] -= f * A[i][j]
      b[r] -= f * b[i]
    }
  }
  return b
}

function inverse(Ain: Mat): Mat {
  const n = Ain.length
  const A = Ain.map((r, i) => r.concat(identity(n)[i]))
  // Gauss-Jordan
  for (let i = 0; i < n; i++) {
    let piv = i
    for (let r = i + 1; r < n; r++) if (Math.abs(A[r][i]) > Math.abs(A[piv][i])) piv = r
    ;[A[i], A[piv]] = [A[piv], A[i]]
    const pivot = A[i][i]
    if (Math.abs(pivot) < 1e-20) throw new Error("Singular matrix")
    for (let j = 0; j < 2 * n; j++) A[i][j] /= pivot
    for (let r = 0; r < n; r++) {
      if (r === i) continue
      const f = A[r][i]
      for (let j = 0; j < 2 * n; j++) A[r][j] -= f * A[i][j]
    }
  }
  return A.map(r => r.slice(n))
}

function powerIteration(A: Mat, iters = 1000, tol = 1e-12): { lambda: number; v: Vec } {
  const n = A.length
  let v = Array.from({ length: n }, () => Math.random())
  let nv = norm2(v)
  v = v.map(x => x / nv)
  let lambdaOld = 0
  for (let k = 0; k < iters; k++) {
    const Av = matVec(A, v)
    const nv2 = norm2(Av)
    if (nv2 < 1e-20) return { lambda: 0, v }
    const v2 = Av.map(x => x / nv2)
    const lambda = dot(v2, matVec(A, v2))
    if (Math.abs(lambda - lambdaOld) < tol) return { lambda, v: v2 }
    v = v2
    lambdaOld = lambda
  }
  const lambda = dot(v, matVec(A, v))
  return { lambda, v }
}

function qrDecomposition(Ain: Mat): { Q: Mat; R: Mat } {
  // Classical Gram-Schmidt
  const A = Ain.map(r => r.slice())
  const [m, n] = matShape(A)
  const Q: Mat = Array.from({ length: m }, () => Array(n).fill(0))
  const R: Mat = Array.from({ length: n }, () => Array(n).fill(0))
  const aCols: Vec[] = Array.from({ length: n }, (_, j) => A.map(r => r[j]))
  const qCols: Vec[] = []
  for (let j = 0; j < n; j++) {
    let v = aCols[j].slice()
    for (let i = 0; i < j; i++) {
      const qi = qCols[i]
      const rij = dot(qi, aCols[j])
      R[i][j] = rij
      for (let k = 0; k < m; k++) v[k] -= rij * qi[k]
    }
    const rjj = norm2(v) || 1e-20
    R[j][j] = rjj
    const qj = v.map(x => x / rjj)
    qCols.push(qj)
  }
  for (let i = 0; i < m; i++) for (let j = 0; j < n; j++) Q[i][j] = qCols[j][i]
  return { Q, R }
}

function qrEigen(Ain: Mat, iters = 100, tol = 1e-12): { eigenvalues: number[]; Ak: Mat } {
  let A = Ain.map(r => r.slice())
  const n = A.length
  for (let k = 0; k < iters; k++) {
    const { Q, R } = qrDecomposition(A)
    A = matMul(R, Q)
  }
  const eigenvalues = Array.from({ length: n }, (_, i) => A[i][i])
  return { eigenvalues, Ak: A }
}

function rank(Ain: Mat): number {
  const A = Ain.map(r => r.slice())
  const [m, n] = matShape(A)
  let r = 0
  let lead = 0
  for (let i = 0; i < m; i++) {
    if (n <= lead) break
    let j = i
    while (Math.abs(A[j][lead]) < 1e-12) {
      j++
      if (j === m) {
        j = i
        lead++
        if (n === lead) return r
      }
    }
    ;[A[i], A[j]] = [A[j], A[i]]
    const div = A[i][lead]
    for (let k = 0; k < n; k++) A[i][k] /= div
    for (let l = 0; l < m; l++) {
      if (l === i) continue
      const mult = A[l][lead]
      for (let k = 0; k < n; k++) A[l][k] -= mult * A[i][k]
    }
    lead++
    r++
  }
  return r
}

function parseMatrix(input: string): Mat {
  // Accepts either JSON-like [[...],[...]] or rows separated by ; and spaces/commas
  let text = input.trim()
  if (text.startsWith("[[")) {
    try {
      const A = JSON.parse(text)
      if (!isMat(A)) throw new Error("Invalid matrix")
      return A
    } catch {
      // fallthrough
    }
  }
  // Row separator ; and numbers separated by spaces or commas
  const rows = text.split(";").map(r => r.trim()).filter(Boolean)
  const A: Mat = rows.map(r => r.split(/[,\s]+/).filter(Boolean).map(Number))
  if (!isMat(A)) throw new Error("Invalid matrix format")
  const n = A[0].length
  if (!A.every(r => r.length === n)) throw new Error("Jagged matrix not allowed")
  return A
}

function parseVector(input: string): Vec {
  let text = input.trim()
  if (text.startsWith("[")) {
    try {
      const v = JSON.parse(text)
      if (!isVec(v)) throw new Error("Invalid vector")
      return v
    } catch {
      // fallthrough
    }
  }
  const v = text.split(/[,\s]+/).filter(Boolean).map(Number)
  if (!isVec(v)) throw new Error("Invalid vector")
  return v
}

// -------------------------------------------------------------------------------------
// Optimization
// -------------------------------------------------------------------------------------
function gradientDescent(
  f: (x: Vec) => number,
  x0: Vec,
  alpha = 1e-2,
  iters = 1000,
  tol = 1e-10
): { x: Vec; fx: number; k: number } {
  const x = x0.slice()
  const n = x.length
  function gradAt(x: Vec): Vec {
    const g = Array(n).fill(0)
    const h = 1e-6
    const fx = f(x)
    for (let i = 0; i < n; i++) {
      const old = x[i]
      x[i] = old + h
      const f1 = f(x)
      g[i] = (f1 - fx) / h
      x[i] = old
    }
    return g
  }
  let fx = f(x)
  for (let k = 0; k < iters; k++) {
    const g = gradAt(x)
    const ng = norm2(g)
    if (ng < tol) return { x: x.slice(), fx, k }
    for (let i = 0; i < n; i++) x[i] -= alpha * g[i]
    fx = f(x)
  }
  return { x: x.slice(), fx, k: iters }
}

function conjugateGradient(A: Mat, b: Vec, tol = 1e-12, maxIter = 1000): { x: Vec; k: number } {
  const n = b.length
  let x = Array(n).fill(0)
  let r = b.map((bi, i) => bi - dot(A[i], x))
  let p = r.slice()
  let rsold = dot(r, r)
  for (let k = 0; k < maxIter; k++) {
    const Ap = matVec(A, p)
    const alpha = rsold / dot(p, Ap)
    x = x.map((xi, i) => xi + alpha * p[i])
    r = r.map((ri, i) => ri - alpha * Ap[i])
    const rsnew = dot(r, r)
    if (Math.sqrt(rsnew) < tol) return { x, k }
    const beta = rsnew / rsold
    p = r.map((ri, i) => ri + beta * p[i])
    rsold = rsnew
  }
  return { x, k: maxIter }
}

// -------------------------------------------------------------------------------------
// Probability & Statistics
// -------------------------------------------------------------------------------------
function mean(v: Vec): number {
  return v.reduce((s, x) => s + x, 0) / v.length
}
function variance(v: Vec, sample = true): number {
  const m = mean(v)
  const denom = sample ? v.length - 1 : v.length
  return v.reduce((s, x) => s + (x - m) ** 2, 0) / denom
}
function stddev(v: Vec, sample = true): number {
  return Math.sqrt(variance(v, sample))
}
function covariance(x: Vec, y: Vec, sample = true): number {
  if (x.length !== y.length) throw new Error("x and y length mismatch")
  const mx = mean(x), my = mean(y)
  const denom = sample ? x.length - 1 : x.length
  let s = 0
  for (let i = 0; i < x.length; i++) s += (x[i] - mx) * (y[i] - my)
  return s / denom
}
function corrcoef(x: Vec, y: Vec): number {
  return covariance(x, y, true) / (stddev(x, true) * stddev(y, true))
}
function linearRegression(x: Vec, y: Vec): { slope: number; intercept: number; r: number } {
  const sx = mean(x), sy = mean(y)
  let sxy = 0, sxx = 0, syy = 0
  for (let i = 0; i < x.length; i++) {
    sxy += (x[i] - sx) * (y[i] - sy)
    sxx += (x[i] - sx) * (x[i] - sx)
    syy += (y[i] - sy) * (y[i] - sy)
  }
  const slope = sxy / sxx
  const intercept = sy - slope * sx
  const r = sxy / Math.sqrt(sxx * syy)
  return { slope, intercept, r }
}
function covarianceMatrix(X: Mat, sample = true): Mat {
  const [n, d] = matShape(X)
  const means = Array(d).fill(0)
  for (let j = 0; j < d; j++) means[j] = mean(X.map(r => r[j]))
  const cov: Mat = Array.from({ length: d }, () => Array(d).fill(0))
  const denom = sample ? n - 1 : n
  for (let i = 0; i < d; i++) {
    for (let j = i; j < d; j++) {
      let s = 0
      for (let k = 0; k < n; k++) s += (X[k][i] - means[i]) * (X[k][j] - means[j])
      const v = s / denom
      cov[i][j] = v
      cov[j][i] = v
    }
  }
  return cov
}

// -------------------------------------------------------------------------------------
// Number Theory
// -------------------------------------------------------------------------------------
function gcd(a: number, b: number): number {
  a = Math.abs(a); b = Math.abs(b)
  while (b !== 0) {
    const t = b
    b = a % b
    a = t
  }
  return a
}
function egcd(a: number, b: number): { g: number; x: number; y: number } {
  if (b === 0) return { g: a, x: 1, y: 0 }
  const { g, x, y } = egcd(b, a % b)
  return { g, x: y, y: x - Math.floor(a / b) * y }
}
function modInv(a: number, m: number): number {
  const { g, x } = egcd(a, m)
  if (g !== 1 && g !== -1) throw new Error("Inverse does not exist")
  return ((x % m) + m) % m
}
function sieve(n: number): number[] {
  const isPrime = Array(n + 1).fill(true)
  isPrime[0] = isPrime[1] = false
  for (let p = 2; p * p <= n; p++) if (isPrime[p]) for (let k = p * p; k <= n; k += p) isPrime[k] = false
  return isPrime.map((v, i) => (v ? i : 0)).filter(Boolean)
}

// -------------------------------------------------------------------------------------
// Polynomials
// -------------------------------------------------------------------------------------
function polyEval(coeffs: Vec, x: number): number {
  // coeffs: [a0, a1, ..., an] for a0 + a1 x + ...
  let y = 0
  for (let i = coeffs.length - 1; i >= 0; i--) y = y * x + coeffs[i]
  return y
}

function polyDeriv(coeffs: Vec): Vec {
  const d: Vec = []
  for (let i = 1; i < coeffs.length; i++) d.push(i * coeffs[i])
  return d
}

function syntheticDivision(coeffs: Vec, r: number): { q: Vec; rem: number } {
  const n = coeffs.length
  const q: Vec = Array(n - 1).fill(0)
  let b = coeffs[n - 1]
  for (let i = n - 2; i >= 0; i--) {
    q[i] = b
    b = coeffs[i] + r * b
  }
  const rem = b
  return { q: q.reverse(), rem }
}

function newtonPolyRoot(coeffs: Vec, x0: number, tol = 1e-12, maxIter = 100): number {
  const f = (x: number) => polyEval(coeffs, x)
  const df = (x: number) => polyEval(polyDeriv(coeffs), x)
  let x = x0
  for (let k = 0; k < maxIter; k++) {
    const fx = f(x), dfx = df(x)
    if (Math.abs(dfx) < 1e-20) break
    const x1 = x - fx / dfx
    if (Math.abs(x1 - x) < tol) return x1
    x = x1
  }
  return x
}

// -------------------------------------------------------------------------------------
// Differential Equations (ODE)
// -------------------------------------------------------------------------------------
function rk4(
  f: (t: number, y: number) => number,
  t0: number,
  y0: number,
  t1: number,
  n = 100
): { t: number[]; y: number[] } {
  const h = (t1 - t0) / n
  const t: number[] = [t0]
  const y: number[] = [y0]
  let ti = t0, yi = y0
  for (let i = 0; i < n; i++) {
    const k1 = f(ti, yi)
    const k2 = f(ti + 0.5 * h, yi + 0.5 * h * k1)
    const k3 = f(ti + 0.5 * h, yi + 0.5 * h * k2)
    const k4 = f(ti + h, yi + h * k3)
    yi = yi + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
    ti = ti + h
    t.push(ti); y.push(yi)
  }
  return { t, y }
}

function rk4System(
  f: (t: number, y: Vec) => Vec,
  t0: number,
  y0: Vec,
  t1: number,
  n = 100
): { t: number[]; y: Vec[] } {
  const h = (t1 - t0) / n
  const t: number[] = [t0]
  const y: Vec[] = [y0.slice()]
  let ti = t0, yi = y0.slice()
  for (let i = 0; i < n; i++) {
    const k1 = f(ti, yi)
    const k2 = f(ti + 0.5 * h, yi.map((yij, j) => yij + 0.5 * h * k1[j]))
    const k3 = f(ti + 0.5 * h, yi.map((yij, j) => yij + 0.5 * h * k2[j]))
    const k4 = f(ti + h, yi.map((yij, j) => yij + h * k3[j]))
    yi = yi.map((yij, j) => yij + (h / 6) * (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]))
    ti = ti + h
    t.push(ti); y.push(yi.slice())
  }
  return { t, y }
}

// -------------------------------------------------------------------------------------
// Transforms (DFT/FFT-lite numeric)
// -------------------------------------------------------------------------------------
function dft(x: Complex[] | number[]): Complex[] {
  const N = x.length
  const out: Complex[] = Array.from({ length: N }, () => ({ re: 0, im: 0 }))
  for (let k = 0; k < N; k++) {
    let sumRe = 0, sumIm = 0
    for (let n = 0; n < N; n++) {
      const xn = typeof x[n] === "number" ? (x[n] as number) : (x[n] as Complex).re
      const angle = (-2 * Math.PI * k * n) / N
      sumRe += xn * Math.cos(angle)
      sumIm += xn * Math.sin(angle)
    }
    out[k] = { re: sumRe, im: sumIm }
  }
  return out
}

function idft(X: Complex[]): number[] {
  const N = X.length
  const out: number[] = Array.from({ length: N }, () => 0)
  for (let n = 0; n < N; n++) {
    let sumRe = 0, sumIm = 0
    for (let k = 0; k < N; k++) {
      const angle = (2 * Math.PI * k * n) / N
      sumRe += X[k].re * Math.cos(angle) - X[k].im * Math.sin(angle)
      sumIm += X[k].re * Math.sin(angle) + X[k].im * Math.cos(angle)
    }
    out[n] = (sumRe / N)
  }
  return out
}

// -------------------------------------------------------------------------------------
// Expression Analysis
// -------------------------------------------------------------------------------------
function expressionTreeDepth(expr: string): number {
  try {
    const node = math.parse(expr)
    function depth(n: MathNode): number {
      if (!("args" in n) || !(n as any).args) return 1
      const args: MathNode[] = (n as any).args || []
      if (!args.length) return 1
      return 1 + Math.max(...args.map(depth))
    }
    return depth(node)
  } catch {
    return 0
  }
}

function countSymbols(expr: string): number {
  try {
    const node = math.parse(expr)
    let count = 0
    node.traverse(() => count++)
    return count
  } catch {
    return 0
  }
}

// -------------------------------------------------------------------------------------
// Rendering Helpers
// -------------------------------------------------------------------------------------
function section(title: string, body: string) {
  return `${mdSub(title)}\n\n${body}\n`
}

function fmtVec(v: Vec, digits = 6) {
  return `[${v.map(x => toFixed(x, digits)).join(", ")}]`
}

function fmtMat(A: Mat, digits = 6) {
  return "[\n" + A.map(r => "  " + fmtVec(r, digits)).join(",\n") + "\n]"
}

function complexToString(z: Complex, digits = 6) {
  const re = toFixed(z.re, digits)
  const im = toFixed(z.im, digits)
  return `${re} ${im >= 0 ? "+" : "-"} ${Math.abs(im)}i`
}

// -------------------------------------------------------------------------------------
// UI Flows
// -------------------------------------------------------------------------------------
async function showResult(title: string, content: string, extras?: { latex?: string; copy?: string }) {
  const latex = extras?.latex
  const copyText = extras?.copy
  const html = md(
    [
      mdHeader(title),
      content,
      latex ? section("LaTeX", mdCode(latex)) : "",
    ].join("\n\n")
  )
  await div(
    {
      html,
      enter: copyText ? "Copy Result" : "Close",
      shortcuts: [
        ...(copyText
          ? [
              {
                name: "Copy",
                key: `${cmd}+c`,
                bar: "right",
                onPress: async () => {
                  await copy(copyText)
                  await toast("Copied")
                },
              },
            ]
          : []),
      ],
    },
    [
      ...(copyText
        ? [
            {
              name: "Copy Result",
              shortcut: `${cmd}+c`,
              onAction: async () => {
                await copy(copyText)
                await toast("Copied")
              },
            },
          ]
        : []),
    ]
  )
}

async function symbolicMenu() {
  const task = await arg("Symbolic Algebra", [
    "Simplify",
    "Expand",
    "Factor",
    "Differentiate",
    "Integrate",
    "Series",
    "Limit",
    "Solve Equation",
  ])
  if (task === "Simplify") {
    const expr = await arg("Enter expression to simplify (e.g., (x^2 - 1)/(x-1))")
    const out = nerdSimplify(expr)
    const body = [section("Input", mdCode(expr)), section("Simplified", mdCode(out.expr))].join("\n")
    await showResult("Symbolic Simplify", body, { latex: out.latex, copy: out.expr })
    return
  }
  if (task === "Expand") {
    const expr = await arg("Enter expression to expand (e.g., (x+1)^5)")
    const out = nerdExpand(expr)
    const body = [section("Input", mdCode(expr)), section("Expanded", mdCode(out.expr))].join("\n")
    await showResult("Symbolic Expand", body, { latex: out.latex, copy: out.expr })
    return
  }
  if (task === "Factor") {
    const expr = await arg("Enter expression to factor (e.g., x^2-1)")
    const out = nerdFactor(expr)
    const body = [section("Input", mdCode(expr)), section("Factored", mdCode(out.expr))].join("\n")
    await showResult("Symbolic Factor", body, { latex: out.latex, copy: out.expr })
    return
  }
  if (task === "Differentiate") {
    const expr = await arg("Enter expression (e.g., sin(x)*e^x)")
    const v = await arg("Differentiate w.r.t variable", ["x", "y", "t"])
    const out = nerdDifferentiate(expr, v)
    const body = [section("Input", mdCode(expr)), section("d/d" + v, mdCode(out.expr))].join("\n")
    await showResult("Symbolic Derivative", body, { latex: out.latex, copy: out.expr })
    return
  }
  if (task === "Integrate") {
    const expr = await arg("Enter expression (e.g., x^2 * sin(x))")
    const v = await arg("Integrate w.r.t variable", ["x", "y", "t"])
    const out = nerdIntegrate(expr, v)
    const body = [section("Input", mdCode(expr)), section("∫ d" + v, mdCode(out.expr))].join("\n")
    await showResult("Symbolic Integral", body, { latex: out.latex, copy: out.expr })
    return
  }
  if (task === "Series") {
    const expr = await arg("Expression (e.g., exp(x))")
    const v = await arg("Variable", ["x", "y", "t"])
    const about = await arg("About (e.g., 0)")
    const order = await arg("Order (e.g., 6)")
    const out = nerdSeries(expr, v, Number(about), Number(order))
    const body = [
      section("Input", mdCode(expr)),
      section(`Series about ${v}=${about} to order ${order}`, mdCode(out.expr)),
    ].join("\n")
    await showResult("Series Expansion", body, { latex: out.latex, copy: out.expr })
    return
  }
  if (task === "Limit") {
    const expr = await arg("Expression (e.g., sin(x)/x)")
    const v = await arg("Variable", ["x", "y", "t"])
    const to = await arg("Approach to (use 'infinity' or a number)", ["infinity", "-infinity", "0", "1"])
    const toVal = to === "infinity" || to === "-infinity" ? to : Number(to)
    const out = nerdLimit(expr, v, toVal)
    const body = [
      section("Input", mdCode(expr)),
      section(`Limit as ${v}->${to}`, mdCode(out.expr)),
    ].join("\n")
    await showResult("Limit", body, { latex: out.latex, copy: out.expr })
    return
  }
  if (task === "Solve Equation") {
    const expr = await arg("Equation (e.g., x^3 - 2x - 5 = 0)")
    const v = await arg("Solve for var", ["x", "y", "t"])
    const normalized = expr.includes("=") ? `solve(${expr.replace("=", "==")}, ${v})` : `solve(${expr}=0, ${v})`
    const { solutions, latex } = nerdSolve(expr.includes("=") ? expr : `${expr}=0`, v)
    const body = [
      section("Input", mdCode(expr)),
      section("Call", mdCode(normalized)),
      section("Solutions", solutions.length ? mdList(solutions) : "No solutions found"),
    ].join("\n")
    await showResult("Solve Equation", body, {
      latex: latex.join("\n"),
      copy: solutions.join(", "),
    })
    return
  }
}

async function calculusMenu() {
  const task = await arg("Calculus (Numeric + Symbolic)", [
    "Symbolic Derivative",
    "Symbolic Integral",
    "Numeric Derivative",
    "Numeric Integral (Trapezoid/Simpson)",
    "Adaptive Integral",
    "Taylor Series (Symbolic)",
  ])
  if (task === "Symbolic Derivative") {
    const expr = await arg("Expression (e.g., tan(x)*ln(x))")
    const v = await arg("Variable", ["x", "y", "t"])
    const out = nerdDifferentiate(expr, v)
    const body = [section("Input", mdCode(expr)), section("Result", mdCode(out.expr))].join("\n")
    await showResult("Symbolic Derivative", body, { latex: out.latex, copy: out.expr })
    return
  }
  if (task === "Symbolic Integral") {
    const expr = await arg("Expression (e.g., (x^2+1)/(x))")
    const v = await arg("Variable", ["x", "y", "t"])
    const out = nerdIntegrate(expr, v)
    const body = [section("Input", mdCode(expr)), section("Result", mdCode(out.expr))].join("\n")
    await showResult("Symbolic Integral", body, { latex: out.latex, copy: out.expr })
    return
  }
  if (task === "Numeric Derivative") {
    const expr = await arg("f(x) (e.g