Provide round_up and round_down for converting double to float

[aprokop]

In some situations, it is helpful when we can control which way the rounding of a number is done. For example, right now, given data in double precision (like, points) we simply do float() conversion to the single precision. This is subject to round off errors and may lead to false negatives during the intersection tests. If we want to be conservative (i.e., only allow false positives but do not allow false negatives), the data should be rounded down for the minimum corner, and rounded up for the maximum corner, as to include the original data.

One of the ways to achieve that is using nextafterf functions. I'm thinking something like this:

#include <cmath>
float round_up(double x) {
  float f = x;
  if (double(f) < x)
    x = std::nextafterf(x, FLT_MAX);
}
float round_down(double x) {
  float f = x;
  if (double(f) > x)
    x = std::nextafterf(x, FLT_MIN);
}

(This may need to be changed slightly to account for nan, inf, etc).

Some papers do it informally (i.e., without using nextafterf) by rounding a positive number up by one ulp through multiplying it with 1 + 2^(-23) and down by multiplying with 1 - 2^(-23) (see pages 73-74 in Woop et al paper).

There are also some other situations where we may want to use nextafterf even for floats to account for losing precision using floating point operations. For example, this may help #560 (see this comment).

[dalg24]

Did you check whether it actually helps?
First thing that comes to mind is std::nextafterf() is not portable and it the function is currently not provided in <Kokkos_MathematicalFunctions.hpp> but if we act fast (I mean NOW) I could add it before we release 3.6
If you want it please open an issue and ask for it.

[aprokop]

Did you check whether it actually helps?

No. I do think that it should be possible to design an example producing false negative. For example, you could set use a double point (e, 0, 0), where e is the first non-zero double. It will be rounded to 0 in floats. If you use intersection with a sphere (0, 1) that does not include the surface itself (<r as opposed to <=r we do know), then this will result in false negative.

First thing that comes to mind is std::nextafterf() is not portable

Is it? I was not sure.

it the function is currently not provided in <Kokkos_MathematicalFunctions.hpp> but if we act fast (I mean NOW) I could add it before we release 3.6
If you want it please open an issue and ask for it.

I'm not sure there is a rush on this. This is simply to start the discussion. I'm not ready to have a detailed analysis of exactly how it may work, and whether it absolutely must be nextafterf, or whether this could be done differently. It's just this is the first time I stumbled upon a way to implement that, as I was not aware of it before.

[dalg24]

It means that if it turns out it is useful we will need workarounds in Arborx until we require Kokkos 3.7

[aprokop]

It means that if it turns out it is useful we will need workarounds in Arborx until we require Kokkos 3.7

I think that's fine. I am not ready to devote enough time to it prior to 3.6 anyway. If anybody wants to take a crack, go ahead.

[dalg24]

Kokkos::nextafter() is available since Kokkos 3.7
Did you want to keep this issue open to track possibly investigating controlling how fp rounding is done?
If so, maybe update the title line. Otherwise I suppose we can close this.

[aprokop]

The main issue issue (implicit conversion of double to float in ArborX::Point) is still there at the moment, but will be gone once we switch to multi-dimensional point.

I still think it is a good idea if ArborX provides these. One use case I can see if someone write an indexable getter that returns some double data, while explicitly specifying that the bounding volume is based on floats.