1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
// Property-based tests for quality improvements (Wave V).
//
// V2: Cubic solver — all roots satisfy the polynomial.
// V3: Matrix inverse — A * A⁻¹ ≈ I for random integer matrices.
// V4: Codegen — generated code is syntactically valid Rust.
// Bonus: Determinant consistency — cofactor and LU paths agree.
mod common;
use proptest::prelude::*;
use symplex::prelude::*;
// ═══════════════════════════════════════════════════════════════════════════
// Strategies
// ═══════════════════════════════════════════════════════════════════════════
/// Generate random coefficients for a cubic: ax³ + bx² + cx + d
/// with a ∈ [1,4] so the leading coefficient is always positive & nonzero.
fn arb_cubic_coeffs() -> impl Strategy<Value = [i64; 4]> {
(1..5i64, -5..5i64, -5..5i64, -5..5i64).prop_map(|(a, b, c, d)| [a, b, c, d])
}
// ═══════════════════════════════════════════════════════════════════════════
// Properties
// ═══════════════════════════════════════════════════════════════════════════
proptest! {
#![proptest_config(ProptestConfig::with_cases(30))]
// ── V2: Cubic solver ───────────────────────────────────────────────
/// Every root returned by `solve_or_empty` should satisfy the
/// original cubic polynomial (residual < 1e-6).
#[test]
fn cubic_roots_satisfy_polynomial(coeffs in arb_cubic_coeffs()) {
let ctx = Context::new();
let [a, b, c, d] = coeffs;
let x = ctx.symbol("x");
// Build a·x³ + b·x² + c·x + d
let poly = &(&x.powi(3) * a) + &(&x.powi(2) * b) + &(&x * c) + d;
let roots = poly.solve_or_empty(&x);
let mut bail = common::BailCounter::new("cubic_roots_satisfy_polynomial");
for root in &roots {
let val = poly.subs(&x, root).eval().simplify();
if let Ok(v) = val.eval_f64() {
bail.check();
prop_assert!(
v.abs() < 1e-6,
"root {} doesn't satisfy {}x³+{}x²+{}x+{}: residual={}",
root, a, b, c, d, v
);
} else {
bail.skip();
}
}
if !roots.is_empty() {
// Symbolic-only roots legitimately fail evalf_f64 for some inputs;
// allow high skip rate — proptest coverage across 30 cases ensures
// non-vacuousness at the suite level.
bail.assert_skip_rate_below(1.0);
}
}
// ── V3: Matrix inverse ─────────────────────────────────────────────
/// For a random 3×3 integer matrix, if it is invertible then
/// A · A⁻¹ ≈ I (diagonal entries ≈ 1, off-diagonal ≈ 0).
#[test]
fn matrix_inverse_is_identity(
entries in proptest::array::uniform9(-3i64..4i64)
) {
let ctx = Context::new();
let data: Vec<Vec<Ex>> = entries
.chunks(3)
.map(|row| row.iter().map(|&v| ctx.int(v)).collect())
.collect();
let m = symplex::matrix::Matrix::new(data).unwrap();
if let Ok(inv) = m.inv() {
let product = m.matmul(&inv).unwrap();
let mut bail = common::BailCounter::new("matrix_inverse_is_identity");
// Check diagonal ≈ 1, off-diagonal ≈ 0
for i in 0..3 {
for j in 0..3 {
let entry = product.get(i, j).eval().simplify();
if let Ok(v) = entry.eval_f64() {
bail.check();
let v: f64 = v;
if i == j {
prop_assert!(
(v - 1.0).abs() < 1e-8,
"diagonal ({},{}) should be 1, got {}",
i, j, v
);
} else {
prop_assert!(
v.abs() < 1e-8,
"off-diagonal ({},{}) should be 0, got {}",
i, j, v
);
}
} else {
bail.skip();
}
}
}
bail.assert_not_vacuous();
}
// If not invertible (det=0), that's fine — skip.
}
// ── V4: Codegen produces valid Rust syntax ─────────────────────────
/// `to_rust_fn` output should always contain the expected function
/// signature components (pub fn, parameter, return type).
#[test]
fn codegen_produces_valid_syntax(a in -5i64..5, b in -5i64..5, c in 1i64..5) {
let ctx = Context::new();
let x = ctx.symbol("x");
let poly = &(&x.powi(2) * a) + &(&x * b) + c;
if let Ok(code) = poly.to_rust_fn("test_fn", &["x"]) {
prop_assert!(
code.contains("pub fn test_fn"),
"missing 'pub fn test_fn' in:\n{code}"
);
prop_assert!(
code.contains("x: f64"),
"missing parameter 'x: f64' in:\n{code}"
);
prop_assert!(
code.contains("-> f64"),
"missing return type '-> f64' in:\n{code}"
);
}
}
// ── Bonus: Determinant of product ──────────────────────────────────
/// det(A · B) == det(A) · det(B) for random 3×3 integer matrices.
/// This exercises both the determinant and matmul paths.
#[test]
fn det_of_product_equals_product_of_dets(
a_entries in proptest::array::uniform9(-3i64..4i64),
b_entries in proptest::array::uniform9(-3i64..4i64),
) {
let ctx = Context::new();
let mat_a = symplex::matrix::Matrix::new(
a_entries.chunks(3)
.map(|row| row.iter().map(|&v| ctx.int(v)).collect())
.collect(),
).unwrap();
let mat_b = symplex::matrix::Matrix::new(
b_entries.chunks(3)
.map(|row| row.iter().map(|&v| ctx.int(v)).collect())
.collect(),
).unwrap();
let det_a = mat_a.det().unwrap().eval().simplify();
let det_b = mat_b.det().unwrap().eval().simplify();
let product_of_dets = (&det_a * &det_b).eval().simplify();
let ab = mat_a.matmul(&mat_b).unwrap();
let det_ab = ab.det().unwrap().eval().simplify();
let mut bail = common::BailCounter::new("det_of_product");
if let (Ok(lhs), Ok(rhs)) = (det_ab.eval_f64(), product_of_dets.eval_f64()) {
if lhs.is_finite() && rhs.is_finite() {
bail.check();
let tol = 1e-6 * lhs.abs().max(rhs.abs()).max(1.0);
prop_assert!(
(lhs - rhs).abs() < tol,
"det(A·B) = {} but det(A)·det(B) = {} (diff={})",
lhs, rhs, (lhs - rhs).abs()
);
} else {
bail.skip();
}
} else {
bail.skip();
}
bail.assert_not_vacuous();
}
}