[BUG]: `DeviceReduce::ReduceByKey` is not run-to-run deterministic for fp64 sums over multi-tile runs

### Is this a duplicate?

- [x] I confirmed there appear to be no [duplicate issues](https://github.com/NVIDIA/cccl/issues) for this bug and that I agree to the [Code of Conduct](CODE_OF_CONDUCT.md)

### Type of Bug

Something else

### Component

CUB

### Describe the bug

## Summary

`cub::DeviceReduce::ReduceByKey` documents:

> - Provides "run-to-run" determinism for pseudo-associative reduction
> (e.g., addition of floating point types) on the same GPU device.
> However, results for pseudo-associative reduction may be inconsistent
> from one device to a another device of a different compute-capability
> because CUB can employ different tile-sizing for different architectures.

([v3.1.4](https://github.com/NVIDIA/cccl/blob/v3.1.4/cub/cub/device/device_reduce.cuh#L1645-L1650),
[main @ 0d919ed](https://github.com/NVIDIA/cccl/blob/0d919ed21f7d4366956d7f56277c2059986dd432/cub/cub/device/device_reduce.cuh#L2390-L2394)
— same sentence in the `ReduceByKey` doc block.)

Empirically this does not hold when a run of equal keys spans several
tiles: repeated invocations of the **same** `ReduceByKey` call on the
**same device** with **byte-identical device input** and the **same
pre-allocated temp storage** return the long run's `float64` aggregate
with **5–6 distinct bit patterns across 200 repetitions**, spread over the
last ~3 ulps. Short (single-tile) runs are bit-stable, and the same
experiment with integer-valued doubles (whose sums are exact under any
association order) is bit-stable — so the variation is a run-to-run
change in the **floating-point combination order** (presumably the
cross-tile partial combination in the decoupled look-back), not memory
corruption or a harness artifact.

Reproduced on two architectures with the same CUB version, so it is a
property of the implementation, not of a particular GPU:

| GPU | SM | Driver | Runtime | `CUB_VERSION` | Result (experiment A) |
|---|---|---|---|---|---|
| NVIDIA T600 Laptop GPU | sm_75 | 13010 | 13010 | 300104 (3.1.4) | 6 distinct bit patterns / 200 reps |
| NVIDIA H100 80GB HBM3 | sm_90 | 13020 | 12080¹ | 300104 (3.1.4) | 5 distinct bit patterns / 200 reps |

¹ the H100 host links a CUDA 12.8 runtime while compiling against a
CCCL 3.1.4 include tree (NVHPC toolchain); the T600 datapoint is a pure
CUDA 13.1 stack, so the finding does not depend on the mixed toolkit.

### How to Reproduce

## Steps to reproduce

Compile the self-contained program below (no dependencies beyond the
toolkit's bundled CCCL) and run it:

```bash
nvcc -O2 -arch=native \
-I$CUDA_HOME/targets/x86_64-linux/include/cccl \
repro.cu -o repro
./repro
```

The program runs three experiments over 65536 `(uint64 key, double value)`
pairs, 200 repetitions each, and prints a census of the bit patterns
observed for the first run's aggregate:

- **A** — one 50 000-item run of key `0` (values: mixed fp64 and zeros),
followed by 64-item runs. **Expected: 1 bit pattern. Actual: 5–6.**
- **B** — identical keys, values replaced by small integers stored as
fp64. Integer sums below 2^53 are exact under any association order, so
any flip here would indicate corruption rather than reordering.
**Actual: 1 bit pattern, equal to the sequential host sum.**
- **C** — the same fp64 values with every run 64 items (single tile).
**Actual: 1 bit pattern.**

Input buffers are uploaded once and checksummed before/after (reported
`UNMODIFIED`); the output buffer is poisoned before every repetition;
temp storage is allocated once and its size printed.

repro.cu

```cpp
// Reproducer: cub::DeviceReduce::ReduceByKey is not run-to-run
// deterministic for float64 addition over a key run spanning many tiles,
// despite the documented guarantee ("Provides 'run-to-run' determinism for
// pseudo-associative reduction (e.g., addition of floating point types) on
// the same GPU device" -- cub/device/device_reduce.cuh, ReduceByKey doc
// block; present in v3.1.4 and on main).
//
// Observed (NVIDIA T600 Laptop GPU, sm_75, CUDA 13.1, CUB 3.1.4 / 300104):
// experiment A returns 5 distinct bit patterns for the long run's aggregate
// across 200 identical invocations, spread over the last ~3 ulps.
//
// Protocol:
// - keys/values are uploaded to the device ONCE and never touched again
// (verified by checksumming the input buffers before and after);
// - temp storage is allocated ONCE (constant size, printed);
// - the SAME ReduceByKey call runs REPS times on the default stream; the
// first (long) run's aggregate is compared bit-for-bit across reps.
//
// Controls (same code path, only the data changes):
// A) fp64 values, one 50k-item key run -> flips if the internal
// combination order varies between invocations;
// B) same keys, values = small integers stored as fp64. Integer sums
// below 2^53 are EXACT in any association order, so if A's flips came
// from memory corruption or a harness bug, B would flip too; if they
// come from reordered floating-point combination, B is bit-stable;
// C) fp64 values but every run fits one tile (64 items) -> bit-stable,
// isolating the nondeterminism to the cross-tile combination.
//
// Build (no dependencies beyond the CUDA toolkit's bundled CCCL):
// nvcc -O2 -arch=native \
// -I$CUDA_HOME/targets/x86_64-linux/include/cccl \
// cub_reducebykey_determinism.cu -o cub_reducebykey_determinism
// ./cub_reducebykey_determinism

#include
#include

#include
#include
#include
#include
#include
#include

#define CK(x) \
if ((x) != cudaSuccess) { \
printf("cuda error %s @%d\n", cudaGetErrorString(x), __LINE__); \
return 1; \
}

struct Sum {
__host__ __device__ double operator()(double a, double b) const {
return a + b;
}
};

static uint64_t fnv1a(const void* p, size_t n) {
const auto* b = static_cast(p);
uint64_t h = 1469598103934665603ULL;
for (size_t i = 0; i < n; ++i) h = (h ^ b[i]) * 1099511628211ULL;
return h;
}

static int experiment(const char* name, const std::vector& keys,
const std::vector& vals, int reps) {
const size_t n = keys.size();
uint64_t *d_keys, *d_ukeys;
double *d_vals, *d_aggs;
size_t* d_nruns;
CK(cudaMalloc(&d_keys, n * 8));
CK(cudaMalloc(&d_ukeys, n * 8));
CK(cudaMalloc(&d_vals, n * 8));
CK(cudaMalloc(&d_aggs, n * 8));
CK(cudaMalloc(&d_nruns, 8));
CK(cudaMemcpy(d_keys, keys.data(), n * 8, cudaMemcpyHostToDevice));
CK(cudaMemcpy(d_vals, vals.data(), n * 8, cudaMemcpyHostToDevice));
const uint64_t in_hash =
fnv1a(keys.data(), n * 8) ^ fnv1a(vals.data(), n * 8);

size_t temp_bytes = 0;
cub::DeviceReduce::ReduceByKey(nullptr, temp_bytes, d_keys, d_ukeys,
d_vals, d_aggs, d_nruns, Sum{}, n);
void* d_temp;
CK(cudaMalloc(&d_temp, temp_bytes));

std::map observed; // aggregate[0] bit pattern -> count
size_t nruns_first = 0;
for (int rep = 0; rep < reps; ++rep) {
CK(cudaMemset(d_aggs, 0xEE, n * 8)); // poison: must be overwritten
cub::DeviceReduce::ReduceByKey(d_temp, temp_bytes, d_keys, d_ukeys,
d_vals, d_aggs, d_nruns, Sum{}, n);
uint64_t bits = 0;
size_t nruns = 0;
CK(cudaMemcpy(&bits, d_aggs, 8, cudaMemcpyDeviceToHost));
CK(cudaMemcpy(&nruns, d_nruns, 8, cudaMemcpyDeviceToHost));
if (rep == 0) nruns_first = nruns;
if (nruns != nruns_first) printf(" !! run count changed\n");
++observed[bits];
}

// Re-download the input to prove it was never modified.
std::vector keys2(n);
std::vector vals2(n);
CK(cudaMemcpy(keys2.data(), d_keys, n * 8, cudaMemcpyDeviceToHost));
CK(cudaMemcpy(vals2.data(), d_vals, n * 8, cudaMemcpyDeviceToHost));
const uint64_t out_hash =
fnv1a(keys2.data(), n * 8) ^ fnv1a(vals2.data(), n * 8);

// Host reference: strict left-to-right sum of run 0.
double seq = 0;
for (size_t i = 0; i < n && keys[i] == keys[0]; ++i) seq += vals[i];

printf("%s: n=%zu runs=%zu temp=%zu B, input %s, %d reps -> %zu "
"distinct bit patterns for aggregate[0]\n",
name, n, nruns_first, temp_bytes,
in_hash == out_hash ? "UNMODIFIED" : "MODIFIED(!)", reps,
observed.size());
for (const auto& [bits, count] : observed) {
double v;
memcpy(&v, &bits, 8);
printf(" %.17g (0x%016llx) x%d%s\n", v,
static_cast(bits), count,
memcmp(&v, &seq, 8) == 0 ? " == sequential host sum" : "");
}
cudaFree(d_keys);
cudaFree(d_ukeys);
cudaFree(d_vals);
cudaFree(d_aggs);
cudaFree(d_nruns);
cudaFree(d_temp);
return 0;
}

int main() {
int dev = 0;
cudaDeviceProp prop{};
CK(cudaGetDeviceProperties(&prop, dev));
int drv = 0, rt = 0;
cudaDriverGetVersion(&drv);
cudaRuntimeGetVersion(&rt);
printf("GPU: %s (sm_%d%d), driver %d, runtime %d, CUB %d\n\n", prop.name,
prop.major, prop.minor, drv, rt, CUB_VERSION);

const size_t n = 65536;
std::mt19937_64 rng(42);
std::uniform_real_distribution dist(-1.0, 1.0);

// A: one 50k-item run of key 0 (fp64 values), then short runs.
std::vector keys(n);
std::vector vals_fp(n), vals_int(n);
for (size_t i = 0; i < n; ++i) {
keys[i] = (i < 50000) ? 0 : 1 + (i - 50000) / 64;
vals_fp[i] = (i % 7 == 0) ? dist(rng) * 1e-2 : 0.0;
vals_int[i] = static_cast(static_cast(rng() % 1000));
}
if (experiment("A) fp64, 50k-item run ", keys, vals_fp, 200)) return 1;

// B: identical keys, integer-valued doubles (any order sums exactly).
if (experiment("B) int-as-fp64, same run", keys, vals_int, 200)) return 1;

// C: fp64 values, all runs short (64 items, single-tile).
std::vector keys_short(n);
for (size_t i = 0; i < n; ++i) keys_short[i] = i / 64;
if (experiment("C) fp64, 64-item runs ", keys_short, vals_fp, 200))
return 1;
return 0;
}
```

### Expected behavior

## Observed output

T600 (sm_75, CUDA 13.1):

```
GPU: NVIDIA T600 Laptop GPU (sm_75), driver 13010, runtime 13010, CUB 300104

A) fp64, 50k-item run : n=65536 runs=244 temp=5119 B, input UNMODIFIED, 200 reps -> 6 distinct bit patterns for aggregate[0]
-0.042204912555673282 (0xbfa59be1de5076d2) x10
-0.042204912555673296 (0xbfa59be1de5076d4) x2
-0.04220491255567331 (0xbfa59be1de5076d6) x53
-0.042204912555673324 (0xbfa59be1de5076d8) x84
-0.042204912555673338 (0xbfa59be1de5076da) x5
-0.042204912555673352 (0xbfa59be1de5076dc) x46
B) int-as-fp64, same run: n=65536 runs=244 temp=5119 B, input UNMODIFIED, 200 reps -> 1 distinct bit patterns for aggregate[0]
24857026 (0x4177b49c20000000) x200 == sequential host sum
C) fp64, 64-item runs : n=65536 runs=1024 temp=5119 B, input UNMODIFIED, 200 reps -> 1 distinct bit patterns for aggregate[0]
0.024264975624164437 (0x3f98d8eaf2f9e148) x200
```

H100 (sm_90):

```
GPU: NVIDIA H100 80GB HBM3 (sm_90), driver 13020, runtime 12080, CUB 300104

A) fp64, 50k-item run : n=65536 runs=244 temp=5119 B, input UNMODIFIED, 200 reps -> 5 distinct bit patterns for aggregate[0]
-0.042204912555673296 (0xbfa59be1de5076d4) x6
-0.04220491255567331 (0xbfa59be1de5076d6) x85
-0.042204912555673324 (0xbfa59be1de5076d8) x26
-0.042204912555673338 (0xbfa59be1de5076da) x11
-0.042204912555673352 (0xbfa59be1de5076dc) x72
B) int-as-fp64, same run: n=65536 runs=244 temp=5119 B, input UNMODIFIED, 200 reps -> 1 distinct bit patterns for aggregate[0]
24857026 (0x4177b49c20000000) x200 == sequential host sum
C) fp64, 64-item runs : n=65536 runs=1024 temp=5119 B, input UNMODIFIED, 200 reps -> 1 distinct bit patterns for aggregate[0]
0.024264975624164437 (0x3f98d8eaf2f9e148) x200
```

## Expected behavior

Per the documentation, identical invocations on the same device should
return bit-identical aggregates, including for floating-point addition
("pseudo-associative reduction").

### Reproduction link

_No response_

### Operating System

_No response_

### nvidia-smi output

_No response_

### NVCC version

_No response_

Run the supplied repro.cu with the documented nvcc command to confirm the varying fp64 bit patterns across repeated ReduceByKey calls. Then inspect cub/device/device_reduce.cuh, focusing on the ReduceByKey implementation and its cross-tile reduction behavior. Done means experiment A is bit-stable across repetitions while the documented same-device determinism guarantee remains valid.

Written by the indexing model from the issue text.