Julia stability vs Rust for Scientific Computing

I asked a coding agent to find correctness bugs in the latest stable Julia release, then let it run for a bit. It found many examples in the standard library where a function, called with a realistic input value, returns the wrong result.

I possess neither Mythos access nor a strong math background, so this was done with Opus 4.7 and some of these reports may be noise. But it’s an experiment I would encourage anyone curious to conduct for themselves.

Reproducer Got Expected Why it matters
quantile!(reshape([1.,2,3,4,5],5,1), [.25,.5,.75]; alpha=0, beta=1) [1.5, 3.0, 4.5] [1.25, 2.5, 3.75] The beta keyword is silently ignored for matrix inputs; you get a different quantile definition.
foldr(*, Iterators.flatten([["a","b"],["c","d"]])) "badc" "abcd" Right-associative reduce over a flattened sequence doesn’t give the same result as over a pre-flattened array.
invmod(2, UInt(3)) 0x00 0x2 Switching one operand from Int to UInt flips the answer to 0, never a valid modular inverse.
rval = [0]; findmin!(rval, [0], [1]; init=false); rval [1] [0] The documented kwarg that says “keep my running result” has no effect; breaks streaming argmin. False positive?
findmin(sparsevec([1], [NaN], 2)) (0.0, 2) (NaN, 1) The universal “missing/invalid” sentinel disappears from sparse findmin/findmax. SparseArrays#714
R = givens(1.,1.,1,2)[1] * givens(1.,1.,2,3)[1]; R*(R'*[1.,0,0]) [0.85, -0.15, -0.5] [1.0, 0.0, 0.0] The most basic identity in linear algebra — a rotation times its inverse — fails.
A = UnitUpperTriangular([1 0; 0 1]); rmul!(A, Diagonal([2,3])); det(A) 1 6 An in-place scale silently does nothing through a UnitTriangular wrapper; det doesn’t change. LinearAlgebra#1600
dot([1.,0.], Symmetric(sparse([1],[2],[1.],2,2)), [0.,1.]) 0.0 1.0 Whenever x or y has a zero at the column of an off-diagonal entry, the contribution from that entry is silently dropped even though it’s mathematically nonzero. SparseArrays#713 SparseArrays#715
pinv(Diagonal([1.0, 1e-100])).diag [1.0, 1e100] [1.0, 0.0] (Not a bug, missing feature; see stevengj below) Pseudoinverse of near-singular data isn’t actually regularized: tiny entries become huge inverses (compare to the same on a dense matrix).
d = [1.,2.]; Diagonal(d) .= Diagonal(view(d, 2:-1:1)); d [2., 2.] [2., 1.] (Not a bug; see stevengj below) A .= B silently corrupts when B is a view of A through a Diagonal/etc. wrapper.

The pinv issue is known since 2022.

Edit: Here’s a link to a zip containing a larger set of LLM-generated bug reports as individual Markdown files.