Rounding issue with the power function on integers

[magnatelee]

@ejmeitz recently reported this weird behavior that the power function on 32-bit integers don't produce correct results:

> LEGATE_TEST=1 LEGATE_CONFIG="--gpus 1 --auto-config 0" python
Python 3.13.7 [...]
>>> import cupynumeric as np
>>> a = np.array([117] * 10, dtype="int32")
>>> b = np.array([1] * 10, dtype="int32")
>>> print(a ** b)
[116 116 116 116 116 116 116 116 116 116]         # Oops

This happens only when the code goes through the GPU code path, and as it turns out, the correctness issue is due to GPU's rounding behavior: because the power function on integer values is implemented using the power function on doubles and type conversions, the float-to-integer conversions can cause correctness issues.

Unfortunately, there's no single correct rounding mode we can consistently use blindly for all values, here's a CUDA program that checks the correctness of roundtrip conversion (int32->double->int32) up to 1024. It'd report only one of two rounding modes (round-up vs. round-down) is correct for some of the values and it's not always the same one that is correct.

I can't think of a reasonable solution here other than implementing a power function on integer values that doesn't do type conversions. Here's a CuPy implementation we can copy from: https://github.com/cupy/cupy/blob/39dcf763dd339afcd49e51f4e5a0311370baa61c/cupy/_core/_routines_math.pyx#L983-L997

[ejmeitz]

I've also observed it with Int64 (on 95 and 147). I only see this through the Julia wrapper I have though which is on v25.05.

I cannot reproduce this on my local install of cupynumeric (Python not Julia). I have 24.11.02 installed with a local CUDA toolkit (v12.4)

[manopapad]

Assigning to @aschaffer

[ejmeitz]

Pretty sure this also affects the FLOOR_DIVIDE binary op.

[aschaffer]

cuda::std::round() fixes the problem, AFAICT.

#include <algorithm>
#include <cmath>
#include <iostream>
#include <iterator>
#include <vector>

#include <thrust/copy.h>
#include <thrust/device_vector.h>
#include <thrust/transform.h>

#include <cuda/std/cmath>

//

template <typename T, typename... Args, template <typename, typename...> class Vector>
void print_v(const Vector<T, Args...>& v, std::ostream& os)
{
  thrust::copy(v.begin(), v.end(), std::ostream_iterator<T>(os, ",")); 
  os << "\n";
}

// $ nvcc -c -std=c++17 --expt-extended-lambda std_pow.cu
// $ nvcc -std=c++17 --expt-extended-lambda std_pow.cu -o spw.exe

int main(void)
{
  int x = 117;
  int e = 1;

  int res = cuda::std::pow(x, e);
  std::cout << x << "^" << e << "=" << res << '\n';

  std::vector<int> h_v{2, 3, 5, 7, 11, 13, 17, 117};
  {
    thrust::device_vector<int> d_v(h_v);

    thrust::transform(
      d_v.begin(), d_v.end(), d_v.begin(), [e] __device__(int x) { return cuda::std::pow(x, e); });

    std::cout << "device vector (1.0):\n";
    print_v(d_v, std::cout);
  }

  {
    thrust::device_vector<int> d_v(h_v);

    thrust::transform(d_v.begin(), d_v.end(), d_v.begin(), [e] __device__(int x) {
      return cuda::std::round(cuda::std::pow(x, e));
    });

    std::cout << "device vector (1.1):\n";
    print_v(d_v, std::cout);
  }

  std::cout << "Done!\n";

  return 0;
}

[magnatelee]

cuda::std::round() fixes the problem, AFAICT.

I don't think that's the right fix though. Rounding can lead to other diverging behaviors for different input values. We should implement a power function for integer values to eliminate the source of inaccuracy.

[aschaffer]

A better approach will soon be cuda::ipow(int, int) but for that we need to migrate to CCCL 3.x or CTK 13.x. Which is in the works.

[manopapad]

We have migrated to CCCL 3.0 in Legate, @aschaffer could you move to the proper fix using ipow?

[aschaffer]

We need CTK 13.1+, CCCL 3.1+ for ipow().

[sbuneti]

Un-assigning old sprint as we wait for CTK 13.1 to make progress on this. Thanks