Reference › Numerics — what the solver will and will not do
kind: reference#numerics#singular#inverse#eigenvalues#decomposition#discretization#zoh#fixed-point#q16.16#hdl#nan#inf#integer#double

Numerics — what the solver will and will not do#

Everything ICoreBlocks computes — a console expression, a block's state update, a discretized model — is done in IEEE double precision, and every matrix operation goes through one linear-algebra layer whose conventions this page states as you meet them: what a singular matrix gives you, which decomposition orderings and signs you get, which discretization spellings are understood, and what precision the HDL exports carry. The equations themselves are on Solver mathematics — the forms, the discretizations and the stepping scheme; the pacing of a run is on Sample time and loops — how the simulator paces a diagram. Every claim below is either read from the named source file or observed in the console runs pasted at the end (binary built 2026-09-25 from approximately commit 1dac24c57).

The rules you meet#

  • Every value is a double. Matrices store double (ICoreMatrix.cpp, std::vector<double> data); a scalar is a 1×1 matrix. The console's # Integer / # Double annotation is a display verdict, not a type: an anonymous result is labelled Integer when it is finite and has no fractional part (ICoreCommandWindow.cpp, formatResultRepr); a named variable is labelled from its text — x = 1e6 says Double because 1e6 does not parse as an integer literal (ICoreVariable.cpp, assignType). A complex result is a pair of doubles and is labelled Complex.
  • A variable is stored at full precision and shown at six digits. The variables space holds text, so what an assignment writes is what every later expression reads. It writes enough digits to name the double exactly (17 significant, max_digits10), while the console still prints six — so x = 1/3 echoes 0.333333 and x * 3 is 1, and P = inv(S) followed by P * S is the identity to rounding (1e-17). Before 2026-09-02 the store wrote six digits too, and both of those came back wrong in the sixth place. The same split applies to a transfer function, a state space, a polynomial and a recorded signal: what is stored round-trips, what is printed is rounded for reading.
  • A matrix operation the console cannot perform either refuses or logs. A wrong shape is refused with an error: line and exit code 1 (det([1 2 3]) → error: det() needs a square matrix; inv, eig and solve the same). A failure inside a computation is reported through the run diagnosis ([ICoreRunDiagnosis] - Error - … on the console) and returns the default matrix, which is a 1×1 zero — printed as 0 # Integer, exit code 0: [1 2; 2 4]^-1 does that. Check the log line, not the value. chol is the pattern to copy: it REFUSES by name, because "the default matrix is a 1×1 zero" and "the answer is a 1×1 zero" are the same value, and reading one as the other is how chol of a non-positive-definite matrix answered a silent 0 for as long as the console had it.
  • inv(A) of a singular matrix answers all Inf, as MATLAB does, and solve(A,b), A \ b, pinv(A), det(A), cond(A), rank(A) do not complain either. inv runs a partial-pivot LU and, when a pivot is exactly zero, returns a matrix of Inf with no message (ICoreMatrix.cpp, inverseOrInf). solve() is a column-pivoted QR and returns some vector for a singular or inconsistent system with no message; pinv() is an SVD pseudo-inverse with MATLAB's automatic tolerance max(rows, cols) · eps(σ_max), where eps(σ_max) is the spacing of doubles at the largest singular value; det() returns whatever the LU produces (a signed zero for a singular matrix); cond() returns Inf when the smallest singular value is exactly 0; rank() counts singular values above the same tolerance as pinv.
  • An ill-conditioned matrix inverts without a warning. inv([1 1; 1 1.0000000001]) returns entries of 1e+10 and says nothing; only cond() (4e+10) tells you. There is no automatic conditioning check anywhere on the inverse path.
  • / by a square matrix is multiplication by its inverse — and that inverse, unlike the console's inv, refuses a singular matrix: [1 2] / [1 2; 2 4] logs Unable to compute matrix inverse. Matrix is singular. and then stops with a dimension error (ICoreExpressionEvaluator.cpp, matrixDiv). / by a non-square matrix is (B' \ A')'. \ is left division: A \ b solves A x = b by the same column-pivoted QR as solve (least squares for a tall A), and refuses an underdetermined (wide) A, pointing you to pinv(A) * b.
  • Division by zero follows IEEE, as in MATLAB. 1/0 is Inf, 0/0 is NaN, and [1 2] ./ 0 and [1 2] ./ [0 0] are both [[Inf, Inf]] — all exit code 0.
  • NaN and Inf otherwise propagate as IEEE says. log(0) is -Inf and 1e308*10 is Inf — # Double, exit code 0, and a variable assigned one prints back as x = NaN # Double. sqrt(-1) leaves the real line instead: it is 0 + 1i # Complex, as in MATLAB.
  • ^ is the matrix power and .^ the element-wise one. 2^62 is 4.61169e+18; A^n for a square A and an integer n is repeated multiplication, a negative n going through the inverse (so a singular A logs and answers the zero placeholder above). mpower(A,n) and pow(A,n) are the same two operations by name.

