This explanation is relatively reductive when it comes to its criticism of computational geometry.
The thing with computational geometry is, that its usually someone else's geometry, i.e you have no control over its quality or intention. In other words, whether two points or planes or lines actually align or align within 1e-4 is no longer really mathematically interesting because its all about the intention of the user: does the user think these planes overlap?.
This is why most geometry kernels (see open cascade) sport things like "fuzzy boolean operations" [0]) that lean into epsilons. These epsilons mask the error-prone supply chain of these meshes that arrive in your program by allowing some tolerance.
Finally, the remark "There are many ways of solving this problem" is also overly reductive, everyone reading here should really understand that this is a topic that is being actively researched right now in 2026, hence there are currently no blessed solutions to this problem, otherwise this research would not be needed. Even more so, to some extent this problem is fundamentally unsolvable depending on what you mean by "solvable", because your input is inexact not all geometrical operations are topologically valid, hence an "exact" or let alone "correct along some dimension" result cannot be achieved for all (combination of) inputs.
> This is why most geometry kernels (see open cascade) sport things like "fuzzy boolean operations" [0]) that lean into epsilons. These epsilons mask the error-prone supply chain of these meshes that arrive in your program by allowing some tolerance.
They don’t just lean into epsilons, the session context tolerance is used for almost every single point classification operation in geometric kernels and many primitives carry their own accumulating error component for downstream math.
Even then the current state of the art (in production kernels) is tolerance expansion where the kernel goes through up to 7 expansion steps retrying point classification until it just gives up. Those edge cases were some of the hardest parts of working on a kernel.
This is a fundamentally unsolvable problem with floating point math (I worked on both Parasolid and ACIS in the 2000s). Even the ray-box intersection example TFA gives is a long standing thorn - raytracing is one of the last fallbacks for nasty point classification problems.
I was working on an SDF-based CAD tool but gave up when I couldn't find a good way to do fillets.
It's very deceptive because the easy way works so well (Use smoothmin instead of min and you get smooth blends for free! You can even use a circular approximation of smoothmin and get proper fillets!). But when you want the user to be able to pick a couple of surfaces and fillet between them, it gets really hard.
It worked by rewriting the expression tree so that the blend arguments become sibling nodes and then applying the blend to the union/intersection that is their parent.
That works every time if you only want 1 single targeted blend, but if you want several of them then you can run into unsatisfiable cases where the same object needs to blended with several others and can't be siblings of all of them.
So I gave up :(. For me, CAD without fillets and chamfers is no CAD at all.
(Also, apropos for this thread: the discontinuity in the chamfer was a floating point precision problem...)
Well, user picking a couple of surfaces is literally an operation on a boundary representation, so of course it's a PITA with fields :)
I think the future is CAD is combined fields and breps. They're literally dual, one is covariant, the other contravariant (breps facilitate pushforwards, fields facilitate pullbacks).
One without the other is necessarily going to be limited in some way.
The distance field tells you the distance to the nearest surface at any point. You can have a "surface id field" tell you the id of the nearest surface to any point, and then when you raymarch to find the intersection of a line with a surface, you can read out the ID from the ID field after finding the intersection point. (Of course the ID field is also implemented as a function mapping points to surfaces).
So when the mouse is hovered or clicked in the 3d view you can easily find the ID of the surface under the pointer, and you can draw that surface in a different colour to show it is selected. No boundary representation needed.
The hard part is, given 2 surface ids, how do you add a fillet between them in the general case?
Another idea I had was to set the fillet radius on every min/max node based on the derived surface id's from the child nodes, but I couldn't find a good way to do this without making the field discontinuous.
That depends on what we mean by surfaces, and in the case of filleting, the user really wants to be picking adjacent faces (as in: an edge between two adjacent faces). That, or even a region to roll a ball along to generate a fillet.
The semantics of fillets even in the simplest case is that it's doing something to the edges, i.e. elements of the boundary representation, so that's a more natural structure for filleting.
>The distance field tells you the distance to the nearest surface at any point.
What you're describing isn't the same. You really are picking solids, not faces.
This wouldn't work even in the simplest case of a cube.
You can define a cube by a distance field:
f(x, y, z) = max(|x|, |y|, |z|) - 1
If the user wants to fillet just one the edges, then what? You only have one surface (the boundary of a cube), and one ID.
The field doesn't know anything about the edges.
OK, OK, we can ignore this edge case (badum-tss), but even if you only allow filleting "two surfaces", those two "surfaces" (really: boundaries of solids) aren't necessarily going to intersect along one edge (which is what the user wants to fillet).
The intersection may have multiple components. Or may not be manifold.
As a concrete example:
f(x, y, z) = z - cos(x)
g(x, y, z) = z - cos(y)
Look ma, no absolute values! Let me smooth-talk a little though:
f(x, y, z) = z - cos(x)cos(y)
g(x, y, z) = z - 0.25
.....and that's before we get to the reality where the user's understanding of "edge" isn't topological (as in, component of intersection between surfaces), but geometric (sharp corners).
B-reps can get away with making no distinction between them... Until you have messy geometry from elsewhere.
Say, an STL file from a scan. Or a mesh that came from an F-rep by marching cubes. Or whatever unholy mess OpenSCAD can cook with CGAL.
It doesn't matter if you use F-rep or convert to one: chisel out a cube as an intersection of half-spaces, then cut with half-spaces that narrowly touch the edges.
It'll look like a cube, and it'll be a cube functionally if you manufacture it.
Good luck with that one.
>If you have good ideas for this I'd love to hear them and resume working on Isoform
Well. The good news is that putting fillets on every edge is kind of easy with fields because you can do it with offsets.
If F(x, y, z) is a distance field that defines a solid, G(x, y, z) = F(x, y, z) + c offsets F inwards by c.
G is not a distance field anymore though, it's giving values that arent distances on the outside of convex corners.
Renormalize G to be a distance field, call it G'.
Now offset G' outwards by c: H = G' - c.
Ta-da! Concave corners aren't touched, convex corners are rounded.
Flip the + and -, and you're filleting concave corners (G = F - c is a field that defines an outwards offset that fails to be a distance field inside the body near concave corners; compute G' — the distance field for G; offset G inwards: H = G' + c).
Now, the "just normalize a field into a distance field" is doing a lot of heavy lifting here.
It works in the case of a cube if you define the cube to be the intersection of 6 half-spaces. There is a video demonstration of it working (partly) on a cube defined this way in the YouTube link in my comment above.
I define a surface to be the region of space where a particular SDF evaluates to 0. You define a solid to be the region of space where that SDF evaluates to <0, but they're broadly the same concept.
It is no problem to ensure that all primitives & all extruded sketches are defined so that each face gets a different surface id, and you would of course want to do this if you want to be able to fillet them.
You're right that there is a difference between an edge and between a pair of surfaces, but finding edges in SDFs is much harder than finding pairs of surfaces. If they intersect along more than one edge then you'll get the fillet along more than one edge. SDFs don't even have concrete "edges" in the general case. I'm not worried about this. Being able to fillet the intersection of 2 surfaces (solids) would satisfy me, but I haven't even got that far.
I'm not trying to find a solution that involves treating edges as "special". That's B-rep thinking. I don't mind if a "fillet" between 2 surfaces that do not touch but are closer together than the fillet radius creates a bridge between them, as long as it is smooth, continuous, and predictable.
It doesn't have to approximate the B-rep way, it just needs to be practically useful for turning sharp edges into rounded ones in a way that lets the user decide where.
> They don’t just lean into epsilons, the session context tolerance is used
for almost every single point classification operation in geometric kernels
and many primitives carry their own accumulating error component for downstream math.
The GP wasn't wrong. To "lean in" means to fully commit to, go all in on, (or, equivalently, go all out on).
I think his point is: rather than "leaning into" it as in, masking through epsilons, he argues that tolerance is fundamental to the problem space, not a way to resolve edge cases.
To me it seems it's used with the intent "They don’t just , they ," implying that Y is a proper superset of X. My point is that X is in fact a superset of Y, making the most charitable reading "They don’t just , they ."
Is there another potential reading of "just" that I'm missing?
im surprised terminology isnt borrowed from mechanical engineering on the type of fit that two pieces are supposed to have. Interference fits vs clearance do a physical job of describing whats happening
Something I've observed as someone who works in the physical sciences and used to work as a software engineer is:
Very few software engineers understand that tolerances are fundamental.
In the physical sciences, strict equality - of actual quantities, not variables in equations - is almost never a thing. Even though an equation might show that two theoretical quantities are exactly equal, the practical fact is that every quantity begins life as a measurement, and measurements have inherent uncertainty.
Even in design work, nothing is exact. It's simply not possible. A resistor's value is subject to manufacturing tolerances, will vary with temperature, and can drift as it ages. A mechanical part will also have manufacturing tolerance, and changes size and shape with temperature, applied forces, and wear. So even if a spec sheet states an exact number, the heading or notes will tell you that this is a nominal value under specific conditions. (Those conditions are also impossible to achieve and maintain exactly for all the same reasons.)
Even the voltages that represent 0 and 1 inside a computer aren't exact. Digital parts like CPUs, GPUs, RAM, etc. specify low and high thresholds, under or over which a voltage is considered a 0 or 1.
Floating-point numbers have uses outside the physical sciences, so there's no one-size-fits-all approach to using them correctly. But if you are writing code that deals with physical quantities, making equality comparisons is almost always going to be wrong even if floating-point numbers had infinite precision and no rounding error. Physical quantities simply can't be used that way.
If there's a way to make it more precise and/or specific and/or faster, or create similar macros with better functionality and/or correctness, that's great.
See the same directory for corresponding assert_* macros for less than, greater than, etc.
Is there any constant more misused in compsci than ieee epsilon? :)
It's defined as the difference between 1.0 and the smallest number larger than 1.0. More usefully, it's the spacing between adjacent representable float numbers in the range 1.0 to 2.0.
Because floats get less precise at every integer power of two, it's impossible for two numbers greater than or equal to 2.0 to be epsilon apart. The spacing between 2.0 and the next larger number is 2*epsilon.
That means abs(a - b) <= epsilon is equivalent to a == b for any a or b greater than or equal to 2.0. And if you use < then the limit will be 1.0 instead.
Epsilon is the wrong tool for the job in 99.9% of cases.
A (perhaps initially) counterintuitive part of the above more explicitly stated: The doubling/halving also means numbers between 0 and 1 actually have _more_ precision than the epsilon would suggest.
Considerably more in many cases. The point of floating point is to have as many distinct values in the range 2-4 as are in the range 1-2 as are between 1/2 and 1, 1/4 and 1/2, 1/8 and 1/4, etc. the smallest representable difference between consecutive floating point numbers down around the size of 1/64 is on the order of epsilon/64
Multiplying epsilon by the largest number you are dealing with is a strategy that makes using epsilons at least somewhat logical.
So I'd probably rewrite that code to first find the ulp of the larger of the abs of a and b and then assert that their difference is less than or equal to that.
Edit: Or maybe the smaller of the abs of the two, I haven't totally thought through the consequences. It might not matter, because the ulps will only differ when the numbers are significantly apart and then it doesn't matter which one you pick. Perhaps you can just always pick the first number and get its ULP.
This is what was done to a raytracer I used. People kept making large-scale scenes with intricate details, think detailed ring placed on table in a room with a huge field in view through the window. For a while one could override the fixed epsilon based on scene scale, but for such high dynamic range scenes a fixed epsilon just didn't cut it.
IIRC it would compute the "dynamic" epsilon value essentially by adding one to the mantissa (treated as an integer) to get the next possible float. Then subtract from that the initial value to get the dynamic epsilon value.
Definitely use library functions if you got 'em though.
Because of the representation of floats, couldn't you just bitwise cast to uints and see if the (abs) difference was less than or equal to one? But practically you probably should check if it's less than or equal to say ten, depending on your tolerance.
It would be very useful to be able to compare the significant directly then. I realize there is a boundary issue when a significant is very close to 0x00..000 or 0xFFF..FFF
Everyone has already made several comments on the incorrect use of EPSILON here, but there's one more thing I want to add that hasn't yet been mentioned:
EPSILON = (1 ulp for numbers in the range [1, 2)). is a lousy choice for tolerance. Every operation whose result is in the range [1, 2) has a mathematical absolute error of ½ ulp. Doing just a few operations in a row has a chance to make the error term larger than your tolerance, simply because of the inherent inaccuracy of floating-point operations. Randomly generate a few doubles in the range [1, 10], then randomize the list and compute the sum of different random orders in the list, and your assertion should fail. I'd guess you haven't run into this issue because either very few people are using this particular assertion, or the people who do happen to be testing it in cases where the result is fully deterministic.
If you look at professional solvers for numerical algorithms, one of the things you'll notice is that not only is the (relative!) tolerance tunable, but there's actually several different tolerance values. The HiGHS linear solver for example uses 5 different tolerance values for its simplex algorithm. Furthermore, the default values for these tolerances tend to be in the region of 10^-6 - 10^-10... about the square root of f64::EPSILON. There's a basic rule of thumb in numerical analysis that you need your internal working precision to be roughly twice the number of digits as your output precision.
Your last comment is essential for numerical analysis, indeed. There is this "surprising" effect that increasing the precision of the input ends up by decreasing that of the output (roughly speaking). So "I shall just use a s very small discretization" can be harmful.
Your assertion code here doesn't make a ton of sense. The epsilon of choice here is the distance between 1 and the next number up, and it's completely separated from the scale of the numbers in question. 1e-50 will compare equal to 2e-50, for example.
I would suggest that "equals" actually is for "exactly equals" as in (a == b). In many pieces of floating point code this is the correct thing to test. Then also add a function for "within range of" so your users can specify an epsilon of interest, using the formula (abs(a - b) < eps). You may also want to support multidimensional quantities by allowing the user to specify a distance metric. You probably also want a relative version of the comparison in addition to an absolute version.
Auto-computing epsilons for an equality check is really hard and depends on the usage, as well as the numerics of the code that is upstream and downstream of the comparison. I don't see how you would do it in an assertion library.
Ignoring the misuse of epsilon, I'd also say that you'd be helping your users more by not providing a general assert_f64_eq macro, but rather force the user to decide the error model. Add a required "precision" parameter as an enum with different modes:
// Precise matching:
assert_f64_eq!(a, 0.1, Steps(2))
// same as: assert!(a == 0.1.next_down().next_down())
// Number of digits (after period) that are matching:
assert_f64_eq!(a, 0.1, Digits(5))
// Relative error:
assert_f64_eq!(a, 0.1, Rel(0.5))
You generally want both relative and absolute tolerances. Relative handles scale, absolute handles values near zero (raw EPSILON isn’t a universal threshold per IEEE 754).
The usual pattern is abs(a - b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol) to avoid both large-value and near-zero pitfalls.
It depends on the use case, but do you consider NaN to be equal to NaN? For an assert macro, I would expect so. Also, your code works differently for very large and very small numbers, eg. 1.0000001, 1.0000002 vs 1e-100, 1.0000002e-100.
For my own soft-floating point math library, I expect the value is off by a some percentage, not just off by epsilon. And so I have my own almostSame method [1] which accounts for that and is quite a bit more complex. Actually multiple such methods. But well, that's just my own use case.
Machine eps provides the maximum rounding error for a single op. Let's say I write:
let y = 2.0;
let x = sqrt(y);
Now is xactually the square root of 2? Of course not - because the digit expansion of sqrt(2) doesn't terminate, the only way to precisely represent it is with symbolics. So what do we actually have? x was either rounded up or down to a number that does have an exact FP representation. So, x / sqrt(2) is in [1 - eps, 1 + eps]. The eps tells you, on a relative scale, the maximum distance to an adjacent FP number for any real number. (Full disclosure, IDK how this interacts with weird stuff like denormals).
Note that in general we can only guarantee hitting this relative error for single ops. More elaborate computations may develop worse error as things compound. But it gets even worse. This error says nothing about errors that don't occur in the machine. For example, say I have a test that takes some experimental data, runs my whiz-bang algorithm, and checks if the result is close to elementary charge of an electron. Now I can't just worry about machine error but also a zillion different kinds of experimental error.
There are also cases where we want to enforce a contract on a number so we stay within acceptable domains. Author alluded to this. For example - if I compute some x s.t. I'm later going to take acos(x), x had better be between [-1, 1]. x >= -1 - EPS && x <= 1 + EPS wouldn't be right because it would include two numbers, -1 - EPS and 1 + EPS, that are outside the acceptable domain.
- "I want to relax exact equality because my computation has errors" -> Make assert_rel_tol and assert_abs_tol.
- "I want to enforce determinism" -> exact equality.
- "I want to enforce a domain" -> exact comparison
Your code here is using eps for controlling absolute error, which is already not great since eps is about relative error. Unfortunately your assertion degenerates to a == b for large numbers but is extremely loose for small numbers.
Apart from what others have commented, IMO an “assertables” crate should not invent new predicates of its own, especially for domains (like math) that are orthogonal to assertability.
EQ should be exactly equal, I think. Although we often (incorrectly) model floats as a real plus some non-deterministic error, there are cases where you can expect an exact bit pattern, and that’s what EQ is for (the obvious example is, you could be writing a library and accept a scaling factor from the user—scaling factors of 1 or 0 allow you to optimize).
You probably also want an isclose and probably want to push most users toward using that.
Numeric comparison implies subtraction or similar a - b sign and zero extraction at some lower level (in an ALU micro-op perhaps), so that's duplicated effort of a - b unless it can be optimized away.
match a - b {
d if d >= 0.0 => d < f64::EPSILON,
d => d >= -f64::EPSILON, /* true if -EPSILON to -0.0 */
}
Um, what? You've linked an RFC for Rust, but the CPP Reference article for C++ So yeah, the Rust RFC documents a proposed change, and the C++ reference documents an implemented feature, but you could equally link the C++ Proposal document and the Rust library docs to make the opposite point if you wanted.
You can also rely on the fact (not promised in C++) that these are actually the IEEE floats and so they have all the resulting properties you can (entirely in safe Rust) just ask for the integers with the same bit pattern, compare integers and because of how IEEE is designed that tells you how far away in some proportional sense, the two values are.
On an actual CPU manufactured this century that's almost free because the type system evaporates during compilation -- for example f32::to_bits is literally zero CPU instructions.
Oh, my research was wrong and the line from the RFC doc...
>Currently it is not possible to answer the question ‘which floating point value comes after x’ in Rust without intimate knowledge of the IEEE 754 standard.
So nevermind on it not being present in Rust I guess I was finding old documentation
> In reality it is a pretty deterministic (modulo compiler options, CPU flags, etc)
IIRC this was not ALWAYS the case, on x86 not too long ago the CPU might choose to put your operation in an 80-bit fp register, and if due to multitasking the CPU state got evicted, it would only be able to store it in a 32-bit slot while it's waiting to be scheduled back in?
It might not be the case now in a modern system if based on load patterns the software decides to schedule some math operations or another on the GPU vs the CPU, or maybe some sort of corner case where you are horizontally load balancing on two different GPUs (one AMD, one Nvidia) -- I'm speculating here.
Think about this. It's silly to use floating point numbers to represent geometry, because it gives coordinates closer to the origin more precision and in most cases the origin is just an arbitrary point.
The thing with floating point numbers is they are meant to work with physical quantities: distances, durations, etc...
Physical quantities involve imprecision: measurement devices, tools, display devices, ADC/DACs etc... They all have some tolerances. And when you are using epsilons, the epsilon value should be chosen based on that physical value. For example, you set the epsilon to 1e-4 because that's 100 microns and you can't display 100 micron details.
That's also the reason why there is not one size fits all solution. If you are working with microscopic objects, 100 microns is huge, and if you are doing a space simulation, 1 km may be negligible. Some operations involve a huge loss of precision, some don't, and sometimes you really want exact numbers and therefore you have to know your fractional powers of 2.
I guess I'm confused. I thought epsilon was the smallest possible value to account for accuracy drift across the range of a floating point representation, not just "1e-4".
Done some reading. Thanks to the article to waking me up to this fact at least. I didn't realize that the epsilon provided by languages tends to be the one that only works around 1.0, and if you want to use episilons globally (which the article would say is generally a bad idea) you need to be more dynamic as your ranges, and potential errors, increase.
The problem with floating point comparison is not that it's nondeterministic, it's that what should be the same number may have different representations, often with different rounding behavior as well, so depending on the exact operations you use to arrive at it, it may not compare as equal, hence the need for the epsilon trick.
If all you're comparing is the result from the same operations, you _may_ be fine using equality, but you should really know that you're never getting a number from an uncontrolled source.
This is mostly about game logic, where I can understand the reliance on floating point numbers. I've also seen these epsilon comparisons in code that had nothing to do with game engines or positions in continuous space, and it has always hurt my eyes.
I think if you want to work with values that might be exactly equal to other values, floating point is simply not the right choice. For money, use BigDecimal or something like that. For lots of purposes, int might be more appropriate. If you do need floating point, maybe compare whether the value is larger than the other value.
I cringed very hard in the slerp example seeing acos(dot(a,b)). Clamping to [-1,1] to avoid NaNs still gives you bad answers and numerical sensitivity around small angles. acos and asin in general lose ~half your sig figs around their singularities[0]. Working around this by introducing a threshold seems like exactly same flavor of issue he's complaining about to begin with.
There are perfectly good solutions for finding the angle using atan and the cross product. Calculating A x B will yield a vector which lies is along their normal with length tan(θ)(A · B) so we can straightforwardly say e.g.:
θ = atan2(norm(A ⅹ B), (A · B))
No thresholds, no branching, only ~1 bit of significance lost. As the original author says himself, when you're introducing arbitrary-feeling thresholds, it's likely you're missing a solution which would improve more than just the weirdness at those thresholds.
> > the myth about exactness is that you can't use strict equality with floating point numbers because they are somehow fuzzy. They are not.
> They are though. All arithmetic operations involve rounding, so e.g. (7.0 / 1234 + 0.5) * 1234 is not equal to 7.0 + 617 (it differs in 1 ULP). On the other hand, (9.0 / 1234 + 0.5) * 1234 is equal to 9.0 + 617, so the end result is sometimes exact and sometimes is not. How can you know beforehand which one is the case in your specific case? Generally, you can't, any arithmetic operation can potentially give you 1 ULP of error, and it can (and likely, will) slowly accumulate.
Also, please don't comment how nobody has a use for "f(x) = (x / 1234 + 0.5) * 1234": there are all kinds of queer computations people do in floating point, and for most of them, figuring out the exactness of the end result requires an absurd amount of applied numerical analysis, doing which would undermine most of the "just let the computer crunch the numbers" point of doing this computation on a computer.
My normal issue with floating-point epsilon shenanigans is that they don't usually pass the sniff test, suggesting something fundamentally wrong with the problem framing or its solution.
It's a classic, so let's take vector normalization as an example. Topologically, you're ripping a hole in the space, and that's causing your issues. It manifests as NaN for length-zero vectors, weird precision issues too close to zero, etc, but no matter what you employ to try to fix it you're never going to have a good time squishing N-D space onto the surface of an N-D sphere if you need it to be continuous.
Some common subroutines where I see this:
1. You want to know the average direction of a bunch of objects and thus have to normalize each vector contributing to that average. Solution 1: That's not what you want almost ever. In any of the sciences, or anything loosely approximating the real world, you want to average the un-normalized vectors 99.999% of the time. Solution 2: Maybe you really do need directions for some reason (e.g., tracking where birds are looking in a game). Then don't rely on vectors for your in-band signaling. Explicitly track direction and magnitude separately and observe the magic of never having direction-related precision errors.
2. You're doing some sort of lighting normalization and need to compute something involving areas of potentially near-degenerate triangles, dividing by those values to weight contributions appropriately. Solution: Same as above, this is kind of like an average of averages problem. It can make fuzzy, intuitive sense, but you'll get better results if you do your summing and averaging in an un-normalized space. If you really do need surface normals, store those explicitly and separate from magnitude.
3. You're doing some sort of ML voodoo to try to get better empirical results via some vague appeal to vanishing gradients or whatever. Solution: The core property you want is a somewhat strange constraint on your layer's Jacobian matrix, and outside of like two papers nobody is willing to put up with the code complexity or runtime costs, even when they recognize it as the right thing to do. Everything you're doing is a hack anyway, so make your normalization term x/(|x|+eps) with eps > 0 rather than equal to zero like normal. Choose eps much smaller than most of the vectors you're normalizing this way and much bigger than zero. Something like 1e-3, 1e-20, and 1e-150 should be fine for f16, f32, and f64. You don't have to tune because it's a pretty weak constraint on the model, and it's able to learn around it.
Unless you use any compiler I write. As a purist, '==' for floats and doubles will be undefined. And by undefined, I mean unrecoverable. And by unrecoverable, I mean they may be able to extract most of the fragments of silicon, aluminum, and dynamic island, and "put you back together again". But you will "think different".
Unless you are comparing references. References are ok.
My preference in tests is a little different than just using IEEE 754 ==,
_Bool equiv(float x, float y) {
return (x <= y && y <= x)
|| (x != x && y != y);
}
which both handles NaNs sensibly (all NaNs are equivalent) and won't warn about using == on floats. I find it also easy to remember how to write when starting a new project.
My one small nitpick is that vector length is usually 2 instructions with SSE4:
dpps xmm0, xmm0, 0x17 ; dot product of 3 lanes, write lane 0
sqrtss xmm0, xmm0
ret
And is considerably faster than the fancy version, mainly because Intel still hasn't given us horizontal-max vector instruction! ARM is a bit better in that regard with their fancy vmaxvq_f32 and vmaxnmvq_f32...
There is another way to compare floats for rough equality that I haven't seen much explored anywhere: bit-cast to integer, strip few least significant bits and then compare for equality.
This is agnostic to magnitude, unlike epsilon which has to be tuned for range of values you expect to get a meaningful result.
If your code may be compiled, to use the Intel x87 numerical coprocessor, an important issue is the so called "excess precision": Different values on chip can collapse after being rounded and stored to their memory locations, invalidating previous comparisons. Spilling can happen unexpectedly. Note that Intel calls the x87 "legacy"
It's a good advice to avoid using floating point at all, if it's not that hard to do so. In one of my hobby projects I have written a simple map application based on OSM data with geometry processing code operation on integers only.
I used floating timestamps as some kind of an identity. If there is ever a conflict, I just increase it by 1 ulp until it doesn't collide with anything. Sorry.
One of the goals of comparing floating points with an epsilon is precisely so that you can apply these types of accuracy increasing (or decreasing) changes to the operations, and still get similar results.
Anything else is basically a nightmare to however has to maintain the code in the future.
Also, good luck with e.g. checking if points are aligned to a grid or the like without introducing a concept of epsilon _somewhere_.
So they say you should use epsilons, then their solution to the first problem is to use an epsilon. There may be some cases when you can get by without using epsilon comparison, but in many cases epsilon comparison is the right thing to do, but you need to choose a good value for it.
This is especially true in cases where the number comes from some kind of input (user controls, sensor reading, etc.) or random number generation.
144 comments
The thing with computational geometry is, that its usually someone else's geometry, i.e you have no control over its quality or intention. In other words, whether two points or planes or lines actually align or align within 1e-4 is no longer really mathematically interesting because its all about the intention of the user: does the user think these planes overlap?.
This is why most geometry kernels (see open cascade) sport things like "fuzzy boolean operations" [0]) that lean into epsilons. These epsilons mask the error-prone supply chain of these meshes that arrive in your program by allowing some tolerance.
Finally, the remark "There are many ways of solving this problem" is also overly reductive, everyone reading here should really understand that this is a topic that is being actively researched right now in 2026, hence there are currently no blessed solutions to this problem, otherwise this research would not be needed. Even more so, to some extent this problem is fundamentally unsolvable depending on what you mean by "solvable", because your input is inexact not all geometrical operations are topologically valid, hence an "exact" or let alone "correct along some dimension" result cannot be achieved for all (combination of) inputs.
[0] https://dev.opencascade.org/content/fuzzy-boolean-operations
They don’t just lean into epsilons, the session context tolerance is used for almost every single point classification operation in geometric kernels and many primitives carry their own accumulating error component for downstream math.
Even then the current state of the art (in production kernels) is tolerance expansion where the kernel goes through up to 7 expansion steps retrying point classification until it just gives up. Those edge cases were some of the hardest parts of working on a kernel.
This is a fundamentally unsolvable problem with floating point math (I worked on both Parasolid and ACIS in the 2000s). Even the ray-box intersection example TFA gives is a long standing thorn - raytracing is one of the last fallbacks for nasty point classification problems.
Could you point to any literature/freely available resource that comes close to the SOTA for these kinds of operations? I would be greatly helped.
> This is a fundamentally unsolvable problem with floating point math
It's a fundamentally unsolvable problem with B-reps! The problem completely disappears with F-reps. (In exchange for some other difficult problems).
Ahhaha.
(I used to work in nTop, and boy is this an understatement when it comes to field based solid modeling)
It's very deceptive because the easy way works so well (Use smoothmin instead of min and you get smooth blends for free! You can even use a circular approximation of smoothmin and get proper fillets!). But when you want the user to be able to pick a couple of surfaces and fillet between them, it gets really hard.
This is the best I got: https://www.youtube.com/watch?v=LOvqdlDbkBs
It worked by rewriting the expression tree so that the blend arguments become sibling nodes and then applying the blend to the union/intersection that is their parent.
That works every time if you only want 1 single targeted blend, but if you want several of them then you can run into unsatisfiable cases where the same object needs to blended with several others and can't be siblings of all of them.
So I gave up :(. For me, CAD without fillets and chamfers is no CAD at all.
(Also, apropos for this thread: the discontinuity in the chamfer was a floating point precision problem...)
I think the future is CAD is combined fields and breps. They're literally dual, one is covariant, the other contravariant (breps facilitate pushforwards, fields facilitate pullbacks).
One without the other is necessarily going to be limited in some way.
The distance field tells you the distance to the nearest surface at any point. You can have a "surface id field" tell you the id of the nearest surface to any point, and then when you raymarch to find the intersection of a line with a surface, you can read out the ID from the ID field after finding the intersection point. (Of course the ID field is also implemented as a function mapping points to surfaces).
So when the mouse is hovered or clicked in the 3d view you can easily find the ID of the surface under the pointer, and you can draw that surface in a different colour to show it is selected. No boundary representation needed.
The hard part is, given 2 surface ids, how do you add a fillet between them in the general case?
Another idea I had was to set the fillet radius on every min/max node based on the derived surface id's from the child nodes, but I couldn't find a good way to do this without making the field discontinuous.
I have more notes in this blog post: https://incoherency.co.uk/blog/stories/frep-cad-building-blo...
If you have good ideas for this I'd love to hear them and resume working on Isoform.
That depends on what we mean by surfaces, and in the case of filleting, the user really wants to be picking adjacent faces (as in: an edge between two adjacent faces). That, or even a region to roll a ball along to generate a fillet.
The semantics of fillets even in the simplest case is that it's doing something to the edges, i.e. elements of the boundary representation, so that's a more natural structure for filleting.
>The distance field tells you the distance to the nearest surface at any point.
What you're describing isn't the same. You really are picking solids, not faces.
This wouldn't work even in the simplest case of a cube.
You can define a cube by a distance field:
If the user wants to fillet just one the edges, then what? You only have one surface (the boundary of a cube), and one ID.The field doesn't know anything about the edges.
OK, OK, we can ignore this edge case (badum-tss), but even if you only allow filleting "two surfaces", those two "surfaces" (really: boundaries of solids) aren't necessarily going to intersect along one edge (which is what the user wants to fillet).
The intersection may have multiple components. Or may not be manifold.
As a concrete example:
Look ma, no absolute values! Let me smooth-talk a little though: .....and that's before we get to the reality where the user's understanding of "edge" isn't topological (as in, component of intersection between surfaces), but geometric (sharp corners).B-reps can get away with making no distinction between them... Until you have messy geometry from elsewhere.
Say, an STL file from a scan. Or a mesh that came from an F-rep by marching cubes. Or whatever unholy mess OpenSCAD can cook with CGAL.
It doesn't matter if you use F-rep or convert to one: chisel out a cube as an intersection of half-spaces, then cut with half-spaces that narrowly touch the edges.
It'll look like a cube, and it'll be a cube functionally if you manufacture it.
Good luck with that one.
>If you have good ideas for this I'd love to hear them and resume working on Isoform
Well. The good news is that putting fillets on every edge is kind of easy with fields because you can do it with offsets.
If F(x, y, z) is a distance field that defines a solid, G(x, y, z) = F(x, y, z) + c offsets F inwards by c.
G is not a distance field anymore though, it's giving values that arent distances on the outside of convex corners.
Renormalize G to be a distance field, call it G'.
Now offset G' outwards by c: H = G' - c.
Ta-da! Concave corners aren't touched, convex corners are rounded.
Flip the + and -, and you're filleting concave corners (G = F - c is a field that defines an outwards offset that fails to be a distance field inside the body near concave corners; compute G' — the distance field for G; offset G inwards: H = G' + c).
Now, the "just normalize a field into a distance field" is doing a lot of heavy lifting here.
But that can give you something to think about.
I define a surface to be the region of space where a particular SDF evaluates to 0. You define a solid to be the region of space where that SDF evaluates to <0, but they're broadly the same concept.
It is no problem to ensure that all primitives & all extruded sketches are defined so that each face gets a different surface id, and you would of course want to do this if you want to be able to fillet them.
You're right that there is a difference between an edge and between a pair of surfaces, but finding edges in SDFs is much harder than finding pairs of surfaces. If they intersect along more than one edge then you'll get the fillet along more than one edge. SDFs don't even have concrete "edges" in the general case. I'm not worried about this. Being able to fillet the intersection of 2 surfaces (solids) would satisfy me, but I haven't even got that far.
I'm not trying to find a solution that involves treating edges as "special". That's B-rep thinking. I don't mind if a "fillet" between 2 surfaces that do not touch but are closer together than the fillet radius creates a bridge between them, as long as it is smooth, continuous, and predictable.
It doesn't have to approximate the B-rep way, it just needs to be practically useful for turning sharp edges into rounded ones in a way that lets the user decide where.
> They don’t just lean into epsilons, the session context tolerance is used for almost every single point classification operation in geometric kernels and many primitives carry their own accumulating error component for downstream math.
The GP wasn't wrong. To "lean in" means to fully commit to, go all in on, (or, equivalently, go all out on).
I'm wondering if people have heard the expression "leaning in" from people who were insincere/lying, and assumed that that was what the phrase means?
Is there another potential reading of "just" that I'm missing?
Very few software engineers understand that tolerances are fundamental.
In the physical sciences, strict equality - of actual quantities, not variables in equations - is almost never a thing. Even though an equation might show that two theoretical quantities are exactly equal, the practical fact is that every quantity begins life as a measurement, and measurements have inherent uncertainty.
Even in design work, nothing is exact. It's simply not possible. A resistor's value is subject to manufacturing tolerances, will vary with temperature, and can drift as it ages. A mechanical part will also have manufacturing tolerance, and changes size and shape with temperature, applied forces, and wear. So even if a spec sheet states an exact number, the heading or notes will tell you that this is a nominal value under specific conditions. (Those conditions are also impossible to achieve and maintain exactly for all the same reasons.)
Even the voltages that represent 0 and 1 inside a computer aren't exact. Digital parts like CPUs, GPUs, RAM, etc. specify low and high thresholds, under or over which a voltage is considered a 0 or 1.
Floating-point numbers have uses outside the physical sciences, so there's no one-size-fits-all approach to using them correctly. But if you are writing code that deals with physical quantities, making equality comparisons is almost always going to be wrong even if floating-point numbers had infinite precision and no rounding error. Physical quantities simply can't be used that way.
The Rust code in the assert_f64_eq macro is:
I'm the author of the Rust assertables crate. It provides floating-point assert macros much as described in the article.https://github.com/SixArm/assertables-rust-crate/blob/main/s...
If there's a way to make it more precise and/or specific and/or faster, or create similar macros with better functionality and/or correctness, that's great.
See the same directory for corresponding assert_* macros for less than, greater than, etc.
It's defined as the difference between 1.0 and the smallest number larger than 1.0. More usefully, it's the spacing between adjacent representable float numbers in the range 1.0 to 2.0.
Because floats get less precise at every integer power of two, it's impossible for two numbers greater than or equal to 2.0 to be epsilon apart. The spacing between 2.0 and the next larger number is 2*epsilon.
That means
abs(a - b) <= epsilonis equivalent toa == bfor any a or b greater than or equal to 2.0. And if you use<then the limit will be 1.0 instead.Epsilon is the wrong tool for the job in 99.9% of cases.
Multiplying epsilon by the largest number you are dealing with is a strategy that makes using epsilons at least somewhat logical.
So I'd probably rewrite that code to first find the ulp of the larger of the abs of a and b and then assert that their difference is less than or equal to that.
Edit: Or maybe the smaller of the abs of the two, I haven't totally thought through the consequences. It might not matter, because the ulps will only differ when the numbers are significantly apart and then it doesn't matter which one you pick. Perhaps you can just always pick the first number and get its ULP.
IIRC it would compute the "dynamic" epsilon value essentially by adding one to the mantissa (treated as an integer) to get the next possible float. Then subtract from that the initial value to get the dynamic epsilon value.
Definitely use library functions if you got 'em though.
epsilons are fine in the case that you actually want to put a static error bound on an equality comparison. numpy's relative errors are better for floats at arbitrary scales (https://numpy.org/doc/stable/reference/generated/numpy.isclo...).
edit: ahh i forgot all about ulps. that is what people often confuse ieee eps with. also, good background material in the necronomicon (https://en.wikipedia.org/wiki/Numerical_Recipes).
EPSILON = (1 ulp for numbers in the range [1, 2)). is a lousy choice for tolerance. Every operation whose result is in the range [1, 2) has a mathematical absolute error of ½ ulp. Doing just a few operations in a row has a chance to make the error term larger than your tolerance, simply because of the inherent inaccuracy of floating-point operations. Randomly generate a few doubles in the range [1, 10], then randomize the list and compute the sum of different random orders in the list, and your assertion should fail. I'd guess you haven't run into this issue because either very few people are using this particular assertion, or the people who do happen to be testing it in cases where the result is fully deterministic.
If you look at professional solvers for numerical algorithms, one of the things you'll notice is that not only is the (relative!) tolerance tunable, but there's actually several different tolerance values. The HiGHS linear solver for example uses 5 different tolerance values for its simplex algorithm. Furthermore, the default values for these tolerances tend to be in the region of 10^-6 - 10^-10... about the square root of f64::EPSILON. There's a basic rule of thumb in numerical analysis that you need your internal working precision to be roughly twice the number of digits as your output precision.
I would suggest that "equals" actually is for "exactly equals" as in (a == b). In many pieces of floating point code this is the correct thing to test. Then also add a function for "within range of" so your users can specify an epsilon of interest, using the formula (abs(a - b) < eps). You may also want to support multidimensional quantities by allowing the user to specify a distance metric. You probably also want a relative version of the comparison in addition to an absolute version.
Auto-computing epsilons for an equality check is really hard and depends on the usage, as well as the numerics of the code that is upstream and downstream of the comparison. I don't see how you would do it in an assertion library.
assert_f64_eqmacro, but rather force the user to decide the error model. Add a required "precision" parameter as an enum with different modes:The usual pattern is abs(a - b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol) to avoid both large-value and near-zero pitfalls.
https://github.com/python/cpython/blob/d61fcf834d197f0113a6a...
For my own soft-floating point math library, I expect the value is off by a some percentage, not just off by epsilon. And so I have my own almostSame method [1] which accounts for that and is quite a bit more complex. Actually multiple such methods. But well, that's just my own use case.
[1] https://github.com/thomasmueller/bau-lang/blob/main/src/test...
xactually the square root of 2? Of course not - because the digit expansion of sqrt(2) doesn't terminate, the only way to precisely represent it is with symbolics. So what do we actually have?xwas either rounded up or down to a number that does have an exact FP representation. So,x/ sqrt(2) is in[1 - eps, 1 + eps]. The eps tells you, on a relative scale, the maximum distance to an adjacent FP number for any real number. (Full disclosure, IDK how this interacts with weird stuff like denormals).Note that in general we can only guarantee hitting this relative error for single ops. More elaborate computations may develop worse error as things compound. But it gets even worse. This error says nothing about errors that don't occur in the machine. For example, say I have a test that takes some experimental data, runs my whiz-bang algorithm, and checks if the result is close to elementary charge of an electron. Now I can't just worry about machine error but also a zillion different kinds of experimental error.
There are also cases where we want to enforce a contract on a number so we stay within acceptable domains. Author alluded to this. For example - if I compute some
xs.t. I'm later going to takeacos(x),xhad better be between[-1, 1].x >= -1 - EPS && x <= 1 + EPSwouldn't be right because it would include two numbers, -1 - EPS and 1 + EPS, that are outside the acceptable domain.- "I want to relax exact equality because my computation has errors" -> Make
assert_rel_tolandassert_abs_tol.- "I want to enforce determinism" -> exact equality.
- "I want to enforce a domain" -> exact comparison
Your code here is using eps for controlling absolute error, which is already not great since eps is about relative error. Unfortunately your assertion degenerates to
a == bfor large numbers but is extremely loose for small numbers.https://numpy.org/doc/stable/reference/generated/numpy.allcl...
[1] https://arxiv.org/html/2403.07492v2
You probably also want an isclose and probably want to push most users toward using that.
C++ implements this https://en.cppreference.com/cpp/numeric/math/nextafter
Rust does not https://rust-lang.github.io/rfcs/3173-float-next-up-down.htm... but people have in various places.
Rust's https://doc.rust-lang.org/std/primitive.f32.html#method.next... https://doc.rust-lang.org/std/primitive.f32.html#method.next... of course exist, they're even actually constant expressions (the C++ functions are constexpr since 2023 but of course you're not promised they actually work as constant expressions because C++ is a stupid language and "constexpr" means almost nothing)
You can also rely on the fact (not promised in C++) that these are actually the IEEE floats and so they have all the resulting properties you can (entirely in safe Rust) just ask for the integers with the same bit pattern, compare integers and because of how IEEE is designed that tells you how far away in some proportional sense, the two values are.
On an actual CPU manufactured this century that's almost free because the type system evaporates during compilation -- for example f32::to_bits is literally zero CPU instructions.
>Currently it is not possible to answer the question ‘which floating point value comes after x’ in Rust without intimate knowledge of the IEEE 754 standard.
So nevermind on it not being present in Rust I guess I was finding old documentation
if a.abs()+b.abs() >= (a-b).abs() * 2f64.powi(48)
It remains accurate for small and for big numbers. 48 is slightly less than 52.
‘a.to_bits() == b.to_bits()’
Alternatively, use ‘partial_eq’ and fall back to bit equality if it returns None.
Fix your precision so it matches.
You only need so many significant digits.
> In reality it is a pretty deterministic (modulo compiler options, CPU flags, etc)
IIRC this was not ALWAYS the case, on x86 not too long ago the CPU might choose to put your operation in an 80-bit fp register, and if due to multitasking the CPU state got evicted, it would only be able to store it in a 32-bit slot while it's waiting to be scheduled back in?
It might not be the case now in a modern system if based on load patterns the software decides to schedule some math operations or another on the GPU vs the CPU, or maybe some sort of corner case where you are horizontally load balancing on two different GPUs (one AMD, one Nvidia) -- I'm speculating here.
Physical quantities involve imprecision: measurement devices, tools, display devices, ADC/DACs etc... They all have some tolerances. And when you are using epsilons, the epsilon value should be chosen based on that physical value. For example, you set the epsilon to 1e-4 because that's 100 microns and you can't display 100 micron details.
That's also the reason why there is not one size fits all solution. If you are working with microscopic objects, 100 microns is huge, and if you are doing a space simulation, 1 km may be negligible. Some operations involve a huge loss of precision, some don't, and sometimes you really want exact numbers and therefore you have to know your fractional powers of 2.
Done some reading. Thanks to the article to waking me up to this fact at least. I didn't realize that the epsilon provided by languages tends to be the one that only works around 1.0, and if you want to use episilons globally (which the article would say is generally a bad idea) you need to be more dynamic as your ranges, and potential errors, increase.
If all you're comparing is the result from the same operations, you _may_ be fine using equality, but you should really know that you're never getting a number from an uncontrolled source.
I think if you want to work with values that might be exactly equal to other values, floating point is simply not the right choice. For money, use BigDecimal or something like that. For lots of purposes, int might be more appropriate. If you do need floating point, maybe compare whether the value is larger than the other value.
acos(dot(a,b)). Clamping to [-1,1] to avoid NaNs still gives you bad answers and numerical sensitivity around small angles. acos and asin in general lose ~half your sig figs around their singularities[0]. Working around this by introducing a threshold seems like exactly same flavor of issue he's complaining about to begin with.There are perfectly good solutions for finding the angle using atan and the cross product. Calculating A x B will yield a vector which lies is along their normal with length tan(θ)(A · B) so we can straightforwardly say e.g.:
θ = atan2(norm(A ⅹ B), (A · B))No thresholds, no branching, only ~1 bit of significance lost. As the original author says himself, when you're introducing arbitrary-feeling thresholds, it's likely you're missing a solution which would improve more than just the weirdness at those thresholds.
[0] A good analysis of this is on Page 46 here, including a more optimized (but less obviously true) alternative to the formula above: https://people.eecs.berkeley.edu/~wkahan/Mindless.pdf
> > the myth about exactness is that you can't use strict equality with floating point numbers because they are somehow fuzzy. They are not.
> They are though. All arithmetic operations involve rounding, so e.g. (7.0 / 1234 + 0.5) * 1234 is not equal to 7.0 + 617 (it differs in 1 ULP). On the other hand, (9.0 / 1234 + 0.5) * 1234 is equal to 9.0 + 617, so the end result is sometimes exact and sometimes is not. How can you know beforehand which one is the case in your specific case? Generally, you can't, any arithmetic operation can potentially give you 1 ULP of error, and it can (and likely, will) slowly accumulate.
Also, please don't comment how nobody has a use for "f(x) = (x / 1234 + 0.5) * 1234": there are all kinds of queer computations people do in floating point, and for most of them, figuring out the exactness of the end result requires an absurd amount of applied numerical analysis, doing which would undermine most of the "just let the computer crunch the numbers" point of doing this computation on a computer.
It's a classic, so let's take vector normalization as an example. Topologically, you're ripping a hole in the space, and that's causing your issues. It manifests as NaN for length-zero vectors, weird precision issues too close to zero, etc, but no matter what you employ to try to fix it you're never going to have a good time squishing N-D space onto the surface of an N-D sphere if you need it to be continuous.
Some common subroutines where I see this:
1. You want to know the average direction of a bunch of objects and thus have to normalize each vector contributing to that average. Solution 1: That's not what you want almost ever. In any of the sciences, or anything loosely approximating the real world, you want to average the un-normalized vectors 99.999% of the time. Solution 2: Maybe you really do need directions for some reason (e.g., tracking where birds are looking in a game). Then don't rely on vectors for your in-band signaling. Explicitly track direction and magnitude separately and observe the magic of never having direction-related precision errors.
2. You're doing some sort of lighting normalization and need to compute something involving areas of potentially near-degenerate triangles, dividing by those values to weight contributions appropriately. Solution: Same as above, this is kind of like an average of averages problem. It can make fuzzy, intuitive sense, but you'll get better results if you do your summing and averaging in an un-normalized space. If you really do need surface normals, store those explicitly and separate from magnitude.
3. You're doing some sort of ML voodoo to try to get better empirical results via some vague appeal to vanishing gradients or whatever. Solution: The core property you want is a somewhat strange constraint on your layer's Jacobian matrix, and outside of like two papers nobody is willing to put up with the code complexity or runtime costs, even when they recognize it as the right thing to do. Everything you're doing is a hack anyway, so make your normalization term x/(|x|+eps) with eps > 0 rather than equal to zero like normal. Choose eps much smaller than most of the vectors you're normalizing this way and much bigger than zero. Something like 1e-3, 1e-20, and 1e-150 should be fine for f16, f32, and f64. You don't have to tune because it's a pretty weak constraint on the model, and it's able to learn around it.
> It's OK to compare floating-points for equality
Unless you use any compiler I write. As a purist, '==' for floats and doubles will be undefined. And by undefined, I mean unrecoverable. And by unrecoverable, I mean they may be able to extract most of the fragments of silicon, aluminum, and dynamic island, and "put you back together again". But you will "think different".
Unless you are comparing references. References are ok.
> The title of this post is an intentional clickbait.
Unfortunately, that's where the honesty ends.
> It's NOT OK to compare floating-points using epsilons.
> So, are epsilons good or bad? Usually bad, but sometimes okay.
So which is it? Emphatically NOT OK, or sometimes okay?
> that's how maths works
Wait is British "maths" a singular noun or is this a typo? I was willing to go along with it if it was plural, but I have to draw the line here.
Anything else is basically a nightmare to however has to maintain the code in the future.
Also, good luck with e.g. checking if points are aligned to a grid or the like without introducing a concept of epsilon _somewhere_.
This is especially true in cases where the number comes from some kind of input (user controls, sensor reading, etc.) or random number generation.