diff --git a/src/check_expr.cpp b/src/check_expr.cpp index 9689d0812..cce251ddf 100644 --- a/src/check_expr.cpp +++ b/src/check_expr.cpp @@ -2367,6 +2367,14 @@ gb_internal void complex_quaternion_element_float_format(Type *type, int *mantis } } +gb_internal bool exact_value_component_overflows_float(ExactValue comp, int mantissa_bits, int ebias) { + if (comp.kind != ExactValue_Integer && comp.kind != ExactValue_Rational) { + return false; + } + ExactValue r = exact_value_round_component_to_float(comp, mantissa_bits, ebias); + return isinf(r.value_float) || isnan(r.value_float); +} + gb_internal bool check_representable_as_constant(CheckerContext *c, ExactValue in_value, Type *type, ExactValue *out_value) { if (in_value.kind == ExactValue_Invalid) { // NOTE(bill): There's already been an error @@ -2628,6 +2636,12 @@ gb_internal bool check_representable_as_constant(CheckerContext *c, ExactValue i imag.kind != ExactValue_Invalid) { int mantissa_bits, ebias; complex_quaternion_element_float_format(type, &mantissa_bits, &ebias); + // A finite component that overflows the element float is not representable (parity with + // scalar floats). Leave out_value unset so the diagnostic reports the source value. + if (exact_value_component_overflows_float(real, mantissa_bits, ebias) || + exact_value_component_overflows_float(imag, mantissa_bits, ebias)) { + return false; + } if (out_value) *out_value = exact_value_complex_ev( exact_value_round_component_to_float(real, mantissa_bits, ebias), exact_value_round_component_to_float(imag, mantissa_bits, ebias)); @@ -2660,6 +2674,14 @@ gb_internal bool check_representable_as_constant(CheckerContext *c, ExactValue i imag.kind != ExactValue_Invalid) { int mantissa_bits, ebias; complex_quaternion_element_float_format(type, &mantissa_bits, &ebias); + // A finite component that overflows the element float is not representable (parity with + // scalar floats). Leave out_value unset so the diagnostic reports the source value. + if (exact_value_component_overflows_float(real, mantissa_bits, ebias) || + exact_value_component_overflows_float(imag, mantissa_bits, ebias) || + exact_value_component_overflows_float(jmag, mantissa_bits, ebias) || + exact_value_component_overflows_float(kmag, mantissa_bits, ebias)) { + return false; + } if (out_value) *out_value = exact_value_quaternion_ev( exact_value_round_component_to_float(real, mantissa_bits, ebias), exact_value_round_component_to_float(imag, mantissa_bits, ebias), diff --git a/src/exact_value.cpp b/src/exact_value.cpp index 05e9f75c5..e1266ff01 100644 --- a/src/exact_value.cpp +++ b/src/exact_value.cpp @@ -1329,6 +1329,22 @@ gb_internal Entity *strip_entity_wrapping(Entity *e); gb_internal gbString write_expr_to_string(gbString str, Ast *node, bool shorthand); +gb_internal gbString write_exact_value_to_string(gbString str, ExactValue const &v, isize string_limit); + +gb_internal gbString write_exact_complex_component_to_string(gbString str, ExactValue comp, isize string_limit) { + f64 f = exact_value_to_f64(comp); + + // The float formatter cannot render a magnitude at or beyond 2**63 (it prints 2**63's digits on a + // loop), and that also catches Inf/NaN since the range test below is false for them. For an exact + // integer/rational component that large, print its exact form (digits, or `num.0/den`) instead. + static f64 const LIMIT = 9223372036854775808.0; // 2**63 + bool formatter_safe = (f >= -LIMIT) && (f <= LIMIT); + if (!formatter_safe && (comp.kind == ExactValue_Integer || comp.kind == ExactValue_Rational)) { + return write_exact_value_to_string(str, comp, string_limit); + } + return gb_string_append_fmt(str, "%.17g", f); +} + gb_internal gbString write_exact_value_to_string(gbString str, ExactValue const &v, isize string_limit=36) { switch (v.kind) { case ExactValue_Invalid: @@ -1395,9 +1411,19 @@ gb_internal gbString write_exact_value_to_string(gbString str, ExactValue const return str; } case ExactValue_Complex: - return gb_string_append_fmt(str, "%.17g+%.17gi", exact_value_to_f64(v.value_complex->real), exact_value_to_f64(v.value_complex->imag)); + str = write_exact_complex_component_to_string(str, v.value_complex->real, string_limit); + str = gb_string_append_fmt(str, "+"); + str = write_exact_complex_component_to_string(str, v.value_complex->imag, string_limit); + return gb_string_append_fmt(str, "i"); case ExactValue_Quaternion: - return gb_string_append_fmt(str, "%.17g+%.17gi+%.17gj+%.17gk", exact_value_to_f64(v.value_quaternion->real), exact_value_to_f64(v.value_quaternion->imag), exact_value_to_f64(v.value_quaternion->jmag), exact_value_to_f64(v.value_quaternion->kmag)); + str = write_exact_complex_component_to_string(str, v.value_quaternion->real, string_limit); + str = gb_string_append_fmt(str, "+"); + str = write_exact_complex_component_to_string(str, v.value_quaternion->imag, string_limit); + str = gb_string_append_fmt(str, "i+"); + str = write_exact_complex_component_to_string(str, v.value_quaternion->jmag, string_limit); + str = gb_string_append_fmt(str, "j+"); + str = write_exact_complex_component_to_string(str, v.value_quaternion->kmag, string_limit); + return gb_string_append_fmt(str, "k"); case ExactValue_Pointer: return str; diff --git a/tests/internal/test_number_literals.odin b/tests/internal/test_number_literals.odin index 8053fe080..8e2ff5b7f 100644 --- a/tests/internal/test_number_literals.odin +++ b/tests/internal/test_number_literals.odin @@ -94,6 +94,47 @@ float_literal_f16_f32_precision :: proc(t: ^testing.T) { testing.expect_value(t, f32(1e38), strconv.parse_f32("1e38") or_else 0) } +@(test) +rational_arithmetic_precision_cap :: proc(t: ^testing.T) { + // Exact-rational constant folding is bounded so a pathological expression cannot grow the + // numerator/denominator without limit. Repeatedly squaring a non-dyadic fraction doubles the + // denominator's bit-length at every step; past the precision cap the fold falls back to a rounded + // f64. This whole chain therefore folds in a few kilobytes; without the cap `X32` alone would need a + // denominator of ~10**(2**32) (gigabytes) to fold exactly, stalling or OOMing the compiler. The test + // passing quickly *is* the regression check — that the guard keeps runaway folding bounded. + X0 :: 0.3 + X1 :: X0*X0 + X2 :: X1*X1 + X3 :: X2*X2 + X4 :: X3*X3 + X5 :: X4*X4 + X6 :: X5*X5 + X7 :: X6*X6 + X8 :: X7*X7 + X9 :: X8*X8 + X10 :: X9*X9 + X11 :: X10*X10 + X12 :: X11*X11 + X13 :: X12*X12 + X14 :: X13*X13 + X15 :: X14*X14 + X16 :: X15*X15 + X17 :: X16*X16 + X18 :: X17*X17 + X19 :: X18*X18 + X20 :: X19*X19 + X24 :: (X20*X20)*(X20*X20) // 0.3 ** 2**24 + X28 :: (X24*X24)*(X24*X24) + X32 :: (X28*X28)*(X28*X28) // 0.3 ** 2**32 + + // 0.3 ** 2**32 is astronomically small, so once folding falls back to f64 it underflows to 0. + testing.expect(t, X32 == 0.0, "rational precision cap: deeply-squared fraction folds to a bounded f64") + testing.expect(t, !(X32 != X32), "capped value is a real number, not NaN") + + // A shallow fold is still exact: 0.3 ** 4 == 81/10000 rounds to the same f64 as the literal 0.0081. + testing.expect(t, X2 == 0.0081, "shallow rational folding stays exact") +} + @(test) float_constant_builtins :: proc(t: ^testing.T) { // Constant-folded builtins on decimal-float (rational) constants must not crash or mis-fold. diff --git a/tests/internal/test_quat_cmplx.odin b/tests/internal/test_quat_cmplx.odin index bf8df6f0d..e7cc5c22d 100644 --- a/tests/internal/test_quat_cmplx.odin +++ b/tests/internal/test_quat_cmplx.odin @@ -1,6 +1,7 @@ package test_internal import "core:testing" +import "core:strconv" // Constant folding of complex and quaternion values, against the answer the backend // produces. Every constant case is paired with the same expression on variables: the @@ -327,3 +328,91 @@ accessors_keep_their_bits_through_transmute :: proc(t: ^testing.T) { testing.expect_value(t, transmute(u64)jmag(q), transmute(u64)j) testing.expect_value(t, transmute(u64)kmag(q), transmute(u64)k) } + +// Untyped complex/quaternion constant components are now held exactly (as big rationals) and the +// arithmetic folds exactly, rounding only once when the value is finally given a concrete type. This +// is the complex/quaternion analogue of `0.1 + 0.2 == 0.3` for untyped floats: before, the components +// were pre-rounded to f64 and every lane folded in binary floating point, coming out a ULP off. + +@(test) +constant_complex_exact_folding :: proc(t: ^testing.T) { + X :: 0.1 + 0.2i + Y :: 0.2 + 0.1i + testing.expect(t, complex128(X + Y) == complex128(0.3 + 0.3i), "exact complex constant addition") + + // (1/10) * 3 == 3/10 exactly, in the real lane + testing.expect(t, complex128((0.1+0i) * 3) == complex128(0.3+0i), "exact complex constant multiplication") + + // (0.3 + 0.6i) / 3 == (0.1 + 0.2i) exactly (true rational division, not fmod) + testing.expect(t, complex128((0.3 + 0.6i) / 3) == complex128(0.1 + 0.2i), "exact complex constant division") + + // The same expression on runtime f64 variables still rounds in binary floating point, so the + // folded and runtime answers genuinely differ here. + xr, xi := 0.1, 0.2 + yr, yi := 0.2, 0.1 + runtime_sum := complex(xr, xi) + complex(yr, yi) + testing.expect(t, runtime_sum != complex128(0.3 + 0.3i), "runtime complex arithmetic still rounds") +} + +@(test) +constant_quaternion_exact_folding :: proc(t: ^testing.T) { + X :: 0.1 + 0.2i + 0.3j + 0.4k + Y :: 0.2 + 0.1i + 0.4j + 0.3k + testing.expect(t, quaternion256(X + Y) == quaternion256(0.3 + 0.3i + 0.7j + 0.7k), + "exact quaternion constant addition") + + // Multiplying by a real scalar keeps every lane exact. + testing.expect(t, quaternion256((0.1 + 0.2i + 0.3j + 0.4k) * 10) == quaternion256(1 + 2i + 3j + 4k), + "exact quaternion constant multiplication") + + q64 := quaternion(w=0.1, x=0.2, y=0.3, z=0.4) + quaternion(w=0.2, x=0.1, y=0.4, z=0.3) + testing.expect(t, q64 != quaternion256(0.3 + 0.3i + 0.7j + 0.7k), "runtime quaternion arithmetic still rounds") +} + +// A finite component that overflows the element float's range is rejected at compile time (parity with +// scalar floats), so `complex128(1.0e400)` is an error rather than a silent +Inf. The rejection itself is +// a compile error and can't be asserted here; these near-max constants guard the other side - that +// representable components are NOT falsely rejected. + +@(test) +constant_complex_representable_bounds :: proc(t: ^testing.T) { + c128 :: complex128(1e308 + 1e308i) // < f64 max (~1.8e308) + testing.expect_value(t, real(c128), 1e308) + testing.expect_value(t, imag(c128), 1e308) + + c64 :: complex64(3e38 + 3e38i) // < f32 max (~3.4e38) + testing.expect(t, real(c64) > 0 && imag(c64) > 0, "complex64 near-max components fold") + + c32 :: complex32(60000 + 60000i) // < f16 max (65504) + testing.expect(t, real(c32) > 0 && imag(c32) > 0, "complex32 near-max components fold") + + q256 :: quaternion256(1e308 + 1e308i + 1e308j + 1e308k) + testing.expect_value(t, jmag(q256), 1e308) + testing.expect_value(t, kmag(q256), 1e308) + + q64 :: quaternion64(60000 + 60000i + 60000j + 60000k) + testing.expect(t, jmag(q64) > 0 && kmag(q64) > 0, "quaternion64 near-max lanes fold") +} + +// complex64/complex32 (and quaternion128/quaternion64) components are rounded exactly once, directly +// from the exact literal to the element format, so a folded constant matches the correctly-rounded +// runtime parse instead of double-rounding through f64. + +@(test) +constant_complex_f32_rounding :: proc(t: ^testing.T) { + c32 :: proc(t: ^testing.T, got: f32, lit: string) { + want, _ := strconv.parse_f32(lit) + testing.expectf(t, transmute(u32)got == transmute(u32)want, + "f32 %s: got %08x, want %08x", lit, transmute(u32)got, transmute(u32)want) + } + rr : f32 : real(complex64(0.1 + 0.2i)) + ii : f32 : imag(complex64(0.1 + 0.2i)) + c32(t, rr, "0.1") + c32(t, ii, "0.2") + c32(t, real(complex64(3.14159265358979323846 + 2.71828182845904523536i)), "3.14159265358979323846") + c32(t, imag(complex64(3.14159265358979323846 + 2.71828182845904523536i)), "2.71828182845904523536") + + // quaternion128 lanes as well + q1 : f32 : jmag(quaternion128(0 + 0i + 0.1j + 0k)) + c32(t, q1, "0.1") +}