From 535d158809b23ac3df179f76b9a8b85bd05a3c43 Mon Sep 17 00:00:00 2001 From: Evgen Byelozorov Date: Wed, 12 Aug 2026 09:46:09 +0200 Subject: [PATCH 1/3] fix(arith): keep integer div exact and UB-free in ray_idiv_fn div routed integer operands through doubles, silently corrupting every result above 2^53 (div 9007199254740993 1 -> 9007199254740992) and tripping UBSan at q == 2^63 (div -9223372036854775807 -1) because the q > (double)INT64_MAX guard can never fire. Integer operands now divide in int64 space with a floor correction, matching the temporal mod path; the double path is kept only for float operands with a tightened guard. --- src/ops/arith.c | 24 ++++++++++++++++++------ test/rfl/integration/math.rfl | 27 +++++++++++++++++++++++++++ 2 files changed, 45 insertions(+), 6 deletions(-) diff --git a/src/ops/arith.c b/src/ops/arith.c index d22c8aba7..829d78b22 100644 --- a/src/ops/arith.c +++ b/src/ops/arith.c @@ -277,13 +277,25 @@ ray_t* ray_idiv_fn(ray_t* a, ray_t* b) { ray_type_name(a->type), ray_type_name(b->type)); if (RAY_ATOM_IS_NULL(a) || RAY_ATOM_IS_NULL(b)) return ray_typed_null(-RAY_I64); - double bv = as_f64(b); - if (bv == 0.0) - return ray_typed_null(-RAY_I64); - double q = floor(as_f64(a) / bv); - if (q < (double)INT64_MIN || q > (double)INT64_MAX) + if (is_float_op(a, b)) { + double bv = as_f64(b); + if (bv == 0.0) + return ray_typed_null(-RAY_I64); + double q = floor(as_f64(a) / bv); + if (q >= 9223372036854775808.0 || q < (double)INT64_MIN) + return ray_typed_null(-RAY_I64); + return make_i64((int64_t)q); + } + /* Integer div: stay in int64 space. The double round-trip above silently + * loses precision for magnitudes > 2^53 and is UB for q == 2^63. */ + int64_t bv = as_i64(b); + if (bv == 0) return ray_typed_null(-RAY_I64); - return make_i64((int64_t)q); + int64_t la = as_i64(a); + int64_t q = la / bv; + if (la % bv != 0 && ((la < 0) != (bv < 0))) + q--; /* floor toward -inf */ + return make_i64(q); } ray_t* ray_mod_fn(ray_t* a, ray_t* b) { diff --git a/test/rfl/integration/math.rfl b/test/rfl/integration/math.rfl index e0b668c7c..cb51fac79 100644 --- a/test/rfl/integration/math.rfl +++ b/test/rfl/integration/math.rfl @@ -2478,3 +2478,30 @@ ;; Test division with type conversion (/ [100 200 300] 3.0) -- [33.33 66.67 100.0] (% [100 200 300] 3.0) -- [1.0 2.0 0.0] + +;; ========== I64 FLOOR DIV EXACTNESS / UB REGRESSIONS ========== +;; Integer div must stay in int64 space: the old double round-trip lost +;; precision above 2^53 and hit UB at q == 2^63 (div -9223372036854775807 -1). +(div 9007199254740993 1) -- 9007199254740993 +(div 9007199254740993 3) -- 3002399751580331 +(div 18014398509481985 1) -- 18014398509481985 +(div 123456789012345678 1) -- 123456789012345678 +(div 9223372036854775807 1) -- 9223372036854775807 +(div 9223372036854775807 -1) -- -9223372036854775807 +(div -9223372036854775807 -1) -- 9223372036854775807 +(div -9223372036854775807 2) -- -4611686018427387904 +(div -9223372036854775807 -2) -- 4611686018427387903 +(div 9223372036854775807 2) -- 4611686018427387903 +;; vectors (1-element vectors route through the scalar path) +(div [9007199254740993 9007199254740993] 3) -- [3002399751580331 3002399751580331] +(div [123456789012345678 9007199254740993] 1) -- [123456789012345678 9007199254740993] +(div [9223372036854775807] -1) -- [-9223372036854775807] +(div [-9223372036854775807] -1) -- [9223372036854775807] +;; out-of-int64-range (float input) still clamps to null +(div 1e308 1) -- 0Nl +(div -1e308 1) -- 0Nl +(div [1e308 9.3e18] 1) -- [0Nl 0Nl] +(div 9223372036854775808.0 1) -- 0Nl +;; floor semantics preserved for int operands +(div [-7 7 -7] [2 -2 -2]) -- [-4 -4 3] +(div [7 8 9] 2) -- [3 4 4] From 8c3957c4cd21e7bd6144a51838cb95a75a71188b Mon Sep 17 00:00:00 2001 From: Evgen Byelozorov Date: Sun, 16 Aug 2026 18:01:09 +0200 Subject: [PATCH 2/3] fix(expr): keep compiled integer div exact --- src/ops/expr.c | 76 +++++++++++++++++++++++++---------- test/rfl/integration/math.rfl | 7 ++++ 2 files changed, 62 insertions(+), 21 deletions(-) diff --git a/src/ops/expr.c b/src/ops/expr.c index 08061f65c..8a3dc5d2a 100644 --- a/src/ops/expr.c +++ b/src/ops/expr.c @@ -42,6 +42,15 @@ static inline uint8_t truthy_f64ish(double v) { return (v == v && v != 0.0) ? 1 : 0; } +static inline bool floor_idiv_i64_checked(int64_t a, int64_t b, int64_t* out) { + if (b == 0 || (a == INT64_MIN && b == -1)) return false; + int64_t q = a / b; + int64_t rem = a % b; + if (rem != 0 && ((a < 0) != (b < 0))) q--; + *out = q; + return true; +} + static bool atom_to_numeric(ray_t* atom, double* out_f, int64_t* out_i, bool* out_is_f64) { if (!atom || !ray_is_atom(atom)) return false; switch (atom->type) { @@ -192,12 +201,10 @@ static bool eval_const_numeric_expr(ray_graph_t* g, ray_op_t* op, case OP_SUB: r = (int64_t)((uint64_t)li - (uint64_t)ri); break; case OP_MUL: r = (int64_t)((uint64_t)li * (uint64_t)ri); break; case OP_DIV: - if (ri==0) return false; - r = li/ri; if ((li^ri)<0 && r*ri!=li) r--; + if (!floor_idiv_i64_checked(li, ri, &r)) return false; break; case OP_IDIV: - if (ri==0) return false; - r = li/ri; if ((li^ri)<0 && r*ri!=li) r--; + if (!floor_idiv_i64_checked(li, ri, &r)) return false; break; case OP_MOD: if (ri==0) return false; @@ -583,11 +590,11 @@ static bool expr_null_capable(uint8_t op, int8_t dt, int8_t t1) { /* Task 7: I64 arithmetic and unary ops. * OP_ADD=20..OP_MOD=24, OP_MIN2=33, OP_MAX2=34 in the binary block. * OP_NEG=10, OP_ABS=11 in the unary block. - * OP_IDIV=49 is NOT contiguous with OP_ADD..OP_MOD and has no I64 plain - * case in exec_elementwise_binary, so it is excluded. */ + * OP_IDIV=49 is NOT contiguous with OP_ADD..OP_MOD, but has explicit I64 + * cases in both fused and fallback kernels. */ if (dt == RAY_I64 && t1 != RAY_F64 && ((op >= OP_ADD && op <= OP_MOD) || op == OP_MIN2 || op == OP_MAX2 || - op == OP_NEG || op == OP_ABS || op == OP_SIGNUM)) + op == OP_IDIV || op == OP_NEG || op == OP_ABS || op == OP_SIGNUM)) return true; if (dt == RAY_I64 && t1 == RAY_F64 && op == OP_SIGNUM) return true; @@ -1053,6 +1060,7 @@ static void expr_exec_binary(uint8_t opcode, uint8_t null_aware, int8_t dt, void case OP_SUB: for (int64_t j = 0; j < n; j++) d[j] = ray_f64_fin(a[j] - b[j]); break; case OP_MUL: for (int64_t j = 0; j < n; j++) d[j] = ray_f64_fin(a[j] * b[j]); break; case OP_DIV: for (int64_t j = 0; j < n; j++) d[j] = b[j] != 0.0 ? ray_f64_fin(a[j] / b[j]) : NULL_F64; break; + case OP_IDIV: for (int64_t j = 0; j < n; j++) d[j] = b[j] != 0.0 ? ray_f64_fin(floor(a[j] / b[j])) : NULL_F64; break; case OP_MOD: for (int64_t j = 0; j < n; j++) { if (b[j] == 0.0) { d[j] = NULL_F64; continue; } double m = fmod(a[j], b[j]); @@ -1081,12 +1089,11 @@ static void expr_exec_binary(uint8_t opcode, uint8_t null_aware, int8_t dt, void * b==0 (non-null) → NULL_I64 mirrors fallback's zero-divisor post-pass. * b==-1 && a==INT64_MIN → overflow: direct NULL_I64 (vs fallback: loop writes 0 * then zero-divisor pass marks null — same observable result). */ - case OP_DIV: for (int64_t j=0;jtype == -RAY_F64 || lhs->type == RAY_F64 || lhs->type == -RAY_F32 || lhs->type == RAY_F32)) ? l_f64 : (double)l_i64) #define RV_READ(i) (rp_f64 ? rp_f64[i] : rp_f32 ? (double)rp_f32[i] : rp_i64 ? (double)rp_i64[i] : rp_i32 ? (double)rp_i32[i] : rp_u32 ? (double)rp_u32[i] : rp_i16 ? (double)rp_i16[i] : rp_bool ? (double)rp_bool[i] : (r_scalar && (rhs->type == -RAY_F64 || rhs->type == RAY_F64 || rhs->type == -RAY_F32 || rhs->type == RAY_F32)) ? r_f64 : (double)r_i64) +#define LV_READ_I64(i) (lp_i64 ? lp_i64[i] : lp_i32 ? (int64_t)lp_i32[i] : lp_u32 ? (int64_t)lp_u32[i] : lp_i16 ? (int64_t)lp_i16[i] : lp_bool ? (int64_t)lp_bool[i] : l_i64) +#define RV_READ_I64(i) (rp_i64 ? rp_i64[i] : rp_i32 ? (int64_t)rp_i32[i] : rp_u32 ? (int64_t)rp_u32[i] : rp_i16 ? (int64_t)rp_i16[i] : rp_bool ? (int64_t)rp_bool[i] : r_i64) /* Compute once: is lhs/rhs integer-family (not float)? Used by BOOL path. */ int l_is_int = !(lp_f64 || lp_f32 || (l_scalar && @@ -3181,7 +3189,13 @@ static void binary_range(ray_op_t* op, int8_t out_type, case OP_SUB: for(int64_t i=0;iri?li:ri;}break; @@ -3195,7 +3209,13 @@ static void binary_range(ray_op_t* op, int8_t out_type, case OP_MUL: for(int64_t i=0;iri?li:ri;}break; @@ -3208,7 +3228,13 @@ static void binary_range(ray_op_t* op, int8_t out_type, case OP_SUB: for(int64_t i=0;iri?li:ri;}break; @@ -3221,7 +3247,13 @@ static void binary_range(ray_op_t* op, int8_t out_type, case OP_SUB: for(int64_t i=0;iri?li:ri;}break; @@ -3262,6 +3294,8 @@ static void binary_range(ray_op_t* op, int8_t out_type, } #undef LV_READ #undef RV_READ +#undef LV_READ_I64 +#undef RV_READ_I64 done: if (lsym_buf) ray_free_raw(lsym_buf); if (rsym_buf) ray_free_raw(rsym_buf); diff --git a/test/rfl/integration/math.rfl b/test/rfl/integration/math.rfl index cb51fac79..4567d38b1 100644 --- a/test/rfl/integration/math.rfl +++ b/test/rfl/integration/math.rfl @@ -2505,3 +2505,10 @@ ;; floor semantics preserved for int operands (div [-7 7 -7] [2 -2 -2]) -- [-4 -4 3] (div [7 8 9] 2) -- [3 4 4] +;; compiled select/update must not round integer columns through double +(set TdivI64 (table [v] (list [9007199254740993 123456789012345678]))) +(at (select {q: (div v 1) from: TdivI64}) 'q) -- [9007199254740993 123456789012345678] +(at (select {q: (div v 3) from: TdivI64}) 'q) -- [3002399751580331 41152263004115226] +(set TdivUpd (table [v q] (list [9007199254740993 123456789012345678] [0 0]))) +(at (update {q: (div v 1) from: TdivUpd}) 'q) -- [9007199254740993 123456789012345678] +(at (update {q: (div v 3) from: TdivUpd}) 'q) -- [3002399751580331 41152263004115226] From c6a9cfd1c4497d2ef67b6a9c85e704d435e09ca6 Mon Sep 17 00:00:00 2001 From: Evgen Byelozorov Date: Sun, 16 Aug 2026 21:39:31 +0200 Subject: [PATCH 3/3] =?UTF-8?q?perf(expr):=20fast-path=20integer=20floor-d?= =?UTF-8?q?iv=20within=20=C2=B12^53=20to=20the=20double=20kernel?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Keeping every integer OP_IDIV on the exact scalar int64 kernel (PR #403) made the hot columnar div path ~2x slower: hardware 64-bit idiv is scalar and replaced the vectorizable divsd+floor. Restore the fast path without losing exactness — per morsel, when both operands are within ±2^53 (a cheap scan, over the in-cache morsel in the fused kernels) delegate to the vectorizable double loop, which is bit-exact there (round(a/b) cannot cross an integer boundary in that range); fall back to the exact int64 kernel only for large magnitudes or nulls. Applied to the fused null-aware and non-null I64 kernels and all four binary_range output arms. Measured (cache-proof, distinct constant divisors, 20M rows ×10 iters, release): exact-only 0.41s -> dual-path 0.24s. Also folds in the review's smaller points: - ray_idiv_fn guards INT64_MIN/-1 in the integer path (symmetry with floor_idiv_i64_checked; unreachable today but removes latent UB); - the exact narrow I32/I16/U8 arms saturate to match their ray_cast_f64_to_iN_null double counterparts instead of wrapping; - binary_range's I64 OP_DIV arm uses floor_idiv_i64_checked, not the open-coded form; - 0x1p63 documents the scalar float-cast bound. Extends test/rfl/integration/math.rfl: fast/exact boundary agreement at 2^53, mixed small+big morsels, div-by-zero through the integer path, and narrow-output floor div. --- src/ops/arith.c | 10 ++-- src/ops/expr.c | 90 ++++++++++++++++++++++++++++------- test/rfl/integration/math.rfl | 18 +++++++ 3 files changed, 98 insertions(+), 20 deletions(-) diff --git a/src/ops/arith.c b/src/ops/arith.c index b8f15a1b7..9be671efd 100644 --- a/src/ops/arith.c +++ b/src/ops/arith.c @@ -291,16 +291,20 @@ ray_t* ray_idiv_fn(ray_t* a, ray_t* b) { if (bv == 0.0) return ray_typed_null(-RAY_I64); double q = floor(as_f64(a) / bv); - if (q >= 9223372036854775808.0 || q < (double)INT64_MIN) + if (q >= 0x1p63 /* 2^63: first double past INT64_MAX */ || q < (double)INT64_MIN) return ray_typed_null(-RAY_I64); return make_i64((int64_t)q); } /* Integer div: stay in int64 space. The double round-trip above silently * loses precision for magnitudes > 2^53 and is UB for q == 2^63. */ int64_t bv = as_i64(b); - if (bv == 0) - return ray_typed_null(-RAY_I64); int64_t la = as_i64(a); + /* bv==0 → null; INT64_MIN/-1 overflow → null (unreachable while INT64_MIN + * is the i64 null sentinel and caught above, but keep the guard so the + * scalar path can never trip signed-overflow UB — symmetric with the + * vector floor_idiv_i64_checked kernel). */ + if (bv == 0 || (la == INT64_MIN && bv == -1)) + return ray_typed_null(-RAY_I64); int64_t q = la / bv; if (la % bv != 0 && ((la < 0) != (bv < 0))) q--; /* floor toward -inf */ diff --git a/src/ops/expr.c b/src/ops/expr.c index 8a3dc5d2a..7fb61d70e 100644 --- a/src/ops/expr.c +++ b/src/ops/expr.c @@ -51,6 +51,24 @@ static inline bool floor_idiv_i64_checked(int64_t a, int64_t b, int64_t* out) { return true; } +/* Integer floor-division can be delegated to the vectorizable double kernel + * only while both operands stay within ±2^53. Such integers are exact + * doubles, and (proof) round(a/b) never crosses an integer boundary: a + * non-integer a/b is at least 1/|b| from any integer, while the rounded + * result's half-ULP is at most 1/|b| — with 2^e ≤ |a/b| we have + * |b|·2^e ≤ |a| ≤ 2^53, so 2^(e-53) ≤ 1/|b|. Thus floor of the rounded + * quotient equals the true integer floor. Outside the range the scalar + * exact kernel (floor_idiv_i64_checked) is required; NULL_I64 (INT64_MIN) + * lies far outside it, so a null in the span also falls to the exact kernel + * (which is where nulls are handled anyway). */ +#define DBL_EXACT_INT_LIM (0x1p53) /* 2^53 as a double bound */ +static inline bool i64_span_dbl_exact(const int64_t* v, int64_t n) { + for (int64_t i = 0; i < n; i++) + if (v[i] > (int64_t)DBL_EXACT_INT_LIM || v[i] < -(int64_t)DBL_EXACT_INT_LIM) + return false; + return true; +} + static bool atom_to_numeric(ray_t* atom, double* out_f, int64_t* out_i, bool* out_is_f64) { if (!atom || !ray_is_atom(atom)) return false; switch (atom->type) { @@ -1090,11 +1108,22 @@ static void expr_exec_binary(uint8_t opcode, uint8_t null_aware, int8_t dt, void * b==-1 && a==INT64_MIN → overflow: direct NULL_I64 (vs fallback: loop writes 0 * then zero-divisor pass marks null — same observable result). */ case OP_DIV: - case OP_IDIV: for (int64_t j=0;jtype == -RAY_F64 || rhs->type == RAY_F64 || rhs->type == -RAY_F32 || rhs->type == RAY_F32))); int src_is_i64_all = l_is_int && r_is_int; + /* Integer floor-div: decide once whether the whole range fits ±2^53, so the + * exact scalar int64 kernel is only paid for genuinely large magnitudes and + * the common small-value case stays on the vectorizable double kernel. */ + bool idiv_i64_small = false; + if (src_is_i64_all && op->opcode == OP_IDIV) { + idiv_i64_small = true; + for (int64_t i = 0; i < n; i++) { + int64_t lv = LV_READ_I64(i), rv = RV_READ_I64(i); + if (lv > (int64_t)DBL_EXACT_INT_LIM || lv < -(int64_t)DBL_EXACT_INT_LIM || + rv > (int64_t)DBL_EXACT_INT_LIM || rv < -(int64_t)DBL_EXACT_INT_LIM) { + idiv_i64_small = false; break; + } + } + } /* Hoist out_type outside the loop. Each branch is a tight per-element kernel. */ if (out_type == RAY_F64) { @@ -3188,9 +3240,9 @@ static void binary_range(ray_op_t* op, int8_t out_type, case OP_ADD: for(int64_t i=0;i=INT32_MAX?INT32_MAX:q<=INT32_MIN?INT32_MIN+1:(int32_t)q):0;} } else { for(int64_t i=0;i=INT16_MAX?INT16_MAX:q<=INT16_MIN?INT16_MIN+1:(int16_t)q):0;} } else { for(int64_t i=0;i=UINT8_MAX?UINT8_MAX:(uint8_t)q):0;} } else { for(int64_t i=0;i2^53 int64) kernels +;; must agree. 2^53 itself takes the fast path; 2^53+1 the exact one. ── +(div 9007199254740992 3) -- 3002399751580330 +(div [9007199254740992 9007199254740993] 1) -- [9007199254740992 9007199254740993] +;; a mixed small+big morsel forces the exact kernel for the whole span +(set TdivMix (table [a b] (list [10 9007199254740993] [3 7]))) +(at (select {q: (div a b) from: TdivMix}) 'q) -- [3 1286742750677284] +;; div-by-zero through the integer path: scalar/vector null, narrow arm 0 +(div 5 0) -- 0Nl +(div [5 6] [0 2]) -- [0Nl 3] +(set TdivZ (table [a b] (list [5 6] [0 2]))) +(at (select {q: (div a b) from: TdivZ}) 'q) -- [0Nl 3] +;; narrow-output floor div (I32 / I16 / U8 arms), incl. div-by-zero -> 0 +(div (as 'I32 [100 55 -7]) (as 'I32 [7 8 2])) -- (as 'I32 [14 6 -4]) +(div (as 'I16 [100 55 -7]) (as 'I16 [7 8 2])) -- (as 'I16 [14 6 -4]) +(set TdivN (table [a b] (list (as 'I16 [100 55 6]) (as 'I16 [7 0 2])))) +(at (select {q: (div a b) from: TdivN}) 'q) -- (as 'I16 [14 0 3])