From d1dd9c6a8d65f4807c3669368ef13de565507da2 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 18 Oct 2026 23:39:40 +0700 Subject: [PATCH] analog: linear_zoh.t27 -- exact zero-order-hold constants for the linear cell, and the h/2 convergence check (Closes #10503) M1c of #10486 (weaknesses W1 and W10 of docs/analog/competitors.md). specs/analog/linear_zoh.t27 `use`s the M1 cell (linear.t27) and adds a second constant set for the same step x' = x + rnd(D x + N u): the exact zero-order hold, M = exp(A_c dt), N = A_c^-1 (M - I) B_c, computed in the tests in f64 (scaling and squaring, Pade (6,6), M1's Gauss-Jordan eliminate for the solves) and quantized to Q.24 like M1. No hardware change. Circuits: M0's RC, the series RLC of M1 (Q = 1), the same with R = 5 Ohm (Q = 20) and R = 0.5 Ohm (Q = 200, the LC tank). Every stored integer is asserted equal to round(value * 2^24); M agrees with the closed form of the series RLC, with exp of the augmented matrix and with det M = exp(tr(A_c) dt). Convergence (h, h/2, h/4 against the analytic step response): backward Euler measures order 0.985/0.993 (RC), 0.975/0.988 (RLC Q 1), 0.915/0.957 (RLC Q 20); the exact scheme is within 1e-12 relative of the analytic solution at every grid point. The Q = 20 circuit leaves the asymptotic regime over 400 steps (0.84), and a test pins that. The three errors are reported apart; the model error is stated as not applicable. Stability: ||M_q^256||_inf < 1 for the RC, the RLC and the Q = 20 circuit; the Q = 200 tank needs K = 4096 (0.83; 1.17 at K = 256), asserted both ways. LC tank after 100 periods: the exact scheme keeps the amplitude to 1.00000000000006 (f64) and 1.0003 (Q.24 constants); backward Euler keeps 1.78e-7 of it. Negative controls are tests (Pade order 1 and 2, wrong dt, negative resistance, dropped integral, BE constants offered as exact); eight mutations of the file make t27c test-report fail, and t27c seal --verify refuses a mutated copy. Co-Authored-By: Claude Sonnet 5.5 --- .trinity/seals/analog_AnalogLinearZoh.json | 33 + specs/analog/linear_zoh.t27 | 1102 ++++++++++++++++++++ 2 files changed, 1135 insertions(+) create mode 100644 .trinity/seals/analog_AnalogLinearZoh.json create mode 100644 specs/analog/linear_zoh.t27 diff --git a/.trinity/seals/analog_AnalogLinearZoh.json b/.trinity/seals/analog_AnalogLinearZoh.json new file mode 100644 index 0000000000..b329f1ff7f --- /dev/null +++ b/.trinity/seals/analog_AnalogLinearZoh.json @@ -0,0 +1,33 @@ +{ + "built_by": "t27c-bootstrap@0.5.2+1d8477c89", + "closure": [ + { + "spec": "specs/analog/linear.t27", + "spec_hash": "sha256:3ddef7be183aba583ed285334ee389b48a65b41335e98914514235db7c4e5d62" + } + ], + "config": "gen=c,rust,verilog,zig;tests=zig", + "gen_hash_c": "sha256:cfe65aeacc84e2ef7c4f413af9af6efcd2e822d94077b346f0cf98d3ac2e79fb", + "gen_hash_rust": "sha256:a8f9123547ada3c60d11b510199af61813a635716162b47e1a6881dcf4cf144b", + "gen_hash_verilog": "sha256:36bd6536290b5b4d43aba8e3486c19cbf0260d34ce12d62f12ac4ebf2012b27c", + "gen_hash_zig": "sha256:2be378a1b5e61863821678240b494a8163b054b58ff4d67af5f0cabb7ad2c5d2", + "module": "AnalogLinearZoh", + "ring": 12, + "seal_schema": 2, + "sealed_at": "2026-10-18T16:38:16Z", + "sealed_by": "t27c-bootstrap@0.5.2", + "spec_hash": "sha256:b1dde6dbbfb691939c803847318166635d0e5b709510f527200b9ad2d796f64a", + "spec_path": "specs/analog/linear_zoh.t27", + "tests": { + "failed": 0, + "failing": [], + "forced": false, + "passed": 15, + "total": 15, + "vacuous": 0 + }, + "toolchain": { + "test_runner": "sha256:1982ae948e24d00c5d221ccec12dfce1eb01204b53e37b0f857f93b16d64629d", + "zig": "0.16.0" + } +} \ No newline at end of file diff --git a/specs/analog/linear_zoh.t27 b/specs/analog/linear_zoh.t27 new file mode 100644 index 0000000000..1ffb14aca4 --- /dev/null +++ b/specs/analog/linear_zoh.t27 @@ -0,0 +1,1102 @@ +// SPDX-License-Identifier: Apache-2.0 +; specs/analog/linear_zoh.t27 -- exact zero-order-hold constants for the linear analog cell, and the h / h/2 convergence check +; M1c of the analog emulation epic (#10486, issue #10503), weaknesses W1 and W10 of docs/analog/competitors.md. Everything +; hand-written is this file. It `use`s specs/analog/linear.t27 (M1, #10487) and copies nothing from it: the step, the number +; format, the constants checks, the stability test and the f64 reference loop are M1's, called here by name. +; - +; THE PROBLEM. Backward Euler is first order and damps resonances: a high-Q LC tank decays faster on the board than in the +; circuit. "Bit-exact versus float64" cannot show that, because both sides run the same scheme. M1's constants are the +; backward Euler ones; the fabric does not care which constants it is given, because the step x' = x + rnd(D x + N u) is +; the same for any linear scheme. So the cure is a second constant set and a check that tells method error from arithmetic +; error. No hardware changes. +; - +; THE EXACT SCHEME. For x' = A_c x + B_c u with u held constant over a step (zero-order hold, which is what the board does: +; one input word per step), the solution over one step is exact: +; M = exp(A_c dt), N = A_c^-1 (M - I) B_c, x(t + dt) = M x(t) + N u. +; The tests compute them in f64: exp by scaling and squaring (halve the matrix until ||A dt||_inf <= 1/2, a Pade (6,6) +; approximant, one linear solve with M1's Gauss-Jordan eliminate, then square back), N by a second solve. They are quantized +; like M1's: D = M - I and N at 2^-24 (sign-magnitude u32 words), and every stored integer is asserted equal to +; round(value * 2^24), with no value within 1e-4 of a rounding tie. The states are M1's scaling, one volt-valued number per +; state: the capacitor voltage v_C and the inductor current times ZS = sqrt(L / C). The node between R and L has no +; capacitor: it is not a state (it is u - R i), so the exact RLC cell has two unknowns where M1's has three. +; - +; THE CIRCUITS (the element values are M1's, read from its tables). (1) the RC of M0 (R = 1 kOhm, C = 1 nF, dt = 50 ns): +; M = exp(-1/20), D = -818234 and N = 818234, where backward Euler has 798915 (1/21). (2) the series RLC (R = 100 Ohm, +; L = 10 uH, C = 1 nF, dt = 5 ns; Q = 1, w0 dt = 0.05). (3) the same with R = 5 Ohm: Q = 20, the circuit of the issue's +; scenarios. (4) the same with R = 0.5 Ohm: Q = 200, the LC tank. Q = 1 / (w0 dt) = 20 is also where backward Euler's +; numerical damping per step equals the physical one, so (3) is where W1 begins to show. +; - +; STABILITY. The same proof as M1: ||M_q^K||_inf < 1 for the quantized matrix, K = 256 (eight squarings), by M1's own +; `stable`. It holds for (1), (2) and (3). For the tank it does not at K = 256 (1.17), and the reason is physical, not a flaw: +; a Q = 200 tank forgets after 1 / (alpha dt) = 8000 steps, so 256 steps cannot show the decay. The same proof at K = 4096 +; (twelve squarings, M1's `power_norm`) gives 0.83. rho^K <= ||M^K|| holds for every K, so any K that goes below 1 proves +; rho < 1. A lossless tank (R = 0) is not certified at either K, and neither is a negative resistance. +; - +; THE CONVERGENCE CHECK (W10, cpphdl's h and h/2). The step response from rest, u = 1 V, against the analytic solution +; v(t) = 1 - e^(-a t) (cos wd t + (a / wd) sin wd t), z(t) = (w0 / wd) e^(-a t) sin wd t, a = R / 2L, w0 = 1 / sqrt(LC), +; at dt, dt/2 and dt/4, over the same time: M1's 200 steps at dt (1 us for the RLC, 10 us for the RC), 400 and 800 steps +; on the finer grids. The error is the largest |scheme - analytic| over both states and all grid points. The measured +; order is log2 of the ratio of two errors. +; RC backward Euler 9.01e-3 4.55e-3 2.29e-3 order 0.985 0.993 +; RLC Q 1 backward Euler 1.89e-2 9.62e-3 4.85e-3 order 0.975 0.988 +; RLC Q 20 backward Euler 1.66e-1 8.78e-2 4.52e-2 order 0.915 0.957 +; and the exact scheme, in every circuit and at every dt, is within 1e-12 (measured: 3e-14 or less) of the analytic solution +; relative to full scale (the issue asks for 1e-9; the Pade (6,6) of ||A dt|| = 0.1 is good to f64 round-off). Its errors at +; dt, dt/2, dt/4 are round-off, so no order is measured for it. The order is a statement about a horizon: over 400 steps +; the Q = 20 circuit measures 0.84 for the first pair, which is below 0.9, and W10's rule applies: at dt = 5 ns that circuit +; is not in the asymptotic regime for three periods, and the report says so (a test pins this). A second-order option +; would be measured against [1.9, 2.1]; none is added here. +; - +; THE THREE ERRORS, SEPARATELY (FR-001, W10). Arithmetic: the 200-step pulse run of M1 (1 V for 100 steps, then 0), the +; integer step against the f64 step of the same scheme, the largest difference in LSB of 2^-24. Discretization: the f64 +; scheme against the analytic solution, above. Model: not applicable here; the reference is the analytic solution of the +; same R, L, C equations, there is no device model, and ngspice joins at M3 (MODEL_ERROR says so). +; arithmetic, LSB RC RLC Q 1 RLC Q 20 tank +; backward Euler 2.65 8.44 (M1) +; exact 4.84 9.61 34.71 120.50 +; The exact scheme's LSB bounds are asserted from both sides: 5, 10, 35, 121. They grow with 1 / (1 - rho) of the slowest +; pole, as M1's do. The backward Euler arithmetic error of the RLC is 5.0e-7 V and its discretization error 1.89e-2 V: +; some 37000 times larger. That is W1 as a number: a bit-exact comparison with float64 could not see it. +; - +; THE TANK, AFTER 100 PERIODS (T = 2 pi / wd = 125.66 steps; 12566 steps; the amplitude is the largest |1 - v_C| over the +; last 127 steps, the ring about the final value of 1 V). Analytic 0.211153. +; exact scheme, f64 0.211153 ratio to analytic 1.00000000000006 +; exact scheme, Q.24 0.211219 ratio 1.0003 (the rounding of D moves the decay rate by 0.02 %) +; backward Euler, f64 3.75e-8 ratio 1.78e-7: it loses 99.99998 % of the amplitude +; Backward Euler's energy factor per step is 1 / (1 + R dt / L + (w0 dt)^2) against exp(-R dt / L): the closed form gives +; 1.54e-7 after 12566 steps, within 15 % of the run. At 125 steps per period the numerical damping of backward Euler is ten +; times the physical damping of the tank (alpha dt = 1.25e-4 against (w0 dt)^2 / 2 = 1.25e-3 per step). +; - +; NEGATIVE CONTROLS, as tests. A Pade (1,1) (the trapezoid rule) changes the stored integers and misses the 1e-9 grid +; test by five orders; a Pade (2,2) rounds to the SAME integers for the RLC of (2) (its error is 6e-9, below half an LSB) and +; is caught by the f64 grid test only; so the f64 grid test is the stricter of the two and both are kept. A wrong dt (2.5 ns, 10 ns), a +; negative resistance, the input term B dt without the integral, and M1's backward Euler constants offered as exact each +; change a stored integer, and the negative resistance also fails the stability test. The same mutations made to a copy of +; this file make `t27c test-report` fail (the PR records them). +; - +; WHAT A PASS SAYS. The integer constants of this file are the Q.24 roundings of the exact zero-order-hold matrices of four +; element sets; M and N agree with the closed form of the series RLC, with exp of the augmented matrix and with the +; determinant identity det M = exp(tr(A_c) dt); the unsigned sign-magnitude step of M1 leaves the DC solution (1 V, 0 A) +; exactly where it is; backward Euler at dt, dt/2, dt/4 converges at order 0.9 to 1.1 on the RC and both RLC circuits; +; the exact scheme meets the analytic solution to 1e-9 relative at every grid point; the exact RC charge at 5 us +; (16664172) is the rounding of 2^24 (1 - e^-5), where M0's backward Euler is 14537 LSB short; the quantized matrices of +; (1), (2), (3) pass the K = 256 stability test and the tank passes it at K = 4096; the fixed-point runs are within the stated +; LSB of the f64 exact scheme, and their folds are pinned so a board running RTL generated from the step can be compared. +; WHAT IT DOES NOT SAY. Nothing about a board: no bitstream, no LUT, FF, DSP or Fmax, no UART line. Nothing about a netlist: +; A_c and B_c are written here for the RC and the series RLC from their element values (the general reader is M2); the +; exact scheme has two unknowns for the RLC and is not a drop-in for M1's three-unknown tables, nor for a circuit whose A_c +; is singular. Nothing about a nonlinear device (M3), a corner, a PDK, ngspice or Symbiotic EDA's evaluator. Nothing about +; inputs that change inside a step (the hold is exact only for a piecewise-constant u), nor about another dt or circuit. +; Nothing about power. The convergence order is stated for the 1 us horizon and the three grids named above. +; ASCII only, English only. +; phi^2 + 1/phi^2 = 3 | TRINITY + +module AnalogLinearZoh; + +use analog::linear; + +pub const KIND : str = "analog-cell-constants"; +pub const ID : str = "analog/linear_zoh"; +pub const REFERENCE : str = "Stanford msdsl exact discretization of linear blocks, mirekez/cpphdl h and h/2 check, M1 specs/analog/linear.t27"; +/// The third error of FR-001, stated rather than left out. +pub const MODEL_ERROR : str = "not applicable: the reference is the analytic solution of the same equations; ngspice joins at M3"; + +; --- the method --------------------------------------------------------------------------------- +/// The order p of the Pade (p, p) approximant. +pub const PADE_ORDER : u32 = 6; +/// 2 pi. +pub const TWO_PI : f64 = 6.283185307179586; +/// The tank rings for this many periods. +pub const PERIODS : u32 = 100; +/// 2^12 steps for the tank's stability proof; M1's K is 2^8. +pub const TANK_SQUARINGS : u32 = 12; +/// The exact scheme meets the analytic solution to this, relative to full scale (the issue's 1e-9). +pub const GRID_TOL : f64 = 0.000000001; + +; --- the circuits: element values are M1's; only the resistance of (3) and (4) is new --------------------- +/// (3) R = 5 Ohm: Q = ZS / R = 20. +pub const R_Q20 : f64 = 5.0; +/// (4) R = 0.5 Ohm: Q = 200, the LC tank. +pub const R_TANK : f64 = 0.5; + +; --- the stored constants: D = M - I and N at 2^-24, sign-magnitude words (as M1) ------------------------ +; (1) the RC of M0, exact. D = -(1 - e^-0.05), N = +(1 - e^-0.05). +pub const ZRC_N : usize = 1; +pub const ZRC_D : [1]u32 = [NEG | 818234]; +/// (`** 1`: a lone literal in a one-element module-level array lowers to Zig's array TYPE, as in M1.) +pub const ZRC_NV : [1]u32 = [818234] ** 1; +pub const ZRC_EXP_S : u64 = 0x22D36DC066F4E558; +pub const ZRC_MAX_ERR_LSB : u32 = 5; + +; (2) the series RLC of M1, exact; states [v_C, i_L * ZS]. +pub const ZRLC_N : usize = 2; +pub const ZRLC_D : [4]u32 = [NEG | 20622, 817894, NEG | 817894, NEG | 838516]; +pub const ZRLC_NV : [2]u32 = [20622, 817894]; +pub const ZRLC_EXP_S : u64 = 0x69134C22123B5111; +pub const ZRLC_MAX_ERR_LSB : u32 = 10; + +; (3) the same with R = 5 Ohm (Q = 20). +pub const ZQ20_N : usize = 2; +pub const ZQ20_D : [4]u32 = [NEG | 20950, 837464, NEG | 837464, NEG | 62823]; +pub const ZQ20_NV : [2]u32 = [20950, 837464]; +pub const ZQ20_EXP_S : u64 = 0x09F569FF9BCFAD1A; +pub const ZQ20_MAX_ERR_LSB : u32 = 35; + +; (4) the same with R = 0.5 Ohm (Q = 200, the tank). +pub const ZTK_N : usize = 2; +pub const ZTK_D : [4]u32 = [NEG | 20965, 838407, NEG | 838407, NEG | 25157]; +pub const ZTK_NV : [2]u32 = [20965, 838407]; +pub const ZTK_EXP_S : u64 = 0xD69BF5F73A7CBE21; +pub const ZTK_MAX_ERR_LSB : u32 = 121; + +; --- exp(A): scaling and squaring with a Pade (p, p) approximant ----------------------------------------- +; Every fn below is f64 and is reached only by the tests, like M1's reference side. + +/// |a - b| <= tol. +pub fn near(a: f64, b: f64, tol: f64) -> bool { + return fabs64(a - b) <= tol; +} + +/// c_k of the Pade (p, p) approximant of exp: c_0 = 1, c_k = c_(k-1) (p - k + 1) / (k (2p - k + 1)), for k <= p. +pub fn pade_coef(p: u32, k: u32) -> f64 { + var c : f64 = 1.0; + var i : u32 = 1; + while (i <= k) { + c = c * ((p - i + 1) as f64) / ((i as f64) * ((2 * p - i + 1) as f64)); + i = i + 1; + } + return c; +} + +/// exp(a) for an n x n matrix (n <= 5, row-major): s halvings until ||a / 2^s||_inf <= 1/2, X = D^-1 N with +/// N = sum c_k Y^k and D = sum (-1)^k c_k Y^k for Y = a / 2^s, k = 0 .. order, then s squarings. +pub fn pade_expm(n: usize, a: []const f64, order: u32) -> [25]f64 { + var s : u32 = 0; + var nrm : f64 = norm_inf(n, a); + while (nrm > 0.5) { + nrm = nrm / 2.0; + s = s + 1; + } + var scale : f64 = 1.0; + var q : u32 = 0; + while (q < s) { + scale = scale / 2.0; + q = q + 1; + } + var y : [25]f64 = [0.0] ** 25; + var pw : [25]f64 = [0.0] ** 25; + var nx : [25]f64 = [0.0] ** 25; + var num : [25]f64 = [0.0] ** 25; + var den : [25]f64 = [0.0] ** 25; + var i : usize = 0; + var j : usize = 0; + while (i < n * n) { + y[i] = a[i] * scale; + i = i + 1; + } + i = 0; + while (i < n) { + pw[i * n + i] = 1.0; + num[i * n + i] = 1.0; + den[i * n + i] = 1.0; + i = i + 1; + } + var sg : f64 = 1.0; + var k : u32 = 1; + while (k <= order) { + const c : f64 = pade_coef(order, k); + sg = 0.0 - sg; + matmul(n, &pw, &y, &nx); + i = 0; + while (i < n * n) { + pw[i] = nx[i]; + num[i] = num[i] + c * pw[i]; + den[i] = den[i] + sg * c * pw[i]; + i = i + 1; + } + k = k + 1; + } + var w : [55]f64 = [0.0] ** 55; + i = 0; + while (i < n) { + j = 0; + while (j < n) { + w[i * WIDE + j] = den[i * n + j]; + w[i * WIDE + n + j] = num[i * n + j]; + j = j + 1; + } + i = i + 1; + } + const sol : [30]f64 = mn_of(eliminate(w, n), n); + var e : [25]f64 = [0.0] ** 25; + i = 0; + while (i < n * n) { + e[i] = sol[i]; + i = i + 1; + } + q = 0; + while (q < s) { + matmul(n, &e, &e, &nx); + i = 0; + while (i < n * n) { + e[i] = nx[i]; + i = i + 1; + } + q = q + 1; + } + return e; +} + +/// The zero-order-hold pair of x' = ac x + bc u over dt, dense as M1's derive returns it: M row-major in [0, n*n), N after: +/// M = exp(ac dt), N = ac^-1 (M - I) bc. +pub fn zoh(n: usize, ac: []const f64, bc: []const f64, dt: f64, order: u32) -> [30]f64 { + var adt : [25]f64 = [0.0] ** 25; + var i : usize = 0; + var j : usize = 0; + while (i < n * n) { + adt[i] = ac[i] * dt; + i = i + 1; + } + const m : [25]f64 = pade_expm(n, &adt, order); + var w : [55]f64 = [0.0] ** 55; + i = 0; + while (i < n) { + var rhs : f64 = 0.0; + j = 0; + while (j < n) { + var mij : f64 = m[i * n + j]; + if (i == j) { mij = mij - 1.0; } + w[i * WIDE + j] = ac[i * n + j]; + rhs = rhs + mij * bc[j]; + j = j + 1; + } + w[i * WIDE + 2 * n] = rhs; + i = i + 1; + } + const sol : [30]f64 = mn_of(eliminate(w, n), n); + var o : [30]f64 = [0.0] ** 30; + i = 0; + while (i < n * n) { + o[i] = m[i]; + i = i + 1; + } + i = 0; + while (i < n) { + o[n * n + i] = sol[n * n + i]; + i = i + 1; + } + return o; +} + +; --- A_c and B_c of the circuits, from M1's element values ------------------------------------------------ + +/// The series RLC of resistance r in the states [v_C, z = i_L * ZS]: dv/dt = z / (ZS C), dz/dt = (ZS / L) (u - v) - (r / L) z. +pub fn rlc_ac(r: f64) -> [4]f64 { + var a : [4]f64 = [0.0] ** 4; + a[1] = 1.0 / (RLC_ZS * RLC_VAL[3]); + a[2] = 0.0 - RLC_ZS / RLC_VAL[2]; + a[3] = 0.0 - r / RLC_VAL[2]; + return a; +} + +/// B_c of the series RLC: the source enters the current's row only. +pub fn rlc_bc() -> [2]f64 { + var b : [2]f64 = [0.0] ** 2; + b[1] = RLC_ZS / RLC_VAL[2]; + return b; +} + +/// The exact pair of the series RLC with resistance r. +pub fn rlc_zoh(r: f64, dt: f64, order: u32) -> [30]f64 { + const a : [4]f64 = rlc_ac(r); + const b : [2]f64 = rlc_bc(); + return zoh(2, &a, &b, dt, order); +} + +/// The exact pair of M0's RC: dv/dt = (u - v) / RC. +pub fn rc_zoh(dt: f64, order: u32) -> [30]f64 { + var a : [1]f64 = [0.0] ** 1; + var b : [1]f64 = [0.0] ** 1; + a[0] = 0.0 - 1.0 / (RC_VAL[1] * RC_VAL[2]); + b[0] = 1.0 / (RC_VAL[1] * RC_VAL[2]); + return zoh(1, &a, &b, dt, order); +} + +/// M1's backward Euler pair for the series RLC with resistance r at step dt: three unknowns [node, v_C, z]. +pub fn be_rlc(r: f64, dt: f64) -> [30]f64 { + var v : [4]f64 = RLC_VAL; + v[1] = r; + return derive(&RLC_KIND, &RLC_A, &RLC_B, &v, RLC_NE, dt, RLC_ZS); +} + +/// M1's backward Euler pair for the RC at step dt. +pub fn be_rc(dt: f64) -> [30]f64 { + return derive(&RC_KIND, &RC_A, &RC_B, &RC_VAL, RC_NE, dt, RC_ZS); +} + +; --- the analytic step response from rest --------------------------------------------------------------- + +/// w0 = 1 / sqrt(L C). +pub fn omega0() -> f64 { + return 1.0 / @sqrt(RLC_VAL[2] * RLC_VAL[3]); +} + +/// The decay rate R / 2L. +pub fn alpha_of(r: f64) -> f64 { + return r / (2.0 * RLC_VAL[2]); +} + +/// The damped frequency sqrt(w0^2 - alpha^2); underdamped circuits only. +pub fn omega_d(r: f64) -> f64 { + const w : f64 = omega0(); + const al : f64 = alpha_of(r); + return @sqrt(w * w - al * al); +} + +/// v_C(t) for u = 1 V from rest. +pub fn an_v(r: f64, t: f64) -> f64 { + const al : f64 = alpha_of(r); + const wd : f64 = omega_d(r); + return 1.0 - @exp(0.0 - al * t) * (@cos(wd * t) + al / wd * @sin(wd * t)); +} + +/// z(t) = i_L(t) * ZS for u = 1 V from rest. +pub fn an_z(r: f64, t: f64) -> f64 { + return omega0() / omega_d(r) * @exp(0.0 - alpha_of(r) * t) * @sin(omega_d(r) * t); +} + +/// v(t) of M0's RC for u = 1 V from rest. +pub fn an_rc(t: f64) -> f64 { + return 1.0 - @exp(0.0 - t / (RC_VAL[1] * RC_VAL[2])); +} + +; --- the error of a scheme against the analytic solution ----------------------------------------------- + +pub struct GridErr { err: f64, scale: f64 } + +/// k steps of x' = M x + N u from rest with u = 1 (the dense pair mn of n states; the capacitor voltage is state iv and +/// z is state iz) at step dt of the series RLC with resistance r. err is the largest |x - analytic| over both states and +/// all grid points t = dt .. k dt; scale is the largest analytic magnitude, so err / scale is relative to full scale. +pub fn grid_err(n: usize, mn: []const f64, iv: usize, iz: usize, r: f64, dt: f64, k: u32) -> GridErr { + var x : [5]f64 = [0.0] ** 5; + var y : [5]f64 = [0.0] ** 5; + var err : f64 = 0.0; + var scale : f64 = 0.0; + var t : u32 = 0; + var i : usize = 0; + var j : usize = 0; + while (t < k) { + i = 0; + while (i < n) { + var acc : f64 = mn[n * n + i]; + j = 0; + while (j < n) { + acc = acc + mn[i * n + j] * x[j]; + j = j + 1; + } + y[i] = acc; + i = i + 1; + } + i = 0; + while (i < n) { + x[i] = y[i]; + i = i + 1; + } + const tt : f64 = ((t + 1) as f64) * dt; + const v : f64 = an_v(r, tt); + const z : f64 = an_z(r, tt); + err = fmax64(err, fmax64(fabs64(x[iv] - v), fabs64(x[iz] - z))); + scale = fmax64(scale, fmax64(fabs64(v), fabs64(z))); + t = t + 1; + } + return GridErr { err: err, scale: scale }; +} + +/// The same for the RC: one state, the analytic 1 - e^(-t / RC). +pub fn grid_err_rc(mn: []const f64, dt: f64, k: u32) -> GridErr { + var x : f64 = 0.0; + var err : f64 = 0.0; + var scale : f64 = 0.0; + var t : u32 = 0; + while (t < k) { + x = mn[0] * x + mn[1]; + const v : f64 = an_rc(((t + 1) as f64) * dt); + err = fmax64(err, fabs64(x - v)); + scale = fmax64(scale, fabs64(v)); + t = t + 1; + } + return GridErr { err: err, scale: scale }; +} + +/// The larger of two numbers. +pub fn fmax64(a: f64, b: f64) -> f64 { + if (a > b) { return a; } + return b; +} + +/// The measured convergence order: log2 of the ratio of the coarse error to the fine error (the step was halved). +pub fn order_of(coarse: f64, fine: f64) -> f64 { + return @log2(coarse / fine); +} + +/// Backward Euler on the series RLC (resistance r): the error at step dt / 2^j over steps * 2^j steps, the same horizon. +pub fn rlc_be_err(r: f64, j: u32, steps: u32) -> GridErr { + const m : u32 = 1 << j; + const dt : f64 = RLC_DT / (m as f64); + const mn : [30]f64 = be_rlc(r, dt); + return grid_err(3, &mn, 1, 2, r, dt, steps * m); +} + +/// The exact scheme on the series RLC, same grid. +pub fn rlc_ex_err(r: f64, j: u32, steps: u32) -> GridErr { + const m : u32 = 1 << j; + const dt : f64 = RLC_DT / (m as f64); + const mn : [30]f64 = rlc_zoh(r, dt, PADE_ORDER); + return grid_err(2, &mn, 0, 1, r, dt, steps * m); +} + +/// Backward Euler on the RC, step RC_DT / 2^j. +pub fn rc_be_err(j: u32, steps: u32) -> GridErr { + const m : u32 = 1 << j; + const dt : f64 = RC_DT / (m as f64); + const mn : [30]f64 = be_rc(dt); + return grid_err_rc(&mn, dt, steps * m); +} + +/// The exact scheme on the RC, same grid. +pub fn rc_ex_err(j: u32, steps: u32) -> GridErr { + const m : u32 = 1 << j; + const dt : f64 = RC_DT / (m as f64); + const mn : [30]f64 = rc_zoh(dt, PADE_ORDER); + return grid_err_rc(&mn, dt, steps * m); +} + +; --- the ring: amplitude after 100 periods ------------------------------------------------------------- + +pub struct Ring { scheme: f64, analytic: f64 } + +/// Steps per damped period of the series RLC of resistance r at step dt. +pub fn period_steps(r: f64, dt: f64) -> f64 { + return TWO_PI / omega_d(r) / dt; +} + +/// The step count of PERIODS periods, rounded to the nearest step. +pub fn ring_k(r: f64, dt: f64) -> u32 { + return ((PERIODS as f64) * period_steps(r, dt) + 0.5) as u32; +} + +/// The window that closes the run: one period and one more step. +pub fn ring_w(r: f64, dt: f64) -> u32 { + return (period_steps(r, dt) + 1.5) as u32; +} + +/// The f64 scheme (dense pair mn of n states, capacitor voltage in state iv) on the step input for PERIODS periods: +/// the largest |1 - v_C| in the closing window, beside the analytic one over the same steps. +pub fn ring_f64(n: usize, mn: []const f64, iv: usize, r: f64, dt: f64) -> Ring { + const k : u32 = ring_k(r, dt); + const w : u32 = ring_w(r, dt); + var x : [5]f64 = [0.0] ** 5; + var y : [5]f64 = [0.0] ** 5; + var ps : f64 = 0.0; + var pa : f64 = 0.0; + var t : u32 = 0; + var i : usize = 0; + var j : usize = 0; + while (t < k) { + i = 0; + while (i < n) { + var acc : f64 = mn[n * n + i]; + j = 0; + while (j < n) { + acc = acc + mn[i * n + j] * x[j]; + j = j + 1; + } + y[i] = acc; + i = i + 1; + } + i = 0; + while (i < n) { + x[i] = y[i]; + i = i + 1; + } + if (t + w >= k) { + ps = fmax64(ps, fabs64(1.0 - x[iv])); + pa = fmax64(pa, fabs64(1.0 - an_v(r, ((t + 1) as f64) * dt))); + } + t = t + 1; + } + return Ring { scheme: ps, analytic: pa }; +} + +/// The same with the integer step of M1 (sign-magnitude words, d and nv the stored constants), u = 1 V throughout. +pub fn ring_fixed(n: usize, d: []const u32, nv: []const u32, iv: usize, r: f64, dt: f64) -> Ring { + const k : u32 = ring_k(r, dt); + const w : u32 = ring_w(r, dt); + var x : [5]u32 = [0] ** 5; + var y : [5]u32 = [0] ** 5; + var ps : f64 = 0.0; + var pa : f64 = 0.0; + var t : u32 = 0; + var i : usize = 0; + while (t < k) { + step(n, d, nv, &x, ONE, &y); + i = 0; + while (i < n) { + x[i] = y[i]; + i = i + 1; + } + if (t + w >= k) { + ps = fmax64(ps, fabs64(1.0 - sm_f64(x[iv]))); + pa = fmax64(pa, fabs64(1.0 - an_v(r, ((t + 1) as f64) * dt))); + } + t = t + 1; + } + return Ring { scheme: ps, analytic: pa }; +} + +; --- the tests ---------------------------------------------------------------------------------- + +test the_pade_exponential_matches_matrices_whose_exponential_is_known { + // the zero matrix + var z : [25]f64 = [0.0] ** 25; + const ez : [25]f64 = pade_expm(2, &z, PADE_ORDER); + assert(near(ez[0], 1.0, 0.000000000000001)); + assert(near(ez[1], 0.0, 0.000000000000001)); + assert(near(ez[2], 0.0, 0.000000000000001)); + assert(near(ez[3], 1.0, 0.000000000000001)); + // a nilpotent matrix: exp([[0, 1], [0, 0]]) = [[1, 1], [0, 1]] + var nl : [25]f64 = [0.0] ** 25; + nl[1] = 1.0; + const en : [25]f64 = pade_expm(2, &nl, PADE_ORDER); + assert(near(en[0], 1.0, 0.00000000000001)); + assert(near(en[1], 1.0, 0.00000000000001)); + assert(near(en[2], 0.0, 0.00000000000001)); + assert(near(en[3], 1.0, 0.00000000000001)); + // a diagonal one that needs four halvings: exp(diag(-5, 1)) + var dg : [25]f64 = [0.0] ** 25; + dg[0] = 0.0 - 5.0; + dg[3] = 1.0; + const ed : [25]f64 = pade_expm(2, &dg, PADE_ORDER); + assert(near(ed[0], 0.006737946999085467, 0.0000000000000001)); + assert(near(ed[3], 2.718281828459045, 0.0000000000001)); + assert(near(ed[1], 0.0, 0.00000000000001)); + assert(near(ed[2], 0.0, 0.00000000000001)); + // a rotation by 3 rad (norm 3, two halvings): exp([[0, 3], [-3, 0]]) = [[cos 3, sin 3], [-sin 3, cos 3]] + var rt : [25]f64 = [0.0] ** 25; + rt[1] = 3.0; + rt[2] = 0.0 - 3.0; + const er : [25]f64 = pade_expm(2, &rt, PADE_ORDER); + assert(near(er[0], @cos(3.0), 0.0000000000001)); + assert(near(er[1], @sin(3.0), 0.0000000000001)); + assert(near(er[2], 0.0 - @sin(3.0), 0.0000000000001)); + assert(near(er[3], @cos(3.0), 0.0000000000001)); + // exp(A) exp(-A) = I for a full matrix + var fa : [25]f64 = [0.0] ** 25; + var fb : [25]f64 = [0.0] ** 25; + fa[0] = 0.3; + fa[1] = 0.0 - 1.2; + fa[2] = 0.7; + fa[3] = 0.1; + fb[0] = 0.0 - 0.3; + fb[1] = 1.2; + fb[2] = 0.0 - 0.7; + fb[3] = 0.0 - 0.1; + const ea : [25]f64 = pade_expm(2, &fa, PADE_ORDER); + const eb : [25]f64 = pade_expm(2, &fb, PADE_ORDER); + var pr : [25]f64 = [0.0] ** 25; + matmul(2, &ea, &eb, &pr); + assert(near(pr[0], 1.0, 0.00000000000001)); + assert(near(pr[1], 0.0, 0.00000000000001)); + assert(near(pr[2], 0.0, 0.00000000000001)); + assert(near(pr[3], 1.0, 0.00000000000001)); + // a three by three, to exercise the stride: exp of a diagonal block and a scalar + var t3 : [25]f64 = [0.0] ** 25; + t3[0] = 0.0 - 1.0; + t3[4] = 0.0 - 2.0; + t3[8] = 0.5; + const e3 : [25]f64 = pade_expm(3, &t3, PADE_ORDER); + assert(near(e3[0], 0.36787944117144233, 0.00000000000001)); + assert(near(e3[4], 0.1353352832366127, 0.00000000000001)); + assert(near(e3[8], 1.6487212707001282, 0.00000000000001)); + assert(near(e3[1], 0.0, 0.00000000000001)); + // the coefficients: Pade (2, 2) has c = 1, 1/2, 1/12 + assert(near(pade_coef(2, 0), 1.0, 0.0000000000000001)); + assert(near(pade_coef(2, 1), 0.5, 0.0000000000000001)); + assert(near(pade_coef(2, 2), 1.0 / 12.0, 0.0000000000000001)); + // a lower order is visibly wrong where the order matters: Pade (1, 1) on diag(-5, 1) is off by percents + const e1 : [25]f64 = pade_expm(2, &dg, 1); + assert(fabs64(e1[0] - 0.006737946999085467) > 0.0001); + assert(fabs64(e1[3] - 2.718281828459045) > 0.0005); +} + +test the_rc_exact_constants_are_the_rounding_of_the_derived_ones { + const mn : [30]f64 = rc_zoh(RC_DT, PADE_ORDER); + assert_eq(mismatches(ZRC_N, &mn, &ZRC_D, &ZRC_NV), 0); + assert(ties_clear(ZRC_N, &mn)); + // the answer is e^(-dt / RC) = e^(-1/20), and N = 1 - M + assert(near(mn[0], @exp(0.0 - 0.05), 0.0000000000000003)); + assert(near(mn[1], 1.0 - @exp(0.0 - 0.05), 0.0000000000000003)); + assert_eq(ZRC_D[0], NEG | 818234); + assert_eq(ZRC_NV[0], 818234); + // M0's backward Euler has 1/21 = 798915 where the exact scheme has 1 - e^-0.05 = 818234 (0.0488 against 0.0476) + assert(RC_NV[0] != ZRC_NV[0]); + assert_eq(ZRC_NV[0] - RC_NV[0], 19319); + // the offered pairs are not each other's: backward Euler constants fail the exact check and the exact fail the BE one + const be : [30]f64 = be_rc(RC_DT); + assert(mismatches(ZRC_N, &be, &ZRC_D, &ZRC_NV) > 0); + assert(mismatches(RC_N, &mn, &RC_D, &RC_NV) > 0); + assert_eq(mismatches(RC_N, &be, &RC_D, &RC_NV), 0); +} + +test the_rlc_exact_constants_are_the_rounding_of_the_derived_ones { + const mn : [30]f64 = rlc_zoh(RLC_VAL[1], RLC_DT, PADE_ORDER); + assert_eq(mismatches(ZRLC_N, &mn, &ZRLC_D, &ZRLC_NV), 0); + assert(ties_clear(ZRLC_N, &mn)); + // the closed form of exp(A_c dt) for the series RLC, from nothing but R, L, C and dt + const al : f64 = alpha_of(RLC_VAL[1]); + const wd : f64 = omega_d(RLC_VAL[1]); + const w0 : f64 = omega0(); + const ee : f64 = @exp(0.0 - al * RLC_DT); + const cd : f64 = @cos(wd * RLC_DT); + const sd : f64 = @sin(wd * RLC_DT); + assert(near(mn[0], ee * (cd + al / wd * sd), 0.00000000000001)); + assert(near(mn[1], ee * w0 / wd * sd, 0.00000000000001)); + assert(near(mn[2], 0.0 - ee * w0 / wd * sd, 0.00000000000001)); + assert(near(mn[3], ee * (cd - al / wd * sd), 0.00000000000001)); + // N follows from the steady state (1 V, 0 A): N = (I - M) (1, 0) + assert(near(mn[4], 1.0 - mn[0], 0.00000000000001)); + assert(near(mn[5], 0.0 - mn[2], 0.00000000000001)); + // Liouville: det exp(A dt) = exp(tr(A) dt) = exp(-R dt / L), the identity that backward Euler misses (M1: 0.9501) + const det : f64 = mn[0] * mn[3] - mn[1] * mn[2]; + assert(near(det, @exp(0.0 - RLC_VAL[1] * RLC_DT / RLC_VAL[2]), 0.00000000000001)); + assert(near(det, 0.951229424500714, 0.000000000000001)); + // N from the matrix exponential of the augmented matrix [[A dt, B dt], [0, 0]], whose last column is N + const a : [4]f64 = rlc_ac(RLC_VAL[1]); + const b : [2]f64 = rlc_bc(); + var g : [25]f64 = [0.0] ** 25; + g[0] = a[0] * RLC_DT; + g[1] = a[1] * RLC_DT; + g[2] = b[0] * RLC_DT; + g[3] = a[2] * RLC_DT; + g[4] = a[3] * RLC_DT; + g[5] = b[1] * RLC_DT; + const eg : [25]f64 = pade_expm(3, &g, PADE_ORDER); + assert(near(eg[2], mn[4], 0.00000000000001)); + assert(near(eg[5], mn[5], 0.00000000000001)); + assert(near(eg[8], 1.0, 0.00000000000001)); + assert(near(eg[0], mn[0], 0.00000000000001)); + assert(near(eg[4], mn[3], 0.00000000000001)); + // ZS = sqrt(L / C) keeps every stored magnitude below 1 (2^24), as in M1 + var i : usize = 0; + var top : u32 = 0; + while (i < 4) { + if ((ZRLC_D[i] & MAG) > top) { top = ZRLC_D[i] & MAG; } + i = i + 1; + } + assert(top < ONE); + assert_eq(top, 838516); +} + +test the_q20_and_tank_exact_constants_are_the_rounding_of_the_derived_ones { + const m20 : [30]f64 = rlc_zoh(R_Q20, RLC_DT, PADE_ORDER); + assert_eq(mismatches(ZQ20_N, &m20, &ZQ20_D, &ZQ20_NV), 0); + assert(ties_clear(ZQ20_N, &m20)); + const mt : [30]f64 = rlc_zoh(R_TANK, RLC_DT, PADE_ORDER); + assert_eq(mismatches(ZTK_N, &mt, &ZTK_D, &ZTK_NV), 0); + assert(ties_clear(ZTK_N, &mt)); + // Q = ZS / R + assert(near(RLC_ZS / R_Q20, 20.0, 0.0000000000001)); + assert(near(RLC_ZS / R_TANK, 200.0, 0.0000000000001)); + // Q = 1 / (w0 dt) = 20 is where backward Euler's numerical damping equals the physical one + assert(near(1.0 / (omega0() * RLC_DT), 20.0, 0.0000000000001)); + // the same determinant identity at both resistances + assert(near(m20[0] * m20[3] - m20[1] * m20[2], @exp(0.0 - R_Q20 * RLC_DT / RLC_VAL[2]), 0.00000000000001)); + assert(near(mt[0] * mt[3] - mt[1] * mt[2], @exp(0.0 - R_TANK * RLC_DT / RLC_VAL[2]), 0.00000000000001)); + // less loss leaves D_zz smaller: the damping term moves from 838516 to 62823 to 25157 + assert((ZRLC_D[3] & MAG) > (ZQ20_D[3] & MAG)); + assert((ZQ20_D[3] & MAG) > (ZTK_D[3] & MAG)); +} + +test the_stored_integers_leave_the_dc_solution_where_it_is { + // u = 1 V holds x = (1 V, 0 A) fixed: D x + N = 0 exactly, so the DC gain is exactly 1 and no sum rounds away from it + var o : [5]u32 = [0] ** 5; + const sc : [1]u32 = [ONE]; + step(1, &ZRC_D, &ZRC_NV, &sc, ONE, &o); + assert_eq(o[0], ONE); + const s2 : [2]u32 = [ONE, 0]; + step(2, &ZRLC_D, &ZRLC_NV, &s2, ONE, &o); + assert_eq(o[0], ONE); + assert_eq(o[1], 0); + step(2, &ZQ20_D, &ZQ20_NV, &s2, ONE, &o); + assert_eq(o[0], ONE); + assert_eq(o[1], 0); + step(2, &ZTK_D, &ZTK_NV, &s2, ONE, &o); + assert_eq(o[0], ONE); + assert_eq(o[1], 0); + // and with u = 0 the state decays toward 0 and a state at rest stays at rest + const zz : [2]u32 = [0, 0]; + step(2, &ZRLC_D, &ZRLC_NV, &zz, 0, &o); + assert_eq(o[0], 0); + assert_eq(o[1], 0); + // the first step from rest is N itself (D x = 0 there) + step(2, &ZRLC_D, &ZRLC_NV, &zz, ONE, &o); + assert_eq(o[0], 20622); + assert_eq(o[1], 817894); +} + +test backward_euler_converges_at_order_one_on_the_rc_and_both_rlc_circuits { + // the RC of M0 + const c0 : f64 = rc_be_err(0, STEPS).err; + const c1 : f64 = rc_be_err(1, STEPS).err; + const c2 : f64 = rc_be_err(2, STEPS).err; + assert(order_of(c0, c1) > 0.9); + assert(order_of(c0, c1) < 1.1); + assert(order_of(c1, c2) > 0.9); + assert(order_of(c1, c2) < 1.1); + assert(c0 > 0.0090); + assert(c0 < 0.0091); + assert(order_of(c0, c1) > 0.98); + assert(order_of(c0, c1) < 0.99); + // the series RLC of M1, Q = 1 + const a0 : f64 = rlc_be_err(RLC_VAL[1], 0, STEPS).err; + const a1 : f64 = rlc_be_err(RLC_VAL[1], 1, STEPS).err; + const a2 : f64 = rlc_be_err(RLC_VAL[1], 2, STEPS).err; + assert(order_of(a0, a1) > 0.9); + assert(order_of(a0, a1) < 1.1); + assert(order_of(a1, a2) > 0.9); + assert(order_of(a1, a2) < 1.1); + assert(a0 > 0.0189); + assert(a0 < 0.0190); + assert(order_of(a0, a1) > 0.97); + assert(order_of(a0, a1) < 0.98); + // the same with R = 5 Ohm, Q = 20 + const b0 : f64 = rlc_be_err(R_Q20, 0, STEPS).err; + const b1 : f64 = rlc_be_err(R_Q20, 1, STEPS).err; + const b2 : f64 = rlc_be_err(R_Q20, 2, STEPS).err; + assert(order_of(b0, b1) > 0.9); + assert(order_of(b0, b1) < 1.1); + assert(order_of(b1, b2) > 0.9); + assert(order_of(b1, b2) < 1.1); + assert(b0 > 0.165); + assert(b0 < 0.166); + assert(order_of(b0, b1) > 0.91); + assert(order_of(b0, b1) < 0.92); + // W10's form of the same fact: error(h/2) / error(h) in [0.4, 0.6] + assert(c1 / c0 > 0.4); + assert(c1 / c0 < 0.6); + assert(a1 / a0 > 0.4); + assert(a1 / a0 < 0.6); + assert(b1 / b0 > 0.4); + assert(b1 / b0 < 0.6); + assert(b2 / b1 > 0.4); + assert(b2 / b1 < 0.6); + // the finer the grid the smaller the error, strictly + assert(c0 > c1); + assert(c1 > c2); + assert(a0 > a1); + assert(a1 > a2); + assert(b0 > b1); + assert(b1 > b2); +} + +test the_q20_circuit_is_not_in_the_asymptotic_regime_over_longer_horizons { + // over 400 steps (2 us, three periods) the first pair measures 0.84: below 0.9, and W10 says to mark that circuit + const e0 : f64 = rlc_be_err(R_Q20, 0, 2 * STEPS).err; + const e1 : f64 = rlc_be_err(R_Q20, 1, 2 * STEPS).err; + const e2 : f64 = rlc_be_err(R_Q20, 2, 2 * STEPS).err; + assert(order_of(e0, e1) < 0.9); + assert(order_of(e0, e1) > 0.8); + // it recovers on the finer pair, as the asymptotic regime is approached from below + assert(order_of(e1, e2) > order_of(e0, e1)); + assert(order_of(e1, e2) < 1.0); + // while the Q = 1 circuit has the same order over this horizon as over the first + const f0 : f64 = rlc_be_err(RLC_VAL[1], 0, 2 * STEPS).err; + const f1 : f64 = rlc_be_err(RLC_VAL[1], 1, 2 * STEPS).err; + assert(order_of(f0, f1) > 0.9); + assert(order_of(f0, f1) < 1.1); +} + +test the_exact_scheme_meets_the_analytic_solution_at_every_grid_point { + var j : u32 = 0; + while (j < 3) { + const e1 : GridErr = rlc_ex_err(RLC_VAL[1], j, STEPS); + assert(e1.err / e1.scale < GRID_TOL); + assert(e1.err / e1.scale < 0.000000000001); + assert(e1.scale > 1.1); + const e2 : GridErr = rlc_ex_err(R_Q20, j, STEPS); + assert(e2.err / e2.scale < GRID_TOL); + assert(e2.err / e2.scale < 0.000000000001); + assert(e2.scale > 1.8); + const e3 : GridErr = rlc_ex_err(R_TANK, j, STEPS); + assert(e3.err / e3.scale < GRID_TOL); + assert(e3.err / e3.scale < 0.000000000001); + assert(e3.scale > 1.9); + const e4 : GridErr = rc_ex_err(j, STEPS); + assert(e4.err / e4.scale < GRID_TOL); + assert(e4.err / e4.scale < 0.000000000001); + assert(e4.scale > 0.99); + j = j + 1; + } + // while backward Euler at the same dt is some ten orders away from it, in the same measure + const b : GridErr = rlc_be_err(RLC_VAL[1], 0, STEPS); + assert(b.err / b.scale > 0.0155); + assert(b.err / b.scale < 0.0170); + const r : GridErr = rc_be_err(0, STEPS); + assert(r.err / r.scale > 0.0090); + assert(r.err / r.scale < 0.0091); +} + +test the_three_errors_are_separate_numbers_and_the_arithmetic_one_cannot_see_the_method { + // backward Euler of the RLC: arithmetic error from M1's run, discretization error from the analytic solution + const mb : [30]f64 = be_rlc(RLC_VAL[1], RLC_DT); + assert_eq(mismatches(RLC_N, &mb, &RLC_D, &RLC_NV), 0); + var xb : [5]u32 = [0] ** 5; + const rb : RunOut = run(RLC_N, &RLC_D, &RLC_NV, &mb, STEPS, &xb); + assert(rb.worst < 9.0); + assert(rb.worst > 8.0); + const db : GridErr = rlc_be_err(RLC_VAL[1], 0, STEPS); + // the arithmetic error in volts is ~5e-7, the discretization error ~2e-2 + const arith_v : f64 = rb.worst / SCALE; + assert(arith_v < 0.0000006); + assert(db.err > 0.0189); + assert(db.err > 30000.0 * arith_v); + // the exact scheme of the same circuit: its own arithmetic error and a discretization error at round-off + const me : [30]f64 = rlc_zoh(RLC_VAL[1], RLC_DT, PADE_ORDER); + var xe : [5]u32 = [0] ** 5; + const re : RunOut = run(ZRLC_N, &ZRLC_D, &ZRLC_NV, &me, STEPS, &xe); + const de : GridErr = rlc_ex_err(RLC_VAL[1], 0, STEPS); + assert(re.worst < (ZRLC_MAX_ERR_LSB as f64)); + assert(re.worst > (ZRLC_MAX_ERR_LSB as f64) - 1.0); + assert(de.err < 0.00000000001); + // so for the exact scheme the arithmetic error is the whole error, and for backward Euler it is not + assert(de.err < re.worst / SCALE / 1000.0); + // the model error has no number here, and says why + assert(MODEL_ERROR != ""); +} + +test the_exact_rc_charges_to_the_analytic_value_where_backward_euler_is_14537_lsb_short { + const mn : [30]f64 = rc_zoh(RC_DT, PADE_ORDER); + var x : [5]u32 = [0] ** 5; + const r : RunOut = run(ZRC_N, &ZRC_D, &ZRC_NV, &mn, CHARGE_STEPS, &x); + // 1 - e^-5 at 5 us, rounded to 2^-24 + assert_eq(x[0], q_sm(1.0 - @exp(0.0 - 5.0))); + assert_eq(x[0], 16664172); + // M0's backward Euler: 16649635, i.e. cpphdl's 0.992396 against 0.993262 + assert_eq(x[0] - 16649635, 14537); + const vv : f64 = sm_f64(x[0]); + assert(near(vv, 0.9932620530009145, 0.0000001)); + assert(r.worst < (ZRC_MAX_ERR_LSB as f64)); + assert(r.peak < HEADROOM); +} + +test the_fixed_point_runs_of_the_exact_scheme_stay_within_the_stated_lsb_of_f64 { + // RC + const m1 : [30]f64 = rc_zoh(RC_DT, PADE_ORDER); + var x1 : [5]u32 = [0] ** 5; + const r1 : RunOut = run(ZRC_N, &ZRC_D, &ZRC_NV, &m1, STEPS, &x1); + assert_eq(r1.sig, ZRC_EXP_S); + assert r1.worst < (ZRC_MAX_ERR_LSB as f64); + assert r1.worst > (ZRC_MAX_ERR_LSB as f64) - 1.0; + assert r1.peak < HEADROOM; + assert_eq(x1[0], 112283); + // RLC, Q = 1 + const m2 : [30]f64 = rlc_zoh(RLC_VAL[1], RLC_DT, PADE_ORDER); + var x2 : [5]u32 = [0] ** 5; + const r2 : RunOut = run(ZRLC_N, &ZRLC_D, &ZRLC_NV, &m2, STEPS, &x2); + assert_eq(r2.sig, ZRLC_EXP_S); + assert r2.worst < (ZRLC_MAX_ERR_LSB as f64); + assert r2.worst > (ZRLC_MAX_ERR_LSB as f64) - 1.0; + assert r2.peak < HEADROOM; + assert_eq(x2[0], NEG | 1215008); + assert_eq(x2[1], 1565780); + var x2a : [5]u32 = [0] ** 5; + const r2a : RunOut = run(ZRLC_N, &ZRLC_D, &ZRLC_NV, &m2, CHARGE_STEPS, &x2a); + assert_eq(x2a[0], 18028634); + assert_eq(x2a[1], NEG | 1475429); + assert(x2a[0] > ONE); + // against the analytic solution at 500 ns: the exact scheme's own error is the arithmetic one, within the bound + const dv : f64 = fabs64(sm_f64(x2a[0]) - an_v(RLC_VAL[1], 100.0 * RLC_DT)) * SCALE; + const dz : f64 = fabs64(sm_f64(x2a[1]) - an_z(RLC_VAL[1], 100.0 * RLC_DT)) * SCALE; + assert(dv < (ZRLC_MAX_ERR_LSB as f64)); + assert(dz < (ZRLC_MAX_ERR_LSB as f64)); + assert(dv > 3.0); + // RLC, Q = 20 + const m3 : [30]f64 = rlc_zoh(R_Q20, RLC_DT, PADE_ORDER); + var x3 : [5]u32 = [0] ** 5; + const r3 : RunOut = run(ZQ20_N, &ZQ20_D, &ZQ20_NV, &m3, STEPS, &x3); + assert_eq(r3.sig, ZQ20_EXP_S); + assert r3.worst < (ZQ20_MAX_ERR_LSB as f64); + assert r3.worst > (ZQ20_MAX_ERR_LSB as f64) - 1.0; + assert r3.peak < HEADROOM; + // the tank + const m4 : [30]f64 = rlc_zoh(R_TANK, RLC_DT, PADE_ORDER); + var x4 : [5]u32 = [0] ** 5; + const r4 : RunOut = run(ZTK_N, &ZTK_D, &ZTK_NV, &m4, STEPS, &x4); + assert_eq(r4.sig, ZTK_EXP_S); + assert r4.worst < (ZTK_MAX_ERR_LSB as f64); + assert r4.worst > (ZTK_MAX_ERR_LSB as f64) - 1.0; + assert r4.peak < HEADROOM; + // the error grows with 1 / (1 - rho) of the slowest pole, as M1's does: rho 0.9512, 0.9753, 0.9988, 0.99988 + assert(r1.worst < r2.worst); + assert(r2.worst < r3.worst); + assert(r3.worst < r4.worst); +} + +test the_k_256_power_test_holds_for_the_exact_constants_and_the_tank_needs_k_4096 { + const q1 : [25]f64 = mq_of(ZRC_N, &ZRC_D); + const q2 : [25]f64 = mq_of(ZRLC_N, &ZRLC_D); + const q3 : [25]f64 = mq_of(ZQ20_N, &ZQ20_D); + const q4 : [25]f64 = mq_of(ZTK_N, &ZTK_D); + assert(stable(ZRC_N, &q1)); + assert(stable(ZRLC_N, &q2)); + assert(stable(ZQ20_N, &q3)); + assert power_norm(ZRC_N, &q1, K_SQUARINGS) < 0.00001; + assert power_norm(ZRLC_N, &q2, K_SQUARINGS) > 0.0030; + assert power_norm(ZRLC_N, &q2, K_SQUARINGS) < 0.0031; + assert power_norm(ZQ20_N, &q3, K_SQUARINGS) > 0.87; + assert power_norm(ZQ20_N, &q3, K_SQUARINGS) < 0.88; + // the plain norm is not enough here either: the underdamped exact RLC has ||M_q||_inf above 1 + assert norm_inf(ZRLC_N, &q2) > 1.0; + assert norm_inf(ZRLC_N, &q2) < 1.1; + // the Q = 200 tank: not certified at K = 256 (1.17), certified at K = 4096 (0.83) + assert(!stable(ZTK_N, &q4)); + assert power_norm(ZTK_N, &q4, K_SQUARINGS) > 1.1; + assert power_norm(ZTK_N, &q4, K_SQUARINGS) < 1.2; + assert power_norm(ZTK_N, &q4, TANK_SQUARINGS) < 1.0; + assert power_norm(ZTK_N, &q4, TANK_SQUARINGS) > 0.8; + // and the same larger K still certifies the others + assert power_norm(ZQ20_N, &q3, TANK_SQUARINGS) < 0.01; + assert_eq(K_SQUARINGS, 8); + assert_eq(TANK_SQUARINGS, 12); +} + +test a_lossless_tank_and_a_negative_resistance_are_not_certified { + // R = 0: exp(A dt) is a rotation. Its Q.24 rounding does not decay (||M^256|| = 1.20, ||M^4096|| = 1.39) + const m0 : [30]f64 = rlc_zoh(0.0, RLC_DT, PADE_ORDER); + var d0 : [4]u32 = [0] ** 4; + var i : usize = 0; + while (i < 4) { + var v : f64 = m0[i]; + if (i == 0 or i == 3) { v = v - 1.0; } + d0[i] = q_sm(v); + i = i + 1; + } + const q0 : [25]f64 = mq_of(2, &d0); + assert(!stable(2, &q0)); + assert(power_norm(2, &q0, TANK_SQUARINGS) > 1.0); + // R = -100: growth, e^0.025 per step + const mn : [30]f64 = rlc_zoh(0.0 - RLC_VAL[1], RLC_DT, PADE_ORDER); + var dn : [4]u32 = [0] ** 4; + i = 0; + while (i < 4) { + var v : f64 = mn[i]; + if (i == 0 or i == 3) { v = v - 1.0; } + dn[i] = q_sm(v); + i = i + 1; + } + const qn : [25]f64 = mq_of(2, &dn); + assert(!stable(2, &qn)); + assert(power_norm(2, &qn, K_SQUARINGS) > 1000.0); + // a sign error in the stored D of the stable RLC is caught the same way + var flip : [4]u32 = ZRLC_D; + flip[0] = flip[0] ^ NEG; + flip[3] = flip[3] ^ NEG; + const qf : [25]f64 = mq_of(2, &flip); + assert(!stable(2, &qf)); +} + +test the_lc_tank_keeps_its_amplitude_after_100_periods_and_backward_euler_does_not { + const dt : f64 = RLC_DT; + // the run: 100 periods of the damped frequency, 125.66 steps each + assert(period_steps(R_TANK, dt) > 125.6); + assert(period_steps(R_TANK, dt) < 125.7); + assert_eq(ring_k(R_TANK, dt), 12566); + assert_eq(ring_w(R_TANK, dt), 127); + // the exact scheme, f64 + const me : [30]f64 = rlc_zoh(R_TANK, dt, PADE_ORDER); + const re : Ring = ring_f64(2, &me, 0, R_TANK, dt); + assert(re.analytic > 0.2111); + assert(re.analytic < 0.2112); + assert(near(re.scheme / re.analytic, 1.0, GRID_TOL)); + assert(near(re.scheme / re.analytic, 1.0, 0.01)); + // the exact scheme, Q.24 constants and the integer step + const rf : Ring = ring_fixed(2, &ZTK_D, &ZTK_NV, 0, R_TANK, dt); + assert(near(rf.analytic, re.analytic, 0.0000000000001)); + assert(near(rf.scheme / rf.analytic, 1.0, 0.01)); + assert(rf.scheme / rf.analytic > 1.0002); + assert(rf.scheme / rf.analytic < 1.0004); + // backward Euler, f64: 3.75e-8 of 0.2111 + const mb : [30]f64 = be_rlc(R_TANK, dt); + const rb : Ring = ring_f64(3, &mb, 1, R_TANK, dt); + assert(near(rb.analytic, re.analytic, 0.0000000000001)); + const loss : f64 = rb.scheme / rb.analytic; + assert(loss > 0.00000017); + assert(loss < 0.00000019); + assert(1.0 - loss > 0.9999998); + // the closed form: (lam2 / exp(-R dt / L)) ^ (K / 2) with lam2 = 1 / (1 + R dt / L + dt^2 / L C) + const lam2 : f64 = 1.0 / (1.0 + R_TANK * dt / RLC_VAL[2] + dt * dt / (RLC_VAL[2] * RLC_VAL[3])); + const ex : f64 = @exp(0.0 - R_TANK * dt / RLC_VAL[2]); + const cf : f64 = @exp(0.5 * (ring_k(R_TANK, dt) as f64) * @log(lam2 / ex)); + assert(cf > 0.00000015); + assert(cf < 0.00000016); + assert(loss / cf > 1.0); + assert(loss / cf < 1.25); + // backward Euler's numerical damping per step against the tank's physical one: ten times + const num : f64 = 0.5 * omega0() * omega0() * dt * dt; + const phys : f64 = alpha_of(R_TANK) * dt; + assert(num / phys > 9.9); + assert(num / phys < 10.1); + // and the Q = 20 circuit, where they are equal + assert(near(num / (alpha_of(R_Q20) * dt), 1.0, 0.0000000001)); +} + +test a_wrong_pade_order_a_wrong_dt_a_wrong_sign_or_a_dropped_integral_is_caught { + // the right constants match (asserted above); each mutation must not + const ok : [30]f64 = rlc_zoh(RLC_VAL[1], RLC_DT, PADE_ORDER); + assert_eq(mismatches(ZRLC_N, &ok, &ZRLC_D, &ZRLC_NV), 0); + // Pade (1, 1) is the trapezoid rule: other integers and a grid error 1.5e-4 of full scale + const p1 : [30]f64 = rlc_zoh(RLC_VAL[1], RLC_DT, 1); + assert(mismatches(ZRLC_N, &p1, &ZRLC_D, &ZRLC_NV) > 0); + const g1 : GridErr = grid_err(2, &p1, 0, 1, RLC_VAL[1], RLC_DT, STEPS); + assert(g1.err / g1.scale > 0.0001); + assert(g1.err / g1.scale > 100000.0 * GRID_TOL); + // Pade (2, 2) rounds to the same integers (its error is below half an LSB) and the f64 grid test alone catches it + const p2 : [30]f64 = rlc_zoh(RLC_VAL[1], RLC_DT, 2); + assert_eq(mismatches(ZRLC_N, &p2, &ZRLC_D, &ZRLC_NV), 0); + const g2 : GridErr = grid_err(2, &p2, 0, 1, RLC_VAL[1], RLC_DT, STEPS); + assert(g2.err / g2.scale > GRID_TOL); + assert(g2.err / g2.scale < 0.00000002); + // Pade (3, 3) is already inside the tolerance + const p3 : [30]f64 = rlc_zoh(RLC_VAL[1], RLC_DT, 3); + const g3 : GridErr = grid_err(2, &p3, 0, 1, RLC_VAL[1], RLC_DT, STEPS); + assert(g3.err / g3.scale < GRID_TOL); + // the RC with Pade (1, 1) + const rp1 : [30]f64 = rc_zoh(RC_DT, 1); + assert(mismatches(ZRC_N, &rp1, &ZRC_D, &ZRC_NV) > 0); + // a wrong dt: 2.5 ns and 10 ns are not 5 ns + const h : [30]f64 = rlc_zoh(RLC_VAL[1], 0.5 * RLC_DT, PADE_ORDER); + assert(mismatches(ZRLC_N, &h, &ZRLC_D, &ZRLC_NV) > 0); + const dbl : [30]f64 = rlc_zoh(RLC_VAL[1], 2.0 * RLC_DT, PADE_ORDER); + assert(mismatches(ZRLC_N, &dbl, &ZRLC_D, &ZRLC_NV) > 0); + const rh : [30]f64 = rc_zoh(0.5 * RC_DT, PADE_ORDER); + assert(mismatches(ZRC_N, &rh, &ZRC_D, &ZRC_NV) > 0); + // a negative resistance + const neg : [30]f64 = rlc_zoh(0.0 - RLC_VAL[1], RLC_DT, PADE_ORDER); + assert(mismatches(ZRLC_N, &neg, &ZRLC_D, &ZRLC_NV) > 0); + // the input term without the integral (N = B dt, forward Euler's): another N, same M + var fe : [30]f64 = ok; + const b : [2]f64 = rlc_bc(); + fe[4] = b[0] * RLC_DT; + fe[5] = b[1] * RLC_DT; + assert(mismatches(ZRLC_N, &fe, &ZRLC_D, &ZRLC_NV) > 0); + // R and L swapped in the stamp: the damping row divides by C instead of L + var sw : [4]f64 = rlc_ac(RLC_VAL[1]); + sw[3] = 0.0 - RLC_VAL[1] / RLC_VAL[3]; + const bb : [2]f64 = rlc_bc(); + const swm : [30]f64 = zoh(2, &sw, &bb, RLC_DT, PADE_ORDER); + assert(mismatches(ZRLC_N, &swm, &ZRLC_D, &ZRLC_NV) > 0); + // a tank constant set used for the Q = 20 circuit + const m20 : [30]f64 = rlc_zoh(R_Q20, RLC_DT, PADE_ORDER); + assert(mismatches(ZQ20_N, &m20, &ZQ20_D, &ZQ20_NV) == 0); + assert(mismatches(ZQ20_N, &m20, &ZTK_D, &ZTK_NV) > 0); +}