[BUG] cudf::round with HALF_UP mode produces non-deterministic output

[ttnghia]

Reproducible code:

using f_wrapper = cudf::test::fixed_width_column_wrapper<double>;
auto const input = f_wrapper{1.95, 2.95, 3.95, 4.95, 5.95, 6.95, 7.95, 8.95, 9.95};

auto const result = cudf::round(input, 1, cudf::rounding_method::HALF_UP);
cudf::test::print(*result);

Output:

2,3,4,5,6,7,8,8.9000000000000004,9.9000000000000004

I digged into the code and see that the rounding operation is indeed doesn't do round up properly:

template <typename T>
struct half_up_positive {
  T n;
  template <typename U = T, std::enable_if_t<cudf::is_floating_point<U>()>* = nullptr>
  __device__ U operator()(U e)
  {
    T integer_part;
    T const fractional_part = generic_modf(e, &integer_part);
    return integer_part + generic_round(fractional_part * n) / n;
  }

When debugging generic_round, given the input I see that:

input:    9.5000000000, output:   10.0000000000
input:    9.5000000000, output:   10.0000000000
input:    9.5000000000, output:   10.0000000000
input:    9.5000000000, output:   10.0000000000
input:    9.5000000000, output:   10.0000000000
input:    9.5000000000, output:   10.0000000000
input:    9.5000000000, output:   10.0000000000
input:    9.5000000000, output:    9.0000000000
input:    9.5000000000, output:    9.0000000000

In other word, generic_round can produce different outputs due to very small round-off error.