From 49f2f4b1fcc86fe7808cc9376adca85773c34153 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 18 Oct 2026 23:07:58 +0700 Subject: [PATCH] analog: mos_table.t27 -- CMOS inverter with a tabulated square-law MOS, backward Euler and fixed Newton, against cpphdl (Closes #10515) M3a of the analog emulation epic (#10486). The device is a 25 x 25 Id(Vgs, Vds) table on a 0.05 V grid (Vdd and Vth are grid lines), read by exact-integer bilinear interpolation; gm and gds are the slopes of the same interpolant. The output node is integrated by backward Euler with K = 3 Newton iterations, fixed. One 2^20-unit cell is 0.05 V, so cell index, fraction and the interpolation divisor are shifts; the only division is Newton's. Tests (13, 0 vacuous): table against the law at the nodes and at all 576 cell midpoints; the fixed-point inverter run against cpphdl's requirement; the three errors kept apart and bounded from both sides (arithmetic 0.61 unit, table 0.428 mV, discretization 27.2 mV); the last residual at most 1 unit at every step; negative controls K = 1 and Vth = 0.30; bit-exact mirror symmetry. fold and fabs64 come from analog/linear (M1) by use; no board top, no PDK. Co-Authored-By: Claude Sonnet 5.5 --- .trinity/seals/analog_AnalogMosTable.json | 28 + specs/analog/mos_table.t27 | 899 ++++++++++++++++++++++ 2 files changed, 927 insertions(+) create mode 100644 .trinity/seals/analog_AnalogMosTable.json create mode 100644 specs/analog/mos_table.t27 diff --git a/.trinity/seals/analog_AnalogMosTable.json b/.trinity/seals/analog_AnalogMosTable.json new file mode 100644 index 0000000000..16550c7676 --- /dev/null +++ b/.trinity/seals/analog_AnalogMosTable.json @@ -0,0 +1,28 @@ +{ + "built_by": "t27c-bootstrap@0.5.2+9b556b7c3", + "closure": [ + { + "spec": "specs/analog/linear.t27", + "spec_hash": "sha256:3ddef7be183aba583ed285334ee389b48a65b41335e98914514235db7c4e5d62" + } + ], + "config": "gen=c,rust,verilog,zig;tests=zig", + "gen_hash_c": "sha256:a4c0ec9ec8298af635a1cb7363b0b9644b99f7f0a1b80bd709fb8b14184a82fb", + "gen_hash_rust": "sha256:d0cbeb4b09d5da7c886eaeefaa2c395a2b699b643a1b25bb7372c62d857f0d87", + "gen_hash_verilog": "sha256:97d7e477baee00d846fc850c5888d0fbbf7cf9c072124ee4f55e4b76eeffc851", + "gen_hash_zig": "sha256:d9752320d155e8ace406d0008627ea48410160f8b25eea9eb424b9e10eb4d540", + "module": "AnalogMosTable", + "ring": 12, + "seal_schema": 2, + "sealed_at": "2026-10-18T16:07:33Z", + "sealed_by": "t27c-bootstrap@0.5.2", + "spec_hash": "sha256:168cd93374cec92ece02ddab1be51f49308d9efd3288ac0427f912d4a7521782", + "spec_path": "specs/analog/mos_table.t27", + "tests": { + "blocked": "does not compile: error: failed to create output directory '/Users/playom/.cache/zig/tmp/efdc11fc0a8df6c4': PermissionDenied" + }, + "toolchain": { + "test_runner": "sha256:1982ae948e24d00c5d221ccec12dfce1eb01204b53e37b0f857f93b16d64629d", + "zig": "0.16.0" + } +} \ No newline at end of file diff --git a/specs/analog/mos_table.t27 b/specs/analog/mos_table.t27 new file mode 100644 index 0000000000..a9f96e0c46 --- /dev/null +++ b/specs/analog/mos_table.t27 @@ -0,0 +1,899 @@ +// SPDX-License-Identifier: Apache-2.0 +; specs/analog/mos_table.t27 -- a nonlinear analog cell: a CMOS inverter whose MOS currents are read from a table +; M3a of the analog emulation epic (#10486, issue #10515). M0 (specs/fpga/rc_filter_ax7203.t27) is the linear RC; this is +; the first cell with a device that is not a resistor. M3b swaps the analytic table below for tables made from the IHP +; SG13G2 PSP103 model, per corner. Everything hand-written is this file. +; - +; THE CIRCUIT. The reference is mirekez/cpphdl, analog/examples/cmos_inverter: Vdd = 1.2 V, an NMOS from Y to ground and a +; PMOS from Y to Vdd, both gates on A, A = 0 V and then 1.2 V from 2 us, Y loaded by 10 pF, dt = 10 ns, stop at 4 us (400 +; steps). cpphdl's MOS law (analog/include/solver.hpp) is the square law with beta = 1 mA/V^2 and Vth = 0.35 V: +; ov = Vgs - Vth; below threshold I = 0; Vds < ov: I = beta (ov Vds - Vds^2 / 2); else I = beta ov^2 / 2. +; The law is symmetric in drain and source and the PMOS is the mirror (Vsg, Vsd for Vgs, Vds); here Vds is never negative, +; so the symmetry is not exercised. cpphdl's README reports Y rising toward 1.2 V, then falling to about 11.5 nV at 4 us, +; and requires |Y(4 us)| <= 1 mV. +; - +; THE DEVICE IS A TABLE. Id(Vgs, Vds) is stored on a uniform 25 x 25 grid over [0, 1.2] V x [0, 1.2] V, 0.05 V a cell, 24 +; cells a side. 0.05 V divides both 1.2 V and 0.35 V, so Vdd and Vth are grid lines and the threshold kink is never +; inside a cell. One table serves both transistors: the NMOS reads it at (A, Y), the PMOS at (Vdd - A, Vdd - Y). The cell +; reads the table by bilinear interpolation. gm and gds are the slopes of that same interpolant (differences of the +; neighbouring entries, weighted), so the value and its derivative are one object. This one-node Newton uses gds only +; (A is an input); gm is built and tested for the cell with a second nonlinear node. The tests build the table from the +; analytic law, node by node, quantized; M3b will build it from PSP103 and change nothing else. +; - +; THE UNITS. A voltage word counts units of 0.05 V / 2^20 = 47.68 nV, so one grid cell is exactly 2^20 units: the cell +; index of a voltage is a shift and the position inside the cell a mask, and the interpolation divides by 2^40, a shift. +; Vdd = 24 * 2^20 = 25165824 and Vth = 7 * 2^20, both exact. (M0 counts 2^-24 V = 59.6 nV; here the grid sets the unit, +; 1.25 times finer, and every division but Newton's becomes a shift.) A table entry is the current times +; dt / C, that is the voltage the node moves in one step, in the same unit. For this circuit dt / C = 10 ns / 10 pF = +; 1000 V/A = 1 V/mA, so an entry is also the current in units of 47.68 pA; a different dt or C is one constant at the +; moment the table is built. The largest entry, Id(1.2 V, 1.2 V) = 0.36125 mA, is 7575962 < 2^23. +; - +; THE STEP. Backward Euler at the one unknown node: C (y' - y0) / dt = Ip(y') - In(y'), which is +; g(y) = y + In(A, y) - y0 - Ip(Vdd - A, Vdd - y) = 0 (volts of charge per step, all terms unsigned). +; Newton solves it with K = 3 iterations, K fixed in the spec, whatever the input: y <- y - g(y) / g'(y), with +; g'(y) = 1 + gds_n + gds_p in units of dt / C, the sum of the two cell slopes. The positive terms (y + In) and the +; negative terms (y0 + Ip) are two unsigned sums, g is their difference with a direction flag, and the update moves y +; toward the root, clamped to [0, Vdd] (the root is inside: g(0) <= 0 <= g(Vdd), and g' >= 1 because gds >= 0). The one +; division is q = round(|g| 2^20 / D) with D = 2^20 + the two cell slopes. The residual |g(y_K)| after the last +; iteration is computed at every step and bounded. +; - +; THE RUN. 400 steps from Y = 0: A = 0 V for steps 0..199 (0 .. 2 us), A = 1.2 V for steps 200..399, as M0 steps its +; source. Y(2 us) is the sample after step 199, Y(4 us) after step 399. S is an order-sensitive fold of all 400 samples. +; - +; THE THREE ERRORS, kept apart (W10): +; arithmetic the fixed-point run against an f64 run of the SAME iteration (same table values, same K, same clamps): +; largest |difference| 0.61 unit (29 nV), bound 0.7, asserted from both sides. +; table the f64 run of the table model against the analytic square law solved to convergence at the same dt: +; largest |difference| 0.428 mV, bound 0.45 mV, from both sides. Newton at K = 3 is already converged on +; the table (K = 3 and K = 12 agree to 1e-3 unit), so this is the table and not the iteration. +; discretization the analytic run at dt against the analytic run at dt / 2: largest |difference| 27.2 mV (at the first +; rising edge, 40 ns in), bound 27.5 mV, from both sides. This is the biggest error here and it is the +; scheme's: backward Euler at dt = 10 ns on a node whose time constant is about 12 ns. +; THE REQUIREMENT. |Y(4 us)| is 0 units (0 V) in the fixed-point run and 8e-54 V in the analytic one; cpphdl requires +; <= 1 mV. Y(2 us) is 1.2 V in both (25165824 units; the analytic run is within 1e-6 V of it, measured 2e-16). Y +; reaches 1.2 V - 1 mV at step 12 and falls below 1 mV 12 steps after the edge. +; THE TABLE AGAINST THE LAW, at all 576 cell midpoints, the worst place for a chord: absolute error at most 3.13e-7 A +; (0.09 % of the 0.36 mA saturation current), which is beta h^2 / 8 = 3.125e-7 A for h = 0.05 V plus the table's own +; rounding (bound 3.2e-7, asserted from both sides). Relative error, with a floor of 1e-9 A on the law: at most 1.0000 +; over all midpoints, and that is a limit of the method, not a defect of this table: in the first cell above Vth, +; where the law is ov^2 / 2 (Vds above the overdrive), a chord reads exactly twice the law at the midpoint, on any +; grid. From the second cell above Vth on the relative error is at most 0.2000 (bound 0.2002, from both sides). The +; slopes: gm and gds differ from the analytic derivatives at the midpoints by at most 0.0125 of dt / C (bound 0.0126). +; - +; THE NEGATIVE CONTROLS, each a test that passes by seeing the failure: +; K = 1 the residual reaches 1199464 units (57 mV) against the bound of 2; K = 2 reaches 29760 (1.4 mV). +; Vth = 0.30 a table built with the wrong threshold is 917504 units off the law at a node, 4.3e-5 A off at a +; midpoint (bound 3.2e-7), 1.89 relative (bound 0.2002), and the inverter run is 46.6 mV from the law +; (bound 0.45 mV). +; And one symmetry: the PMOS reads the table that the NMOS reads, mirrored, so the circuit run with A -> Vdd - A from +; Y = Vdd gives Vdd - Y at every step, bit for bit (a test). +; - +; THE HARDWARE. Everything from `gate_at` to `board_step`, and `fold`, takes and returns bool and unsigned words and has +; no float: it is what a board would run (`t27b verilog` refuses signed types). `fold` and `fabs64` are not written +; again: they come from M1 (`use analog::linear`, specs/analog/linear.t27), whose other names this file avoids (`step`, +; `residual`, `STEPS` became `board_step`, `g_resid`, `RUN_STEPS`). The table is passed as a slice. The f64 side (the law, the +; table builder, the f64 model of the iteration, the analytic run) is reached only by the tests. No board top and no +; clock are in this file; the cost of the division and of the table memory is measured in M3b. +; - +; WHAT A PASS SAYS. The cell reads a 25 x 25 table by exact-integer bilinear interpolation, solves backward Euler at +; the inverter's output node with a fixed three Newton iterations whose last residual is at most 1 unit (47.68 nV) at +; every one of 400 steps, ends the run at 0 V where cpphdl requires at most 1 mV, stays within 0.61 unit of the f64 +; model of the same iteration, within 0.428 mV of the analytic square law solved to convergence, and differs from its +; own dt / 2 run by at most 27.2 mV; the 400 samples fold to the S pinned below. +; WHAT IT DOES NOT SAY. Nothing about a transistor: the table is the square law, not PSP103 (M3b), and nothing here is a +; PDK corner. Nothing about a board: no bitstream, no clock, no LUT, DSP or BRAM, no division cost, no UART line. +; Nothing about the continuous circuit: the 27 mV is the scheme at dt = 10 ns, and the table error adds to it. Nothing +; about a circuit with a second nonlinear node (the Newton here is scalar), a netlist (the reader is M2), a device +; with Vds < 0, or an input other than the two levels. Nothing about power or about speed against SPICE. cpphdl's +; own number (11.5 nV) was not re-run: it is quoted from its README. +; ASCII only, English only. +; phi^2 + 1/phi^2 = 3 | TRINITY + +module AnalogMosTable; + +use analog::linear::{fold, fabs64}; + +pub const KIND : str = "analog-cell"; +pub const ID : str = "analog/mos_table"; +pub const REFERENCE : str = "mirekez/cpphdl analog/examples/cmos_inverter, analog/include/solver.hpp"; + +; --- the circuit -------------------------------------------------------------------------------- +/// beta of both transistors, A/V^2. +pub const BETA : f64 = 0.001; +/// Threshold of both transistors, V. +pub const VTH_V : f64 = 0.35; +/// The deliberately wrong threshold of the negative control, V. +pub const VTH_WRONG_V : f64 = 0.30; +pub const VDD_V : f64 = 1.2; +pub const DT_NS : u32 = 10; +pub const C_PF : u32 = 10; +/// dt / C in V per A: 10 ns / 10 pF. +pub const STEP_V_PER_A : f64 = 1000.0; + +; --- the grid and the units --------------------------------------------------------------------- +pub const CELL_BITS : u32 = 20; +/// One grid cell, 0.05 V, in voltage units: 2^20. +pub const CELL : u32 = 1048576; +pub const NCELL : u32 = 24; +pub const NODES : u32 = 25; +pub const TBL_LEN : usize = 625; +/// 24 cells = 1.2 V. +pub const VDD : u32 = 25165824; +/// 2^39: half of the interpolation's divisor 2^40, for rounding half up. +pub const HALF_Q : u64 = 549755813888; +/// Voltage units per volt: 20 cells * 2^20. +pub const UNITS_PER_V : f64 = 20971520.0; +pub const CELL_F : f64 = 1048576.0; +pub const VDD_F : f64 = 25165824.0; +/// The largest table entry, Id(1.2 V, 1.2 V) in table units. +pub const TBL_MAX : u32 = 7575962; + +; --- the run ------------------------------------------------------------------------------------ +pub const RUN_STEPS : u32 = 400; +/// The first step with A = Vdd: 2 us. +pub const EDGE_STEP : u32 = 200; +/// Newton iterations per step, fixed. +pub const K_NEWTON : u32 = 3; +/// The most the last residual may be, in units (95 nV). +pub const RES_BOUND : u32 = 2; +/// Iterations of the analytic reference's Newton (it needs fewer than 10). +pub const ANA_ITERS : u32 = 20; + +; --- what the run must give, derived by the tests below ----------------------------------------------- +pub const EXP_Y0 : u32 = 7570207; +pub const EXP_Y3 : u32 = 21437831; +pub const EXP_Y200 : u32 = 17595617; +pub const EXP_S : u64 = 0x703F6836B1F5B275; + +; --- the hardware side: unsigned words only ----------------------------------------------------------- + +/// The gate voltage at step n: 0 V before 2 us, Vdd from then on. +pub fn gate_at(n: u32) -> u32 { + if (n < EDGE_STEP) { return 0; } + return VDD; +} + +/// A voltage held to the table's range. +pub fn lim(v: u32) -> u32 { + if (v > VDD) { return VDD; } + return v; +} + +/// The grid cell of a voltage; the top node belongs to the last cell. +pub fn cell_of(v: u32) -> u32 { + const c : u32 = v >> CELL_BITS; + if (c >= NCELL) { return NCELL - 1; } + return c; +} + +/// The table entry at node (cx, cy), cx along Vgs and cy along Vds, as a wide word. +pub fn node_at(t: []const u32, cx: u32, cy: u32) -> u64 { + return (t[(cx * NODES + cy) as usize]) as u64; +} + +/// a - b, or 0 when b is larger; a table of a law that does not decrease never needs the 0. +pub fn dpos(a: u64, b: u64) -> u64 { + if (a >= b) { return a - b; } + return 0; +} + +/// The bilinear blend of four nodes: t00 at (0, 0), t10 at (1, 0), t01 at (0, 1), t11 at (1, 1), with fx, fy in [0, CELL] +/// the position inside the cell. In table units, rounded half up. The four weights add up to 2^40 exactly, so the sum +/// is at most TBL_MAX * 2^40 < 2^64. Scalars only: this is the arithmetic of one table read, apart from the memory. +pub fn blend(t00: u64, t10: u64, t01: u64, t11: u64, fx: u32, fy: u32) -> u32 { + const wx1 : u64 = fx as u64; + const wx0 : u64 = (CELL - fx) as u64; + const lo : u64 = t00 * wx0 + t10 * wx1; + const hi : u64 = t01 * wx0 + t11 * wx1; + return ((lo * ((CELL - fy) as u64) + hi * (fy as u64) + HALF_Q) >> 40) as u32; +} + +/// The linear blend of two cell differences d0 (at 0) and d1 (at 1) with the position f in [0, CELL], rounded half up. +pub fn blend1(d0: u64, d1: u64, f: u32) -> u32 { + return ((d0 * ((CELL - f) as u64) + d1 * (f as u64) + ((CELL >> 1) as u64)) >> CELL_BITS) as u32; +} + +/// The position of voltage v inside its cell, 0..CELL (CELL only at the top node). +pub fn frac_of(v: u32) -> u32 { + return v - (cell_of(v) << CELL_BITS); +} + +/// Id(vg, vd) by bilinear interpolation, in table units. +pub fn id_of(t: []const u32, vg: u32, vd: u32) -> u32 { + const gv : u32 = lim(vg); + const dv : u32 = lim(vd); + const cx : u32 = cell_of(gv); + const cy : u32 = cell_of(dv); + return blend(node_at(t, cx, cy), node_at(t, cx + 1, cy), node_at(t, cx, cy + 1), node_at(t, cx + 1, cy + 1), frac_of(gv), frac_of(dv)); +} + +/// gds: how much Id moves across one cell of Vds, at (vg, vd), in table units per cell. Divided by CELL it is the slope +/// of the interpolant, d(Id dt/C) / dVds, a pure number. +pub fn gds_cell(t: []const u32, vg: u32, vd: u32) -> u32 { + const gv : u32 = lim(vg); + const dv : u32 = lim(vd); + const cx : u32 = cell_of(gv); + const cy : u32 = cell_of(dv); + const d0 : u64 = dpos(node_at(t, cx, cy + 1), node_at(t, cx, cy)); + const d1 : u64 = dpos(node_at(t, cx + 1, cy + 1), node_at(t, cx + 1, cy)); + return blend1(d0, d1, frac_of(gv)); +} + +/// gm: how much Id moves across one cell of Vgs, at (vg, vd), in table units per cell. +pub fn gm_cell(t: []const u32, vg: u32, vd: u32) -> u32 { + const gv : u32 = lim(vg); + const dv : u32 = lim(vd); + const cx : u32 = cell_of(gv); + const cy : u32 = cell_of(dv); + const d0 : u64 = dpos(node_at(t, cx + 1, cy), node_at(t, cx, cy)); + const d1 : u64 = dpos(node_at(t, cx + 1, cy + 1), node_at(t, cx, cy + 1)); + return blend1(d0, d1, frac_of(dv)); +} + +/// The positive terms of g at y: y + In(A, y). +pub fn g_pos(t: []const u32, y: u32, a: u32) -> u32 { + return y + id_of(t, a, y); +} + +/// The negative terms of g at y: y0 + Ip(Vdd - A, Vdd - y). The PMOS reads the same table, mirrored. +pub fn g_neg(t: []const u32, y: u32, y0: u32, a: u32) -> u32 { + return y0 + id_of(t, VDD - a, VDD - y); +} + +/// |p - n|. +pub fn mag_diff(p: u32, n: u32) -> u32 { + if (p >= n) { return p - n; } + return n - p; +} + +/// |g(y)|, the residual of the backward Euler equation at y, in units. +pub fn g_resid(t: []const u32, y: u32, y0: u32, a: u32) -> u32 { + return mag_diff(g_pos(t, y, a), g_neg(t, y, y0, a)); +} + +/// round(g * 2^20 / d): the Newton move for a residual of magnitude g and a derivative of d / 2^20. d >= 2^20. +pub fn newton_move(g: u32, d: u64) -> u32 { + return ((((g as u64) << CELL_BITS) + (d >> 1)) / d) as u32; +} + +/// One Newton iteration of backward Euler at the output node: from y toward the root of g, clamped to [0, Vdd]. +pub fn newton_iter(t: []const u32, y: u32, y0: u32, a: u32) -> u32 { + const p : u32 = g_pos(t, y, a); + const n : u32 = g_neg(t, y, y0, a); + const d : u64 = (CELL as u64) + (gds_cell(t, a, y) as u64) + (gds_cell(t, VDD - a, VDD - y) as u64); + if (p >= n) { + const down : u32 = newton_move(p - n, d); + if (down >= y) { return 0; } + return y - down; + } + const up : u32 = newton_move(n - p, d); + if (up >= VDD - y) { return VDD; } + return y + up; +} + +/// One backward Euler step from y0 with the gate at a: k Newton iterations. +pub fn inv_step(t: []const u32, y0: u32, a: u32, k: u32) -> u32 { + var y : u32 = y0; + var i : u32 = 0; + while (i < k) { + y = newton_iter(t, y, y0, a); + i = i + 1; + } + return y; +} + +/// The step the board runs: K_NEWTON iterations. The trip count is a constant, so the loop unrolls (`t27b verilog` takes +/// a `while` only with a static bound); inv_step above takes k as an argument so the tests can run K = 1 and K = 2. +pub fn board_step(t: []const u32, y0: u32, a: u32) -> u32 { + var y : u32 = y0; + var i : u32 = 0; + while (i < K_NEWTON) { + y = newton_iter(t, y, y0, a); + i = i + 1; + } + return y; +} + +; --- the f64 side: references that only the tests reach ------------------------------------------------- +; A module-level f64 fn is carried to Verilog only when the hardware reaches it; nothing below is reached by `step`. + +/// The larger of two reals. +pub fn fmax64(a: f64, b: f64) -> f64 { + if (a > b) { return a; } + return b; +} + +/// cpphdl's MOS law: the drain current in A for Vgs, Vds in V, with threshold vth. +pub fn law_amps(vgs: f64, vds: f64, vth: f64) -> f64 { + const ov : f64 = vgs - vth; + if (ov <= 0.0 or vds <= 0.0) { return 0.0; } + if (vds < ov) { return BETA * (ov * vds - vds * vds / 2.0); } + return BETA * ov * ov / 2.0; +} + +/// dId/dVds of the law, A/V. +pub fn law_gds(vgs: f64, vds: f64, vth: f64) -> f64 { + const ov : f64 = vgs - vth; + if (ov <= 0.0 or vds <= 0.0) { return 0.0; } + if (vds < ov) { return BETA * (ov - vds); } + return 0.0; +} + +/// dId/dVgs of the law, A/V. +pub fn law_gm(vgs: f64, vds: f64, vth: f64) -> f64 { + const ov : f64 = vgs - vth; + if (ov <= 0.0 or vds <= 0.0) { return 0.0; } + if (vds < ov) { return BETA * vds; } + return BETA * ov; +} + +/// Fill the table from the law with threshold vth: entry (ix, iy) = round(Id(ix * 0.05 V, iy * 0.05 V) * dt/C), in units. +pub fn build_table(t: []u32, vth: f64) { + var ix : u32 = 0; + while (ix < NODES) { + var iy : u32 = 0; + while (iy < NODES) { + const i : f64 = law_amps((ix as f64) / 20.0, (iy as f64) / 20.0, vth); + t[(ix * NODES + iy) as usize] = (i * STEP_V_PER_A * UNITS_PER_V + 0.5) as u32; + iy = iy + 1; + } + ix = ix + 1; + } +} + +/// The largest |entry - law| over the 625 nodes, in table units, against the law with the true threshold. +pub fn node_err(t: []const u32) -> f64 { + var worst : f64 = 0.0; + var ix : u32 = 0; + while (ix < NODES) { + var iy : u32 = 0; + while (iy < NODES) { + const want : f64 = law_amps((ix as f64) / 20.0, (iy as f64) / 20.0, VTH_V) * STEP_V_PER_A * UNITS_PER_V; + worst = fmax64(worst, fabs64((t[(ix * NODES + iy) as usize] as f64) - want)); + iy = iy + 1; + } + ix = ix + 1; + } + return worst; +} + +pub struct MidStats { abs_a: f64, rel_all: f64, rel_hi: f64, gds_dev: f64, gm_dev: f64 } + +/// The table against the law at the 576 cell midpoints: the largest absolute error in A; the largest relative error with +/// a floor of 1e-9 A on the law, over all midpoints and over the cells from the second above Vth on (cx >= 8); the +/// largest deviation of gds and gm from the law's derivatives, in units of dt/C. +pub fn mid_stats(t: []const u32) -> MidStats { + var abs_a : f64 = 0.0; + var rel_all : f64 = 0.0; + var rel_hi : f64 = 0.0; + var gds_dev : f64 = 0.0; + var gm_dev : f64 = 0.0; + var cx : u32 = 0; + while (cx < NCELL) { + var cy : u32 = 0; + while (cy < NCELL) { + const vg : u32 = cx * CELL + (CELL >> 1); + const vd : u32 = cy * CELL + (CELL >> 1); + const gf : f64 = (vg as f64) / UNITS_PER_V; + const df : f64 = (vd as f64) / UNITS_PER_V; + const law : f64 = law_amps(gf, df, VTH_V); + const tab : f64 = (id_of(t, vg, vd) as f64) / (STEP_V_PER_A * UNITS_PER_V); + const e : f64 = fabs64(tab - law); + const r : f64 = e / fmax64(law, 0.000000001); + abs_a = fmax64(abs_a, e); + rel_all = fmax64(rel_all, r); + if (cx >= 8) { rel_hi = fmax64(rel_hi, r); } + gds_dev = fmax64(gds_dev, fabs64(((gds_cell(t, vg, vd) as f64) / CELL_F) - law_gds(gf, df, VTH_V) * STEP_V_PER_A)); + gm_dev = fmax64(gm_dev, fabs64(((gm_cell(t, vg, vd) as f64) / CELL_F) - law_gm(gf, df, VTH_V) * STEP_V_PER_A)); + cy = cy + 1; + } + cx = cx + 1; + } + return MidStats { abs_a: abs_a, rel_all: rel_all, rel_hi: rel_hi, gds_dev: gds_dev, gm_dev: gm_dev }; +} + +/// Is the table nondecreasing along both axes (a law with gm >= 0 and gds >= 0 gives one)? +pub fn monotone(t: []const u32) -> bool { + var ok : bool = true; + var ix : u32 = 0; + while (ix < NODES - 1) { + var iy : u32 = 0; + while (iy < NODES - 1) { + const here : u32 = t[(ix * NODES + iy) as usize]; + if (t[((ix + 1) * NODES + iy) as usize] < here) { ok = false; } + if (t[(ix * NODES + iy + 1) as usize] < here) { ok = false; } + iy = iy + 1; + } + ix = ix + 1; + } + return ok; +} + +/// The fixed-point run: 400 steps of `inv_step` with k iterations from Y = y_init, the gate mirrored (Vdd - A) when +/// mirror is set. out gets every Y, res the last residual of every step. Returns the fold of the 400 samples. +pub fn run_fix(t: []const u32, k: u32, y_init: u32, mirror: bool, out: []u32, res: []u32) -> u64 { + var y : u32 = y_init; + var s : u64 = 0; + var n : u32 = 0; + while (n < RUN_STEPS) { + var a : u32 = gate_at(n); + if (mirror) { a = VDD - a; } + const y0 : u32 = y; + y = inv_step(t, y0, a, k); + res[n as usize] = g_resid(t, y, y0, a); + out[n as usize] = y; + s = fold(s, y); + n = n + 1; + } + return s; +} + +/// The largest residual of a run, in units. +pub fn max_of(r: []const u32) -> u32 { + var m : u32 = 0; + var n : u32 = 0; + while (n < RUN_STEPS) { + if (r[n as usize] > m) { m = r[n as usize]; } + n = n + 1; + } + return m; +} + +/// The cell of a real voltage in units (the f64 model's cell_of). +pub fn cellf(v: f64) -> u32 { + var c : u32 = (v / CELL_F) as u32; + if (c >= NCELL) { c = NCELL - 1; } + return c; +} + +pub fn nodef(t: []const u32, cx: u32, cy: u32) -> f64 { + return (t[(cx * NODES + cy) as usize]) as f64; +} + +/// id_of in real arithmetic on the same table values: no rounding. +pub fn idf(t: []const u32, vg: f64, vd: f64) -> f64 { + const cx : u32 = cellf(vg); + const cy : u32 = cellf(vd); + const fx : f64 = vg - (cx as f64) * CELL_F; + const fy : f64 = vd - (cy as f64) * CELL_F; + const lo : f64 = nodef(t, cx, cy) * (CELL_F - fx) + nodef(t, cx + 1, cy) * fx; + const hi : f64 = nodef(t, cx, cy + 1) * (CELL_F - fx) + nodef(t, cx + 1, cy + 1) * fx; + return (lo * (CELL_F - fy) + hi * fy) / (CELL_F * CELL_F); +} + +/// gds_cell in real arithmetic. +pub fn gdsf(t: []const u32, vg: f64, vd: f64) -> f64 { + const cx : u32 = cellf(vg); + const cy : u32 = cellf(vd); + const fx : f64 = vg - (cx as f64) * CELL_F; + const d0 : f64 = fmax64(nodef(t, cx, cy + 1) - nodef(t, cx, cy), 0.0); + const d1 : f64 = fmax64(nodef(t, cx + 1, cy + 1) - nodef(t, cx + 1, cy), 0.0); + return (d0 * (CELL_F - fx) + d1 * fx) / CELL_F; +} + +/// newton_iter in real arithmetic: the same equation, the same slopes, the same clamp, no rounding. +pub fn newton_iter_f(t: []const u32, y: f64, y0: f64, a: f64) -> f64 { + const p : f64 = y + idf(t, a, y); + const n : f64 = y0 + idf(t, VDD_F - a, VDD_F - y); + const d : f64 = CELL_F + gdsf(t, a, y) + gdsf(t, VDD_F - a, VDD_F - y); + var z : f64 = y - (p - n) * CELL_F / d; + if (z < 0.0) { z = 0.0; } + if (z > VDD_F) { z = VDD_F; } + return z; +} + +/// The f64 model of the table run: 400 steps of k real Newton iterations, out in units. +pub fn run_tabf(t: []const u32, k: u32, out: []f64) { + var y : f64 = 0.0; + var n : u32 = 0; + while (n < RUN_STEPS) { + const a : f64 = (gate_at(n)) as f64; + const y0 : f64 = y; + var i : u32 = 0; + while (i < k) { + y = newton_iter_f(t, y, y0, a); + i = i + 1; + } + out[n as usize] = y; + n = n + 1; + } +} + +/// The analytic reference: backward Euler of the inverter with cpphdl's law, Newton to convergence (ANA_ITERS). f scales +/// dt (1.0 is 10 ns, 0.5 is 5 ns); the gate steps at step edge. out gets Y in volts after each of nsteps steps. +pub fn run_ana(f: f64, nsteps: u32, edge: u32, out: []f64) { + const k : f64 = STEP_V_PER_A * f; + var y : f64 = 0.0; + var n : u32 = 0; + while (n < nsteps) { + var a : f64 = 0.0; + if (n >= edge) { a = VDD_V; } + const y0 : f64 = y; + var i : u32 = 0; + while (i < ANA_ITERS) { + const g : f64 = y - y0 - k * law_amps(VDD_V - a, VDD_V - y, VTH_V) + k * law_amps(a, y, VTH_V); + const gp : f64 = 1.0 + k * law_gds(VDD_V - a, VDD_V - y, VTH_V) + k * law_gds(a, y, VTH_V); + y = y - g / gp; + if (y < 0.0) { y = 0.0; } + if (y > VDD_V) { y = VDD_V; } + i = i + 1; + } + out[n as usize] = y; + n = n + 1; + } +} + +; --- the tests ---------------------------------------------------------------------------------- + +test the_grid_and_the_units_are_exact { + // Variables, not constants: a constant assert is folded by the compiler and counts as no assert at all. + var cell : u32 = CELL; + var bits : u32 = CELL_BITS; + var ncell : u32 = NCELL; + var nodes : u32 = NODES; + var vdd : u32 = VDD; + var half : u64 = HALF_Q; + var tmax : u32 = TBL_MAX; + var upv : f64 = UNITS_PER_V; + var cellf_ : f64 = CELL_F; + var vddf : f64 = VDD_F; + assert_eq(cell, 1 << bits); + assert_eq(vdd, ncell * cell); + assert_eq(nodes, ncell + 1); + assert_eq(nodes * nodes, TBL_LEN as u32); + assert_eq(half, 1 << 39); + assert_eq(upv, 20.0 * cellf_); + assert_eq(vddf, (vdd as f64)); + assert_eq(cellf_, (cell as f64)); + // 1.2 V is 24 cells and 0.35 V is 7, so both are grid lines and nothing is rounded. + assert_eq(vddf / upv, VDD_V); + assert_eq(7.0 * cellf_ / upv, VTH_V); + // dt / C is 10 ns / 10 pF = 1000 V/A. + const dtc : f64 = (DT_NS as f64) * 0.000000001 / ((C_PF as f64) * 0.000000000001) - STEP_V_PER_A; + assert dtc * dtc < 0.000001 * 0.000001; + // The widest product of the interpolation: TBL_MAX * 2^40 + the rounding half fits in a u64. + assert tmax < 8388608; + assert (tmax as u64) * 1099511627776 + half < 9223372036854775808; +} + +test the_table_is_the_law_at_the_nodes { + var t : [625]u32 = [0] ** 625; + build_table(&t, VTH_V); + assert node_err(&t) <= 0.5 + 0.0000001; + assert node_err(&t) > 0.3; + assert(monotone(&t)); + assert_eq(t[624], TBL_MAX); + assert_eq(t[0], 0); + // Vgs on the threshold line (node 7) and below: no current. Vds = 0: no current. + var ix : u32 = 0; + var iy : u32 = 0; + while (ix <= 7) { + iy = 0; + while (iy < NODES) { + assert_eq(t[(ix * NODES + iy) as usize], 0); + iy = iy + 1; + } + ix = ix + 1; + } + ix = 0; + while (ix < NODES) { + assert_eq(t[(ix * NODES) as usize], 0); + ix = ix + 1; + } + // Saturation is flat in Vds: Vgs = 0.55 V (node 11) has ov = 0.2 V, so Vds >= 0.2 V (node 4) gives one value. + assert_eq(t[(11 * NODES + 4) as usize], t[(11 * NODES + 24) as usize]); + assert t[(11 * NODES + 3) as usize] < t[(11 * NODES + 4) as usize]; +} + +test the_interpolation_is_exact_at_the_nodes_and_bilinear_inside { + var t : [625]u32 = [0] ** 625; + build_table(&t, VTH_V); + var ix : u32 = 0; + var iy : u32 = 0; + while (ix < NODES) { + iy = 0; + while (iy < NODES) { + assert_eq(id_of(&t, ix * CELL, iy * CELL), t[(ix * NODES + iy) as usize]); + iy = iy + 1; + } + ix = ix + 1; + } + // A cell centre is the mean of the four corners, rounded half up: cell (12, 4). + const c4 : u32 = t[(12 * NODES + 4) as usize] + t[(13 * NODES + 4) as usize] + t[(12 * NODES + 5) as usize] + t[(13 * NODES + 5) as usize]; + assert_eq(id_of(&t, 12 * CELL + (CELL >> 1), 4 * CELL + (CELL >> 1)), (c4 + 2) / 4); + // Along an edge it is linear in the position: a quarter of the way along Vds at Vgs = node 12. + const e0 : u32 = t[(12 * NODES + 8) as usize]; + const e1 : u32 = t[(12 * NODES + 9) as usize]; + assert_eq(id_of(&t, 12 * CELL, 8 * CELL + (CELL >> 2)), (3 * e0 + e1 + 2) / 4); + // The scalar arithmetic alone: the corners at the corners, the mean at the centre, one rounding in all. + assert_eq(blend(100, 200, 300, 400, 0, 0), 100); + assert_eq(blend(100, 200, 300, 400, CELL, 0), 200); + assert_eq(blend(100, 200, 300, 400, 0, CELL), 300); + assert_eq(blend(100, 200, 300, 400, CELL, CELL), 400); + assert_eq(blend(100, 200, 300, 400, CELL >> 1, CELL >> 1), 250); + assert_eq(blend(1, 0, 0, 0, CELL >> 1, CELL >> 1), 0); + assert_eq(blend(2, 0, 0, 0, CELL >> 1, CELL >> 1), 1); + assert_eq(blend1(10, 30, 0), 10); + assert_eq(blend1(10, 30, CELL), 30); + assert_eq(blend1(10, 30, CELL >> 1), 20); + assert_eq(frac_of(VDD), CELL); + assert_eq(frac_of(3 * CELL + 7), 7); + // The step the board runs is inv_step at K_NEWTON, and the loop count is the only difference. + assert_eq(board_step(&t, 4000000, 0), inv_step(&t, 4000000, 0, K_NEWTON)); + assert_eq(board_step(&t, 4000000, VDD), inv_step(&t, 4000000, VDD, K_NEWTON)); + assert_eq(board_step(&t, VDD, 0), VDD); + assert_eq(board_step(&t, 0, VDD), 0); + // Inputs past Vdd are held at the last node. + assert_eq(id_of(&t, VDD + 5, VDD + 9), TBL_MAX); + assert_eq(cell_of(VDD), NCELL - 1); + assert_eq(cell_of(VDD - 1), NCELL - 1); + assert_eq(cell_of(CELL - 1), 0); + assert_eq(cell_of(CELL), 1); +} + +test gm_and_gds_are_the_slopes_of_the_interpolant { + var t : [625]u32 = [0] ** 625; + build_table(&t, VTH_V); + // Across a whole cell of Vds the interpolant moves by exactly gds_cell (the two reads and the slope each round). + var cy : u32 = 3; + var dy : u32 = 0; + var cx : u32 = 8; + var dx : u32 = 0; + var fx : u32 = 0; + while (cx < 24) { + fx = (cx * 977) % CELL; + cy = 3; + while (cy < 23) { + dy = id_of(&t, cx * CELL + fx, (cy + 1) * CELL) - id_of(&t, cx * CELL + fx, cy * CELL); + assert mag_diff(dy, gds_cell(&t, cx * CELL + fx, cy * CELL + 12345)) <= 2; + dx = id_of(&t, (cx + 1) * CELL, cy * CELL + fx) - id_of(&t, cx * CELL, cy * CELL + fx); + assert mag_diff(dx, gm_cell(&t, cx * CELL + 12345, cy * CELL + fx)) <= 2; + cy = cy + 1; + } + cx = cx + 1; + } + // Spot values: saturation is flat in Vds (gds 0), gm is positive there, and below threshold both are 0. + assert_eq(gds_cell(&t, 11 * CELL, 10 * CELL), 0); + assert gm_cell(&t, 11 * CELL, 10 * CELL) > 0; + assert_eq(gm_cell(&t, 3 * CELL, 10 * CELL), 0); + assert_eq(gds_cell(&t, 3 * CELL, 10 * CELL), 0); + // A table that decreased would be clamped at 0, not wrapped: dpos. + assert_eq(dpos(5, 9), 0); + assert_eq(dpos(9, 5), 4); +} + +test the_table_against_the_law_at_the_cell_midpoints { + var t : [625]u32 = [0] ** 625; + build_table(&t, VTH_V); + const m : MidStats = mid_stats(&t); + // Absolute: beta h^2 / 8 = 3.125e-7 A, from both sides. + assert m.abs_a < 0.00000032; + assert m.abs_a > 0.0000003; + // Relative, floor 1e-9 A: all midpoints, a chord against a square reads twice the law in the first cell above Vth. + assert m.rel_all < 1.0001; + assert m.rel_all > 0.999; + // From the second cell above Vth on. + assert m.rel_hi < 0.2002; + assert m.rel_hi > 0.19; + // The slopes of the interpolant against the law's derivatives, in units of dt / C. + assert m.gds_dev < 0.0126; + assert m.gds_dev > 0.012; + assert m.gm_dev < 0.0126; + assert m.gm_dev > 0.012; +} + +test the_first_steps_are_pinned_and_the_inverter_ends_at_zero { + var t : [625]u32 = [0] ** 625; + build_table(&t, VTH_V); + var out : [400]u32 = [0] ** 400; + var res : [400]u32 = [0] ** 400; + const s : u64 = run_fix(&t, K_NEWTON, 0, false, &out, &res); + assert_eq(out[0], EXP_Y0); + assert_eq(out[3], EXP_Y3); + assert_eq(out[199], VDD); + assert_eq(out[200], EXP_Y200); + assert_eq(out[399], 0); + assert_eq(s, EXP_S); + assert_eq(board_step(&t, 0, 0), EXP_Y0); + // cpphdl's requirement: |Y(4 us)| <= 1 mV, and Y(2 us) is 1.2 V. + const y4 : f64 = (out[399] as f64) / UNITS_PER_V; + assert y4 * y4 <= 0.001 * 0.001; + const y2 : f64 = (out[199] as f64) / UNITS_PER_V; + assert_eq(y2, VDD_V); + // Y stays in [0, Vdd], rises while the input is low, falls once it is high. + var n : u32 = 0; + while (n < RUN_STEPS) { + assert out[n as usize] <= VDD; + if (n > 0 and n < EDGE_STEP) { assert out[n as usize] >= out[(n - 1) as usize]; } + if (n > EDGE_STEP) { assert out[n as usize] <= out[(n - 1) as usize]; } + n = n + 1; + } + // Y is within 1 mV of 1.2 V by step 12 and within 1 mV of 0 by step 212. + const mv : f64 = 0.001 * UNITS_PER_V; + assert ((VDD - out[11]) as f64) > mv; + assert ((VDD - out[12]) as f64) <= mv; + assert (out[211] as f64) > mv; + assert (out[212] as f64) <= mv; +} + +test the_last_residual_is_bounded_at_every_step { + var t : [625]u32 = [0] ** 625; + build_table(&t, VTH_V); + var out : [400]u32 = [0] ** 400; + var res : [400]u32 = [0] ** 400; + const s : u64 = run_fix(&t, K_NEWTON, 0, false, &out, &res); + var n : u32 = 0; + var ones : u32 = 0; + while (n < RUN_STEPS) { + assert res[n as usize] <= RES_BOUND; + if (res[n as usize] == 1) { ones = ones + 1; } + n = n + 1; + } + // From both sides: the worst is exactly one unit (the rounding of the last move), so the bound is within a factor 2. + assert_eq(max_of(&res), 1); + assert_eq(ones, 22); + assert_eq(s, EXP_S); + // The residual is |g|: recompute it for the first step by hand from the pieces. + const y1 : u32 = out[0]; + assert_eq(g_resid(&t, y1, 0, 0), mag_diff(y1 + id_of(&t, 0, y1), id_of(&t, VDD, VDD - y1))); +} + +test k1_is_the_negative_control_it_breaks_the_residual_bound { + var t : [625]u32 = [0] ** 625; + build_table(&t, VTH_V); + var out : [400]u32 = [0] ** 400; + var res : [400]u32 = [0] ** 400; + var res2 : [400]u32 = [0] ** 400; + var out2 : [400]u32 = [0] ** 400; + const s1 : u64 = run_fix(&t, 1, 0, false, &out, &res); + // One Newton iteration leaves 1199464 units (57 mV) of residual, 600000 times the bound. + assert max_of(&res) > RES_BOUND; + assert max_of(&res) > 1000000; + assert max_of(&res) < 1300000; + assert s1 != EXP_S; + const s2 : u64 = run_fix(&t, 2, 0, false, &out2, &res2); + assert max_of(&res2) > RES_BOUND; + assert_eq(max_of(&res2), 29760); + assert s2 != EXP_S; + // And the bound itself is not met by K = 1 anywhere near: the first rise step. + assert res[1] > 1000000; +} + +test the_mirrored_circuit_gives_the_mirrored_samples_bit_for_bit { + var t : [625]u32 = [0] ** 625; + build_table(&t, VTH_V); + var out : [400]u32 = [0] ** 400; + var res : [400]u32 = [0] ** 400; + var outm : [400]u32 = [0] ** 400; + var resm : [400]u32 = [0] ** 400; + const s : u64 = run_fix(&t, K_NEWTON, 0, false, &out, &res); + const sm : u64 = run_fix(&t, K_NEWTON, VDD, true, &outm, &resm); + var n : u32 = 0; + while (n < RUN_STEPS) { + assert_eq(outm[n as usize], VDD - out[n as usize]); + assert_eq(resm[n as usize], res[n as usize]); + n = n + 1; + } + assert s != sm; + // The mirror of a single step, off the run: arbitrary y0 and gate. + assert_eq(inv_step(&t, VDD - 9000000, VDD - 5000000, 3), VDD - inv_step(&t, 9000000, 5000000, 3)); + assert_eq(inv_step(&t, VDD - 123456, VDD, 3), VDD - inv_step(&t, 123456, 0, 3)); +} + +test the_three_errors_arithmetic_table_discretization { + var t : [625]u32 = [0] ** 625; + build_table(&t, VTH_V); + var out : [400]u32 = [0] ** 400; + var res : [400]u32 = [0] ** 400; + var tabf : [400]f64 = [0.0] ** 400; + var tab12 : [400]f64 = [0.0] ** 400; + var ana : [400]f64 = [0.0] ** 400; + var half : [800]f64 = [0.0] ** 800; + const s : u64 = run_fix(&t, K_NEWTON, 0, false, &out, &res); + run_tabf(&t, K_NEWTON, &tabf); + run_tabf(&t, 12, &tab12); + run_ana(1.0, 400, 200, &ana); + run_ana(0.5, 800, 400, &half); + var arith : f64 = 0.0; + var table : f64 = 0.0; + var disc : f64 = 0.0; + var lag : f64 = 0.0; + var n : u32 = 0; + while (n < RUN_STEPS) { + const i : usize = n as usize; + arith = fmax64(arith, fabs64((out[i] as f64) - tabf[i])); + table = fmax64(table, fabs64(tabf[i] / UNITS_PER_V - ana[i])); + disc = fmax64(disc, fabs64(ana[i] - half[2 * i + 1])); + lag = fmax64(lag, fabs64(tabf[i] - tab12[i])); + n = n + 1; + } + // Arithmetic: fixed point against the f64 model of the same iteration, in units (47.68 nV). Measured 0.610. + assert arith < 0.7; + assert arith > 0.5; + // Table: the f64 table model against the analytic law solved to convergence, in volts. Measured 0.428 mV. + assert table < 0.00045; + assert table > 0.0004; + // Newton at K = 3 is converged on the table: it equals K = 12 to far below a unit. + assert lag < 0.001; + // Discretization: dt against dt/2 on the analytic law, in volts. Measured 27.2 mV. + assert disc < 0.0275; + assert disc > 0.027; + // The issue's scenario: Y at 2 us within 1 uV of the analytic run, Y at 4 us within 1 uV of it and under 1 mV. + const d2 : f64 = (out[199] as f64) / UNITS_PER_V - ana[199]; + assert d2 * d2 < 0.000001 * 0.000001; + const d4 : f64 = (out[399] as f64) / UNITS_PER_V - ana[399]; + assert d4 * d4 < 0.000001 * 0.000001; + assert_eq(s, EXP_S); +} + +test the_analytic_reference_is_a_converged_backward_euler { + var ana : [400]f64 = [0.0] ** 400; + run_ana(1.0, 400, 200, &ana); + // Y(2 us) is Vdd, Y(4 us) is ~0 (cpphdl: 11.5 nV; this solver is converged, so it is below that). + const d2 : f64 = ana[199] - VDD_V; + assert d2 * d2 < 0.000000000001 * 0.000000000001; + assert ana[399] * ana[399] < 0.001 * 0.001; + // The first rise step by hand: Y1 solves y = k Ip(1.2 V, 1.2 V - y) in the saturation/linear mix. + const k : f64 = STEP_V_PER_A; + const y1 : f64 = ana[0]; + const g1 : f64 = y1 - k * law_amps(VDD_V, VDD_V - y1, VTH_V); + assert g1 * g1 < 0.0000000001 * 0.0000000001; + // The law's own checks: zero below threshold, the two regions meet at Vds = ov, symmetric drain/source not needed. + assert_eq(law_amps(0.35, 1.0, VTH_V), 0.0); + assert_eq(law_amps(1.2, 0.0, VTH_V), 0.0); + const sat : f64 = law_amps(1.2, 0.85, VTH_V); + const lin : f64 = law_amps(1.2, 0.8499999999, VTH_V); + assert (sat - lin) * (sat - lin) < 0.000000001 * 0.000000001; + const peak : f64 = law_amps(1.2, 1.2, VTH_V) - 0.00036125; + assert peak * peak < 0.000000000000001 * 0.000000000000001; +} + +test vth_030_is_the_negative_control_the_law_check_fails { + var good : [625]u32 = [0] ** 625; + var bad : [625]u32 = [0] ** 625; + build_table(&good, VTH_V); + build_table(&bad, VTH_WRONG_V); + // At the nodes: the right table is within half a unit of the law, the wrong one is nowhere near. + assert node_err(&good) <= 0.5 + 0.0000001; + assert node_err(&bad) > 900000.0; + // At the midpoints: past the absolute bound by a factor 100, past the relative bound by a factor 9. + const mg : MidStats = mid_stats(&good); + const mb : MidStats = mid_stats(&bad); + assert mg.abs_a < 0.00000032; + assert mb.abs_a > 0.00001; + assert mg.rel_hi < 0.2002; + assert mb.rel_hi > 1.0; + // In the circuit: the inverter run on the wrong table leaves the analytic run by tens of mV (bound 0.45 mV). + var outb : [400]u32 = [0] ** 400; + var resb : [400]u32 = [0] ** 400; + var ana : [400]f64 = [0.0] ** 400; + const sb : u64 = run_fix(&bad, K_NEWTON, 0, false, &outb, &resb); + run_ana(1.0, 400, 200, &ana); + var worst : f64 = 0.0; + var n : u32 = 0; + while (n < RUN_STEPS) { + worst = fmax64(worst, fabs64((outb[n as usize] as f64) / UNITS_PER_V - ana[n as usize])); + n = n + 1; + } + assert worst > 0.04; + assert worst < 0.05; + assert worst > 90.0 * 0.00045; + assert sb != EXP_S; + // The circuit still ends where cpphdl needs it: the wrong threshold is caught by the law check, not by the endpoint. + assert ((outb[399]) as f64) <= 0.001 * UNITS_PER_V; +} + +test a_wrong_reading_of_the_table_is_caught { + var t : [625]u32 = [0] ** 625; + build_table(&t, VTH_V); + // The table is not symmetric: reading it with the axes swapped gives another number at an off-diagonal point. + assert id_of(&t, 12 * CELL, 5 * CELL) != id_of(&t, 5 * CELL, 12 * CELL); + // Truncating the interpolation instead of rounding it moves a cell centre by the rounding half. + const c : u32 = t[(12 * NODES + 4) as usize] + t[(13 * NODES + 4) as usize] + t[(12 * NODES + 5) as usize] + t[(13 * NODES + 5) as usize]; + assert (c % 4) >= 2; + assert_eq(id_of(&t, 12 * CELL + (CELL >> 1), 4 * CELL + (CELL >> 1)), (c + 2) / 4); + assert id_of(&t, 12 * CELL + (CELL >> 1), 4 * CELL + (CELL >> 1)) != c / 4; + // The fold knows the order. + assert fold(fold(0, 1), 2) != fold(fold(0, 2), 1); + assert_eq(fold(0x8000000000000000, 0), 1); +}