Fix numerical cancellation in RectLattice::distance for large pitch values

[MatteoZammataro]

This PR fixes a segmentation fault in RectLattice::distance that was triggered by large lattice pitch values (see #3852). This bug occurred due to floating-point cancellation in the boundary-crossing detection logic.

As a simple fix, I suggest computing each axis distance (dx, dy, dz) independently and updating lattice_trans as the minimum is determined. This should eliminate the risk of floating-point cancellation.

Credits to @MatteoF29 for the help.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[GuySten]

This problem is very similar to #3792.
I suggest we refactor that solution into an isclose function and use it here and there.

[MatteoZammataro]

Good idea, thanks for the suggestion! We could refactor the logic into a small utility like:

inline bool isclose(double a, double b,
                    double rel_tol = 1e-15,
                    double abs_tol = FP_PRECISION)
{
  return std::abs(a - b) <=
         std::max(rel_tol * std::max(std::abs(a), std::abs(b)), abs_tol);
}

The boundary-crossing detection in lattice.cpp would then become:

double dx = (u.x != 0.0) ? (x0 - x) / u.x : INFTY;
double dy = (u.y != 0.0) ? (y0 - y) / u.y : INFTY;
double dz = (u.z != 0.0) ? (z0 - z) / u.z : INFTY;

double d = std::min({dx, dy, dz});
d = std::max(d, FP_PRECISION);

if (isclose(d, dx)) lattice_trans[0] = std::copysign(1, u.x);
if (isclose(d, dy)) lattice_trans[1] = std::copysign(1, u.y);
if (is_3d_ && isclose(d, dz)) lattice_trans[2] = std::copysign(1, u.z);

and similarly in mesh.cpp, e.g.

if (isclose(r0, 0.0))
  return INFTY;

if (isclose(denominator, 0.0))
  return INFTY;

if (isclose(R, r0, RADIAL_MESH_TOL))
  return INFTY;

Does this look close to what you had in mind? If so, do you have any suggestion on where the isclose function should be placed?

[GuySten]

It is close to what I had in mind.
I thought the natural place for that function should be in src/math_functions.cpp.

[MatteoZammataro]

Done! I have a doubt about the relative tolerance that should be used in the lattice. At the moment, I used 1E-12, but it is a guess value. We could eventually consider adding a constant variable named something like "RECT_LATTICE_TOL" to be placed in constants.h.