Julia stability vs Rust for Scientific Computing

What annoys me a bit about @yurivish blogpost is not that its wrong or anything, but the framing that it’s somehow unique to Julia, and that people keep bringing it up as a “well, its proven that Julia is a mess”. Of course if they were too annoying for @yurivish to do their everyday work, its a very fair point, and if they found something that works better, great!

But there are millions of happy Julia users and in general it’s pretty normal to find bugs in software if you look for it.
Secondly, if you compare a 20 year old industry standard with millions in funding to a 3 year old library written by some Phd, the latter one is bound to fail on many metrics.
To be fair, that might be the point of the blogpost, but it ignores the dimension of how things will be in the future and that it’s not necessarily an architectural problem and also not across the whole ecosystem.

I think there’s also a mindset split, some people just like to have things more strict and avoid bugs by having their compiler proof everything, and others like more freedom and are fine with occasional mishaps. I’m pretty deep in the latter group while wanting performance, therefore Julia is pretty perfect for me (btw, I haven’t had a single Julia bug in the recent 6 years or more, and most of the dependencies I use are really stable).
If you’re in the former group, I don’t think Julia is as good (compared to e.g. Rust), although I still think its much better in that regard than Python architecturally.
Regarding Rust, I’ve just recently read a pretty interesting blogpost (sadly cant find it right now) on how the Rust compiler only helps you avoid ~5-20% of classic correctness and security bugs, and the rest still needs having lots of testing and hardening - and part of the current image of security for Rust has mainly been because of security by obscurity, and now that its used more, people start finding pretty big bugs.

Just for the fun of it, I put claude on Python, and it also found some eye watering correctness issues (to be fair, I haven’t taken the time to verify and judge them, but it seems like that’s a similar situation for the Julia version):

Library Reproducer Got Expected Why it matters
stdlib random.choices(['a','b','c'], weights=[-1,5,1], k=10000) Counter({'b': ~8000, 'c': ~2000}) error or proportional sampling Negative weight on 'a' silently shifts
mass to the next bucket. Validation only catches the case where the total ≤ 0.
stdlib random.choices(['a','b','c'], cum_weights=[5,2,7], k=10000) Counter({'a': ~2800, 'c': ~7200}) error on non-monotone Non-monotone cum_weights makes 'b'
unselectable. No validation.
stdlib statistics.fmean([1,2,3], weights=[-1,1,1]) 4.0 error or value in [1,3] “Mean” of three values in [1,3] returns 4 — outside the convex hull.
stdlib json.dumps({1: 'a', '1': 'b'}) '{"1": "a", "1": "b"}' error or merge Produces invalid JSON with duplicate keys; round-trip silently drops one entry. Same for
{True:'a','true':'b'}.
stdlib json.dumps({1: 'a', '1': 'b'}, sort_keys=True) TypeError: '<' not supported between str and int succeed Internal sort runs before the documented int→str key
coercion.
stdlib urlparse('http://example.com/?').geturl() 'http://example.com/' unchanged Trailing ? (empty query) and # (empty fragment) silently stripped — breaks
signing/canonicalization round-trips.
stdlib a=dt(2026,11,1,1,30,tz=NY,fold=0); b=a.replace(fold=1); len({a,b}) 1 2 Two different absolute instants compare equal and dedupe in set/dict even though
a.timestamp() != b.timestamp().
stdlib dt(2026,3,8,1,30,tz=NY) + timedelta(hours=1) 2026-03-08 02:30 EST (a wall time the timezone says doesn’t exist) 03:30 EDT timedelta arithmetic on aware
datetimes is wall-clock, not absolute.
numpy x = np.array([1,2,3], dtype=np.int8); np.where(x>1, x, 1000) array([-24, 2, 3], dtype=int8) error or upcast 1000 silently truncates against the array dtype.
x[0]=1000, x.fill(1000), np.full(3,1000,np.int8) all raise — np.where is the silent outlier.
numpy np.array([np.nan, np.inf, 1e20]).astype(np.int64) array([INT64_MIN, INT64_MIN, INT64_MIN]) (RuntimeWarning only) ValueError np.array([nan], dtype=np.int64)
raises; astype produces garbage. Same call, two contracts.
numpy np.histogram([1, 2, 3, np.nan, 5]) ValueError: autodetected range of [nan, nan] is not finite histogram of the four finite values One NaN among four finite values
poisons auto-range.
numpy np.unique([1, np.nan, np.nan, 2]) vs np.unique([[1,np.nan],[1,np.nan]], axis=0) first collapses NaNs → [1,2,nan]; second does NOT collapse same semantics flat vs
axis Same operation, two answers.
numpy np.intersect1d([1, np.nan], [2, np.nan]) vs np.union1d([1, np.nan], [2, np.nan]) [] vs [1, 2, nan] consistent NaN semantics Union collapses NaNs, intersection
treats them as unequal.
numpy np.isin([1, np.nan, 2], [np.nan]) [False, False, False] matches unique’s view A third NaN semantic in the same family of set ops.
numpy x = np.random.randn(63); np.fft.irfft(np.fft.rfft(x)) length 62 array, max err ~3.1 against x round-trip identity rfft discards parity; irfft defaults
to even length. No warning when input was odd; silent data corruption.
numpy np.average([1, 2, 3], weights=[-1, 1, 1]) 4.0 error or value in [1, 3] Same negative-weight pattern as Python’s random.choices / statistics.fmean.
numpy np.argmin([1.0, np.nan, 0.5, 2.0]); np.argmax([1.0, np.nan, 0.5, 2.0]) 1 and 1 2 and 3 The NaN slot wins both argmin AND argmax.
numpy np.median([1, np.nan, 2]) nan 1.5 (or hard error) Silent NaN propagation; need np.nanmedian. statistics.median raises instead.
numpy np.maximum(1, np.nan) vs np.fmax(1, np.nan) nan vs 1.0 one consistent semantic Two functions, same name shape, opposite NaN behavior.
numpy np.searchsorted([1, 2, 3, 4, np.nan], np.nan) 4 (i.e. before the existing NaN) well-defined Insertion point for NaN lands inside the NaN run.
numpy np.array([1, 'two', 3.0]).dtype dtype('<U32') error or object Silent stringification; arr[0] == 1 is now False.
numpy np.nansum([np.nan, np.nan]) 0.0 nan or empty-input error “Sum of these values” is 0 when every value is NaN.
numpy np.outer(np.eye(2), np.eye(2)).shape (4, 4) (2, 2, 2, 2) Kronecker-like np.outer silently flattens 2-D inputs to 1-D.
numpy A=np.zeros((2,3,4)); B=np.zeros((2,4,5)); np.dot(A,B).shape (2, 3, 2, 5) (2, 3, 5) (batched matmul) np.dot does tensor contraction, not batched matmul. A @ B
does the latter; they diverge at 3-D.