fix(arange): correct CPU length calc + GPU step scaling (#1076)#1077
Conversation
…extent (#1076) Problem: the arange(start, end, step) front-end computed the element count as integer truncation of (end-start)/step plus a `fmod(...) > 1.0e-14` bump. That was neither a reliable ceil nor a stable half-open [start, end) test: `fmod(1.0, 0.1) ~= 0.1` (not 0) so familiar decimal steps included the endpoint, and the fixed 1e-14 absolute threshold miscounted at small scales (arange(0, 1e-15, 2e-15) computed 0 elements and threw instead of yielding the single value 0). Empty and direction-mismatched ranges threw, even though zero-extent tensors are supported. The same count feeds both CPU and GPU kernels. Fix (per @ianmccul's #1071 review, folded into #1076): compute the half-open count as `Nelem = ceil((end - start) / step)`, clamped to 0. A non-positive count is an empty/direction-mismatched range and returns a zero-extent tensor rather than throwing. Added guards for zero step, non-finite start/end/step, and overflow of the double->uint64 count. The arange(Nelem) wrapper is made consistent: Nelem == 0 yields an empty tensor, negative still errors. Behavior change (flagged for review): endpoint handling under the half-open contract, small-scale counts, and empty range -> zero-extent instead of throw. Testing: new tests/arange_test.cpp with independent hand-computed values -- half-open integer/fractional/negative-step, the 1e-15 small-scale case, ULP boundaries via nextafter(end, +/-inf), zero-extent empties, and the step/finite guards. Bug-targeting cases were confirmed to fail on the pre-fix code. Updated the one existing test that asserted the old throw-on-direction-mismatch to expect a zero-extent result. Full CPU suite: 1801 passed, 0 failed. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…1070/#1076) Problem: the generic real-dtype cuSetArange_kernel dropped the parentheses in the value expression -- `start + step * blockIdx.x * blockDim.x + threadIdx.x` -- so the intra-block offset threadIdx.x was added un-scaled by step, and for a single block the values collapsed to `start + threadIdx.x`. arange(10, 40, 10) on Device.cuda produced [10, 11, 12] instead of [10, 20, 30] (#1070). The complex overloads already had the parentheses. Separately, the linear index `blockIdx.x * blockDim.x + threadIdx.x` was 32-bit and overflows past 2^32 elements. Fix: compute the linear index once as 64-bit (uint64_t) in all real and complex kernels, and use `start + step * idx` so the full index is scaled by step. Testing: new tests/gpu/arange_test.cpp with independent hand-computed values (not compared against the CPU path): non-unit step across Int64/Int32/Double/ Float, a >512-element multi-block case exercising the block-boundary index (value[i] == 3*i), fractional step, and a zero-extent empty range. GoogleTest names avoid underscores. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
![gemini-code-assist[bot]](https://avatars.githubusercontent.com/in/956858?s=60&v=4){kind=small}
There was a problem hiding this comment.
Code Review
This pull request refactors the arange implementation for both CPU and GPU to support zero-extent tensors (empty ranges) and fixes several bugs, including 32-bit index overflow in GPU kernels and incorrect step scaling. It also adds comprehensive unit tests for arange on CPU and GPU. The review feedback highlights critical issues regarding floating-point precision: first, std::ceil can include the exclusive endpoint due to roundoff errors (e.g., arange(0.5, 0.8, 0.1)), which can be resolved by subtracting a small epsilon; second, casting std::numeric_limits<cytnx_uint64>::max() to double rounds up to 2^64, leading to potential undefined behavior when checking for overflow, which can be fixed by using an explicit threshold. Additionally, a test case for non-binary fractional steps is suggested to verify floating-point boundary behavior.
Important
The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.
| if (count > 0) { | ||
| const cytnx_double nelem = std::ceil(count); |
There was a problem hiding this comment.
Using pure std::ceil(count) can lead to an extra element (including the exclusive endpoint) due to floating-point roundoff errors. For example, with arange(0.5, 0.8, 0.1), count is computed as 3.0000000000000004 instead of exactly 3.0. Taking std::ceil of this rounds it up to 4.0, causing the exclusive endpoint 0.8 to be incorrectly included in the output.
To prevent this, we should subtract a small epsilon that scales with count before applying std::ceil.
if (count > 0) {
const cytnx_double eps = count * std::numeric_limits<cytnx_double>::epsilon() * 100.0;
const cytnx_double nelem = std::ceil(count - eps);There was a problem hiding this comment.
not sure, I don't think numpy implementation does that adjustment.
$ python3
Python 3.12.3 (main, Jun 19 2026, 12:46:00) [GCC 13.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy as np
>>> print(np.arange(0.5, 0.8, 0.1))
[0.5 0.6 0.7 0.8]
| cytnx_error_msg( | ||
| nelem > static_cast<cytnx_double>(std::numeric_limits<cytnx_uint64>::max()), | ||
| "[ERROR] arange(start=%f,end=%f,step=%f): the requested number of elements exceeds the " | ||
| "maximum representable size.\n", | ||
| start, end, step); |
There was a problem hiding this comment.
std::numeric_limits<cytnx_uint64>::max() is 2^64 - 1. When cast to double (cytnx_double), it rounds up to 2^64 (18446744073709551616.0) because double only has 53 bits of precision.
As a result, the condition nelem > static_cast<cytnx_double>(max) is equivalent to nelem > 2^64. If nelem is exactly 2^64, this condition evaluates to false, bypassing the error check. However, casting 2^64 to uint64_t on the next line is undefined behavior in C++ because it exceeds the maximum representable value of uint64_t (2^64 - 1).
Using nelem >= 18446744073709551616.0 (which is exactly 2^64) avoids this undefined behavior.
| cytnx_error_msg( | |
| nelem > static_cast<cytnx_double>(std::numeric_limits<cytnx_uint64>::max()), | |
| "[ERROR] arange(start=%f,end=%f,step=%f): the requested number of elements exceeds the " | |
| "maximum representable size.\n", | |
| start, end, step); | |
| cytnx_error_msg( | |
| nelem >= 18446744073709551616.0, | |
| "[ERROR] arange(start=%f,end=%f,step=%f): the requested number of elements exceeds the " | |
| "maximum representable size.\n", | |
| start, end, step); |
There was a problem hiding this comment.
If nelem >= 18446744073709551616 I think you've got bigger problems than overflowing a uint64_t...
| EXPECT_DOUBLE_EQ(t.storage().at<cytnx::cytnx_double>(1), 0.25); | ||
| EXPECT_DOUBLE_EQ(t.storage().at<cytnx::cytnx_double>(2), 0.5); | ||
| EXPECT_DOUBLE_EQ(t.storage().at<cytnx::cytnx_double>(3), 0.75); | ||
| } |
There was a problem hiding this comment.
Adding a test case for non-binary fractional steps (like arange(0.5, 0.8, 0.1)) ensures that floating-point inaccuracies near integer boundaries do not cause the exclusive endpoint to be incorrectly included.
}
TEST(Arange, FractionalStepNonBinary) {
// [0.5, 0.6, 0.7]; 0.8 is the exclusive endpoint.
// This tests that floating-point inaccuracies near integer boundaries
// do not cause the exclusive endpoint to be incorrectly included.
auto t = cytnx::arange(0.5, 0.8, 0.1, Type.Double);
ASSERT_EQ(t.shape()[0], 3u);
EXPECT_DOUBLE_EQ(t.storage().at<cytnx::cytnx_double>(0), 0.5);
EXPECT_DOUBLE_EQ(t.storage().at<cytnx::cytnx_double>(1), 0.6);
EXPECT_DOUBLE_EQ(t.storage().at<cytnx::cytnx_double>(2), 0.7);
}
![chatgpt-codex-connector[bot]](https://avatars.githubusercontent.com/in/1144995?s=60&v=4){kind=small}
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: a1b3be3faf
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| const cytnx_double count = (end - start) / step; | ||
| cytnx_uint64 Nelem = 0; | ||
| if (count > 0) { | ||
| const cytnx_double nelem = std::ceil(count); |
There was a problem hiding this comment.
For floating ranges that should land exactly on the endpoint but whose (end - start) / step rounds slightly above an integer, this ceil adds an extra element. For example, arange(0.2, 0.8, 0.2) computes a count of about 3.0000000000000004, so it returns four values and includes 0.8, violating the documented half-open [start, end) contract and changing the tensor shape; this can break reshapes or downstream numerical code. Consider tolerating near-integer counts or validating the generated last value against the endpoint.
Useful? React with 👍 / 👎.
There was a problem hiding this comment.
True, but numpy does this as well. https://numpy.org/doc/stable/reference/generated/numpy.arange.html
consistency with numpy is probably more useful here.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #1077 +/- ##
=======================================
Coverage 72.67% 72.67%
=======================================
Files 226 226
Lines 28177 28177
Branches 71 71
=======================================
Hits 20479 20479
Misses 7677 7677
Partials 21 21
Flags with carried forward coverage won't be shown. Click here to find out more.
Continue to review full report in Codecov by Harness.
🚀 New features to boost your workflow:
|
ianmccul
left a comment
ianmccul
left a comment
There was a problem hiding this comment.
Add some pytests that compare against numpy, and document that we follow the numpy implementation, in particular that the count of elements might include (or even slightly beyond?) the final value in some cases, contrary to the specification but consistent with https://numpy.org/doc/stable/reference/generated/numpy.arange.html
|
Maybe we should remove gemini-code from reviewing .. |
2 participants
Problem
arangecorrectness needed a coordinated CPU + GPU fix (#1076). Two independent bugs, sharing one broken length calculation:src/Generator.cpp::arange) — computed the element count as integer truncation of(end-start)/stepplus afmod(...) > 1.0e-14bump. Per @ianmccul's fix(utils): honor non-unit step in GPU arange for real dtypes (#1070) #1071 review this was "neither a reliable ceil nor a stable half-open test":fmod(1.0, 0.1) ≈ 0.1(not 0) so decimal steps included the endpoint; the fixed1e-14absolute threshold miscounted at small scales (arange(0, 1e-15, 2e-15)computed 0 elements and threw instead of yielding[0]); and empty/direction-mismatched ranges threw despite zero-extent tensors being supported.cuSetArange_gpu.cu) — the generic real-dtype kernel dropped the parentheses (start + step * blk * dim + thread), so the per-thread offset was added un-scaled bystep;arange(10, 40, 10)onDevice.cudagave[10, 11, 12]instead of[10, 20, 30](GPU arange ignores non-unit step for real dtypes (missing parentheses in cuSetArange_kernel) #1070). The linear index was also 32-bit and overflows past 2^32 elements.The front-end computes the count once and feeds it to both backends, so the CPU length bug propagated to the GPU path too.
Fix
Nelem = ceil((end - start) / step), clamped to 0 for an empty/direction-mismatched range (returns a zero-extent tensor). Guards added for zero step, non-finitestart/end/step, and overflow of thedouble → uint64count.arange(Nelem)made consistent (0→ empty, negative → error).start + step * idxin all real and complex kernels.This changes user-visible numeric behavior, as #1076 intended:
[start, end)endpoint handling now via ceil-of-count;arange(0, 1e-15, 2e-15)→[0]);One existing test (
DenseUniTensorTest.arange_step_error) asserted the old throw-on-direction-mismatch; it's updated to expect a zero-extent result and renamed.Testing
tests/arange_test.cpp, new): independent hand-computed values — half-open integer/fractional/negative-step, the1e-15small-scale case, ULP boundaries vianextafter(end, ±inf), zero-extent empties, and the step/finite guards. Bug-targeting cases were confirmed to fail on the pre-fix code. Full CPU suite: 1801 passed, 0 failed.tests/gpu/arange_test.cpp, new): independent literal values (not the CPU path as oracle, per Ian) — non-unit step across Int64/Int32/Double/Float, a >512-element multi-block case exercising the block-boundary index (value[i] == 3*i), fractional step, zero-extent empty. Run locally on an RTX 4070 Ti SUPER (CUDA 13, sm_89). GPU tests don't run in CI (no GPU runner).Not in scope: the remaining #1076 items (e.g. the
<cuda/std/complex>install/PUBLIC-propagation for GPU consumers) belong to #1067.Closes #1076. Fixes the GPU half of #1070.