diff --git a/lib/scikit/linear.flow b/lib/scikit/linear.flow
index 7194d98..feda9e5 100644
--- a/lib/scikit/linear.flow
+++ b/lib/scikit/linear.flow
@@ -2087,6 +2087,328 @@ export function logistic_regression_free(model: LogisticRegression) -> void {
     array_free_f32(model.classes)
 }
 
+# ============================================================================
+# Logistic coefficient inference (Issue #458)
+#
+# The deferred half of #353. PR #446 shipped OLS inference and left logistic
+# models out; PR #443's Newton solver forms the matrix this needs, the observed
+# information at the fitted coefficients.
+#
+# CONVENTION. For a binary logit the asymptotic covariance of the maximum
+# likelihood coefficients is the inverse of the observed information
+#
+#     I(theta) = sum_i p_i (1 - p_i) xt_i xt_i^T,        cov = I^-1
+#
+# over the augmented rows xt_i = [x_i, 1], with p_i the fitted probability of
+# the positive class. Written with the MEAN Hessian H = I / m that
+# lin357_newton_fit assembles, the same quantity is H^-1 / m. The two are equal
+# and this code computes the unscaled form directly. It is what statsmodels
+# reports as Logit.cov_params() and what R's glm reports as vcov(), so these
+# numbers are checkable against either.
+#
+# The information matrix is never formed. sqrt(p_i (1 - p_i)) scales row i of
+# the augmented design into A, and then I = A^T A, so the QR of A gives the
+# triangular factor R with R^T R = I and cov = R^-1 R^-T. That is exactly the
+# route ols_inference_fit takes, reusing _solve_lstsq_qr and
+# lin353_invert_upper unchanged, and it avoids squaring the condition number
+# the way an explicit Hessian and a Gaussian solve would.
+#
+# INDEX CONVENTION follows ols_inference_fit: 0..n_features-1 are the slopes in
+# feature order and index n_features is the intercept. There is no p value, for
+# the same reason #446 shipped none: it needs a normal CDF accurate in the
+# tail, and a hand-rolled one is worse than none.
+#
+# OWNERSHIP. Nothing is retained on the fit. LogisticRegression gains no field,
+# logistic_regression_fit is untouched, and this call recomputes the
+# probabilities at the fitted coefficients into buffers it owns and
+# logistic_inference_free releases. That is #446's answer to the ownership
+# question #443 left open on #353, applied unchanged.
+#
+# WHAT IS REFUSED, each with fitted = false and zeroed buffers.
+#
+#   Penalized fits. The Newton Hessian carries the alpha * I term, so its
+#   inverse is a PENALIZED covariance and not the sampling covariance a reader
+#   of a standard error expects. ols_inference_fit refuses the same case; this
+#   follows it. After #430 the canonical benchmark always fits penalized, so
+#   this surface is for callers who fit with penalty_none().
+#
+#   More than two classes. Above two classes the fit is the joint multinomial
+#   softmax of #433, whose information matrix is a (K-1)d x (K-1)d block system
+#   with cross-class blocks. That is a substantially bigger job than the binary
+#   case and is deliberately out of scope, as #458 scopes it.
+#
+#   A separated fixture. Complete separation has no maximum likelihood
+#   estimate: the coefficients diverge and any standard error read off the
+#   information matrix records where the optimizer stopped rather than a
+#   sampling variance. It is reported through separable = true rather than as a
+#   very large finite number. The test is the fitted hyperplane classifying
+#   every point correctly, which is exactly what separation means.
+#
+#   A singular weighted design, reported through rank_deficient = true, on the
+#   same relative cutoff (1e-11 * max|R_ii|, floored at 1e-15) that #446 and
+#   #450 made standard in this file.
+# ============================================================================
+
+export struct LogisticInference {
+    cov: ptr<f32>,
+    std_errors: ptr<f32>,
+    z_stats: ptr<f32>,
+    n_features: i32,
+    n_params: i32,
+    n_samples: i32,
+    rank: i32,
+    rank_deficient: bool,
+    separable: bool,
+    fitted: bool
+}
+
+# The binomial variance p(1-p) at a logit, as e / (1+e)^2 with e = exp(-|z|).
+# Computing p first and then p * (1 - p) loses the whole quantity to
+# cancellation once p rounds to 1; this form stays accurate into the tails and
+# underflows smoothly to zero.
+#
+# exp() lives in a helper rather than at the call site because the transpiler
+# miscompiles it inside a nested while loop (compiler bug #421). The name
+# carries the module and issue number because non-exported functions sharing a
+# name across modules collide in the generated C (compiler bug #465).
+function lin458_variance_f64(z: f64) -> f64 {
+    let mut az: f64 = z
+    if az < 0.0 { az = 0.0 - z }
+    let e: f64 = exp(0.0 - az)
+    let s: f64 = 1.0 + e
+    return e / (s * s)
+}
+
+# A refused result: zeroed buffers, fitted = false, and the flag that says why.
+function lin458_empty(n: i32, d: i32, rank: i32, deficient: bool, separable: bool) -> LogisticInference {
+    let cov: ptr<f32> = array_new_f32(d * d + 16)
+    let se: ptr<f32> = array_new_f32(d + 8)
+    let zs: ptr<f32> = array_new_f32(d + 8)
+    return LogisticInference {
+        cov: cov,
+        std_errors: se,
+        z_stats: zs,
+        n_features: n,
+        n_params: d,
+        n_samples: 0,
+        rank: rank,
+        rank_deficient: deficient,
+        separable: separable,
+        fitted: false
+    }
+}
+
+# Coefficient inference for an unpenalized binary logistic fit.
+#
+# model must be the result of fitting the same X and y that are passed here,
+# with penalty_none(). The coefficients described are the ones
+# logistic_decision_function reports, that is the model for classes[1]. Class
+# order comes from distinct_f32, which lists labels in order of first
+# appearance in y, so classes[1] is the second distinct label seen and not
+# necessarily the larger one.
+export function logistic_inference_fit(model: LogisticRegression, X: Matrix, y: ptr<f32>) -> LogisticInference {
+    let n: i32 = X.cols
+    let m: i32 = X.rows
+    let d: i32 = n + 1
+
+    if n <= 0 || m <= 0 {
+        let bad0: LogisticInference = lin458_empty(n, d, 0, false, false)
+        return bad0
+    }
+    if model.fitted == false || model.n_features != n {
+        let bad1: LogisticInference = lin458_empty(n, d, 0, false, false)
+        return bad1
+    }
+    # Penalized: refused, see the header.
+    if model.penalty.kind != PENALTY_NONE {
+        let bad2: LogisticInference = lin458_empty(n, d, 0, false, false)
+        return bad2
+    }
+    # Multiclass: refused, out of scope. n_classes < 2 means y carried a single
+    # label, which has no identified coefficient either.
+    if model.n_classes != 2 {
+        let bad3: LogisticInference = lin458_empty(n, d, 0, false, false)
+        return bad3
+    }
+    # Fewer rows than parameters cannot support a full rank information matrix.
+    if m < d {
+        let bad4: LogisticInference = lin458_empty(n, d, 0, true, false)
+        return bad4
+    }
+
+    let w: ptr<f32> = model.weights + n
+    let bias: f32 = model.biases[1]
+    let positive: f32 = model.classes[1]
+
+    let xdata: ptr<f32> = X.data
+    let stride: i32 = X.cols
+
+    # Logits at the fitted coefficients. No exp() in this loop: it is nested.
+    let zbuf: ptr<f64> = array_new_f64(m + 8)
+    for i in 0 to m {
+        let base: i32 = i * stride
+        let mut z: f64 = bias as f64
+        for j in 0 to n {
+            z = z + (w[j] as f64) * (xdata[base + j] as f64)
+        }
+        zbuf[i] = z
+    }
+
+    # Row weights, and the separation test.
+    #
+    # Complete separation IS the condition that some hyperplane classifies
+    # every point correctly. When one exists, scaling that hyperplane up drives
+    # the log likelihood towards zero without ever attaining it, the maximum
+    # likelihood estimate does not exist, and no covariance describes a
+    # coefficient that is not there. So the test is a property of the data, and
+    # the fitted hyperplane is the witness: if every z_i has the sign of its
+    # label, the data are separable.
+    #
+    # This is exact where a threshold is not. An earlier draft asked whether
+    # every fitted probability was within 1e-4 of its label and whether the
+    # total binomial variance had collapsed. On a six point separable fixture
+    # the Newton fit tripped both, and the LBFGS fit of the same data stopped
+    # at a slope of 8.3 and tripped neither, so the same data got standard
+    # errors from one solver and a refusal from the other. The sign test
+    # catches both fits, because both found a separating hyperplane.
+    #
+    # A point with z_i exactly zero sits on the boundary and is counted as not
+    # separated, which is the conservative direction: it answers rather than
+    # refuses.
+    let sw: ptr<f64> = array_new_f64(m + 8)
+    let mut separable: bool = true
+    for i in 0 to m {
+        let zi: f64 = zbuf[i]
+        sw[i] = sqrt(lin458_variance_f64(zi))
+        if y[i] == positive {
+            if zi <= 0.0 { separable = false }
+        } else {
+            if zi >= 0.0 { separable = false }
+        }
+    }
+
+    if separable {
+        array_free_f64(zbuf)
+        array_free_f64(sw)
+        let bad5: LogisticInference = lin458_empty(n, d, 0, true, true)
+        return bad5
+    }
+
+    # A = diag(sqrt(v)) * [X | 1], so that A^T A is the observed information.
+    let A: ptr<f64> = malloc((m as i64) * (d as i64) * 8 + 128) as ptr<f64>
+    let rhs: ptr<f64> = array_new_f64(m + 8)
+    for i in 0 to m {
+        let base: i32 = i * stride
+        let dst: i32 = i * d
+        let s: f64 = sw[i]
+        for j in 0 to n {
+            A[dst + j] = s * (xdata[base + j] as f64)
+        }
+        A[dst + n] = s
+        rhs[i] = 0.0
+    }
+
+    # Leaves R, the triangular factor of A, in the upper triangle of A. rhs is
+    # the zero vector and its solution is discarded; only R is wanted.
+    _solve_lstsq_qr(A, rhs, m, d)
+
+    let mut maxdiag: f64 = 0.0
+    for j in 0 to d {
+        let dj: f64 = fabs(A[j * d + j])
+        if dj > maxdiag { maxdiag = dj }
+    }
+    let mut tol: f64 = 0.00000000001 * maxdiag
+    if tol < 0.000000000000001 { tol = 0.000000000000001 }
+
+    let Z: ptr<f64> = array_new_f64(d * d + 16)
+    let rank: i32 = lin353_invert_upper(A, Z, d, tol)
+
+    if rank < d {
+        array_free_f64(zbuf)
+        array_free_f64(sw)
+        free(A as ptr<void>)
+        array_free_f64(rhs)
+        array_free_f64(Z)
+        let bad6: LogisticInference = lin458_empty(n, d, rank, true, false)
+        return bad6
+    }
+
+    # cov = (A^T A)^-1 = Z Z^T, symmetric.
+    let cov: ptr<f32> = array_new_f32(d * d + 16)
+    let se: ptr<f32> = array_new_f32(d + 8)
+    let zs: ptr<f32> = array_new_f32(d + 8)
+    for a in 0 to d {
+        for b in 0 to d {
+            let mut s: f64 = 0.0
+            for k in 0 to d {
+                s = s + Z[a * d + k] * Z[b * d + k]
+            }
+            cov[a * d + b] = s as f32
+        }
+    }
+
+    for j in 0 to d {
+        let v: f32 = cov[j * d + j]
+        let mut sj: f32 = 0.0
+        if v > 0.0 { sj = sqrt((v as f64)) as f32 }
+        se[j] = sj
+        let mut coef: f32 = bias
+        if j < n { coef = w[j] }
+        if sj > 0.0 {
+            zs[j] = coef / sj
+        } else {
+            zs[j] = 0.0
+        }
+    }
+
+    array_free_f64(zbuf)
+    array_free_f64(sw)
+    free(A as ptr<void>)
+    array_free_f64(rhs)
+    array_free_f64(Z)
+
+    return LogisticInference {
+        cov: cov,
+        std_errors: se,
+        z_stats: zs,
+        n_features: n,
+        n_params: d,
+        n_samples: m,
+        rank: rank,
+        rank_deficient: false,
+        separable: false,
+        fitted: true
+    }
+}
+
+# Standard error of coefficient j. j in 0..n_features-1 is a slope, j ==
+# n_features is the intercept.
+export function logistic_inference_std_error(inf: LogisticInference, j: i32) -> f32 {
+    if j < 0 || j >= inf.n_params { return 0.0 }
+    return inf.std_errors[j]
+}
+
+# z statistic of coefficient j against the null that it is zero. The logistic
+# analogue of ols_inference_t_stat: the reference distribution is normal, not
+# t, because the variance is not estimated separately from the coefficients.
+export function logistic_inference_z_stat(inf: LogisticInference, j: i32) -> f32 {
+    if j < 0 || j >= inf.n_params { return 0.0 }
+    return inf.z_stats[j]
+}
+
+# Entry (a, b) of the coefficient covariance matrix.
+export function logistic_inference_cov(inf: LogisticInference, a: i32, b: i32) -> f32 {
+    if a < 0 || a >= inf.n_params { return 0.0 }
+    if b < 0 || b >= inf.n_params { return 0.0 }
+    return inf.cov[a * inf.n_params + b]
+}
+
+export function logistic_inference_free(inf: LogisticInference) -> void {
+    array_free_f32(inf.cov)
+    array_free_f32(inf.std_errors)
+    array_free_f32(inf.z_stats)
+}
+
 # ============================================================================
 # Ridge Regression (L2-penalized linear regression with closed-form-ish GD)
 # ============================================================================
