数值增量与证书实验源码
研究总览 · 2026-09-20
附录 B：原始隔离探针，非产品执行器。
附录 B：数值实验的完整复现源码
Section titled “附录 B：数值实验的完整复现源码”以下是 §12.3–§12.4 实际执行的隔离 Node 探针，仅重新排版；原始 UTF-8 文件 SHA-256 为 531b0b5730df17655f24f295ce1b4a751967266745e3c73ee35cd75b6cb0bdcc。保存为临时 .cjs 后以 Node 运行；源码中的 folder 改为相同锁定包的本机绝对路径，不安装或修改任何产品依赖。默认复现 240/480 组；把 n=240,m=480 改为 n=24,m=48 复现小规模组。
这是二进制能力和机制比较，包含内存中的研究专用绑定注入，不是生产桥或独立建模执行器。负环中对 ray 的整体符号处理仅适用于这个已知结构；一般产品采用 §6 的双边行映射和精确核验。输出中的 maxViolation 是 binary64 数值检查，不能当任意精确可行性证明。presolve=off 组从 choose 组结束的模型继续同样的修改规则；组内三种模式输入一致，不声称两种 presolve 组逐样本输入完全相同。
const fs = require("node:fs"), Module = require("node:module"), path = require("node:path"), os = require("node:os"), crypto = require("node:crypto");const folder = "<workspace-root>/packages/aira-geometry-kernel/node_modules/@casadi/casadi-wasm";function compile(name, edit, override) { const fn = path.join(folder, name), m = new Module(fn, module); m.filename = fn; m.paths = Module._nodeModulePaths(folder); if (override) m.require = override; m._compile(edit(fs.readFileSync(fn, "utf8")), fn); return m.exports;}const wf = compile("casadi_wasm.js", (s) => s.replace(";return moduleRtn", ";Module.__research={LDSO};return moduleRtn"),);const cf = compile( "casadi.js", (s) => s.replace(" return __m;", " __m.__researchModule=M;return __m;"), (n) => (n === "./casadi_wasm.js" ? wf : require(n)),);const median = (a) => [...a].sort((x, y) => x - y)[Math.floor(a.length / 2)];const stats = (a) => ({ n: a.length, median_ms: median(a), p95_ms: [...a].sort((x, y) => x - y)[Math.min(a.length - 1, Math.ceil(a.length * 0.95) - 1)], min_ms: Math.min(...a), max_ms: Math.max(...a),});(async () => { const ca = await cf(); await ca.load_conic("highs"); const M = ca.__researchModule, H = M.__research.LDSO.loadedLibsByName["libcasadi_conic_highs.so"].exports; const ptrs = new Set(); function alloc(arr, T = Float64Array) { const p = M._malloc(arr.length * T.BYTES_PER_ELEMENT); ptrs.add(p); new T(M.HEAPF64.buffer, p, arr.length).set(arr); return p; } const f = (a) => alloc(a), i = (a) => alloc(a, Int32Array); function str(s) { const bytes = Buffer.from(s + "\0"); return alloc(bytes, Uint8Array); } const si = str("simplex_iteration_count"), bv = str("basis_validity"), o = str("output_flag"), sol = str("solver"), simplex = str("simplex"), pre = str("presolve"), off = str("off"); const inf = 1e30; const check = (s, label) => { if (s !== 0) throw new Error(label + ": " + s); }; const readi = (p) => new Int32Array(M.HEAPF64.buffer, p, 1)[0]; const info = (Hh) => { const v = i([0]); check(H.Highs_getIntInfoValue(Hh, si, v), "info"); const it = readi(v); check(H.Highs_getIntInfoValue(Hh, bv, v), "basis"); const basis = readi(v); M._free(v); ptrs.delete(v); return { it, basis }; }; function create(presolveOff) { const h = H.Highs_create(); check(H.Highs_setBoolOptionValue(h, o, 0), "quiet"); check(H.Highs_setStringOptionValue(h, sol, simplex), "simplex"); if (presolveOff) check(H.Highs_setStringOptionValue(h, pre, off), "presolve"); return h; } function matrix(n, rows, c, cl, cu) { const m = rows.length, cols = Array.from({ length: n }, () => []); rows.forEach((r, k) => r.v.forEach(([j, a]) => cols[j].push([k, a]))); const st = [], ix = [], va = []; cols.forEach((co) => { st.push(ix.length); co.forEach(([k, a]) => { ix.push(k); va.push(a); }); }); st.push(ix.length); return { n, m, nz: ix.length, c, cl, cu, rows, pc: f(c), pl: f(cl), pu: f(cu), prl: f(rows.map((r) => r.l ?? -inf)), pru: f(rows.map((r) => r.u ?? inf)), ps: i(st), pi: i(ix), pv: f(va), }; } function pass(h, p) { check( H.Highs_passLp(h, p.n, p.m, p.nz, 1, 1, 0, p.pc, p.pl, p.pu, p.prl, p.pru, p.ps, p.pi, p.pv), "passLp", ); } function run(h) { check(H.Highs_run(h), "run"); return H.Highs_getModelStatus(h); } function solution(h, n) { const p = f(Array(n).fill(0)); check(H.Highs_getSolution(h, p, 0, 0, 0), "solution"); const x = Array.from(new Float64Array(M.HEAPF64.buffer, p, n)); M._free(p); ptrs.delete(p); return x; } function valid(p, x) { let violation = 0; for (let j = 0; j < p.n; j++) violation = Math.max(violation, p.cl[j] - x[j], x[j] - p.cu[j]); for (const r of p.rows) { const v = r.v.reduce((a, [j, c]) => a + c * x[j], 0); violation = Math.max(violation, (r.l ?? -inf) - v, v - (r.u ?? inf)); } return { maxViolation: violation, objective: x.reduce((a, v, j) => a + v * p.c[j], 0) }; } let seed = 1831565813; function rand() { seed ^= seed << 13; seed ^= seed >>> 17; seed ^= seed << 5; return (seed >>> 0) / 4294967296; } const n = 240, m = 480; const rows = Array.from({ length: m }, () => { const seen = new Set(), v = []; while (v.length < 6) { const j = Math.floor(rand() * n); if (seen.has(j)) continue; seen.add(j); v.push([j, 0.2 + rand() * 1.8]); } return { v, u: v.reduce((a, [, c]) => a + c * 0.5, 0) + 0.1 + rand() * 0.4 }; }); const cost = Array.from({ length: n }, () => -0.5 - rand()), p = matrix(n, rows, cost, Array(n).fill(0), Array(n).fill(1)); const report = { environment: { node: process.version, platform: process.platform, arch: process.arch, cpu: os.cpus()[0].model, coreVersion: M.UTF8ToString(H.Highs_version()), coreGit: M.UTF8ToString(H.Highs_githash()), pluginSHA256: crypto .createHash("sha256") .update(fs.readFileSync(path.join(folder, "libcasadi_conic_highs.so"))) .digest("hex"), seed: 1831565813, }, problem: { n, m, nnz: p.nz, repetitions: 35, warmup: 5, edit: "one row upper bound and one column cost per solve", }, }; report.lp = []; for (const presolveOff of [false, true]) { const hs = { rebuild: create(presolveOff), incremental: create(presolveOff), rebuild_setBasis: create(presolveOff), }; for (const h of Object.values(hs)) { pass(h, p); if (run(h) !== 7) throw new Error("initial solve"); } const cb = i(Array(n).fill(0)), rb = i(Array(m).fill(0)); const times = { rebuild: [], incremental: [], rebuild_setBasis: [] }, iters = { rebuild: [], incremental: [], rebuild_setBasis: [] }; let maxDiff = 0, maxViolation = 0; for (let k = 0; k < 40; k++) { const row = (k * 17) % m, col = (k * 13) % n; p.rows[row].u += k % 2 ? 0.015 : -0.015; p.c[col] += k % 2 ? 0.003 : -0.003; new Float64Array(M.HEAPF64.buffer, p.pru, m)[row] = p.rows[row].u; new Float64Array(M.HEAPF64.buffer, p.pc, n)[col] = p.c[col]; const outcomes = {}; const modes = Object.keys(hs); for (const mode of k % 2 ? modes : [...modes].reverse()) { const h = hs[mode]; const t = performance.now(); if (mode === "rebuild_setBasis") { check(H.Highs_getBasis(h, cb, rb), "getBasis"); pass(h, p); check(H.Highs_setBasis(h, cb, rb), "setBasis"); } else if (mode === "rebuild") pass(h, p); else { check(H.Highs_changeRowBounds(h, row, -inf, p.rows[row].u), "changeRow"); check(H.Highs_changeColCost(h, col, p.c[col]), "changeCost"); } if (run(h) !== 7) throw new Error("solve not optimal"); const ms = performance.now() - t, ii = info(h), vv = valid(p, solution(h, n)); if (ii.basis !== 1) throw new Error("invalid basis"); maxViolation = Math.max(maxViolation, vv.maxViolation); outcomes[mode] = vv.objective; if (k >= 5) { times[mode].push(ms); iters[mode].push(ii.it); } } maxDiff = Math.max( maxDiff, Math.abs(outcomes.rebuild - outcomes.incremental), Math.abs(outcomes.rebuild - outcomes.rebuild_setBasis), ); } report.lp.push({ presolve: presolveOff ? "off" : "choose", maxObjectiveDifference: maxDiff, maxViolation, modes: Object.fromEntries( Object.keys(hs).map((mode) => [ mode, { ...stats(times[mode]), medianIterations: median(iters[mode]), maxIterations: Math.max(...iters[mode]), }, ]), ), }); Object.values(hs).forEach((h) => H.Highs_destroy(h)); } const rn = 120, rr = Array.from({ length: rn }, (_, j) => ({ v: [ [j, -1], [(j + 1) % rn, 1], ], u: -1, })), rp = matrix(rn, rr, Array(rn).fill(0), Array(rn).fill(-inf), Array(rn).fill(inf)); const ar = Array.from({ length: rn }, (_, j) => ({ v: [ [j, -1], [(j + rn - 1) % rn, 1], ], l: 0, u: 0, })); ar.push({ v: Array.from({ length: rn }, (_, j) => [j, 1]), l: 1, u: 1 }); const ap = matrix(rn, ar, Array(rn).fill(-1), Array(rn).fill(0), Array(rn).fill(inf)); const has = i([0]), ray = f(Array(rn).fill(0)), rt = [], at = []; let rayMin = inf, rayDot = 0, rayResidual = 0, auxMin = inf, auxDot = 0, auxResidual = 0; for (let k = 0; k < 40; k++) { const h = create(false); pass(h, rp); if (run(h) !== 8) throw new Error("cycle should infeasible"); const t = performance.now(); check(H.Highs_getDualRay(h, has, ray), "getDualRay"); const ms = performance.now() - t; if (!readi(has)) throw new Error("no ray"); let y = Array.from(new Float64Array(M.HEAPF64.buffer, ray, rn)); if (y.reduce((a, b) => a + b, 0) < 0) y = y.map((v) => -v); rayMin = Math.min(...y); rayDot = -y.reduce((a, b) => a + b, 0); rayResidual = Math.max(...y.map((v, j) => Math.abs(y[(j + rn - 1) % rn] - v))); H.Highs_destroy(h); const ta = performance.now(), ah = create(false); pass(ah, ap); if (run(ah) !== 7) throw new Error("aux fail"); const yy = solution(ah, rn); H.Highs_destroy(ah); const ams = performance.now() - ta; auxMin = Math.min(...yy); auxDot = -yy.reduce((a, b) => a + b, 0); auxResidual = Math.max(...yy.map((v, j) => Math.abs(yy[(j + rn - 1) % rn] - v))); if (k >= 5) { rt.push(ms); at.push(ams); } } report.certificate = { model: "120-row negative difference cycle, b=-1, free variables", ray: stats(rt), auxiliaryLP: stats(at), rayCheck: { min: rayMin, bDot: rayDot, maxAtResidual: rayResidual }, auxCheck: { min: auxMin, bDot: auxDot, maxAtResidual: auxResidual }, timingScope: "ray after primary infeasibility solve; auxiliary includes create/pass/solve/extract/destroy, prebuilt JS/wasm matrices excluded", }; console.log(JSON.stringify(report, null, 2)); for (const p of ptrs) M._free(p);})();