Decomposition conventions — what you get, and in what order#

All decompositions delegate to Eigen; the conventions below are what the delegation makes observable (ICoreMatrix.cpp, functions named).

FunctionConvention observed
eig(A)Eigenvalues in the order Eigen's general solver produces them — not sorted. Real spectrum → N×1 column; any complex value → a complex N×1 column (eig([0 -1; 1 0]) = [0 + 1i; 0 - 1i] # Matrix of Complex). Run: eig([2 0; 0 1]) = [[2], [1]]; eig([2 1; 1 2]) = [[3], [1]].
eigvec(A)Columns are eigenvectors, real part only. A symmetric matrix takes the self-adjoint solver, whose eigenvalues are ascending — so for a symmetric A the column order of eigvec(A) is the reverse of eig(A)'s row order (run: eigvec([2 1; 1 2]) = [[-0.707107, 0.707107], [0.707107, 0.707107]], the first column belonging to eigenvalue 1). A non-symmetric matrix keeps eig's order. A complex eigenpair logs eigenvectors(): matrix has complex eigenpairs; returning real part only. and hands back the real parts.
qrq(A), qrr(A)Thin Householder QR, A = Q·R. Signs are Householder's: the diagonal of R may be negative (qrr([1 2; 3 4]) = [[-3.16228, -4.42719], [0, -0.632456]]).
lul(A), luu(A), lup(A)Partial-pivot LU, P·A = L·U, unit-diagonal L. lup([1 2; 3 4]) = [[0, 1], [1, 0]]. Square only.
svdu(A), svd(A), svdv(A)Thin Jacobi SVD, A = U·diag(S)·Vᵀ. svd(A) is a column of non-negative singular values in descending order; the sign of a singular triple is not fixed (svdu([3 0; 0 -2]) = [[1, -0], [0, -1]], svdv = identity).
svd(A), qr(A), lu(A), ldl(A) with several outputs[U, S, V] = svd(A) is the FULL factorization (S is m×n, not the vector), svd(A, "econ") the economy one; [Q, R] = qr(A) is full and [Q, R, P] = qr(A) is the column-pivoted A·P = Q·R; [L, U] = lu(A) folds the permutation into L so L·U = A, and [L, U, P] = lu(A) keeps it so L·U = P·A; [L, D, P] = ldl(A) states the permutation MATLAB's way, Pᵀ·A·P = L·D·Lᵀ. One output is what MATLAB gives for one: qr(A) is R, lu(A) writes both factors into one matrix, ldl(A) is L.
chol(A)Upper factor R, A = Rᵀ·R, positive diagonal (chol([4 2; 2 3]) = [[2, 1], [0, 1.41421]]); chol(A, "lower") is L = Rᵀ. A non-SPD input is refused by name, and [R, p] = chol(A) is the form that does not fail: p is 0 for a positive-definite A and otherwise the pivot the factorization stopped at, with R the leading block that succeeded.
ldltl/ldltd/ldltp(A)Pivoted LDLᵀ, PᵀLDLᵀP = A; ldltd is the diagonal as a column. The console's D is diagonal where MATLAB's may hold 2×2 blocks.
null(A)Orthonormal basis for the kernel, as columns, from the SVD — the columns of V whose singular values fall at or below max(size(A))·eps(largest). null(A, "r") is the rational basis read off the reduced row echelon form, whose entries are the small whole numbers a hand calculation gives.

