Uint16 x Uint16 multiply is signed-overflow UB in Kron/Outer/Mul (CPU and GPU)

[yingjerkao]

Multiplying two Uint16 tensors is undefined behaviour for operands above ~46341, on both CPU and GPU. The results are currently correct on GCC/nvcc — the wrap happens to land on the right value — so this is a latent correctness hazard rather than a live wrong answer, and nothing in CI would detect it.

Surfaced while reviewing #1108 (Codex raised it; the numbers below are measured).

The construct

Three sites, same shape:

// src/backend/linalg_internal_cpu/Kron_internal.hpp:49
out[i] = Lin[x] * Rin[y];

// src/backend/linalg_internal_gpu/cuKron_internal.cuh
out[...] = TO(Lin[x]) * TO(Rin[y]);

// src/backend/linalg_internal_gpu/cuArithmeticDispatch.cuh:33   (op_code 1 = Mul)
return static_cast<TO>(lhs) * static_cast<TO>(rhs);

For Uint16 (x) Uint16, Type_class::type_promote gives TO = cytnx_uint16, so the casts change nothing. uint16_t * uint16_t undergoes integral promotion — int represents every uint16_t value, so both operands become int and the multiplication is int * int ([conv.prom]). Signed overflow is undefined ([expr.pre]/4). The threshold is sqrt(INT_MAX) ~ 46341.

The third site matters most: cuArithmeticDispatch.cuh is the shared typed dispatch that #1003 is consolidating GPU elementwise arithmetic onto, so this is on the path that work is actively touching, and it reaches Mul generally rather than just Kron.

Which dtype pairs are affected

Only Uint16 (x) Uint16. Everything else in the dtype set is safe, some of it narrowly:

operands	promoted TO	worst product	verdict
Uint16 (x) Uint16	uint16_t, promotes to int	65535^2 = 4,294,836,225	overflows INT_MAX
Int16 (x) Int16	int16_t, promotes to int	32768^2 = 1,073,741,824	safe
Int16 (x) Uint16	int32_t	32767 x 65535 = 2,147,385,345	safe, by 98,302
Uint32 (x) Uint32	uint32_t, no promotion to int	—	defined wraparound

Reproduction

Tensor a = zeros({1}, Type.Uint16), b = zeros({1}, Type.Uint16);
a.at<cytnx_uint16>({0}) = 65350;
b.at<cytnx_uint16>({0}) = 62109;          // 65350 * 62109 = 4,058,823,150 > INT_MAX

linalg::Mul(a, b);                         // CPU  -> 47598
linalg::Mul(a.to(Device.cuda), b.to(Device.cuda));   // GPU -> 47598
linalg::Kron(a, b);                        // CPU  -> 47598
linalg::Kron(a.to(Device.cuda), b.to(Device.cuda));  // GPU -> 47598

47598 is (65350 * 62109) mod 65536, i.e. the mathematically correct uint16 result. All four paths agree with it today: two's-complement wrap into int followed by truncation to uint16 coincides with the right answer.

So there is no wrong result to point at right now. What is missing is the guarantee — UB licenses the optimizer to assume the overflow cannot happen, and the behaviour is not contractual across compilers, flags, or a future vectorization of these loops.

Why CI cannot see it

-fsanitize=address is the only sanitizer configured (CMakeLists.txt:338); there is no -fsanitize=undefined in any preset. ASan does not detect signed-integer overflow, so no existing job would flag this.

The debug-* presets would report it immediately if UBSan were added — worth considering separately, since this class of bug is otherwise invisible to us.

Proposed fix

Compute through an unsigned type wide enough to hold the product, where overflow is defined wraparound. The truncated result is bit-identical to what is produced today, so this is a no-op on observable behaviour for every input:

// only where TO is a narrow integer that promotes to int
return static_cast<TO>(static_cast<unsigned int>(lhs) * static_cast<unsigned int>(rhs));

Applies to the Mul branch of ApplyGpuArithOp, the CPU Kron_internal multiply, and the GPU cuKron_general multiply. Add/Sub need nothing: uint16 + uint16 peaks at 131,070, comfortably inside int.

This changes arithmetic-kernel code, so per CLAUDE.md it wants explicit sign-off rather than being folded into an unrelated PR — hence this issue instead of a drive-by patch.

Already done

#1108 bounds Uint16 sweep operands to [0, 1000] (commit 8a059fe0), matching the range Uint32/Uint64 already use. That stops the test from depending on UB; it does not fix the kernels.

Worth noting how the test reached those magnitudes at all: GetRandRange in tests/gpu/gpu_test_tools.cpp gives the narrow types their full numeric_limits range (Uint16 -> [0, 65535], Int16 -> [-32768, 32767]) while the wider unsigned types use [0, 1000]. With the fixed seeds, the Kron/Outer sweeps generated overflowing pairs on every shape:

sweep case	max a	max b	worst product	overflowing pairs
Kron {2,3}x{3,2}	65472	62109	4,066,400,448	6 / 36
Kron {4,1}x{2,5}	65472	62109	4,066,400,448	12 / 40
Kron {1,1}x{3,3}	65350	62109	4,058,823,150	3 / 9
Outer {5}x{7}	55046	59022	3,248,925,012	5 / 35

Whether GetRandRange should give Uint16 the same [0, 1000] as the other unsigned dtypes is a separate question; it is shared by every GPU suite, so changing it shifts other tests' inputs and tolerance margins.

Refs #1003, #1108.

[IvanaGyro]

I doubt that we really have to fix the overflow problem. We can fix uint16_t * unit16_t, but not uint64_t * uint64_t. Maybe UB can just be a feature.