diff --git a/.trinity/seals/analog_AnalogLinear.json b/.trinity/seals/analog_AnalogLinear.json new file mode 100644 index 0000000000..7f9373ef54 --- /dev/null +++ b/.trinity/seals/analog_AnalogLinear.json @@ -0,0 +1,28 @@ +{ + "built_by": "t27c-bootstrap@0.5.2+9b556b7c3", + "closure": [], + "config": "gen=c,rust,verilog,zig;tests=zig", + "gen_hash_c": "sha256:294ec2cd61e90873b7907ef44f25207a48967d6d323fee13d94612804e385e34", + "gen_hash_rust": "sha256:41c8a5d22739a6cf3bbf2fa7d2db29ca461a8a994d5b05d7838c14c21a781681", + "gen_hash_verilog": "sha256:94464f3fcd304291b93a6729872dfcaa643fbc67e05eaa5d00b0e2e213df25e4", + "gen_hash_zig": "sha256:65082e5f7bc98fd7874fe89d07131844820a32b360dba73116e62c87b4cd889c", + "module": "AnalogLinear", + "ring": 12, + "seal_schema": 2, + "sealed_at": "2026-10-18T15:37:30Z", + "sealed_by": "t27c-bootstrap@0.5.2", + "spec_hash": "sha256:3ddef7be183aba583ed285334ee389b48a65b41335e98914514235db7c4e5d62", + "spec_path": "specs/analog/linear.t27", + "tests": { + "failed": 0, + "failing": [], + "forced": false, + "passed": 14, + "total": 14, + "vacuous": 0 + }, + "toolchain": { + "test_runner": "sha256:1982ae948e24d00c5d221ccec12dfce1eb01204b53e37b0f857f93b16d64629d", + "zig": "0.16.0" + } +} \ No newline at end of file diff --git a/specs/analog/linear.t27 b/specs/analog/linear.t27 new file mode 100644 index 0000000000..3c72e30f8c --- /dev/null +++ b/specs/analog/linear.t27 @@ -0,0 +1,1086 @@ +// SPDX-License-Identifier: Apache-2.0 +; specs/analog/linear.t27 -- a linear analog cell: any small R / C / L circuit with one voltage source, fixed step, fixed point +; M1 of the analog emulation epic (#10486, issue #10487). M0 (specs/fpga/rc_filter_ax7203.t27) is the 1x1 case of this cell +; and this cell reproduces it bit for bit. Everything hand-written is this file. +; - +; THE METHOD. Modified nodal analysis with backward Euler and a fixed step dt gives (G + C/dt) x' = (C/dt) x + B u. With dt +; fixed, A = G + C/dt is a constant, so it is folded here, before any hardware exists, to x' = M x + N u with M = A^-1 C/dt +; and N = A^-1 B. The fabric never solves and never divides. The unknowns x are the free node voltages (node 1, 2, ...), +; then one entry per inductor, its current stored as i * ZS, which is volts (ZS is the circuit's characteristic impedance, +; and it keeps every entry of M near 1). The one driven node is NODE_IN (15), the V element's plus terminal, held at u. +; - +; THE STEP. D = M - I and N are stored as integers at 2^-24. One step is the delta form with ONE rounding per state: +; x_i' = x_i + rnd( sum_j D_ij x_j + N_i u ) +; the sum is exact in a wide integer and is rounded once, half away from zero. For the M0 circuit D = -798915 and +; N = 798915 (1/21 at 2^-24), so the sum is 798915 * (vin - v) and the rounding is M0's: the 200 samples, A, B and S of +; M0 are reproduced by construction. Rounding M itself would give x' = rnd(Mx + Nu) and a different, worse result. +; - +; THE NUMBERS. Every state, every constant and every input is an UNSIGNED 32-bit word in sign-magnitude: bit 31 is the +; direction flag (1 = negative, as M0's vup), bits 0..30 are the magnitude in units of 2^-24 (1 V = 16777216). Zero is +; never negative. No signed type is used, because `t27b verilog` refuses one. A product of two magnitudes is a u64. The +; positive terms and the negative terms add into two u64 registers (pos_after, neg_after); one subtraction gives the +; magnitude and the direction of the sum, round_back rounds it, add_sm adds it to the old state. Magnitudes stay below +; 2^26 in every run here (asserted), so the 31 bits and the u64 sums have wide headroom. +; - +; THE MACHINE (M1b, #10496, not in this file). One multiply-accumulate per clock over the n * (n + 1) constants of a circuit: +; pos_after and neg_after are the two accumulator registers, delta_of the one subtraction and the one rounding, next_x the +; write back. A constant that is a short sum of shifts needs no multiplier (M0's 1/21 is four shifted adds); the others +; take a DSP48. The tables below are arrays and the f64 references below are reached only by the tests. +; - +; THE CONSTANTS ARE DERIVED, NOT TRUSTED. The tests rebuild A, C/dt and B in f64 from each circuit's element table (the +; stamps of R, C, L and V), solve A [C/dt | B] = [M | N] by Gauss-Jordan elimination with partial pivoting, and assert +; that every stored integer equals round(value * 2^24), that no value is within 1e-4 of a rounding tie, and that A X - R is +; zero to 1e-12. A conductance stamped on the wrong node, a capacitor on the wrong node, an inductor with its terminals +; swapped each change a derived value and fail (the negative controls). +; - +; STABILITY IS PROVED, NOT SAMPLED. For the quantized M_q = I + D / 2^24 the tests assert ||M_q^256||_inf < 1 in f64 (eight +; squarings). Since rho^256 <= ||M^256||, that bounds the spectral radius below 1. ||M_q||_inf < 1 alone is not enough: +; the RLC has ||M_q||_inf = 1.045 and is stable. An M_q that is deliberately unstable (the sign of D flipped, a growing +; pair, the identity) fails the same check. +; - +; THE RUN. 200 steps: u = 1 V for steps 0..99, 0 V for steps 100..199. A = the state after step 99, B = after step 199, +; S = an order-sensitive fold of every state of every step (step by step, state by state), as M0 folds its one state. +; - +; THE CIRCUITS. (1) RC, R = 1 kOhm, C = 1 nF, dt = 50 ns: M0. (2) series RLC driven by a step, R = 100 Ohm, L = 10 uH, +; C = 1 nF, dt = 5 ns; the period is 628 ns and the damping ratio 0.5, so the response rings; the unknowns are the node +; between R and L (algebraic: no capacitor, its D_aa is exactly -1), the capacitor voltage, and the inductor current * +; 100 Ohm. (3) a 3-node and (4) a 5-node RC ladder, R = 1 kOhm and C = 1 nF per section, dt = 50 ns. +; - +; THE ERROR. Each circuit's 200 samples against the f64 backward Euler run of the same method (x' = M x + N u with the +; unrounded M and N), the largest difference in LSB of 2^-24 over every state and step: +; RC 2.65 (bound 3), RLC 8.44 (9), ladder-3 9.45 (10), ladder-5 40.03 (41). +; It grows with 1 / (1 - rho) of the slowest pole (rho 0.9524, 0.9747, 0.9902, 0.9960): a rounding error in D is amplified +; by the number of steps the circuit remembers. The bounds are asserted from both sides, so a change that makes the cell +; better or worse must change the number written here. +; - +; THE RLC'S NUMERICAL DAMPING. Backward Euler is not exact. For the RLC the energy factor per step is 1 / (1 + R dt / L + +; dt^2 / LC) = 0.95012, the real circuit's is exp(-R dt / L) = 0.95123: the ringing decays faster than the physical circuit, +; and over the 100 steps of the charge phase its amplitude is 0.943 of the true one. +; - +; WHAT A PASS SAYS. The integer constants of this file are the Q.24 roundings of the backward Euler matrices of the four +; element tables; the unsigned sign-magnitude step runs 200 steps within the stated LSB of the f64 solution of the same +; method; the RC circuit gives M0's A = 16649635, B = 126610 and S = 0x77EDFF0295C3FE9E; the quantized matrix of every +; circuit is stable by the K = 256 power test. The fold S of each circuit is pinned, so a board running RTL generated from +; the step can be compared with it. +; WHAT IT DOES NOT SAY. Nothing about a board: no bitstream, no LUT, FF, DSP or Fmax, no UART line (the board tops and +; their receipts are M1b, #10496). Nothing about a netlist: circuits are element tables written by hand (the reader is M2). +; Nothing about a nonlinear device (M3), a corner, a PDK, ngspice or Symbiotic EDA's evaluator. Nothing about the continuous +; circuit: only about backward Euler at this dt. Nothing about power. Nothing about another dt, another circuit or a driven +; node that also has a capacitor or an inductor on it (table_ok refuses those). Linear, fixed step, one input. +; ASCII only, English only. +; phi^2 + 1/phi^2 = 3 | TRINITY + +module AnalogLinear; + +pub const KIND : str = "analog-cell"; +pub const ID : str = "analog/linear"; +pub const REFERENCE : str = "mirekez/cpphdl analog/examples/rc_filter, M0 specs/fpga/rc_filter_ax7203.t27"; + +; --- the numbers -------------------------------------------------------------------------------- +/// Bit 31 of a state, a constant or an input: set when negative. +pub const NEG : u32 = 0x80000000; +/// Bits 0..30: the magnitude. +pub const MAG : u32 = 0x7FFFFFFF; +/// 1 V in units of 2^-24 V. +pub const ONE : u32 = 16777216; +/// Half an LSB of the product's 2^-48 scale, for rounding half away from zero. +pub const HALF_LSB : u64 = 8388608; +/// Raw units per volt, as f64. +pub const SCALE : f64 = 16777216.0; +/// A magnitude at or above this leaves less than 5 bits of headroom; no run here reaches it. +pub const HEADROOM : u32 = 67108864; + +; --- the sizes ---------------------------------------------------------------------------------- +/// The most unknowns (free nodes plus inductors) a table may have. +pub const NMAX : usize = 5; +/// Row stride of the work matrix [A | C/dt | B]: NMAX + NMAX + 1. +pub const WIDE : usize = 11; +/// NMAX * WIDE. +pub const AUG : usize = 55; +/// M (NMAX * NMAX) then N (NMAX), the dense result of derive. +pub const MN_LEN : usize = 30; + +; --- the run ------------------------------------------------------------------------------------ +pub const STEPS : u32 = 200; +pub const CHARGE_STEPS : u32 = 100; +/// ||M^K|| is taken at K = 2^8 = 256. +pub const K_SQUARINGS : u32 = 8; +/// A derived value this close to a rounding tie (0.5 LSB) could round the other way on another machine. +pub const TIE_GAP : f64 = 0.0001; + +; --- the elements ------------------------------------------------------------------------------- +pub const EL_R : u8 = 1; +pub const EL_C : u8 = 2; +pub const EL_L : u8 = 3; +/// The one voltage source: plus terminal NODE_IN, minus terminal NODE_GND; its value is the run's input u. +pub const EL_V : u8 = 4; +pub const NODE_GND : u8 = 0; +/// The driven node (the V element's plus terminal). Free nodes are 1, 2, 3 ... +pub const NODE_IN : u8 = 15; + +; --- circuit 1: the RC low-pass of M0 ------------------------------------------------------------ +pub const RC_N : usize = 1; +pub const RC_NE : usize = 3; +pub const RC_KIND : [3]u8 = [EL_V, EL_R, EL_C]; +pub const RC_A : [3]u8 = [NODE_IN, NODE_IN, 1]; +pub const RC_B : [3]u8 = [NODE_GND, 1, NODE_GND]; +pub const RC_VAL : [3]f64 = [0.0, 1000.0, 1.0e-9]; +pub const RC_DT : f64 = 5.0e-8; +pub const RC_ZS : f64 = 1.0; +/// D = M - I at 2^-24: -1/21. +pub const RC_D : [1]u32 = [NEG | 798915]; +/// N at 2^-24: 1/21. (`** 1`: a lone literal in a one-element module-level array lowers to Zig's array TYPE `[798915]`.) +pub const RC_NV : [1]u32 = [798915] ** 1; +pub const RC_EXP_S : u64 = 0x77EDFF0295C3FE9E; +pub const RC_MAX_ERR_LSB : u32 = 3; + +; --- circuit 2: series RLC, unknowns [node between R and L, capacitor voltage, inductor current * ZS] ---- +pub const RLC_N : usize = 3; +pub const RLC_NE : usize = 4; +pub const RLC_KIND : [4]u8 = [EL_V, EL_R, EL_L, EL_C]; +pub const RLC_A : [4]u8 = [NODE_IN, NODE_IN, 1, 2]; +pub const RLC_B : [4]u8 = [NODE_GND, 1, 2, NODE_GND]; +pub const RLC_VAL : [4]f64 = [0.0, 100.0, 1.0e-5, 1.0e-9]; +pub const RLC_DT : f64 = 5.0e-9; +/// sqrt(L / C) = 100 Ohm. +pub const RLC_ZS : f64 = 100.0; +pub const RLC_D : [9]u32 = [NEG | 16777216, 797017, NEG | 15940348, 0, NEG | 39851, 797017, 0, NEG | 797017, NEG | 836868]; +pub const RLC_NV : [3]u32 = [15980199, 39851, 797017]; +pub const RLC_EXP_S : u64 = 0x84BE2C7F9115D08A; +pub const RLC_MAX_ERR_LSB : u32 = 9; + +; --- circuit 3: 3-node RC ladder ------------------------------------------------------------------- +pub const LAD3_N : usize = 3; +pub const LAD3_NE : usize = 7; +pub const LAD3_KIND : [7]u8 = [EL_V, EL_R, EL_R, EL_R, EL_C, EL_C, EL_C]; +pub const LAD3_A : [7]u8 = [NODE_IN, NODE_IN, 1, 2, 1, 2, 3]; +pub const LAD3_B : [7]u8 = [NODE_GND, 1, 2, 3, NODE_GND, NODE_GND, NODE_GND]; +pub const LAD3_VAL : [7]f64 = [0.0, 1000.0, 1000.0, 1000.0, 1.0e-9, 1.0e-9, 1.0e-9]; +pub const LAD3_DT : f64 = 5.0e-8; +pub const LAD3_ZS : f64 = 1.0; +pub const LAD3_D : [9]u32 = [NEG | 1493555, 696219, 33153, 696219, NEG | 1460402, 729372, 33153, 729372, NEG | 764183]; +pub const LAD3_NV : [3]u32 = [764183, 34811, 1658]; +pub const LAD3_EXP_S : u64 = 0xD345ED15E5B270C7; +pub const LAD3_MAX_ERR_LSB : u32 = 10; + +; --- circuit 4: 5-node RC ladder ------------------------------------------------------------------- +pub const LAD5_N : usize = 5; +pub const LAD5_NE : usize = 11; +pub const LAD5_KIND : [11]u8 = [EL_V, EL_R, EL_R, EL_R, EL_R, EL_R, EL_C, EL_C, EL_C, EL_C, EL_C]; +pub const LAD5_A : [11]u8 = [NODE_IN, NODE_IN, 1, 2, 3, 4, 1, 2, 3, 4, 5]; +pub const LAD5_B : [11]u8 = [NODE_GND, 1, 2, 3, 4, 5, NODE_GND, NODE_GND, NODE_GND, NODE_GND, NODE_GND]; +pub const LAD5_VAL : [11]f64 = [0.0, 1000.0, 1000.0, 1000.0, 1000.0, 1000.0, 1.0e-9, 1.0e-9, 1.0e-9, 1.0e-9, 1.0e-9]; +pub const LAD5_DT : f64 = 5.0e-8; +pub const LAD5_ZS : f64 = 1.0; +pub const LAD5_D : [25]u32 = [ + NEG | 1493558, 696153, 31709, 1444, 69, + 696153, NEG | 1461849, 697597, 31778, 1513, + 31709, 697597, NEG | 1461780, 697666, 33222, + 1444, 31778, 697666, NEG | 1460336, 729375, + 69, 1513, 33222, 729375, NEG | 764183 +]; +pub const LAD5_NV : [5]u32 = [764183, 34808, 1585, 72, 3]; +pub const LAD5_EXP_S : u64 = 0x29439851A3F5D9EE; +pub const LAD5_MAX_ERR_LSB : u32 = 41; + +; --- the step: unsigned, sign-magnitude ---------------------------------------------------------- +; The scalar fns from vin_at to fold take and return bool and unsigned words only: they are what a board top reaches +; (`t27b verilog` refuses signed types and lowers only the fns the hardware reaches). step_row and step run them over +; the tables as slices, which is how the tests and the f64 comparison drive them. + +/// The input source at step n: 1 V while charging, then 0 V. +pub fn vin_at(n: u32) -> u32 { + if (n < CHARGE_STEPS) { return ONE; } + return 0; +} + +/// The word for magnitude m and direction neg; zero is never negative. +pub fn pack(m: u32, neg: bool) -> u32 { + if (m == 0) { return 0; } + if (neg) { return m | NEG; } + return m; +} + +/// Bit 31: is the word negative? +pub fn is_neg(p: u32) -> bool { + return (p & NEG) != 0; +} + +/// |c| * |x| as a u64. +pub fn prod(c: u32, x: u32) -> u64 { + return ((c & MAG) as u64) * ((x & MAG) as u64); +} + +/// Does c * x have a negative sign (the two directions differ)? +pub fn opposite(c: u32, x: u32) -> bool { + return ((c ^ x) & NEG) != 0; +} + +/// The positive-terms register after one more term c * x. +pub fn pos_after(pos: u64, c: u32, x: u32) -> u64 { + if (opposite(c, x)) { return pos; } + return pos + prod(c, x); +} + +/// The negative-terms register after one more term c * x. +pub fn neg_after(neg: u64, c: u32, x: u32) -> u64 { + if (opposite(c, x)) { return neg + prod(c, x); } + return neg; +} + +/// A product-scale magnitude (2^-48) back to 2^-24, rounded half away from zero. +pub fn round_back(t: u64) -> u32 { + return ((t + HALF_LSB) >> 24) as u32; +} + +/// The rounded sum pos - neg as a word: one subtraction, one rounding. +pub fn delta_of(pos: u64, neg: u64) -> u32 { + if (pos >= neg) { return pack(round_back(pos - neg), false); } + return pack(round_back(neg - pos), true); +} + +/// x + dl, both words. +pub fn add_sm(x: u32, dl: u32) -> u32 { + const mx : u32 = x & MAG; + const md : u32 = dl & MAG; + if (!opposite(x, dl)) { return pack(mx + md, is_neg(x)); } + if (mx >= md) { return pack(mx - md, is_neg(x)); } + return pack(md - mx, is_neg(dl)); +} + +/// The new state from the old one and the two registers. +pub fn next_x(x: u32, pos: u64, neg: u64) -> u32 { + return add_sm(x, delta_of(pos, neg)); +} + +/// The fold after one more sample: rotate left by one, xor the sample in. Order matters. M0's fold, on words. +pub fn fold(sum: u64, x: u32) -> u64 { + return ((sum << 1) | (sum >> 63)) ^ (x as u64); +} + +/// New state i of an n-state circuit: x_i + rnd(sum_j D_ij x_j + N_i u). d is n * n row-major. +pub fn step_row(i: usize, n: usize, d: []const u32, nv: []const u32, x: []const u32, u: u32) -> u32 { + var pos : u64 = 0; + var neg : u64 = 0; + var j : usize = 0; + while (j < n) { + pos = pos_after(pos, d[i * n + j], x[j]); + neg = neg_after(neg, d[i * n + j], x[j]); + j = j + 1; + } + pos = pos_after(pos, nv[i], u); + neg = neg_after(neg, nv[i], u); + return next_x(x[i], pos, neg); +} + +/// All n new states, from the old ones, into out. +pub fn step(n: usize, d: []const u32, nv: []const u32, x: []const u32, u: u32, out: []u32) { + var i : usize = 0; + while (i < n) { + out[i] = step_row(i, n, d, nv, x, u); + i = i + 1; + } +} + +; --- the f64 side: references that only the tests reach ---------------------------------------------- +; A module-level f64 fn is carried to Verilog only when the hardware reaches it (verilog_plan.t27, "the use +; closure"); nothing below is reached by a step, so a board top that uses this file takes the section above only. + +/// |x|. +pub fn fabs64(x: f64) -> f64 { + if (x < 0.0) { return 0.0 - x; } + return x; +} + +/// A word as a real number of volts. +pub fn sm_f64(p: u32) -> f64 { + const v : f64 = ((p & MAG) as f64) / SCALE; + if (is_neg(p)) { return 0.0 - v; } + return v; +} + +/// round(v * 2^24) as a word. +pub fn q_sm(v: f64) -> u32 { + return pack((fabs64(v) * SCALE + 0.5) as u32, v < 0.0); +} + +/// How far v * 2^24 is from a rounding tie (0.5 above an integer); 0.5 is the farthest. +pub fn tie_gap(v: f64) -> f64 { + const x : f64 = fabs64(v) * SCALE; + const f : f64 = x - ((x as u64) as f64); + return fabs64(f - 0.5); +} + +; --- the element table: nodes, unknowns, well-formedness --------------------------------------------- + +/// Is the node a free node (neither ground nor the driven node)? +pub fn is_free(node: u8) -> bool { + return node != NODE_GND and node != NODE_IN; +} + +/// The row of a free node. +pub fn ix(node: u8) -> usize { + return (node as usize) - 1; +} + +/// The highest free node id in the table. +pub fn node_count(a: []const u8, b: []const u8, ne: usize) -> usize { + var top : usize = 0; + var e : usize = 0; + while (e < ne) { + if (is_free(a[e]) and (a[e] as usize) > top) { top = a[e] as usize; } + if (is_free(b[e]) and (b[e] as usize) > top) { top = b[e] as usize; } + e = e + 1; + } + return top; +} + +/// How many inductors the table has. +pub fn ind_count(kind: []const u8, ne: usize) -> usize { + var cnt : usize = 0; + var e : usize = 0; + while (e < ne) { + if (kind[e] == EL_L) { cnt = cnt + 1; } + e = e + 1; + } + return cnt; +} + +/// Free nodes plus inductors. +pub fn unknowns(kind: []const u8, a: []const u8, b: []const u8, ne: usize) -> usize { + return node_count(a, b, ne) + ind_count(kind, ne); +} + +/// A table this cell can stamp: one V from NODE_IN to ground; R, C, L with a positive value and two different nodes; +/// no C or L on the driven node; every node 1..top touched; at most NMAX unknowns. +pub fn table_ok(kind: []const u8, a: []const u8, b: []const u8, val: []const f64, ne: usize) -> bool { + var sources : usize = 0; + var ok : bool = true; + var e : usize = 0; + while (e < ne) { + const k : u8 = kind[e]; + if (k == EL_V) { + sources = sources + 1; + if (a[e] != NODE_IN or b[e] != NODE_GND) { ok = false; } + } else if (k == EL_R or k == EL_C or k == EL_L) { + if (val[e] <= 0.0) { ok = false; } + if (a[e] == b[e]) { ok = false; } + if (k != EL_R and (a[e] == NODE_IN or b[e] == NODE_IN)) { ok = false; } + } else { + ok = false; + } + e = e + 1; + } + const top : usize = node_count(a, b, ne); + var id : usize = 1; + while (id <= top) { + var seen : bool = false; + e = 0; + while (e < ne) { + if ((a[e] as usize) == id or (b[e] as usize) == id) { seen = true; } + e = e + 1; + } + if (!seen) { ok = false; } + id = id + 1; + } + return ok and sources == 1 and unknowns(kind, a, b, ne) <= NMAX; +} + +; --- stamping: A, C/dt and B of the work matrix [A | C/dt | B] ---------------------------------------- + +/// w[idx] += y. +pub fn bump(w: []f64, idx: usize, y: f64) { + w[idx] = w[idx] + y; +} + +/// An admittance y between nodes a and b into the n columns of the work matrix that start at off. +pub fn stamp_pair(w: []f64, off: usize, a: u8, b: u8, y: f64) { + if (is_free(a)) { bump(w, ix(a) * WIDE + off + ix(a), y); } + if (is_free(b)) { bump(w, ix(b) * WIDE + off + ix(b), y); } + if (is_free(a) and is_free(b)) { + bump(w, ix(a) * WIDE + off + ix(b), 0.0 - y); + bump(w, ix(b) * WIDE + off + ix(a), 0.0 - y); + } +} + +/// The driven node's side of an admittance y: its column moves to the right-hand side (column col) as y * u. +pub fn stamp_src(w: []f64, col: usize, a: u8, b: u8, y: f64) { + if (is_free(a) and b == NODE_IN) { bump(w, ix(a) * WIDE + col, y); } + if (is_free(b) and a == NODE_IN) { bump(w, ix(b) * WIDE + col, y); } +} + +/// The work matrix [A | C/dt | B] of a table: R stamps G, C stamps C/dt into both A and the middle block, L adds a +/// branch current (stored times ZS) with KCL +- i at its nodes and the row (L / (dt ZS)) x_k' - v_a' + v_b' = (L / (dt ZS)) x_k. +pub fn stamp(kind: []const u8, a: []const u8, b: []const u8, val: []const f64, ne: usize, dt: f64, zs: f64) -> [55]f64 { + var w : [55]f64 = [0.0] ** 55; + const nn : usize = node_count(a, b, ne); + const n : usize = nn + ind_count(kind, ne); + var l : usize = 0; + var e : usize = 0; + while (e < ne) { + const k : u8 = kind[e]; + if (k == EL_R) { + stamp_pair(&w, 0, a[e], b[e], 1.0 / val[e]); + stamp_src(&w, 2 * n, a[e], b[e], 1.0 / val[e]); + } else if (k == EL_C) { + stamp_pair(&w, 0, a[e], b[e], val[e] / dt); + stamp_pair(&w, n, a[e], b[e], val[e] / dt); + } else if (k == EL_L) { + const ki : usize = nn + l; + const yl : f64 = val[e] / (dt * zs); + if (is_free(a[e])) { + bump(&w, ix(a[e]) * WIDE + ki, 1.0 / zs); + bump(&w, ki * WIDE + ix(a[e]), 0.0 - 1.0); + } + if (is_free(b[e])) { + bump(&w, ix(b[e]) * WIDE + ki, 0.0 - 1.0 / zs); + bump(&w, ki * WIDE + ix(b[e]), 1.0); + } + bump(&w, ki * WIDE + ki, yl); + bump(&w, ki * WIDE + n + ki, yl); + if (a[e] == NODE_IN) { bump(&w, ki * WIDE + 2 * n, 1.0); } + if (b[e] == NODE_IN) { bump(&w, ki * WIDE + 2 * n, 0.0 - 1.0); } + l = l + 1; + } + e = e + 1; + } + return w; +} + +; --- solving: Gauss-Jordan elimination with partial pivoting ------------------------------------------ + +/// Reduce the first n columns of the work matrix to the identity; the rest of each row becomes A^-1 [C/dt | B]. +pub fn eliminate(aug: [55]f64, n: usize) -> [55]f64 { + var w : [55]f64 = aug; + var c : usize = 0; + var r : usize = 0; + var p : usize = 0; + var k : usize = 0; + while (c < n) { + p = c; + r = c + 1; + while (r < n) { + if (fabs64(w[r * WIDE + c]) > fabs64(w[p * WIDE + c])) { p = r; } + r = r + 1; + } + if (p != c) { + k = 0; + while (k < WIDE) { + const t : f64 = w[c * WIDE + k]; + w[c * WIDE + k] = w[p * WIDE + k]; + w[p * WIDE + k] = t; + k = k + 1; + } + } + const piv : f64 = w[c * WIDE + c]; + k = 0; + while (k < WIDE) { + w[c * WIDE + k] = w[c * WIDE + k] / piv; + k = k + 1; + } + r = 0; + while (r < n) { + if (r != c) { + const f : f64 = w[r * WIDE + c]; + k = 0; + while (k < WIDE) { + w[r * WIDE + k] = w[r * WIDE + k] - f * w[c * WIDE + k]; + k = k + 1; + } + } + r = r + 1; + } + c = c + 1; + } + return w; +} + +/// The dense result of a solved work matrix: M row-major in [0, n*n), N in [n*n, n*n + n). +pub fn mn_of(w: [55]f64, n: usize) -> [30]f64 { + var o : [30]f64 = [0.0] ** 30; + var i : usize = 0; + var j : usize = 0; + while (i < n) { + j = 0; + while (j < n) { + o[i * n + j] = w[i * WIDE + n + j]; + j = j + 1; + } + o[n * n + i] = w[i * WIDE + 2 * n]; + i = i + 1; + } + return o; +} + +/// M and N of a table, in f64, from nothing but the table and dt. +pub fn derive(kind: []const u8, a: []const u8, b: []const u8, val: []const f64, ne: usize, dt: f64, zs: f64) -> [30]f64 { + const n : usize = unknowns(kind, a, b, ne); + return mn_of(eliminate(stamp(kind, a, b, val, ne, dt, zs), n), n); +} + +/// The largest |A X - R| over the n + 1 solved columns of the work matrix: a solve that did not solve shows here. +pub fn residual(a0: [55]f64, sol: [55]f64, n: usize) -> f64 { + var worst : f64 = 0.0; + var i : usize = 0; + var c : usize = 0; + var k : usize = 0; + while (i < n) { + c = 0; + while (c <= n) { + var acc : f64 = 0.0 - a0[i * WIDE + n + c]; + k = 0; + while (k < n) { + acc = acc + a0[i * WIDE + k] * sol[k * WIDE + n + c]; + k = k + 1; + } + if (fabs64(acc) > worst) { worst = fabs64(acc); } + c = c + 1; + } + i = i + 1; + } + return worst; +} + +; --- the stored integers against the derived reals ------------------------------------------------------ + +/// How many stored constants differ from round(derived value * 2^24): D_ij against M_ij - delta_ij, N_i against N_i. +pub fn mismatches(n: usize, mn: []const f64, d: []const u32, nv: []const u32) -> usize { + var bad : usize = 0; + var i : usize = 0; + var j : usize = 0; + while (i < n) { + j = 0; + while (j < n) { + var m : f64 = mn[i * n + j]; + if (i == j) { m = m - 1.0; } + if (q_sm(m) != d[i * n + j]) { bad = bad + 1; } + j = j + 1; + } + if (q_sm(mn[n * n + i]) != nv[i]) { bad = bad + 1; } + i = i + 1; + } + return bad; +} + +/// Is every derived value at least TIE_GAP away from a rounding tie? +pub fn ties_clear(n: usize, mn: []const f64) -> bool { + var ok : bool = true; + var i : usize = 0; + var j : usize = 0; + while (i < n) { + j = 0; + while (j < n) { + var m : f64 = mn[i * n + j]; + if (i == j) { m = m - 1.0; } + if (tie_gap(m) <= TIE_GAP) { ok = false; } + j = j + 1; + } + if (tie_gap(mn[n * n + i]) <= TIE_GAP) { ok = false; } + i = i + 1; + } + return ok; +} + +; --- the run against the f64 solution of the same method ------------------------------------------------ + +pub struct RunOut { sig: u64, worst: f64, peak: u32 } + +/// k steps from rest with the integer step and, beside it, the f64 step x' = M x + N u. The final integer state lands in +/// x. sig folds every state of every step; worst is the largest |integer - f64| in LSB; peak is the largest magnitude. +pub fn run(n: usize, d: []const u32, nv: []const u32, mn: []const f64, k: u32, x: []u32) -> RunOut { + var y : [5]u32 = [0] ** 5; + var f : [5]f64 = [0.0] ** 5; + var g : [5]f64 = [0.0] ** 5; + var s : u64 = 0; + var worst : f64 = 0.0; + var peak : u32 = 0; + var t : u32 = 0; + var i : usize = 0; + var j : usize = 0; + while (i < n) { + x[i] = 0; + i = i + 1; + } + while (t < k) { + const uu : u32 = vin_at(t); + const uf : f64 = (uu as f64) / SCALE; + step(n, d, nv, x, uu, &y); + i = 0; + while (i < n) { + var acc : f64 = mn[n * n + i] * uf; + j = 0; + while (j < n) { + acc = acc + mn[i * n + j] * f[j]; + j = j + 1; + } + g[i] = acc; + i = i + 1; + } + i = 0; + while (i < n) { + x[i] = y[i]; + f[i] = g[i]; + s = fold(s, y[i]); + const err : f64 = fabs64(sm_f64(y[i]) - g[i]) * SCALE; + if (err > worst) { worst = err; } + if ((y[i] & MAG) > peak) { peak = y[i] & MAG; } + i = i + 1; + } + t = t + 1; + } + return RunOut { sig: s, worst: worst, peak: peak }; +} + +; --- stability: ||M_q^256||_inf < 1 --------------------------------------------------------------------- + +/// M_q = I + D / 2^24 as an n * n row-major f64 matrix; D / 2^24 is exact in f64. +pub fn mq_of(n: usize, d: []const u32) -> [25]f64 { + var o : [25]f64 = [0.0] ** 25; + var i : usize = 0; + var j : usize = 0; + while (i < n) { + j = 0; + while (j < n) { + o[i * n + j] = sm_f64(d[i * n + j]); + if (i == j) { o[i * n + j] = o[i * n + j] + 1.0; } + j = j + 1; + } + i = i + 1; + } + return o; +} + +/// out = a * b, n * n row-major; out must not be a or b. +pub fn matmul(n: usize, a: []const f64, b: []const f64, out: []f64) { + var i : usize = 0; + var j : usize = 0; + var k : usize = 0; + while (i < n) { + j = 0; + while (j < n) { + var acc : f64 = 0.0; + k = 0; + while (k < n) { + acc = acc + a[i * n + k] * b[k * n + j]; + k = k + 1; + } + out[i * n + j] = acc; + j = j + 1; + } + i = i + 1; + } +} + +/// The infinity norm: the largest row sum of absolute values. +pub fn norm_inf(n: usize, a: []const f64) -> f64 { + var best : f64 = 0.0; + var i : usize = 0; + var j : usize = 0; + while (i < n) { + var row : f64 = 0.0; + j = 0; + while (j < n) { + row = row + fabs64(a[i * n + j]); + j = j + 1; + } + if (row > best) { best = row; } + i = i + 1; + } + return best; +} + +/// ||m^(2^squarings)||_inf by repeated squaring. +pub fn power_norm(n: usize, m: []const f64, squarings: u32) -> f64 { + var p : [25]f64 = [0.0] ** 25; + var q : [25]f64 = [0.0] ** 25; + var i : usize = 0; + var s : u32 = 0; + while (i < n * n) { + p[i] = m[i]; + i = i + 1; + } + while (s < squarings) { + matmul(n, &p, &p, &q); + i = 0; + while (i < n * n) { + p[i] = q[i]; + i = i + 1; + } + s = s + 1; + } + return norm_inf(n, &p); +} + +/// The spectral radius of m is below 1 when ||m^256||_inf is. +pub fn stable(n: usize, m: []const f64) -> bool { + return power_norm(n, m, K_SQUARINGS) < 1.0; +} + +; --- the tests ---------------------------------------------------------------------------------- + +test the_solver_pivots_and_inverts { + // a system with a zero on the diagonal: A = [[0, 1], [1, 0]], C/dt = I, B = [5, 7] + var w : [55]f64 = [0.0] ** 55; + w[1] = 1.0; + w[WIDE] = 1.0; + w[2] = 1.0; + w[WIDE + 3] = 1.0; + w[4] = 5.0; + w[WIDE + 4] = 7.0; + const sol : [55]f64 = eliminate(w, 2); + const mn : [30]f64 = mn_of(sol, 2); + assert mn[0] == 0.0; + assert mn[1] == 1.0; + assert mn[2] == 1.0; + assert mn[3] == 0.0; + assert mn[4] == 7.0; + assert mn[5] == 5.0; + assert residual(w, sol, 2) < 0.000000000001; + assert_eq(AUG, NMAX * WIDE); + assert_eq(WIDE, 2 * NMAX + 1); + assert_eq(MN_LEN, NMAX * NMAX + NMAX); +} + +test the_tables_are_well_formed_and_stamp_to_a_solvable_system { + assert(table_ok(&RC_KIND, &RC_A, &RC_B, &RC_VAL, RC_NE)); + assert(table_ok(&RLC_KIND, &RLC_A, &RLC_B, &RLC_VAL, RLC_NE)); + assert(table_ok(&LAD3_KIND, &LAD3_A, &LAD3_B, &LAD3_VAL, LAD3_NE)); + assert(table_ok(&LAD5_KIND, &LAD5_A, &LAD5_B, &LAD5_VAL, LAD5_NE)); + assert_eq(unknowns(&RC_KIND, &RC_A, &RC_B, RC_NE), RC_N); + assert_eq(unknowns(&RLC_KIND, &RLC_A, &RLC_B, RLC_NE), RLC_N); + assert_eq(unknowns(&LAD3_KIND, &LAD3_A, &LAD3_B, LAD3_NE), LAD3_N); + assert_eq(unknowns(&LAD5_KIND, &LAD5_A, &LAD5_B, LAD5_NE), LAD5_N); + assert_eq(ind_count(&RLC_KIND, RLC_NE), 1); + assert_eq(ind_count(&LAD5_KIND, LAD5_NE), 0); + const w1 : [55]f64 = stamp(&RC_KIND, &RC_A, &RC_B, &RC_VAL, RC_NE, RC_DT, RC_ZS); + const s1 : [55]f64 = eliminate(w1, RC_N); + assert residual(w1, s1, RC_N) < 0.000000000001; + const w2 : [55]f64 = stamp(&RLC_KIND, &RLC_A, &RLC_B, &RLC_VAL, RLC_NE, RLC_DT, RLC_ZS); + const s2 : [55]f64 = eliminate(w2, RLC_N); + assert residual(w2, s2, RLC_N) < 0.000000000001; + const w3 : [55]f64 = stamp(&LAD3_KIND, &LAD3_A, &LAD3_B, &LAD3_VAL, LAD3_NE, LAD3_DT, LAD3_ZS); + const s3 : [55]f64 = eliminate(w3, LAD3_N); + assert residual(w3, s3, LAD3_N) < 0.000000000001; + const w5 : [55]f64 = stamp(&LAD5_KIND, &LAD5_A, &LAD5_B, &LAD5_VAL, LAD5_NE, LAD5_DT, LAD5_ZS); + const s5 : [55]f64 = eliminate(w5, LAD5_N); + assert residual(w5, s5, LAD5_N) < 0.000000000001; + // a table with a capacitor on the driven node, or two sources, or a gap in the node numbers, is refused + const bad_k : [4]u8 = [EL_V, EL_R, EL_C, EL_C]; + const bad_a : [4]u8 = [NODE_IN, NODE_IN, 1, NODE_IN]; + const bad_b : [4]u8 = [NODE_GND, 1, NODE_GND, 1]; + const bad_v : [4]f64 = [0.0, 1000.0, 1.0e-9, 1.0e-9]; + assert(!table_ok(&bad_k, &bad_a, &bad_b, &bad_v, 4)); + const two_k : [4]u8 = [EL_V, EL_V, EL_R, EL_C]; + const two_a : [4]u8 = [NODE_IN, NODE_IN, NODE_IN, 1]; + const two_b : [4]u8 = [NODE_GND, NODE_GND, 1, NODE_GND]; + assert(!table_ok(&two_k, &two_a, &two_b, &bad_v, 4)); + const gap_a : [3]u8 = [NODE_IN, NODE_IN, 3]; + assert(!table_ok(&RC_KIND, &gap_a, &RC_B, &RC_VAL, 3)); +} + +test the_rc_constants_are_the_rounding_of_the_derived_ones { + const mn : [30]f64 = derive(&RC_KIND, &RC_A, &RC_B, &RC_VAL, RC_NE, RC_DT, RC_ZS); + assert_eq(mismatches(RC_N, &mn, &RC_D, &RC_NV), 0); + assert(ties_clear(RC_N, &mn)); + // M0's coefficient: g = dt / (tau + dt) = 1/21, so M = 20/21 and N = 1/21 + const dm : f64 = mn[0] - 20.0 / 21.0; + const dn : f64 = mn[1] - 1.0 / 21.0; + assert dm * dm < 0.000000000000001 * 0.000000000000001; + assert dn * dn < 0.000000000000001 * 0.000000000000001; + assert_eq(RC_D[0], NEG | 798915); + assert_eq(RC_NV[0], 798915); + assert_eq(798915 * 21, ONE - 1); +} + +test the_rlc_constants_are_the_rounding_of_the_derived_ones { + const mn : [30]f64 = derive(&RLC_KIND, &RLC_A, &RLC_B, &RLC_VAL, RLC_NE, RLC_DT, RLC_ZS); + assert_eq(mismatches(RLC_N, &mn, &RLC_D, &RLC_NV), 0); + assert(ties_clear(RLC_N, &mn)); + // the node between R and L has no capacitor: its column of M is zero, its D_aa is exactly -1 + assert_eq(RLC_D[0], NEG | ONE); + assert mn[0] == 0.0; + assert_eq(RLC_D[3], 0); + assert_eq(RLC_D[6], 0); + // ZS is the characteristic impedance sqrt(L / C): ZS^2 C / L = 1 + const zz : f64 = RLC_ZS * RLC_ZS * RLC_VAL[3] / RLC_VAL[2] - 1.0; + assert zz * zz < 0.000000000001; + // and it keeps the current's rows near the voltages': no stored magnitude above 1 (2^24) + var i : usize = 0; + var top : u32 = 0; + while (i < 9) { + if ((RLC_D[i] & MAG) > top) { top = RLC_D[i] & MAG; } + i = i + 1; + } + i = 0; + while (i < 3) { + if ((RLC_NV[i] & MAG) > top) { top = RLC_NV[i] & MAG; } + i = i + 1; + } + assert_eq(top, ONE); +} + +test the_ladder3_constants_are_the_rounding_of_the_derived_ones { + const mn : [30]f64 = derive(&LAD3_KIND, &LAD3_A, &LAD3_B, &LAD3_VAL, LAD3_NE, LAD3_DT, LAD3_ZS); + assert_eq(mismatches(LAD3_N, &mn, &LAD3_D, &LAD3_NV), 0); + assert(ties_clear(LAD3_N, &mn)); + // the ladder is symmetric: D_ij = D_ji + assert_eq(LAD3_D[1], LAD3_D[3]); + assert_eq(LAD3_D[2], LAD3_D[6]); + assert_eq(LAD3_D[5], LAD3_D[7]); +} + +test the_ladder5_constants_are_the_rounding_of_the_derived_ones { + const mn : [30]f64 = derive(&LAD5_KIND, &LAD5_A, &LAD5_B, &LAD5_VAL, LAD5_NE, LAD5_DT, LAD5_ZS); + assert_eq(mismatches(LAD5_N, &mn, &LAD5_D, &LAD5_NV), 0); + assert(ties_clear(LAD5_N, &mn)); + assert_eq(LAD5_D[1], LAD5_D[5]); + assert_eq(LAD5_D[3], LAD5_D[15]); + assert_eq(LAD5_D[9], LAD5_D[21]); + assert_eq(LAD5_D[24], NEG | 764183); + // the last section of both ladders sees one neighbour, the same as the 3-node ladder's: N_n and D_nn agree + assert_eq(LAD5_D[24], LAD3_D[8]); +} + +test the_step_rounds_half_away_from_zero_and_zero_is_never_negative { + assert_eq(delta_of(8388608, 0), 1); + assert_eq(delta_of(8388607, 0), 0); + assert_eq(delta_of(0, 8388608), NEG | 1); + assert_eq(delta_of(0, 8388607), 0); + assert_eq(delta_of(5, 5), 0); + assert_eq(delta_of(25165824, 8388608), 1); + assert_eq(add_sm(3, 4), 7); + assert_eq(add_sm(NEG | 3, NEG | 4), NEG | 7); + assert_eq(add_sm(3, NEG | 5), NEG | 2); + assert_eq(add_sm(NEG | 3, 5), 2); + assert_eq(add_sm(5, NEG | 5), 0); + assert_eq(add_sm(NEG | 5, 5), 0); + assert_eq(pack(0, true), 0); + assert_eq(pos_after(0, NEG | 2, NEG | 3), 6); + assert_eq(pos_after(0, 2, 3), 6); + assert_eq(pos_after(0, NEG | 2, 3), 0); + assert_eq(neg_after(0, NEG | 2, 3), 6); + assert_eq(neg_after(0, 2, NEG | 3), 6); + assert_eq(neg_after(0, 2, 3), 0); + assert_eq(prod(MAG, MAG), 4611686014132420609); + assert(is_neg(NEG | 1)); + assert(!is_neg(MAG)); + assert_eq(fold(0x8000000000000000, 0), 1); + assert fold(fold(0, 1), 2) != fold(fold(0, 2), 1); +} + +test the_rc_circuit_reproduces_m0_bit_for_bit { + // M0's pinned numbers (specs/fpga/rc_filter_ax7203.t27: EXP_A, EXP_B, EXP_S and its first samples) + const m0_a : u32 = 16649635; + const m0_b : u32 = 126610; + const m0_s : u64 = 0x77EDFF0295C3FE9E; + const mn : [30]f64 = derive(&RC_KIND, &RC_A, &RC_B, &RC_VAL, RC_NE, RC_DT, RC_ZS); + var x : [5]u32 = [0] ** 5; + const first : RunOut = run(RC_N, &RC_D, &RC_NV, &mn, 1, &x); + assert_eq(x[0], 798915); + var y : [5]u32 = [0] ** 5; + const second : RunOut = run(RC_N, &RC_D, &RC_NV, &mn, 2, &y); + assert_eq(y[0], 1559786); + var xa : [5]u32 = [0] ** 5; + const ra : RunOut = run(RC_N, &RC_D, &RC_NV, &mn, CHARGE_STEPS, &xa); + assert_eq(xa[0], m0_a); + var xb : [5]u32 = [0] ** 5; + const rb : RunOut = run(RC_N, &RC_D, &RC_NV, &mn, STEPS, &xb); + assert_eq(xb[0], m0_b); + assert_eq(rb.sig, m0_s); + assert_eq(rb.sig, RC_EXP_S); + // discharge from 1 V: M0 gives ONE - 798915 for one step + const sx : [1]u32 = [ONE]; + var so : [1]u32 = [0]; + step(1, &RC_D, &RC_NV, &sx, 0, &so); + assert_eq(so[0], ONE - 798915); + step(1, &RC_D, &RC_NV, &sx, ONE, &so); + assert_eq(so[0], ONE); + // a step with no input and no charge stays at zero + const zx : [1]u32 = [0]; + step(1, &RC_D, &RC_NV, &zx, 0, &so); + assert_eq(so[0], 0); + assert first.worst > 0.0; + assert second.worst > 0.0; + assert ra.peak < HEADROOM; +} + +test the_rc_charge_at_5_us_is_cpphdls_0_992396 { + const mn : [30]f64 = derive(&RC_KIND, &RC_A, &RC_B, &RC_VAL, RC_NE, RC_DT, RC_ZS); + var x : [5]u32 = [0] ** 5; + const r : RunOut = run(RC_N, &RC_D, &RC_NV, &mn, CHARGE_STEPS, &x); + const vv : f64 = sm_f64(x[0]); + // cpphdl's README: 1 - (1/1.05)^100 = 0.992396 at 5 us, against 1 - exp(-5) = 0.993262 + const dc : f64 = vv - 0.992396; + assert dc * dc < 0.0000005 * 0.0000005; + const da : f64 = vv - 0.9932620530009145; + assert da * da > 0.0008 * 0.0008; + assert da * da < 0.0009 * 0.0009; + assert r.worst < 3.0; +} + +test each_circuit_stays_within_its_stated_error_of_f64_backward_euler { + // RC: the bound is M0's 3 LSB + const m1 : [30]f64 = derive(&RC_KIND, &RC_A, &RC_B, &RC_VAL, RC_NE, RC_DT, RC_ZS); + var x1 : [5]u32 = [0] ** 5; + const r1 : RunOut = run(RC_N, &RC_D, &RC_NV, &m1, STEPS, &x1); + assert_eq(r1.sig, RC_EXP_S); + assert r1.worst < (RC_MAX_ERR_LSB as f64); + assert r1.worst > (RC_MAX_ERR_LSB as f64) - 1.0; + assert r1.peak < HEADROOM; + // RLC: node voltages and the inductor current * ZS; the current changes sign + const m2 : [30]f64 = derive(&RLC_KIND, &RLC_A, &RLC_B, &RLC_VAL, RLC_NE, RLC_DT, RLC_ZS); + var x2 : [5]u32 = [0] ** 5; + const r2 : RunOut = run(RLC_N, &RLC_D, &RLC_NV, &m2, STEPS, &x2); + assert_eq(r2.sig, RLC_EXP_S); + assert r2.worst < (RLC_MAX_ERR_LSB as f64); + assert r2.worst > (RLC_MAX_ERR_LSB as f64) - 1.0; + assert r2.peak < HEADROOM; + assert_eq(x2[0], NEG | 1419706); + assert_eq(x2[1], NEG | 1265689); + assert_eq(x2[2], 1419706); + var x2a : [5]u32 = [0] ** 5; + const r2a : RunOut = run(RLC_N, &RLC_D, &RLC_NV, &m2, CHARGE_STEPS, &x2a); + assert_eq(x2a[0], 18100421); + assert_eq(x2a[1], 18050612); + assert_eq(x2a[2], NEG | 1323205); + assert(x2a[1] > ONE); + assert r2a.worst < (RLC_MAX_ERR_LSB as f64); + // ladder-3 + const m3 : [30]f64 = derive(&LAD3_KIND, &LAD3_A, &LAD3_B, &LAD3_VAL, LAD3_NE, LAD3_DT, LAD3_ZS); + var x3 : [5]u32 = [0] ** 5; + const r3 : RunOut = run(LAD3_N, &LAD3_D, &LAD3_NV, &m3, STEPS, &x3); + assert_eq(r3.sig, LAD3_EXP_S); + assert r3.worst < (LAD3_MAX_ERR_LSB as f64); + assert r3.worst > (LAD3_MAX_ERR_LSB as f64) - 1.0; + assert r3.peak < HEADROOM; + // ladder-5 + const m5 : [30]f64 = derive(&LAD5_KIND, &LAD5_A, &LAD5_B, &LAD5_VAL, LAD5_NE, LAD5_DT, LAD5_ZS); + var x5 : [5]u32 = [0] ** 5; + const r5 : RunOut = run(LAD5_N, &LAD5_D, &LAD5_NV, &m5, STEPS, &x5); + assert_eq(r5.sig, LAD5_EXP_S); + assert r5.worst < (LAD5_MAX_ERR_LSB as f64); + assert r5.worst > (LAD5_MAX_ERR_LSB as f64) - 1.0; + assert r5.peak < HEADROOM; +} + +test the_k_256_power_test_holds_for_every_circuit_and_fails_for_unstable_ones { + const q1 : [25]f64 = mq_of(RC_N, &RC_D); + const q2 : [25]f64 = mq_of(RLC_N, &RLC_D); + const q3 : [25]f64 = mq_of(LAD3_N, &LAD3_D); + const q5 : [25]f64 = mq_of(LAD5_N, &LAD5_D); + assert(stable(RC_N, &q1)); + assert(stable(RLC_N, &q2)); + assert(stable(LAD3_N, &q3)); + assert(stable(LAD5_N, &q5)); + assert power_norm(RC_N, &q1, K_SQUARINGS) < 0.00001; + assert power_norm(RLC_N, &q2, K_SQUARINGS) < 0.003; + assert power_norm(LAD3_N, &q3, K_SQUARINGS) < 0.1; + assert power_norm(LAD5_N, &q5, K_SQUARINGS) < 0.45; + // the plain norm is not enough: the underdamped RLC has ||M_q||_inf above 1 and is stable + assert norm_inf(RLC_N, &q2) > 1.0; + assert norm_inf(RLC_N, &q2) < 1.05; + // a sign error in D (M = 1 + 1/21) grows + const wrong : [1]u32 = [798915]; + const qw : [25]f64 = mq_of(1, &wrong); + assert(!stable(1, &qw)); + // M = I does not decay: a circuit with no loss is not certified + const ident : [1]u32 = [0]; + const qi : [25]f64 = mq_of(1, &ident); + assert(!stable(1, &qi)); + // a coupled pair with eigenvalue 1.4 (M = [[0.5, 0.9], [0.9, 0.5]]) grows + const grow : [4]u32 = [NEG | 8388608, 15099494, 15099494, NEG | 8388608]; + const qg : [25]f64 = mq_of(2, &grow); + assert(!stable(2, &qg)); + assert_eq(K_SQUARINGS, 8); +} + +test the_rlc_backward_euler_damps_the_ringing_faster_than_the_circuit { + const r : f64 = RLC_VAL[1]; + const l : f64 = RLC_VAL[2]; + const c : f64 = RLC_VAL[3]; + const rdl : f64 = r * RLC_DT / l; + const dr : f64 = rdl - 0.05; + assert dr * dr < 0.000000000001; + // the energy factor per step of backward Euler: 1 / (1 + R dt / L + dt^2 / L C) + const lam2 : f64 = 1.0 / (1.0 + rdl + RLC_DT * RLC_DT / (l * c)); + const exact : f64 = 0.951229424500714; + assert lam2 > 0.9501; + assert lam2 < 0.9502; + assert lam2 < exact; + // the same number from the derived matrix: its first column is zero, so the other two eigenvalues are those of the + // lower-right 2 x 2 block, whose determinant is |lambda|^2 for a complex pair + const mn : [30]f64 = derive(&RLC_KIND, &RLC_A, &RLC_B, &RLC_VAL, RLC_NE, RLC_DT, RLC_ZS); + const det : f64 = mn[4] * mn[8] - mn[5] * mn[7]; + const tr : f64 = mn[4] + mn[8]; + const dd : f64 = det - lam2; + assert dd * dd < 0.000000000001 * 0.000000000001; + assert tr * tr < 4.0 * det; + // after the 100 charging steps the amplitude is this fraction of the true one + var q : f64 = 1.0; + var k : u32 = 0; + while (k < 50) { + q = q * lam2 / exact; + k = k + 1; + } + assert q > 0.94; + assert q < 0.945; +} + +test a_conductance_on_the_wrong_node_is_caught_and_so_is_a_wrong_capacitor_or_inductor { + // the right tables derive the stored constants (zero mismatches, asserted above); a mis-stamp must not + const g3 : [30]f64 = derive(&LAD3_KIND, &LAD3_A, &LAD3_B, &LAD3_VAL, LAD3_NE, LAD3_DT, LAD3_ZS); + assert_eq(mismatches(LAD3_N, &g3, &LAD3_D, &LAD3_NV), 0); + // ladder-3: the R from node 2 to 3 stamped from node 1 to 3 + var wa : [7]u8 = LAD3_A; + wa[3] = 1; + assert(table_ok(&LAD3_KIND, &wa, &LAD3_B, &LAD3_VAL, LAD3_NE)); + const bad_r : [30]f64 = derive(&LAD3_KIND, &wa, &LAD3_B, &LAD3_VAL, LAD3_NE, LAD3_DT, LAD3_ZS); + assert(mismatches(LAD3_N, &bad_r, &LAD3_D, &LAD3_NV) > 0); + // ladder-3: the capacitor of node 3 on node 2 + var wc : [7]u8 = LAD3_A; + wc[6] = 2; + assert(table_ok(&LAD3_KIND, &wc, &LAD3_B, &LAD3_VAL, LAD3_NE)); + const bad_c : [30]f64 = derive(&LAD3_KIND, &wc, &LAD3_B, &LAD3_VAL, LAD3_NE, LAD3_DT, LAD3_ZS); + assert(mismatches(LAD3_N, &bad_c, &LAD3_D, &LAD3_NV) > 0); + // ladder-5: the R from node 3 to 4 stamped from node 2 to 4 + var w5 : [11]u8 = LAD5_A; + w5[4] = 2; + assert(table_ok(&LAD5_KIND, &w5, &LAD5_B, &LAD5_VAL, LAD5_NE)); + const bad_5 : [30]f64 = derive(&LAD5_KIND, &w5, &LAD5_B, &LAD5_VAL, LAD5_NE, LAD5_DT, LAD5_ZS); + assert(mismatches(LAD5_N, &bad_5, &LAD5_D, &LAD5_NV) > 0); + // RLC: the inductor with its terminals swapped + var la : [4]u8 = RLC_A; + var lb : [4]u8 = RLC_B; + la[2] = 2; + lb[2] = 1; + assert(table_ok(&RLC_KIND, &la, &lb, &RLC_VAL, RLC_NE)); + const bad_l : [30]f64 = derive(&RLC_KIND, &la, &lb, &RLC_VAL, RLC_NE, RLC_DT, RLC_ZS); + assert(mismatches(RLC_N, &bad_l, &RLC_D, &RLC_NV) > 0); + // RLC: the series R on the capacitor node instead of the inductor node + var rb : [4]u8 = RLC_B; + rb[1] = 2; + assert(table_ok(&RLC_KIND, &RLC_A, &rb, &RLC_VAL, RLC_NE)); + const bad_s : [30]f64 = derive(&RLC_KIND, &RLC_A, &rb, &RLC_VAL, RLC_NE, RLC_DT, RLC_ZS); + assert(mismatches(RLC_N, &bad_s, &RLC_D, &RLC_NV) > 0); + // a wrong value is as visible: R of the RC ladder-3 section doubled + var wv : [7]f64 = LAD3_VAL; + wv[2] = 2000.0; + const bad_v : [30]f64 = derive(&LAD3_KIND, &LAD3_A, &LAD3_B, &wv, LAD3_NE, LAD3_DT, LAD3_ZS); + assert(mismatches(LAD3_N, &bad_v, &LAD3_D, &LAD3_NV) > 0); +} + +test a_truncating_step_is_not_the_cell_and_a_wrong_dt_is_caught { + // M0's wrong rounding: dropping the half LSB lands 11 LSB short after 100 steps (16649624, not 16649635) + const g : u64 = 798915; + var v : u32 = 0; + var k : u32 = 0; + while (k < CHARGE_STEPS) { + v = v + ((g * ((ONE - v) as u64)) >> 24) as u32; + k = k + 1; + } + assert v != 16649635; + assert_eq(v, 16649624); + // a wrong dt (25 ns, not 50 ns) is a different coefficient, so the stored RC constants no longer match + const m_half : [30]f64 = derive(&RC_KIND, &RC_A, &RC_B, &RC_VAL, RC_NE, 2.5e-8, RC_ZS); + assert(mismatches(RC_N, &m_half, &RC_D, &RC_NV) > 0); +}