Because these are conventions and not unique answers, the console regression corpus pins them as golden cases; a different sign in your own derivation is not a defect in either.

Discretization — which spellings are understood, and what a typo does#

  • In the app, the method is a fixed choice, not free text. The seven options are exactly Zero-order Hold, First-order Hold, Impulse, Tustin, Matched, Backward Euler, Forward Euler (ICoreModelConfigurator.cpp, DISCRETE_METHOD_*); the model-configuration panel offers them in a combo box, and the console's setModelConfig discretizationMethod <value> accepts an option case-insensitively, by any |-separated segment, or by a unique case-insensitive substring — and rejects anything else with the list (ICoreModelConfigCommands.cpp, resolveChoice). Run: setModelConfig discretizationMethod tustn → no option matches 'tustn' — available: …; tustin → discretizationMethod = Tustin. Note zoh is not a substring of Zero-order Hold and is rejected; type Zero-order Hold (spaces are fine) or the unique substring zero; hold is reported ambiguous.
  • Below the option list, the string is matched loosely and an unknown string silently becomes ZOH. The conversion routine lower-cases the method and recognises foh / first-order hold / firstorderhold / first order hold, impulse / imp, tustin, matched / match, backwardeuler / backward-euler / backward euler, forwardeuler / forward-euler / forward euler; everything else, including zoh and any misspelling, takes the Zero-order Hold branch with no diagnostic (ICoreStateSpaceDiscretization.cpp, discretizeStateSpace). The de-discretization routine recognises only tustin, matched/match, forwardeuler variants, and defaults to ZOH inversion.
  • Where the loose match is reachable by a user: the project file. solver.ini's discretizationMethod key is read back through a setter that does not validate against the option list (ICoreStudioSerialization.cpp reads it; ICoreModelConfigurator::setDiscretizationMethod assigns unconditionally), and the per-block dispatch then finds no matching isDiscretizationMethod_* and takes the // Default ZOH branch (ICoreBlockSolverEnvironment.cpp, discretize). A hand-edited or foreign-version project with discretizationMethod=Tustn therefore simulates as Zero-order Hold with no message. Blocks with their own per-block method (Discrete Nonlinear State Space, Discrete Nonlinear State Space — Control Systems/Discrete) use the same seven option strings and the same silent ZOH default. Which model each method produces is on Solver mathematics — the forms, the discretizations and the stepping scheme.