diff --git a/tests/test_logistic_inference.flow b/tests/test_logistic_inference.flow
new file mode 100644
index 0000000..7a9edd5
--- /dev/null
+++ b/tests/test_logistic_inference.flow
@@ -0,0 +1,637 @@
+# Tests for binary logistic coefficient inference (Issue #458).
+# Run: flow run tests/test_logistic_inference.flow
+#
+# The convention under test is cov = I^-1 with I = sum_i p_i (1 - p_i) xt_i
+# xt_i^T, the unscaled observed information at the fitted coefficients. Two
+# independent references pin it:
+#
+#   Fixture A is checked against statsmodels 0.14.6,
+#   Logit(y, add_constant(X)).fit(method=newton, tol=1e-13).cov_params().
+#
+#   Fixture B is the saturated 2x2 table, whose maximum likelihood estimate and
+#   covariance are closed form and worked out by hand in the comment above the
+#   test. statsmodels agrees with that hand computation to every printed digit.
+#
+# Index convention, the same one ols_inference_fit uses: 0..n_features-1 are
+# the slopes in feature order and index n_features is the intercept.
+
+import "lib/scikit/scikit.flow"
+
+function t458_close(got: f32, want: f32, tol: f32) -> bool {
+    let mut d: f32 = got - want
+    if d < 0.0 { d = 0.0 - d }
+    return d <= tol
+}
+
+function t458_flag(label: string, v: bool) -> void {
+    print("  ")
+    print(label)
+    if v {
+        println(" yes")
+    } else {
+        println(" no")
+    }
+}
+
+function t458_report(label: string, got: f32, want: f32) -> void {
+    print("  ")
+    print(label)
+    print(" got ")
+    printf("%.6f", got)
+    print(" want ")
+    printf("%.6f", want)
+    println("")
+}
+
+# Fixture A, shared by the reference test and by the penalty test below.
+function t458_fill_fixture_a(X: Matrix, y: ptr<f32>) -> void {
+
+    matrix_set(X, 0, 0, 0.5)
+    matrix_set(X, 0, 1, 1.0)
+    matrix_set(X, 1, 0, 1.0)
+    matrix_set(X, 1, 1, 0.5)
+    matrix_set(X, 2, 0, 1.5)
+    matrix_set(X, 2, 1, 2.0)
+    matrix_set(X, 3, 0, 2.0)
+    matrix_set(X, 3, 1, 1.0)
+    matrix_set(X, 4, 0, 2.5)
+    matrix_set(X, 4, 1, 3.0)
+    matrix_set(X, 5, 0, 3.0)
+    matrix_set(X, 5, 1, 2.0)
+    matrix_set(X, 6, 0, 3.5)
+    matrix_set(X, 6, 1, 1.5)
+    matrix_set(X, 7, 0, 4.0)
+    matrix_set(X, 7, 1, 3.5)
+    matrix_set(X, 8, 0, 4.5)
+    matrix_set(X, 8, 1, 2.5)
+    matrix_set(X, 9, 0, 5.0)
+    matrix_set(X, 9, 1, 4.0)
+    matrix_set(X, 10, 0, 0.75)
+    matrix_set(X, 10, 1, 2.5)
+    matrix_set(X, 11, 0, 1.25)
+    matrix_set(X, 11, 1, 3.0)
+    matrix_set(X, 12, 0, 1.75)
+    matrix_set(X, 12, 1, 0.5)
+    matrix_set(X, 13, 0, 2.25)
+    matrix_set(X, 13, 1, 4.0)
+    matrix_set(X, 14, 0, 2.75)
+    matrix_set(X, 14, 1, 1.25)
+    matrix_set(X, 15, 0, 3.25)
+    matrix_set(X, 15, 1, 4.5)
+    matrix_set(X, 16, 0, 3.75)
+    matrix_set(X, 16, 1, 0.75)
+    matrix_set(X, 17, 0, 4.25)
+    matrix_set(X, 17, 1, 4.75)
+    matrix_set(X, 18, 0, 4.75)
+    matrix_set(X, 18, 1, 1.75)
+    matrix_set(X, 19, 0, 5.5)
+    matrix_set(X, 19, 1, 3.25)
+    matrix_set(X, 20, 0, 2.0)
+    matrix_set(X, 20, 1, 2.0)
+    matrix_set(X, 21, 0, 3.0)
+    matrix_set(X, 21, 1, 3.0)
+    matrix_set(X, 22, 0, 1.0)
+    matrix_set(X, 22, 1, 4.0)
+    matrix_set(X, 23, 0, 4.0)
+    matrix_set(X, 23, 1, 1.0)
+
+    y[0] = 0.0
+    y[1] = 1.0
+    y[2] = 0.0
+    y[3] = 0.0
+    y[4] = 1.0
+    y[5] = 0.0
+    y[6] = 1.0
+    y[7] = 1.0
+    y[8] = 0.0
+    y[9] = 1.0
+    y[10] = 0.0
+    y[11] = 1.0
+    y[12] = 0.0
+    y[13] = 0.0
+    y[14] = 1.0
+    y[15] = 1.0
+    y[16] = 0.0
+    y[17] = 1.0
+    y[18] = 0.0
+    y[19] = 1.0
+    y[20] = 0.0
+    y[21] = 1.0
+    y[22] = 1.0
+    y[23] = 0.0
+}
+
+# Fixture A: 24 samples, 2 features, labels deliberately overlapping so the
+# maximum likelihood estimate exists and the information matrix is well
+# conditioned (condition number 140, fitted probabilities between 0.134 and
+# 0.909).
+#
+# statsmodels reports
+#   coefficients   0.0992060632   0.9025423465  -2.4127187328
+#   std errors     0.3402656264   0.4265496942   1.3457563286
+#   z              0.2915547606   2.1159137113  -1.7928347663
+#   cov  [[ 0.1157806965 -0.0237681631 -0.2695920839]
+#         [-0.0237681631  0.1819446416 -0.3572172811]
+#         [-0.2695920839 -0.3572172811  1.8110600959]]
+function test_binary_inference_matches_statsmodels() -> i32 {
+    println("Test: binary logistic inference against statsmodels")
+    let X: Matrix = matrix_new(24, 2)
+    let y: ptr<f32> = array_new_f32(24)
+    t458_fill_fixture_a(X, y)
+
+    # Solver 1 is Newton / IRLS (issue #357). The literal is used rather than
+    # LOGISTIC_SOLVER_NEWTON because an export const is not transitively
+    # visible through the scikit.flow umbrella.
+    let model: LogisticRegression = logistic_regression_fit_solver(X, y, 2, 100, 0.1, penalty_none(), 1)
+    let inf: LogisticInference = logistic_inference_fit(model, X, y)
+
+    let mut fails: i32 = 0
+
+    if inf.fitted == false {
+        println("  FAIL: inference should be available for an unpenalized binary fit")
+        fails = fails + 1
+    }
+    if inf.rank_deficient {
+        println("  FAIL: a well conditioned design was flagged rank deficient")
+        fails = fails + 1
+    }
+    if inf.separable {
+        println("  FAIL: an overlapping fixture was flagged separable")
+        fails = fails + 1
+    }
+    if inf.n_params != 3 {
+        print("  FAIL: n_params should be 3, got ")
+        print(inf.n_params)
+        println("")
+        fails = fails + 1
+    }
+
+    # The covariance is evaluated at the fitted coefficients, so the fit has to
+    # be on the optimum before the standard errors can match the reference.
+    let c0: f32 = model.weights[2]
+    let c1: f32 = model.weights[3]
+    let cb: f32 = model.biases[1]
+    t458_report("coef(x0)    ", c0, 0.0992061)
+    if t458_close(c0, 0.0992061, 0.0005) == false { fails = fails + 1 }
+    t458_report("coef(x1)    ", c1, 0.9025423)
+    if t458_close(c1, 0.9025423, 0.0005) == false { fails = fails + 1 }
+    t458_report("intercept   ", cb, -2.4127187)
+    if t458_close(cb, -2.4127187, 0.0015) == false { fails = fails + 1 }
+
+    let se0: f32 = logistic_inference_std_error(inf, 0)
+    t458_report("se(x0)      ", se0, 0.3402656)
+    if t458_close(se0, 0.3402656, 0.0004) == false { fails = fails + 1 }
+
+    let se1: f32 = logistic_inference_std_error(inf, 1)
+    t458_report("se(x1)      ", se1, 0.4265497)
+    if t458_close(se1, 0.4265497, 0.0005) == false { fails = fails + 1 }
+
+    let se2: f32 = logistic_inference_std_error(inf, 2)
+    t458_report("se(intercpt)", se2, 1.3457563)
+    if t458_close(se2, 1.3457563, 0.0015) == false { fails = fails + 1 }
+
+    let z0: f32 = logistic_inference_z_stat(inf, 0)
+    t458_report("z(x0)       ", z0, 0.2915548)
+    if t458_close(z0, 0.2915548, 0.002) == false { fails = fails + 1 }
+
+    let z1: f32 = logistic_inference_z_stat(inf, 1)
+    t458_report("z(x1)       ", z1, 2.1159137)
+    if t458_close(z1, 2.1159137, 0.003) == false { fails = fails + 1 }
+
+    let z2: f32 = logistic_inference_z_stat(inf, 2)
+    t458_report("z(intercpt) ", z2, -1.7928348)
+    if t458_close(z2, -1.7928348, 0.003) == false { fails = fails + 1 }
+
+    let v00: f32 = logistic_inference_cov(inf, 0, 0)
+    t458_report("cov(0,0)    ", v00, 0.1157807)
+    if t458_close(v00, 0.1157807, 0.0003) == false { fails = fails + 1 }
+
+    let v01: f32 = logistic_inference_cov(inf, 0, 1)
+    t458_report("cov(0,1)    ", v01, -0.0237682)
+    if t458_close(v01, -0.0237682, 0.0003) == false { fails = fails + 1 }
+
+    let v02: f32 = logistic_inference_cov(inf, 0, 2)
+    t458_report("cov(0,2)    ", v02, -0.2695921)
+    if t458_close(v02, -0.2695921, 0.001) == false { fails = fails + 1 }
+
+    let v11: f32 = logistic_inference_cov(inf, 1, 1)
+    t458_report("cov(1,1)    ", v11, 0.1819446)
+    if t458_close(v11, 0.1819446, 0.0005) == false { fails = fails + 1 }
+
+    let v12: f32 = logistic_inference_cov(inf, 1, 2)
+    t458_report("cov(1,2)    ", v12, -0.3572173)
+    if t458_close(v12, -0.3572173, 0.0012) == false { fails = fails + 1 }
+
+    let v22: f32 = logistic_inference_cov(inf, 2, 2)
+    t458_report("cov(2,2)    ", v22, 1.8110601)
+    if t458_close(v22, 1.8110601, 0.004) == false { fails = fails + 1 }
+
+    let v10: f32 = logistic_inference_cov(inf, 1, 0)
+    if t458_close(v10, v01, 0.0000001) == false {
+        println("  FAIL: covariance is not symmetric")
+        fails = fails + 1
+    }
+
+    logistic_inference_free(inf)
+    logistic_regression_free(model)
+    array_free_f32(y)
+    matrix_free(X)
+
+    if fails == 0 {
+        println("  OK: standard errors, z statistics and covariance match statsmodels")
+        return 0
+    }
+    return 1
+}
+
+# Fixture B: one binary feature, the saturated 2x2 table. Eight rows at x = 0
+# of which two are positive, eight rows at x = 1 of which six are positive. The
+# first row carries the negative label on purpose: the class whose coefficients
+# are reported is classes[1], the second distinct label in order of appearance
+# in y, which is the convention logistic_decision_function already uses.
+#
+# The maximum likelihood estimate is the table itself:
+#   p0 = 0.25, p1 = 0.75
+#   intercept = log(p0 / (1 - p0)) = log(1/3)  = -1.0986123
+#   slope     = log(p1 / (1 - p1)) - intercept =  2.1972246
+#
+# The information matrix is block diagonal in the two cell counts, so the
+# covariance is closed form with n0 p0 q0 = n1 p1 q1 = 8 * 0.1875 = 1.5:
+#   var(intercept) = 1 / 1.5                 = 0.6666667  se = 0.8164966
+#   var(slope)     = 1 / 1.5 + 1 / 1.5       = 1.3333333  se = 1.1547005
+#   cov(slope, intercept) = -1 / 1.5         = -0.6666667
+#   z(slope) = 1.9028523     z(intercept) = -1.3455198
+#
+# statsmodels reproduces every one of those digits, so this test is a hand
+# computation that a library happens to agree with rather than the other way
+# round.
+function test_saturated_table_closed_form() -> i32 {
+    println("Test: saturated 2x2 table against the closed form")
+    let X: Matrix = matrix_new(16, 1)
+    let y: ptr<f32> = array_new_f32(16)
+    for i in 0 to 8 {
+        matrix_set(X, i, 0, 0.0)
+        matrix_set(X, i + 8, 0, 1.0)
+        y[i] = 0.0
+        y[i + 8] = 1.0
+    }
+    y[2] = 1.0
+    y[3] = 1.0
+    y[14] = 0.0
+    y[15] = 0.0
+
+    let model: LogisticRegression = logistic_regression_fit_solver(X, y, 2, 100, 0.1, penalty_none(), 1)
+    let inf: LogisticInference = logistic_inference_fit(model, X, y)
+
+    let mut fails: i32 = 0
+    if inf.fitted == false {
+        println("  FAIL: inference should be available for this fit")
+        fails = fails + 1
+        logistic_inference_free(inf)
+        logistic_regression_free(model)
+        array_free_f32(y)
+        matrix_free(X)
+        return 1
+    }
+
+    let slope: f32 = model.weights[1]
+    let intercept: f32 = model.biases[1]
+    t458_report("slope       ", slope, 2.1972246)
+    if t458_close(slope, 2.1972246, 0.002) == false { fails = fails + 1 }
+    t458_report("intercept   ", intercept, -1.0986123)
+    if t458_close(intercept, -1.0986123, 0.002) == false { fails = fails + 1 }
+
+    let se_slope: f32 = logistic_inference_std_error(inf, 0)
+    t458_report("se(slope)   ", se_slope, 1.1547005)
+    if t458_close(se_slope, 1.1547005, 0.0012) == false { fails = fails + 1 }
+
+    let se_int: f32 = logistic_inference_std_error(inf, 1)
+    t458_report("se(intercpt)", se_int, 0.8164966)
+    if t458_close(se_int, 0.8164966, 0.0009) == false { fails = fails + 1 }
+
+    let cov01: f32 = logistic_inference_cov(inf, 0, 1)
+    t458_report("cov(0,1)    ", cov01, -0.6666667)
+    if t458_close(cov01, -0.6666667, 0.0015) == false { fails = fails + 1 }
+
+    let zs: f32 = logistic_inference_z_stat(inf, 0)
+    t458_report("z(slope)    ", zs, 1.9028523)
+    if t458_close(zs, 1.9028523, 0.003) == false { fails = fails + 1 }
+
+    let zi: f32 = logistic_inference_z_stat(inf, 1)
+    t458_report("z(intercpt) ", zi, -1.3455198)
+    if t458_close(zi, -1.3455198, 0.003) == false { fails = fails + 1 }
+
+    logistic_inference_free(inf)
+    logistic_regression_free(model)
+    array_free_f32(y)
+    matrix_free(X)
+
+    if fails == 0 {
+        println("  OK: closed form covariance for the saturated table reproduced")
+        return 0
+    }
+    return 1
+}
+
+# A perfectly separable fixture has no maximum likelihood estimate. The
+# coefficients diverge, every fitted probability saturates and the information
+# matrix goes to zero, so a standard error read off it says only where the
+# optimizer stopped. The result must be refused and flagged, not returned as a
+# large finite number.
+function test_separable_fixture_is_refused() -> i32 {
+    println("Test: separable fixture is refused")
+    let X: Matrix = matrix_new(6, 1)
+    let y: ptr<f32> = array_new_f32(6)
+    matrix_set(X, 0, 0, -3.0)
+    matrix_set(X, 1, 0, -2.0)
+    matrix_set(X, 2, 0, -1.0)
+    matrix_set(X, 3, 0, 1.0)
+    matrix_set(X, 4, 0, 2.0)
+    matrix_set(X, 5, 0, 3.0)
+    y[0] = 0.0
+    y[1] = 0.0
+    y[2] = 0.0
+    y[3] = 1.0
+    y[4] = 1.0
+    y[5] = 1.0
+
+    # Both solvers, because they stop in very different places on separable
+    # data: Newton runs the slope out to 13.4, LBFGS stops at 8.3. Separation
+    # is a property of the data, so both have to be refused.
+    let newton: LogisticRegression = logistic_regression_fit_solver(X, y, 2, 100, 0.1, penalty_none(), 1)
+    let lbfgs: LogisticRegression = logistic_regression_fit(X, y, 2, 100, 0.1, penalty_none())
+    let inf_n: LogisticInference = logistic_inference_fit(newton, X, y)
+    let inf_l: LogisticInference = logistic_inference_fit(lbfgs, X, y)
+
+    let mut fails: i32 = 0
+    print("  newton slope ")
+    printf("%.4f", newton.weights[1])
+    print("   lbfgs slope ")
+    printf("%.4f", lbfgs.weights[1])
+    println("")
+    t458_flag("newton fitted    ", inf_n.fitted)
+    t458_flag("newton separable ", inf_n.separable)
+    t458_flag("lbfgs  fitted    ", inf_l.fitted)
+    t458_flag("lbfgs  separable ", inf_l.separable)
+
+    if inf_n.fitted || inf_l.fitted {
+        println("  FAIL: a separated fit must not report standard errors")
+        fails = fails + 1
+    }
+    if inf_n.separable == false || inf_l.separable == false {
+        println("  FAIL: separation was not detected")
+        fails = fails + 1
+    }
+    if logistic_inference_std_error(inf_n, 0) != 0.0 {
+        println("  FAIL: refused inference must return a zero standard error")
+        fails = fails + 1
+    }
+    if logistic_inference_std_error(inf_l, 0) != 0.0 {
+        println("  FAIL: refused inference must return a zero standard error")
+        fails = fails + 1
+    }
+
+    logistic_inference_free(inf_n)
+    logistic_inference_free(inf_l)
+    logistic_regression_free(newton)
+    logistic_regression_free(lbfgs)
+    array_free_f32(y)
+    matrix_free(X)
+
+    if fails == 0 {
+        println("  OK: separation refused rather than reported")
+        return 0
+    }
+    return 1
+}
+
+# A duplicated column makes the weighted design rank deficient. The refusal is
+# the same one ols_inference_fit gives, on the same relative pivot cutoff.
+function test_collinear_design_is_refused() -> i32 {
+    println("Test: collinear design is refused")
+    let X: Matrix = matrix_new(10, 2)
+    let y: ptr<f32> = array_new_f32(10)
+    for i in 0 to 10 {
+        let v: f32 = (i as f32) * 0.5 + 1.0
+        matrix_set(X, i, 0, v)
+        matrix_set(X, i, 1, v)
+        y[i] = 0.0
+    }
+    y[1] = 1.0
+    y[3] = 1.0
+    y[4] = 1.0
+    y[6] = 1.0
+    y[8] = 1.0
+    y[9] = 1.0
+
+    let model: LogisticRegression = logistic_regression_fit_solver(X, y, 2, 100, 0.1, penalty_none(), 1)
+    let inf: LogisticInference = logistic_inference_fit(model, X, y)
+
+    let mut fails: i32 = 0
+    t458_flag("fitted     ", inf.fitted)
+    print("  rank        ")
+    print(inf.rank)
+    print(" of ")
+    print(inf.n_params)
+    println("")
+
+    if inf.fitted {
+        println("  FAIL: a singular information matrix must not report standard errors")
+        fails = fails + 1
+    }
+    if inf.rank_deficient == false {
+        println("  FAIL: rank deficiency was not flagged")
+        fails = fails + 1
+    }
+
+    logistic_inference_free(inf)
+    logistic_regression_free(model)
+    array_free_f32(y)
+    matrix_free(X)
+
+    if fails == 0 {
+        println("  OK: singular information matrix refused")
+        return 0
+    }
+    return 1
+}
+
+# The Newton Hessian of a penalized fit carries the alpha * I term, so its
+# inverse is a penalized covariance and not a sampling covariance. Refused, as
+# ols_inference_fit refuses the same case.
+function test_penalized_fit_is_refused() -> i32 {
+    println("Test: penalized fit is refused")
+    # Fixture A again, so the only thing that changes between the refusal and
+    # the answer below is the penalty.
+    let X: Matrix = matrix_new(24, 2)
+    let y: ptr<f32> = array_new_f32(24)
+    t458_fill_fixture_a(X, y)
+
+    let model: LogisticRegression = logistic_regression_fit_solver(X, y, 2, 100, 0.1, penalty_l2(0.5), 1)
+    let inf: LogisticInference = logistic_inference_fit(model, X, y)
+
+    let mut fails: i32 = 0
+    if inf.fitted {
+        println("  FAIL: a penalized fit must be refused")
+        fails = fails + 1
+    }
+    if logistic_inference_cov(inf, 0, 0) != 0.0 {
+        println("  FAIL: refused inference must return a zero covariance")
+        fails = fails + 1
+    }
+
+    logistic_inference_free(inf)
+    logistic_regression_free(model)
+
+    # The same design fitted unpenalized is answered, which is what makes the
+    # refusal above about the penalty and not about the data.
+    let free_model: LogisticRegression = logistic_regression_fit_solver(X, y, 2, 100, 0.1, penalty_none(), 1)
+    let free_inf: LogisticInference = logistic_inference_fit(free_model, X, y)
+    if free_inf.fitted == false {
+        println("  FAIL: the same design fitted unpenalized should be answered")
+        fails = fails + 1
+    }
+    if logistic_inference_std_error(free_inf, 0) <= 0.0 {
+        println("  FAIL: unpenalized fit should report a positive standard error")
+        fails = fails + 1
+    }
+    logistic_inference_free(free_inf)
+    logistic_regression_free(free_model)
+
+    array_free_f32(y)
+    matrix_free(X)
+
+    if fails == 0 {
+        println("  OK: penalized refused, unpenalized answered on the same design")
+        return 0
+    }
+    return 1
+}
+
+# Above two classes the fit is the joint multinomial softmax of #433, whose
+# information matrix is a block system this surface does not build. Refused.
+function test_multiclass_fit_is_refused() -> i32 {
+    println("Test: multiclass fit is refused")
+    let X: Matrix = matrix_new(12, 2)
+    let y: ptr<f32> = array_new_f32(12)
+    for i in 0 to 12 {
+        let a: f32 = (i as f32) * 0.5
+        matrix_set(X, i, 0, a)
+        matrix_set(X, i, 1, 3.0 - a * 0.2)
+        y[i] = (i % 3) as f32
+    }
+
+    let model: LogisticRegression = logistic_regression_fit_solver(X, y, 3, 60, 0.1, penalty_none(), 1)
+    let inf: LogisticInference = logistic_inference_fit(model, X, y)
+
+    let mut fails: i32 = 0
+    print("  n_classes   ")
+    print(model.n_classes)
+    println("")
+    t458_flag("fitted     ", inf.fitted)
+
+    if inf.fitted {
+        println("  FAIL: a multiclass fit must be refused")
+        fails = fails + 1
+    }
+    if logistic_inference_z_stat(inf, 0) != 0.0 {
+        println("  FAIL: refused inference must return a zero z statistic")
+        fails = fails + 1
+    }
+
+    logistic_inference_free(inf)
+    logistic_regression_free(model)
+    array_free_f32(y)
+    matrix_free(X)
+
+    if fails == 0 {
+        println("  OK: multiclass refused")
+        return 0
+    }
+    return 1
+}
+
+# The default LBFGS solver optimizes the same objective, so it lands on the
+# same optimum and the inference has to agree with the Newton fit to the
+# accuracy of the two solvers. This is what pins the claim that the inference
+# is a property of the fitted coefficients and not of the optimizer.
+function test_lbfgs_and_newton_agree() -> i32 {
+    println("Test: LBFGS and Newton fits give the same standard errors")
+    let X: Matrix = matrix_new(16, 1)
+    let y: ptr<f32> = array_new_f32(16)
+    for i in 0 to 8 {
+        matrix_set(X, i, 0, 0.0)
+        matrix_set(X, i + 8, 0, 1.0)
+        y[i] = 0.0
+        y[i + 8] = 1.0
+    }
+    y[2] = 1.0
+    y[3] = 1.0
+    y[14] = 0.0
+    y[15] = 0.0
+
+    let newton: LogisticRegression = logistic_regression_fit_solver(X, y, 2, 100, 0.1, penalty_none(), 1)
+    let lbfgs: LogisticRegression = logistic_regression_fit(X, y, 2, 100, 0.1, penalty_none())
+    let inf_n: LogisticInference = logistic_inference_fit(newton, X, y)
+    let inf_l: LogisticInference = logistic_inference_fit(lbfgs, X, y)
+
+    let mut fails: i32 = 0
+    if inf_n.fitted == false || inf_l.fitted == false {
+        println("  FAIL: both fits should be answered")
+        fails = fails + 1
+    }
+
+    let sn: f32 = logistic_inference_std_error(inf_n, 0)
+    let sl: f32 = logistic_inference_std_error(inf_l, 0)
+    t458_report("se newton   ", sn, 1.1547005)
+    t458_report("se lbfgs    ", sl, 1.1547005)
+    if t458_close(sn, sl, 0.01) == false {
+        println("  FAIL: the two solvers disagree on the standard error")
+        fails = fails + 1
+    }
+
+    logistic_inference_free(inf_n)
+    logistic_inference_free(inf_l)
+    logistic_regression_free(newton)
+    logistic_regression_free(lbfgs)
+    array_free_f32(y)
+    matrix_free(X)
+
+    if fails == 0 {
+        println("  OK: both solvers give the same inference")
+        return 0
+    }
+    return 1
+}
+
+function main() -> i32 {
+    println("Running binary logistic inference tests...")
+    println("=========================================")
+    println("")
+
+    let mut failures: i32 = 0
+
+    if test_binary_inference_matches_statsmodels() != 0 { failures = failures + 1 }
+    println("")
+    if test_saturated_table_closed_form() != 0 { failures = failures + 1 }
+    println("")
+    if test_separable_fixture_is_refused() != 0 { failures = failures + 1 }
+    println("")
+    if test_collinear_design_is_refused() != 0 { failures = failures + 1 }
+    println("")
+    if test_penalized_fit_is_refused() != 0 { failures = failures + 1 }
+    println("")
+    if test_multiclass_fit_is_refused() != 0 { failures = failures + 1 }
+    println("")
+    if test_lbfgs_and_newton_agree() != 0 { failures = failures + 1 }
+
+    println("")
+    if failures == 0 {
+        println("All logistic inference tests passed!")
+        return 0
+    }
+    print("FAILED: ")
+    print(failures)
+    println(" tests")
+    return 1
+}