[ANN] MathChecker.jl: Find problems in floating-point calculations

I’m always surprised this isn’t a bigger topic of conversation among Julia users: how to find where problems are coming from in numerics. The gap between our pretty equations and algorithms on one hand, and practical code implementations on the other can simultaneously be frightening and invisible. Most importantly, we need to know:

  1. Where do NaNs, Infs, and subnormal numbers come from, or go to?
  2. Where do we lose precision — due to cancellation or absorption, inexact rounding, or simply mixing number types with unintended precision?
  3. Did we use uninitialized values? Maybe we forgot to initialize a recurrence relation?

Because Julia is generic, we can answer all of these questions by creating our own number type to wrap any given float type. This is what MathChecker.jl — and specifically, the MathChecker.Checked type — does. As of this writing, it’s not in the General registry yet, but it’s on its way.

The type is parameterized by the underlying float type, and flags for whether or not to check each of the items in the list above. Using that type is as easy as wrapping input values or arrays in the constructor:

x = Checked(1.0)

Then, just feed these into your functions, and most should work just fine — until a problem with your numerics is found.

By default, only NaN and Inf or mismatched precisions are considered problems. Features can be turned on or off (or adjusted) using keyword arguments, as in

y = Checked(3.14159265; inf=false, cancellation=true, absorption=true)

If problems are found, the default behavior is to raise an error, showing exactly where the problematic code is. But it’s also possible to handle problems using custom behaviors, such as logging, so we can see where all the problems are in a program.

I first got into this when trying to figure out where a NaN was coming from in some massive pile of code. Basically, I was looking for a signaling NaN. Then, I came across this discourse post by @brianguenter, which was basically the inspiration for this approach. This package has significant overlap (but not precise equivalence) with the IEEE 754 floating-point exception flags.

This project has been on my to-do list for a long time, and I’ve frequently found myself slapping together worse versions of this code for quick tests. But now I had some free time, and a lot of help from Claude, and I’ve managed to put something together that might be useful to more people than just myself.

Do you have a link to the code?

:man_facepalming:

That’s what I get for posting last thing before bed. Thanks for pointing that out. I’ve added the link prominently above.

Hey, just to let you know, I made a package that is very similar to this one a few months ago: PrecisionCarriers.jl

It is already registered and written entirely by me and without LLMs. With how closely your concept matches, I do wonder if your LLM trained on it. Also with how many — your post uses, I wonder if the text is written by you or LLM generated as well, which to my knowledge, is not allowed on here.

Your comment suggests you didn’t read beyond the first two sentences about my package. There is very little similarity whatsoever, beyond the fact that they both worry about numerics and use the standard Julia approach of wrapping a standard float. And given the rather insulting claims you raise, I would think you might at least do a little due diligence before making such unfounded accusations. The design is completely different, the capabilities are different, and the behavior is different. Plus, your timeline is shaky.

PrecisionCarriers stores a BigFloat alongside the native value and uses divergence between the two computations as its diagnostic. MathChecker performs a variety of selectable tests on the native floating-point result itself at each operation, without maintaining a higher-precision shadow computation. MathChecker can show you the exact spot in the code where a particular type of problem is happening without changing anything in the code, while PrecisionCarriers just gives you the end result compared to the BigFloat computation.

My post above and the documentation clearly credit Brian Guenter’s post for its origins and the inspiration for its design. That is a design that I’ve been using for a couple years — even before your package existed. I’ve basically just added a few different checks to Brian’s concept. I also prominently note related packages, all of which are more closely related to my package than your package is. So I’m not trying to hide anything about lineage or influences, or even claim that my contribution is particularly original (though I have helped earlier ideas evolve).

Also, everything I’ve written in this thread was written 100% by me, without any LLM contribution. I have a liberal arts education, and I use em dashes. I’ve been using em dashes since the 90s. As you can see from my very prominent notes in the README and docs, I am also very happy to give credit to Claude where credit is due — and it’s not here.

Also, everything I’ve written in this thread was written 100% by me, without any LLM contribution. I have a liberal arts education, and I use em dashes. I’ve been using em dashes since the 90s. As you can see from my very prominent notes in the README and docs, I am also very happy to give credit to Claude where credit is due — and it’s not here.

That is fair and in this case I’m sorry. I wish this whole thing wasn’t a problem but em dashes still throw me off. For what it’s worth, I didn’t flag the post because I wasn’t sure enough, but I wanted to ask.

I saw your related package mentions in your README. They are mostly packages I’m aware of. I thought I mentioned them myself in my README, but it seems I only did that in a poster I made about it for JuliaCon and a presentation I gave about it earlier and I confused that.
The fact that your related package mentions didn’t include my package was what made me think you (and/or Claude) were probably not aware of it, so I replied.

I still think the packages are very related. The way both are used is to wrap initial values in a custom type and then both packages propagate that value through an entire calculation by overloading all of the relevant basic arithmetic functions. What is different is what happens in the overloads, where you emit some sort of warning/log and I don’t do anything until the value falls out the other side for inspection. This explicitly also includes handling of NaN/Inf values (which seem to be your main concern from what I can tell), which I took care get propagated properly.

I’ve had the idea for an “inspect” functionality in my package too, but haven’t had a good idea of how to visualize it and make it usable. I never opened an issue about it though, so you couldn’t know that of course. I just didn’t want to have more or less arbitrary thresholds for what constitutes “too much” error, so I was imagining some sort of flame graph where you can see which part loses how much precision, or something similar.

Thank you for apologizing. I got a little snippy because I take attribution very seriously, and I do not take accusations of failing to do so lightly.

Okay, but in fairness, a lot of packages do that — probably most famously the Dual of ForwardDiff.jl, and I wouldn’t say that’s very closely related.

That’s the big difference. Again, the SphericalFunctions.jl package I noted above was the driving force for me. The problem was that I had this enormous and complicated code, I knew that the final results were inaccurate, but I needed to see exactly where the problems were coming from.

That would be very cool. I guess the obvious problems are that (1) unlike the profilers, precision loss is not strictly monotonically increasing, and (2) even scalar-to-scalar code can involve intermediate steps with many values, each of which has its own precision problems. Still, I bet with all the existing infrastructure and AD experience in this community something truly awesome could be built.