Fixed point on the HDL path — what VHDL / Verilog / SystemVerilog give you#

  • One format, Q16.16, for every signal. VHDL: subtype Fx is sfixed(ICORE_INT_BITS-1 downto -ICORE_FRAC_BITS) with both constants 16 in the generated icore_pkg.vhd (ICoreVHDLParser.cpp, buildSupportPackage). Verilog / SystemVerilog: ` ICORE_INT_BITS 16 `, ICORE_FRAC_BITS 16 `, ICORE_WIDTH ` = 32 in icore_defs.vh (ICoreVerilogParser.cpp, buildDefs`). That is 16 integer bits including sign — a range of −32768 to +32767.99998 — with a quantum of 2⁻¹⁶ ≈ 1.53e-5. Every signal port is a fixed-size bus of these; there is no per-block or per-signal word length.
  • Rounding at the boundary. Constants and testbench stimuli enter through to_fx: VHDL's to_sfixed rounds; the two Verilogs add ±0.5 before $rtoi, i.e. round-to-nearest, ties away from zero (the comment in fixedPointFunctions records why bare truncation was replaced). Values leave through to_real = value / 2^16. Fixed-point → count conversions (fx_to_int, used for indices, delays, periods) floor by arithmetic shift on all three targets, to agree with the software targets' floor().
  • Inside a block, products are formed at double width and shifted back once. The generated core declares reg signed [2*ICORE_WIDTH-1:0] acc (ICoreVerilogParser.cpp, core builder); a synthesizable block body multiplies into acc (Q32.32) and assigns acc >>> ICORE_FRAC_BITS to the 32-bit signal (read in, e.g., Gear_Train's generateBodyCode_Verilog). The >>> floors (toward −∞), while VHDL bodies use resize(a*b, acc), which rounds — so the two HDL families can differ by one quantum on a value that sits exactly off the lattice, and an antisymmetric operation like q(-v) vs -q(v) shows it. Overflow is not saturated in the Verilog bodies read: the shifted product is assigned to a narrower bus and the high bits are dropped (wrap-around). Whether VHDL's resize saturates was not verified from the generated text and is not claimed here.
  • "Synthesizable" versus "simulation-only". A block whose HDL body is written in Q16.16 integer arithmetic is synthesizable. A block that needs sin/cos, a random generator or other real-valued math has HDL bodies that carry the arithmetic in the simulator's real type (VHDL ieee.math_real, Verilog real), quantizing only at the port boundary — the generated file says SIMULATION-ONLY real arithmetic in its banner, and the block's page says so (e.g. Planar Arm Forward Kinematics — Robotics/Planar Kinematics; Constant — Control Systems/Sources is the synthesizable case). Such output runs in an HDL simulator and passes verification, but a synthesis tool will not build it. 81 block sources carry the simulation-only marker on 2026-08-17. Every software target (C, C++, Python, Java, Rust, MATLAB, PLC ST) uses double / LREAL throughout, so the HDL family is the only place quantization exists. How to export is on Exporting code — the ten targets, what each produces, and what verification proves.

Things that surprise people#

  • "inv printed Inf and my model kept going." A singular matrix inverts to all Inf with no message, and anything downstream computes with it; a negative matrix power of one logs an [ICoreRunDiagnosis] - Error line and hands back a zero instead. In a script check det/cond/rank first, or use pinv.
  • "solve gave me an answer for a singular system." It always does — column-pivoted QR returns a vector even when none or infinitely many exist (solve([1 2; 2 4], [1; 1]) = [[0], [0.3]], and [1 2; 2 4] \ [1; 1] the same). Check the residual A*x - b yourself.
  • "eig and eigvec disagree on the order." Only for symmetric matrices, and always: the symmetric path sorts ascending, eig does not sort. Pair them by re-checking A*v = λ*v, or read eig off the diagonal of the same solver by asking for eigvec and multiplying.
  • "My project silently simulates as ZOH." Its discretizationMethod text is not one of the seven exact spellings — pick the method again in the panel or with setModelConfig.
  • "VHDL and Verilog verify to different residuals on the same diagram." Round versus floor on the product shift; a one-quantum (1.5e-5) disagreement on values off the Q16.16 lattice is the datapath, not the block.
  • "My HDL export runs in the simulator but will not synthesize." The block is simulation-only in HDL; its page says so.
  • "det of an integer matrix came back with a fraction on the end." The determinant goes through a pivoted LU, not the schoolbook cofactor formula, so an integer matrix does not guarantee an integer answer: det([7 3; 2 5]) is 29 + 3.55271e-15 and det([2 3 1; 4 7 2; 6 18 5]) is 4 - 8.88178e-16. Both print as 29 and 4 — the console's Integer/Double tag is what gives it away. Round if you need an exact integer, and do not test a determinant for equality. Measured 2026-09-26.

Real runs (2026-09-26)#

All runs: HOME=<scratch> ICoreBlocks.app/Contents/MacOS/ICoreBlocks --console "<line>", binary built 2026-09-25 (≈ commit 1dac24c57), startup lines stripped, everything else verbatim.

--console "inv([1 2; 2 4])"        →  [[Inf, Inf], [Inf, Inf]]  # Matrix of Double
--console "x = inv([1 2;2 4])"     →  x = [[Inf, Inf], [Inf, Inf]]  # Matrix of Double
--console "det([1 2; 2 4])"        →  -0  # Integer
--console "solve([1 2; 2 4], [1; 2])"  →  [[0], [0.5]]  # Matrix of Double
--console "solve([1 2; 2 4], [1; 1])"  →  [[0], [0.3]]  # Matrix of Double
--console "[1 2; 2 4] \ [1; 1]"    →  [[0], [0.3]]  # Matrix of Double
--console "pinv([1 2; 2 4])"       →  [[0.04, 0.08], [0.08, 0.16]]  # Matrix of Double
--console "cond([1 2; 2 4])"       →  Inf  # Double
--console "rank([1 2; 2 4])"       →  1  # Integer
--console "inv([1 1; 1 1.0000000001])"  →  [[1e+10, -1e+10], [-1e+10, 1e+10]]  # Matrix of Double
--console "cond([1 1; 1 1.0000000001])" →  4e+10  # Double
--console "[1 2] / [1 2; 2 4]"
[ICoreRunDiagnosis] - Error - Unable to compute matrix inverse. Matrix is singular.
error: '/' inner dimensions disagree: 1x2 / 2x2      (exit 1)
--console "[1 2; 2 4]^-1"
[ICoreRunDiagnosis] - Error - Unable to compute matrix inverse. Matrix is singular.
[ICoreRunDiagnosis] - Error - matmul dimension mismatch: inner dims 2 != 1
0  # Integer                              (exit 0)
--console "[1 2 3; 4 5 6] \ [1; 2]"
error: '\' on an underdetermined system (2x3) has no single answer: MATLAB returns a basic solution with at most rank(A) non-zeros -- write pinv(A) * b for the minimum-norm one      (exit 1)
--console "det([1 2 3])"           →  error: det() needs a square matrix   (exit 1)

--console "1/0"          →  Inf  # Double
--console "0/0"          →  NaN  # Double
--console "[1 2] ./ 0"   →  [[Inf, Inf]]  # Matrix of Double
--console "[1 2] ./ [0 0]" →  [[Inf, Inf]]  # Matrix of Double
--console "sqrt(-1)"     →  0 + 1i  # Complex
--console "log(0)"       →  -Inf  # Double
--console "1e308*10"     →  Inf  # Double
--console "x = 0/0"      →  x = NaN  # Double
--console "2/1"          →  2  # Integer
--console "5/2"          →  2.5  # Double
--console "0.1+0.2"      →  0.3  # Double
--console "1e6"          →  1e+06  # Integer
--console "x = 1e6"      →  x = 1e+06  # Double
--console "2^62"         →  4.61169e+18  # Integer
--console "[1 2; 3 4]^2" →  [[7, 10], [15, 22]]  # Matrix of Double
--console "x = 1/3; x * 3"  →  1  # Integer

--console "eig([2 0; 0 1])"        →  [[2], [1]]  # Matrix of Double
--console "eigvec([2 0; 0 1])"     →  [[0, 1], [1, 0]]  # Matrix of Double
--console "eig([2 1; 1 2])"        →  [[3], [1]]  # Matrix of Double
--console "eigvec([2 1; 1 2])"     →  [[-0.707107, 0.707107], [0.707107, 0.707107]]  # Matrix of Double
--console "eig([1 2; 0 3])"        →  [[1], [3]]  # Matrix of Double
--console "eigvec([1 2; 0 3])"     →  [[1, 0.707107], [0, 0.707107]]  # Matrix of Double
--console "eig([0 -1; 1 0])"       →  [0 + 1i; 0 - 1i]  # Matrix of Complex
--console "eigvec([0 -1; 1 0])"
[ICoreRunDiagnosis] - Warning - eigenvectors(): matrix has complex eigenpairs; returning real part only.
[[-0.707107, -0.707107], [0, 0]]  # Matrix of Double
--console "eig([1 2 3; 4 5 6; 7 8 10])" →  [[16.7075], [-0.90574], [0.198247]]  # Matrix of Double

--console "qrq([1 2; 3 4])"   →  [[-0.316228, -0.948683], [-0.948683, 0.316228]]  # Matrix of Double
--console "qrr([1 2; 3 4])"   →  [[-3.16228, -4.42719], [0, -0.632456]]  # Matrix of Double
--console "lul([1 2; 3 4])"   →  [[1, 0], [0.333333, 1]]  # Matrix of Double
--console "luu([1 2; 3 4])"   →  [[3, 4], [0, 0.666667]]  # Matrix of Double
--console "lup([1 2; 3 4])"   →  [[0, 1], [1, 0]]  # Matrix of Double
--console "svd([3 0; 0 -2])"  →  [[3], [2]]  # Matrix of Double
--console "svdu([3 0; 0 -2])" →  [[1, -0], [0, -1]]  # Matrix of Double
--console "svdv([3 0; 0 -2])" →  [[1, 0], [0, 1]]  # Matrix of Double
--console "chol([4 2; 2 3])"  →  [[2, 1], [0, 1.41421]]  # Matrix of Double
--console "chol([4 2; 2 3], \"lower\")" →  [[2, 0], [1, 1.41421]]  # Matrix of Double
--console "chol([1 2; 2 1])"
error: chol(A): the matrix is not positive definite (the factorization stops at pivot 2) -- [R, p] = chol(A) answers where it failed instead of refusing
--console "null([1 1; 1 1])"          →  [[-0.707107], [0.707107]]  # Matrix of Double
--console "null([1 1; 1 1], \"r\")"    →  [[-1], [1]]  # Matrix of Double
--console "qr([1 2; 3 4])"            →  [[-3.16228, -4.42719], [0, -0.632456]]  # Matrix of Double
--console "lu([1 2; 3 4])"            →  [[3, 4], [0.333333, 0.666667]]  # Matrix of Double
--console "ldl([4 2; 2 3])"           →  [[1, 0], [0.5, 1]]  # Matrix of Double
--console "ldltd([4 2; 2 3])"         →  [[4], [2]]  # Matrix of Double
--console "svds([3 0; 0 -2])"
error: svds() is MATLAB's sparse partial SVD, which the console does not have; svd(A) is every singular value as a column
--console "det([7 3; 2 5]) - 29"              →  3.55271e-15  # Double
--console "det([2 3 1; 4 7 2; 6 18 5]) - 4"   →  -8.88178e-16  # Double

--console "setModelConfig discretizationMethod tustn"
setModelConfig discretizationMethod: no option matches 'tustn' — available:
  Zero-order Hold
  First-order Hold
  Impulse
  Tustin
  Matched
  Backward Euler
  Forward Euler
--console "setModelConfig discretizationMethod tustin"           →  discretizationMethod = Tustin
--console "setModelConfig discretizationMethod zoh"              →  (rejected, same list)
--console "setModelConfig discretizationMethod Zero-order Hold"  →  discretizationMethod = Zero-order Hold
--console "setModelConfig discretizationMethod zero"             →  discretizationMethod = Zero-order Hold
--console "setModelConfig discretizationMethod hold"
setModelConfig discretizationMethod: 'hold' is ambiguous — matches:
  Zero-order Hold
  First-order Hold

For contributors#

The module page for the numerics zone holds the invariants behind this page (the Eigen wrapper rule, the recognised-method list, the golden decomposition corpus, the ICoreTimeResponse discrete-only trap); the simulator page holds the discretize-on-build behaviour that reads the method; the code-generation page holds the HDL parser map. All three are internal pages listed in this page's related: front-matter.