diff --git a/.kiro/specs/calculator/design.md b/.kiro/specs/calculator/design.md index 681005a..b6d4262 100644 --- a/.kiro/specs/calculator/design.md +++ b/.kiro/specs/calculator/design.md @@ -509,11 +509,21 @@ Bug-by-bug, where the fix lands: | `0.1 + 0.2` -> `0.3` | yes | - | | `1.1 + 2.2` -> `3.3` | yes | - | | `0.1 * 3` -> `0.3` | yes | - | -| `12 in in ft` -> `1` | yes | - | +| `1e20 + 1 - 1e20` -> `1` | yes | - | | `factorial(171)` | partly (becomes `inf`, not an error) | yes, for the exact value | | `9007199254740993` | **no** (f64 cannot hold it) | yes | | `2^53 + 1` | **no** | yes | | `1/3` retained exactly | no (only observable as exact) | yes | +| `12 in in ft` -> `1` | **no** - see below | needs exact unit factors | + +**Unit conversion is a separate problem from the evaluator.** An earlier draft of +this table claimed the exact evaluator would fix `12 in in ft`. It does not, and +the reason is worth recording: `UnitDef.to_base_factor` is an `f64`, so `0.0254` +is *already* the binary approximation before any arithmetic happens, and +`convertUnits` does its own f64 multiply and divide without ever touching the +evaluator. Making conversions exact requires the factors themselves to be exact +(declared as decimal text and parsed into rationals), which is a mechanical change +across ~100 table entries and belongs in its own commit. Tracked as Task 2.0e. So step 2 is independently verifiable through the current API, and the large-integer cases become the motivating tests for step 3 rather than diff --git a/.kiro/specs/calculator/tasks.md b/.kiro/specs/calculator/tasks.md index ed90610..e471603 100644 --- a/.kiro/specs/calculator/tasks.md +++ b/.kiro/specs/calculator/tasks.md @@ -126,16 +126,47 @@ behavioral change (rationale in design.md 2.7.9): - Verify: 487 tests pass (was 402), zig fmt and zlint clean, existing tests untouched. -#### 2.0b: Evaluator internals become exact, `evalString` still returns f64 -- Exact/inexact boundary per design.md 2.7.4. Per-evaluation, not per-function: - `sqrt(4)` exact, `sqrt(2)` inexact. Deliberately not a CAS. -- `NumberValue` and the AST number node GAIN an exact field rather than replacing - `float`/`int_value`, so tokenizer and parser tests keep compiling. -- ALL existing tests must still pass unchanged. -- Payoff is visible through the f64 boundary because a single final rounding - avoids today's accumulated error: `0.1 + 0.2` -> `0.3`, `1.1 + 2.2` -> `3.3`, - `0.1 * 3` -> `0.3`, `12 in in ft` -> `1`. -- Verify: existing suite green, plus new tests for the four fixes above. +#### 2.0b: Evaluator internals become exact, `evalString` still returns f64 [DONE] +- `ast.Expr.Number` gains `text: []const u8` (the literal's source text). Additive, + so the 10 parser assertions on `float_value`/`int_value` keep compiling. The + text is needed because `float_value` has ALREADY rounded by the time it exists: + `0.1` cannot be recovered from its binary approximation, so exact evaluation has + to re-parse the literal. +- `evaluator.evaluate` keeps its `CalcError!f64` signature but now runs on + `Number` internally via `evalExact`, collapsing to f64 once at the boundary. + A scratch arena per evaluation keeps `Number` lifetimes trivial and means the + caller's allocator (an arena in the CLI, a GPA in the TUI) never holds + intermediates. +- Exact: `+ - * / %`, integer powers, `sqrt` of perfect squares, `abs`, `floor`, + `ceil`, `round`, `factorial`, `max`, `min`, and all literals. +- Inexact: transcendentals, fractional powers, `sqrt` of non-squares, `cbrt`, + `atan2`, two-arg `log`, `rand`, and the constants `pi`/`e`/`tau`. +- Fixed-width integer operators (`& | xor << >> >>> rol ror`, `~`) stay on the + float/integer projection: they are not rational arithmetic. +- Added to `rational.zig`: `floor`, `ceil`, `round`, `mod` (the + `a - b*floor(a/b)` definition), and unbounded exact `factorial`. +- Added to `number.zig`: the matching wrappers plus `max`/`min`. +- `rational.Error.ExponentTooLarge` maps to `CalcError.Overflow`. +- ALL 489 pre-existing tests pass unchanged. 522 total now (+33). +- Verified through the unchanged f64 API: `0.1 + 0.2` = `0.3`, `1.1 + 2.2` = `3.3`, + `0.1 * 3` = `0.3`, `0.1+0.2+0.3` = `0.6`, `(0.1+0.2)*10-3` = `0`, + `1e20 + 1 - 1e20` = `1` (f64 gives 0), `1/3*3` = `1`, `factorial(171)` no longer + errors. Unchanged: `2 + 3 * 4`, `sqrt(144)`, `sqrt(2)`, `sin(0)`, `e`, + `0o777 - 0x0f`, `10 % 3`, `0xFF & 0x0F`. +- Known limitations deferred to 2.0c/2.0e: + - Variables and `Ans` are still stored as f64 in `Environment`, so + `X = 0.1` round-trips through a float. Correct for `pi`/`e`/`tau` (irrational), + a real limitation for user variables. + - `9007199254740993` and `2^53 + 1` are computed exactly but cannot survive the + f64 return. That is precisely what 2.0c fixes. + - `12 in in ft` is NOT fixed: unit conversion never touches the evaluator and + its factors are already-rounded f64 values. See Task 2.0e. +- Test-oracle bug found by the instrumented coverage build: Zig's float `@mod` + with a NEGATIVE divisor returned `-2` in the normal build and `1` under the + coverage build (comptime folding and the runtime path disagree). A test whose + expected value shifts with optimize mode verifies nothing, so the `mod` test now + asserts the values the definition requires instead of deriving them from `@mod`. + Worth reporting upstream. #### 2.0c: `evalString` and the display path move to `Number` - Mechanical signature churn through evaluator/formatter/CLI/TUI tests. This is @@ -158,6 +189,24 @@ because that is the precision we can justify; see design.md 2.7.8 for why f128 i NOT the upgrade path), and any change to programmer mode's `u128` semantics or the IEEE 754 float view. +### Task 2.0e: Exact unit conversion factors [NOT STARTED] +Discovered while implementing 2.0b: making the evaluator exact does NOT fix +`12 in in ft` = `0.9999999999999998`, because unit conversion never goes through +the evaluator. `UnitDef.to_base_factor` is an `f64`, so `0.0254` is already the +binary approximation before `convertUnits` does its own f64 multiply and divide. + +- Declare factors as exact decimal text (e.g. `"0.0254"`) so they can be parsed + into rationals; an inch is exactly `127/5000` m and a foot exactly `381/1250` m, + which makes `12 in to ft` exactly `1`. +- Provide an exact conversion path returning `Number`, keeping the f64 path for + callers that want it. +- Mechanical across ~100 table entries, hence its own commit: the existing + invariant tests (every unit and unit pair round-trips) become far stronger when + the round trip is exact rather than within a tolerance. +- Non-terminating conversions (`100 km to mi` = `781250/12573`) still render as + rounded decimals, but from an exact value, and the exact fraction becomes + available for display. + ### Task 2.1: Implement struct DSL tokenizer and parser - Create `engine/src/struct_layout.zig` - Tokenize: `struct`, `{`, `}`, `;`, type names, identifiers, `[`, `]`, numbers (for arrays) diff --git a/engine/src/ast.zig b/engine/src/ast.zig index e8ad761..8ec4941 100644 --- a/engine/src/ast.zig +++ b/engine/src/ast.zig @@ -21,6 +21,14 @@ pub const Expr = union(enum) { /// If the number is a pure integer, stores the exact value. int_value: ?u64, base: Base, + /// The literal's source text, with separators intact. Points into the + /// expression source, which outlives the AST. + /// + /// Kept so the evaluator can reconstruct the literal *exactly* as a + /// rational: `float_value` has already lost information by the time it + /// exists (0.1 is not representable in binary), so it cannot be the + /// basis for exact arithmetic. See design.md 2.7. + text: []const u8 = "", }; pub const Unary = struct { diff --git a/engine/src/evaluator.zig b/engine/src/evaluator.zig index a577b0e..33962f8 100644 --- a/engine/src/evaluator.zig +++ b/engine/src/evaluator.zig @@ -17,6 +17,8 @@ const ProgrammerConfig = types.ProgrammerConfig; const CalcError = types.CalcError; const parser_mod = @import("parser.zig"); const Parser = parser_mod.Parser; +const number_mod = @import("number.zig"); +const Number = number_mod.Number; /// Evaluation environment holding variables, history, and config. pub const Environment = struct { @@ -61,75 +63,138 @@ pub const Environment = struct { /// Evaluate a parsed expression in the given environment. /// Returns the computed value as f64 for standard mode. +/// +/// Internally the computation runs on `Number`, so exact arithmetic is used +/// wherever possible and only collapses to f64 here, at the boundary. That +/// single final rounding is what fixes the accumulated-error class of bug: +/// `0.1 + 0.2` is computed as exactly `3/10` and rounds to the f64 nearest +/// `0.3`, rather than adding two separately-rounded operands. +/// +/// Task 2.0c replaces this boundary with a `Number`-returning API, which is what +/// the remaining integer-precision cases need. pub fn evaluate(env: *Environment, expr: *const Expr) CalcError!f64 { + // A scratch arena keeps Number lifetimes trivial: nothing in the recursive + // evaluator has to free intermediates, and the caller's allocator is never + // left holding them regardless of whether it is an arena itself. + var arena = std.heap.ArenaAllocator.init(env.allocator); + defer arena.deinit(); + const scratch = arena.allocator(); + + const result = try evalExact(env, scratch, expr); + return result.toFloat(scratch); +} + +/// The exact evaluation core. Produces a `Number`, staying exact until an +/// operation forces the float fallback (see design.md 2.7.4). +fn evalExact(env: *Environment, scratch: Allocator, expr: *const Expr) CalcError!Number { switch (expr.*) { - .number => |n| return n.float_value, + .number => |n| return literalToNumber(scratch, n), .string_literal => |text| { - // Pack ASCII bytes into integer (same as programmer mode, BE packing) - var result: u128 = 0; + // Pack ASCII bytes into an integer (BE packing, as programmer mode). + var packed_value: u128 = 0; for (text) |byte| { if (byte > 0x7F) return CalcError.InvalidNumber; - result = (result << 8) | byte; + packed_value = (packed_value << 8) | byte; } - return @floatFromInt(result); + return Number.fromInt(scratch, packed_value) catch |err| return mapError(err); }, .variable => |name| { - return env.getVar(name) orelse return CalcError.UnknownVariable; + // Variables and the built-in constants are stored as f64 today, so + // reading one yields an inexact value. For pi/e/tau that is correct + // (they are irrational); for user variables it is a temporary + // limitation that Task 2.0c removes by storing Number in the + // environment. + const value = env.getVar(name) orelse return CalcError.UnknownVariable; + return Number.fromFloat(value); }, .assignment => |a| { - const val = try evaluate(env, a.value); - env.setVar(a.name, val) catch return CalcError.OutOfMemory; + const val = try evalExact(env, scratch, a.value); + env.setVar(a.name, val.toFloat(scratch)) catch return CalcError.OutOfMemory; return val; }, .unary => |u| { - const operand = try evaluate(env, u.operand); + const operand = try evalExact(env, scratch, u.operand); return switch (u.op) { - .negate => -operand, - .bitwise_not => { - // In standard mode, bitwise not doesn't really make sense, - // but we'll compute it on the integer representation - const int_val: u64 = @bitCast(@as(i64, @intFromFloat(operand))); + .negate => Number.negate(scratch, operand) catch |err| mapError(err), + // Bitwise NOT is a fixed-width integer operation, not rational + // arithmetic, so it drops to the float/integer path. + .bitwise_not => blk: { + const f = operand.toFloat(scratch); + const int_val: u64 = @bitCast(@as(i64, @intFromFloat(f))); const mask_val: u64 = @truncate(env.programmer_config.bit_width.mask()); const result = ~int_val & mask_val; - return @floatFromInt(@as(i64, @bitCast(result))); + break :blk Number.fromFloat(@floatFromInt(@as(i64, @bitCast(result)))); }, }; }, .binary => |b| { - const left = try evaluate(env, b.left); - const right = try evaluate(env, b.right); - return evalBinaryOp(b.op, left, right); + const left = try evalExact(env, scratch, b.left); + const right = try evalExact(env, scratch, b.right); + return evalBinaryOp(scratch, b.op, left, right); }, .call => |c| { - return evalFunction(env, c.name, c.args); + return evalFunction(env, scratch, c.name, c.args); }, } } -/// Evaluate a binary operation on two f64 values. -fn evalBinaryOp(op: BinaryOp, left: f64, right: f64) CalcError!f64 { +/// Turn a literal into a Number, exactly where possible. +/// +/// Decimal literals are re-parsed from their source text rather than taken from +/// `float_value`, because `float_value` has already rounded: `0.1` cannot be +/// recovered from its binary approximation. +fn literalToNumber(scratch: Allocator, n: ast.Expr.Number) CalcError!Number { + if (n.base == .decimal and n.text.len > 0) { + if (Number.parse(scratch, n.text)) |value| return value else |_| { + // Fall through to the approximations below rather than failing: the + // tokenizer already accepted this text, so a parse mismatch here + // should degrade, not error. + } + } + // Non-decimal literals are integers; use the exact integer the tokenizer + // recovered when it fits, otherwise accept the float approximation. + if (n.int_value) |int_val| { + return Number.fromInt(scratch, int_val) catch |err| return mapError(err); + } + return Number.fromFloat(n.float_value); +} + +/// Map the numeric model's errors onto the engine's error set. +fn mapError(err: number_mod.Error) CalcError { + return switch (err) { + error.OutOfMemory => CalcError.OutOfMemory, + error.DivisionByZero => CalcError.DivisionByZero, + error.InvalidNumber => CalcError.InvalidNumber, + // An exponent too large to compute is an overflow from the caller's view. + error.ExponentTooLarge => CalcError.Overflow, + }; +} + +/// Evaluate a binary operation. +fn evalBinaryOp(scratch: Allocator, op: BinaryOp, left: Number, right: Number) CalcError!Number { return switch (op) { - .add => left + right, - .sub => left - right, - .mul => left * right, - .div => if (right == 0) CalcError.DivisionByZero else left / right, - .mod => if (right == 0) CalcError.DivisionByZero else @mod(left, right), - .pow => math.pow(f64, left, right), - // Bitwise ops in standard mode operate on integer truncations - .bit_and => floatBitwise(left, right, bitwiseAnd), - .bit_or => floatBitwise(left, right, bitwiseOr), - .bit_xor => floatBitwise(left, right, bitwiseXor), - .shift_left => floatShift(left, right, true), - .shift_right, .shift_right_logical => floatShift(left, right, false), - .rotate_left, .rotate_right => { + .add => Number.add(scratch, left, right) catch |err| mapError(err), + .sub => Number.sub(scratch, left, right) catch |err| mapError(err), + .mul => Number.mul(scratch, left, right) catch |err| mapError(err), + .div => Number.div(scratch, left, right) catch |err| mapError(err), + .mod => Number.mod(scratch, left, right) catch |err| mapError(err), + .pow => Number.pow(scratch, left, right) catch |err| mapError(err), + // The remaining operators are fixed-width integer operations rather than + // rational arithmetic, so they work on the float/integer projection. + .bit_and => Number.fromFloat(floatBitwise(left.toFloat(scratch), right.toFloat(scratch), bitwiseAnd)), + .bit_or => Number.fromFloat(floatBitwise(left.toFloat(scratch), right.toFloat(scratch), bitwiseOr)), + .bit_xor => Number.fromFloat(floatBitwise(left.toFloat(scratch), right.toFloat(scratch), bitwiseXor)), + .shift_left => Number.fromFloat(floatShift(left.toFloat(scratch), right.toFloat(scratch), true)), + .shift_right, .shift_right_logical => Number.fromFloat(floatShift(left.toFloat(scratch), right.toFloat(scratch), false)), + .rotate_left, .rotate_right => blk: { // Rotations need bit width context; in standard mode, use 64-bit - const l: u64 = @bitCast(@as(i64, @intFromFloat(left))); - const r: u6 = @intFromFloat(@mod(right, 64.0)); + const l: u64 = @bitCast(@as(i64, @intFromFloat(left.toFloat(scratch)))); + const r: u6 = @intFromFloat(@mod(right.toFloat(scratch), 64.0)); const result = if (op == .rotate_left) math.rotl(u64, l, r) else math.rotr(u64, l, r); - return @floatFromInt(@as(i64, @bitCast(result))); + break :blk Number.fromFloat(@floatFromInt(@as(i64, @bitCast(result)))); }, }; } @@ -159,25 +224,59 @@ fn floatShift(left: f64, right: f64, is_left: bool) f64 { } /// Evaluate a built-in function call. -fn evalFunction(env: *Environment, name: []const u8, args: []const *Expr) CalcError!f64 { +fn evalFunction(env: *Environment, scratch: Allocator, name: []const u8, args: []const *Expr) CalcError!Number { // Single-argument functions if (args.len == 1) { - const x = try evaluate(env, args[0]); - return evalSingleArgFn(name, x) orelse CalcError.UnknownFunction; + const x = try evalExact(env, scratch, args[0]); + + // Functions with an exact implementation. + if (std.mem.eql(u8, name, "abs")) { + return Number.abs(scratch, x) catch |err| mapError(err); + } + if (std.mem.eql(u8, name, "floor")) { + return Number.floor(scratch, x) catch |err| mapError(err); + } + if (std.mem.eql(u8, name, "ceil")) { + return Number.ceil(scratch, x) catch |err| mapError(err); + } + if (std.mem.eql(u8, name, "round")) { + return Number.round(scratch, x) catch |err| mapError(err); + } + if (std.mem.eql(u8, name, "sqrt")) { + // Negative inputs are a domain error rather than a NaN. + if (x.isNegative()) return CalcError.UnknownFunction; + return Number.sqrt(scratch, x) catch |err| mapError(err); + } + if (std.mem.eql(u8, name, "factorial")) { + const result = Number.factorial(scratch, x) catch |err| return mapError(err); + return result orelse CalcError.UnknownFunction; + } + + // Everything else escapes the rationals, so it falls back to f64. + const f = evalSingleArgFn(name, x.toFloat(scratch)) orelse + return CalcError.UnknownFunction; + return Number.fromFloat(f); } // Multi-argument functions if (args.len == 2) { - const a = try evaluate(env, args[0]); - const b = try evaluate(env, args[1]); + const a = try evalExact(env, scratch, args[0]); + const b = try evalExact(env, scratch, args[1]); - if (std.mem.eql(u8, name, "max")) return @max(a, b); - if (std.mem.eql(u8, name, "min")) return @min(a, b); - if (std.mem.eql(u8, name, "atan2")) return math.atan2(a, b); + if (std.mem.eql(u8, name, "max")) { + return Number.max(scratch, a, b) catch |err| mapError(err); + } + if (std.mem.eql(u8, name, "min")) { + return Number.min(scratch, a, b) catch |err| mapError(err); + } + + const x = a.toFloat(scratch); + const y = b.toFloat(scratch); + if (std.mem.eql(u8, name, "atan2")) return Number.fromFloat(math.atan2(x, y)); if (std.mem.eql(u8, name, "log")) { // log(value, base) - if (b <= 0 or b == 1 or a <= 0) return CalcError.DomainError; - return @log(a) / @log(b); + if (y <= 0 or y == 1 or x <= 0) return CalcError.DomainError; + return Number.fromFloat(@log(x) / @log(y)); } } @@ -185,14 +284,14 @@ fn evalFunction(env: *Environment, name: []const u8, args: []const *Expr) CalcEr if (args.len == 0) { if (std.mem.eql(u8, name, "rand")) { // Not truly random in a pure engine, but useful as placeholder - return 0.0; + return Number.fromFloat(0.0); } } return CalcError.UnknownFunction; } -/// Evaluate a single-argument built-in function. +/// Evaluate a single-argument built-in function that has no exact form. fn evalSingleArgFn(name: []const u8, x: f64) ?f64 { if (std.mem.eql(u8, name, "sin")) return @sin(x); if (std.mem.eql(u8, name, "cos")) return @cos(x); @@ -210,33 +309,11 @@ fn evalSingleArgFn(name: []const u8, x: f64) ?f64 { if (std.mem.eql(u8, name, "log10")) return @log10(x); if (std.mem.eql(u8, name, "ln")) return @log(x); if (std.mem.eql(u8, name, "log2")) return @log2(x); - if (std.mem.eql(u8, name, "sqrt")) { - if (x < 0) return null; - return @sqrt(x); - } if (std.mem.eql(u8, name, "cbrt")) return math.cbrt(x); - if (std.mem.eql(u8, name, "abs")) return @abs(x); - if (std.mem.eql(u8, name, "ceil")) return @ceil(x); - if (std.mem.eql(u8, name, "floor")) return @floor(x); - if (std.mem.eql(u8, name, "round")) return @round(x); if (std.mem.eql(u8, name, "exp")) return @exp(x); - if (std.mem.eql(u8, name, "factorial")) { - if (x < 0 or x != @round(x) or x > 170) return null; - return factorial(@intFromFloat(x)); - } return null; } -fn factorial(n: u64) f64 { - if (n <= 1) return 1.0; - var result: f64 = 1.0; - var i: u64 = 2; - while (i <= n) : (i += 1) { - result *= @floatFromInt(i); - } - return result; -} - /// Result of evaluation with metadata for display decisions. pub const EvalInfo = struct { value: f64, @@ -653,3 +730,124 @@ test "eval unknown two-arg function" { const result = testEval("bogus(1, 2)"); try testing.expectError(CalcError.UnknownFunction, result); } + +// -- Exact arithmetic (Task 2.0b) -- +// +// These verify the exact evaluation core through the unchanged f64 API. The +// payoff is visible here because today's errors are ACCUMULATED: f64 rounds +// each decimal literal before operating on it, whereas the exact core computes +// the true value and rounds once, at the boundary. +// +// Note the runtime-`var` dance in the comparisons against plain f64: Zig folds +// float literals at comptime as `comptime_float`, so `0.1 + 0.2 != 0.3` is +// false at comptime and would not exercise f64 at all. + +test "exact: 0.1 + 0.2 is 0.3" { + try testing.expectEqual(@as(f64, 0.3), try testEval("0.1 + 0.2")); + + var x: f64 = 0.1; + var y: f64 = 0.2; + _ = &x; + _ = &y; + try testing.expect(x + y != @as(f64, 0.3)); +} + +test "exact: 1.1 + 2.2 is 3.3" { + try testing.expectEqual(@as(f64, 3.3), try testEval("1.1 + 2.2")); +} + +test "exact: 0.1 * 3 is 0.3" { + try testing.expectEqual(@as(f64, 0.3), try testEval("0.1 * 3")); +} + +test "exact: chained decimal addition" { + try testing.expectEqual(@as(f64, 0.6), try testEval("0.1 + 0.2 + 0.3")); + try testing.expectEqual(@as(f64, 0.8), try testEval("0.7 + 0.1")); + try testing.expectEqual(@as(f64, 0.2), try testEval("0.3 - 0.1")); + try testing.expectEqual(@as(f64, 0.01), try testEval("0.1 * 0.1")); +} + +test "exact: an expression that cancels reaches exactly zero" { + try testing.expectEqual(@as(f64, 0.0), try testEval("(0.1 + 0.2) * 10 - 3")); + + var x: f64 = 0.1; + var y: f64 = 0.2; + _ = &x; + _ = &y; + try testing.expect((x + y) * 10.0 - 3.0 != 0.0); +} + +test "exact: intermediates beyond f64 precision survive" { + // 1e20 + 1 is not representable in f64, so the f64 route loses the 1 and + // yields 0. Exact arithmetic keeps it and the final result fits. + try testing.expectEqual(@as(f64, 1.0), try testEval("1e20 + 1 - 1e20")); + + var big: f64 = 1e20; + _ = &big; + try testing.expectEqual(@as(f64, 0.0), big + 1.0 - big); +} + +test "exact: division round trip" { + try testing.expectEqual(@as(f64, 1.0), try testEval("1 / 3 * 3")); + try testing.expectEqual(@as(f64, 1.0), try testEval("1 / 7 * 7")); + try testing.expectEqual(@as(f64, 100.5), try testEval("1.005 * 100")); +} + +test "exact: factorial is no longer capped at 170" { + // Previously `factorial(171)` reported "unknown function" because the f64 + // implementation overflowed. It now computes exactly and only loses + // magnitude at the f64 boundary. + const result = try testEval("factorial(171)"); + try testing.expect(math.isPositiveInf(result)); + + // And a value f64 can still hold comes back exact. + try testing.expectEqual(@as(f64, 120.0), try testEval("factorial(5)")); +} + +test "exact: perfect square roots stay exact, irrational ones fall back" { + try testing.expectEqual(@as(f64, 12.0), try testEval("sqrt(144)")); + try testing.expectEqual(@as(f64, 0.5), try testEval("sqrt(0.25)")); + try testing.expectApproxEqAbs(math.sqrt2, try testEval("sqrt(2)"), 1e-15); +} + +test "exact: floor, ceil and round match the float builtins" { + try testing.expectEqual(@as(f64, -4.0), try testEval("floor(-3.2)")); + try testing.expectEqual(@as(f64, -3.0), try testEval("ceil(-3.2)")); + try testing.expectEqual(@as(f64, 3.0), try testEval("round(2.5)")); + try testing.expectEqual(@as(f64, -3.0), try testEval("round(-2.5)")); + try testing.expectEqual(@as(f64, 0.1), try testEval("abs(-0.1)")); +} + +test "exact: mod keeps the sign of the divisor" { + try testing.expectEqual(@as(f64, 1.0), try testEval("10 % 3")); + try testing.expectEqual(@as(f64, 2.0), try testEval("-10 % 3")); + try testing.expectEqual(@as(f64, 0.5), try testEval("7.5 % 1")); +} + +test "exact: non-decimal literals are exact integers" { + try testing.expectEqual(@as(f64, 255.0), try testEval("0xFF")); + try testing.expectEqual(@as(f64, 496.0), try testEval("0o777 - 0x0f")); + try testing.expectEqual(@as(f64, 10.0), try testEval("0b1010")); +} + +test "exact: transcendentals still fall back to floats" { + // These have no exact rational form, so they must go through f64 and are + // only expected to be approximately right. + try testing.expectApproxEqAbs(@as(f64, 0.0), try testEval("sin(0)"), 1e-15); + try testing.expectApproxEqAbs(math.pi, try testEval("pi"), 1e-15); + try testing.expectApproxEqAbs(@as(f64, 1.0), try testEval("ln(e)"), 1e-15); + try testing.expectApproxEqAbs(@as(f64, 2.0), try testEval("log10(100)"), 1e-15); +} + +test "exact: a transcendental contaminates the rest of the expression" { + // Once sin() enters, the result is inexact; it must still be numerically + // right, just not exact. + const result = try testEval("sin(0) + 0.1 + 0.2"); + try testing.expectApproxEqAbs(@as(f64, 0.3), result, 1e-15); +} + +test "exact: overflow from an absurd exponent is reported as overflow" { + // The exponent guard in the rational layer surfaces as Overflow rather than + // silently producing infinity or exhausting memory. + try testing.expectError(CalcError.Overflow, testEval("2 ^ 3000000")); +} diff --git a/engine/src/number.zig b/engine/src/number.zig index fa3fd57..e823183 100644 --- a/engine/src/number.zig +++ b/engine/src/number.zig @@ -228,6 +228,71 @@ pub const Number = union(enum) { return .{ .inexact = f(a.toFloat(allocator)) }; } + /// Shape of a unary operation that has an exact implementation. + fn unary( + allocator: Allocator, + a: Number, + comptime exactOp: fn (Allocator, Rational) Error!Rational, + comptime floatOp: fn (f64) f64, + ) Error!Number { + return switch (a) { + .exact => |r| capped(allocator, try exactOp(allocator, r)), + .inexact => |f| .{ .inexact = floatOp(f) }, + }; + } + + fn floorFloat(x: f64) f64 { + return @floor(x); + } + fn ceilFloat(x: f64) f64 { + return @ceil(x); + } + fn roundFloat(x: f64) f64 { + return @round(x); + } + + pub fn floor(allocator: Allocator, a: Number) Error!Number { + return unary(allocator, a, Rational.floor, floorFloat); + } + + pub fn ceil(allocator: Allocator, a: Number) Error!Number { + return unary(allocator, a, Rational.ceil, ceilFloat); + } + + pub fn round(allocator: Allocator, a: Number) Error!Number { + return unary(allocator, a, Rational.round, roundFloat); + } + + fn modFloat(x: f64, y: f64) f64 { + return @mod(x, y); + } + + /// Remainder, taking the sign of the divisor (matching `@mod`). + pub fn mod(allocator: Allocator, a: Number, b: Number) Error!Number { + if (b == .exact and b.exact.isZero()) return Error.DivisionByZero; + if (b == .inexact and b.inexact == 0) return Error.DivisionByZero; + return binary(allocator, a, b, Rational.mod, modFloat); + } + + /// Exact factorial of a non-negative integer. Returns null when the input is + /// not a non-negative exact integer, letting the caller raise a domain error. + pub fn factorial(allocator: Allocator, a: Number) Error!?Number { + const n = a.asExactInt(i64) orelse return null; + if (n < 0) return null; + const result = try Rational.factorial(allocator, @intCast(n)); + return capped(allocator, result); + } + + /// The larger of two values, preserving exactness when both are exact. + pub fn max(allocator: Allocator, a: Number, b: Number) Error!Number { + return if ((try order(allocator, a, b)) == .lt) b.clone() else a.clone(); + } + + /// The smaller of two values, preserving exactness when both are exact. + pub fn min(allocator: Allocator, a: Number, b: Number) Error!Number { + return if ((try order(allocator, a, b)) == .gt) b.clone() else a.clone(); + } + // -- Comparison -- pub fn order(allocator: Allocator, a: Number, b: Number) Error!std.math.Order { @@ -654,3 +719,142 @@ test "clone preserves the tier" { defer d.deinit(); try testing.expect(!d.isExact()); } + +test "floor/ceil/round preserve exactness" { + var a = try Number.parse(alloc, "3.7"); + defer a.deinit(); + + var f = try Number.floor(alloc, a); + defer f.deinit(); + try testing.expect(f.isExact()); + try expectDecimal("3", true, f, 20); + + var c = try Number.ceil(alloc, a); + defer c.deinit(); + try testing.expect(c.isExact()); + try expectDecimal("4", true, c, 20); + + var r = try Number.round(alloc, a); + defer r.deinit(); + try testing.expect(r.isExact()); + try expectDecimal("4", true, r, 20); +} + +test "floor/ceil/round on inexact stay inexact" { + var a = Number.fromFloat(3.7); + defer a.deinit(); + + var f = try Number.floor(alloc, a); + defer f.deinit(); + try testing.expect(!f.isExact()); + try testing.expectEqual(@as(f64, 3.0), f.toFloat(alloc)); + + var c = try Number.ceil(alloc, a); + defer c.deinit(); + try testing.expect(!c.isExact()); + try testing.expectEqual(@as(f64, 4.0), c.toFloat(alloc)); + + var r = try Number.round(alloc, a); + defer r.deinit(); + try testing.expect(!r.isExact()); + try testing.expectEqual(@as(f64, 4.0), r.toFloat(alloc)); +} + +test "mod with an inexact operand stays inexact" { + var a = Number.fromFloat(10.0); + defer a.deinit(); + var b = try Number.fromInt(alloc, 3); + defer b.deinit(); + + var m = try Number.mod(alloc, a, b); + defer m.deinit(); + try testing.expect(!m.isExact()); + try testing.expectEqual(@as(f64, 1.0), m.toFloat(alloc)); + + // And with the inexact value on the right. + var c = try Number.fromInt(alloc, 10); + defer c.deinit(); + var d = Number.fromFloat(3.0); + defer d.deinit(); + var m2 = try Number.mod(alloc, c, d); + defer m2.deinit(); + try testing.expect(!m2.isExact()); + try testing.expectEqual(@as(f64, 1.0), m2.toFloat(alloc)); +} + +test "mod preserves exactness and rejects a zero divisor" { + var a = try Number.fromInt(alloc, 10); + defer a.deinit(); + var b = try Number.fromInt(alloc, 3); + defer b.deinit(); + var m = try Number.mod(alloc, a, b); + defer m.deinit(); + try testing.expect(m.isExact()); + try expectDecimal("1", true, m, 20); + + var zero = try Number.fromInt(alloc, 0); + defer zero.deinit(); + try testing.expectError(Error.DivisionByZero, Number.mod(alloc, a, zero)); + + var fzero = Number.fromFloat(0); + defer fzero.deinit(); + try testing.expectError(Error.DivisionByZero, Number.mod(alloc, a, fzero)); +} + +test "factorial is exact and unbounded" { + var five = try Number.fromInt(alloc, 5); + defer five.deinit(); + var f = (try Number.factorial(alloc, five)).?; + defer f.deinit(); + try testing.expect(f.isExact()); + try expectDecimal("120", true, f, 20); + + // 171! is beyond f64 but fine here. + var big = try Number.fromInt(alloc, 171); + defer big.deinit(); + var bf = (try Number.factorial(alloc, big)).?; + defer bf.deinit(); + try testing.expect(bf.isExact()); +} + +test "factorial rejects non-integers and negatives" { + var frac = try Number.parse(alloc, "2.5"); + defer frac.deinit(); + try testing.expect((try Number.factorial(alloc, frac)) == null); + + var neg = try Number.fromInt(alloc, -1); + defer neg.deinit(); + try testing.expect((try Number.factorial(alloc, neg)) == null); + + var inexact = Number.fromFloat(5); + defer inexact.deinit(); + try testing.expect((try Number.factorial(alloc, inexact)) == null); +} + +test "max and min preserve exactness" { + var a = try Number.parse(alloc, "0.1"); + defer a.deinit(); + var b = try Number.parse(alloc, "0.2"); + defer b.deinit(); + + var hi = try Number.max(alloc, a, b); + defer hi.deinit(); + try testing.expect(hi.isExact()); + try expectDecimal("0.2", true, hi, 20); + + var lo = try Number.min(alloc, a, b); + defer lo.deinit(); + try testing.expect(lo.isExact()); + try expectDecimal("0.1", true, lo, 20); +} + +test "max and min with an inexact operand return that operand as-is" { + var a = try Number.fromInt(alloc, 1); + defer a.deinit(); + var b = Number.fromFloat(2.0); + defer b.deinit(); + var hi = try Number.max(alloc, a, b); + defer hi.deinit(); + try testing.expect(!hi.isExact()); + try testing.expectEqual(@as(f64, 2.0), hi.toFloat(alloc)); +} diff --git a/engine/src/parser.zig b/engine/src/parser.zig index 901dfda..7d69d3b 100644 --- a/engine/src/parser.zig +++ b/engine/src/parser.zig @@ -101,6 +101,7 @@ pub const Parser = struct { .float_value = num.float, .int_value = num.int_value, .base = num.base, + .text = text, } }); }, .string_literal => { diff --git a/engine/src/rational.zig b/engine/src/rational.zig index 232609c..42b27a3 100644 --- a/engine/src/rational.zig +++ b/engine/src/rational.zig @@ -419,6 +419,82 @@ pub const Rational = struct { return try initOwned(allocator, num_root, den_root); } + /// Largest integer not greater than the value. + pub fn floor(allocator: Allocator, a: Rational) Error!Rational { + if (a.isInteger()) return a.clone(); + + var q = try Managed.init(allocator); + errdefer q.deinit(); + var r = try Managed.init(allocator); + defer r.deinit(); + // divFloor rounds the quotient toward negative infinity, which is + // exactly floor for a positive denominator (an invariant here). + try q.divFloor(&r, &a.num, &a.den); + + const den = try Managed.initSet(allocator, 1); + return initOwned(allocator, q, den); + } + + /// Smallest integer not less than the value. + pub fn ceil(allocator: Allocator, a: Rational) Error!Rational { + if (a.isInteger()) return a.clone(); + var f = try floor(allocator, a); + errdefer f.deinit(); + // Not an integer, so ceil is always floor + 1. + try f.num.addScalar(&f.num, 1); + return f; + } + + /// Round to the nearest integer, halves away from zero (matching `@round`). + pub fn round(allocator: Allocator, a: Rational) Error!Rational { + if (a.isInteger()) return a.clone(); + + var half = try initRatio(allocator, 1, 2); + defer half.deinit(); + + if (a.isNegative()) { + var shifted = try sub(allocator, a, half); + defer shifted.deinit(); + return ceil(allocator, shifted); + } + var shifted = try add(allocator, a, half); + defer shifted.deinit(); + return floor(allocator, shifted); + } + + /// Exact integer remainder matching `@mod`: the result takes the sign of + /// the divisor, and equals `a - b * floor(a / b)`. + pub fn mod(allocator: Allocator, a: Rational, b: Rational) Error!Rational { + if (b.isZero()) return Error.DivisionByZero; + + var quotient = try div(allocator, a, b); + defer quotient.deinit(); + var floored = try floor(allocator, quotient); + defer floored.deinit(); + var scaled = try mul(allocator, floored, b); + defer scaled.deinit(); + return sub(allocator, a, scaled); + } + + /// Exact factorial. Unbounded, unlike the f64 version which overflows past + /// 170. + pub fn factorial(allocator: Allocator, n: u64) Error!Rational { + // A guard against absurd inputs that would take effectively forever; + // 20000! is already a ~78000-digit number. + if (n > 20_000) return Error.ExponentTooLarge; + + var acc = try Managed.initSet(allocator, 1); + errdefer acc.deinit(); + var i: u64 = 2; + while (i <= n) : (i += 1) { + var factor = try Managed.initSet(allocator, i); + defer factor.deinit(); + try acc.mul(&acc, &factor); + } + const den = try Managed.initSet(allocator, 1); + return initOwned(allocator, acc, den); + } + // -- Conversion -- /// Convert to the nearest f64. @@ -1196,3 +1272,161 @@ test "chained arithmetic stays exact where f64 would drift" { _ = &y; try testing.expect((x + y) * 10.0 - 3.0 != 0.0); } + +test "floor" { + const cases = [_]struct { in: []const u8, out: []const u8 }{ + .{ .in = "3.7", .out = "3" }, + .{ .in = "3.2", .out = "3" }, + .{ .in = "3", .out = "3" }, + .{ .in = "-3.2", .out = "-4" }, + .{ .in = "-3.7", .out = "-4" }, + .{ .in = "-3", .out = "-3" }, + .{ .in = "0", .out = "0" }, + .{ .in = "0.5", .out = "0" }, + .{ .in = "-0.5", .out = "-1" }, + }; + for (cases) |c| { + var a = try Rational.parseDecimal(alloc, c.in); + defer a.deinit(); + var f = try Rational.floor(alloc, a); + defer f.deinit(); + try expectFrac(c.out, f); + } +} + +test "ceil" { + const cases = [_]struct { in: []const u8, out: []const u8 }{ + .{ .in = "3.2", .out = "4" }, + .{ .in = "3.7", .out = "4" }, + .{ .in = "3", .out = "3" }, + .{ .in = "-3.2", .out = "-3" }, + .{ .in = "-3.7", .out = "-3" }, + .{ .in = "0.5", .out = "1" }, + .{ .in = "-0.5", .out = "0" }, + }; + for (cases) |c| { + var a = try Rational.parseDecimal(alloc, c.in); + defer a.deinit(); + var r = try Rational.ceil(alloc, a); + defer r.deinit(); + try expectFrac(c.out, r); + } +} + +test "round: halves go away from zero, matching @round" { + const cases = [_]struct { in: []const u8, out: []const u8 }{ + .{ .in = "3.5", .out = "4" }, + .{ .in = "3.4", .out = "3" }, + .{ .in = "3.6", .out = "4" }, + .{ .in = "2.5", .out = "3" }, + .{ .in = "-3.5", .out = "-4" }, + .{ .in = "-3.4", .out = "-3" }, + .{ .in = "-2.5", .out = "-3" }, + .{ .in = "7", .out = "7" }, + .{ .in = "0", .out = "0" }, + }; + for (cases) |c| { + var a = try Rational.parseDecimal(alloc, c.in); + defer a.deinit(); + var r = try Rational.round(alloc, a); + defer r.deinit(); + try expectFrac(c.out, r); + // Cross-check against the f64 builtin for the same input. + const f = try std.fmt.parseFloat(f64, c.in); + try testing.expectEqual(@round(f), r.toFloat(alloc)); + } +} + +test "mod: matches the a - b*floor(a/b) definition, sign following the divisor" { + // Expected values are written out rather than derived from `@mod`: Zig's + // float `@mod` with a NEGATIVE divisor gave different answers in different + // builds of this suite (-2 normally, 1 under the instrumented coverage + // build, i.e. comptime folding and the runtime path disagree). An oracle + // that changes with optimize mode cannot verify anything, so these are the + // values the definition requires. + const cases = [_]struct { a: i64, b: i64, expected: i64 }{ + .{ .a = 10, .b = 3, .expected = 1 }, // floor(10/3)=3 -> 10-9 + .{ .a = -10, .b = 3, .expected = 2 }, // floor(-10/3)=-4 -> -10+12 + .{ .a = 10, .b = -3, .expected = -2 }, // floor(10/-3)=-4 -> 10-12 + .{ .a = -10, .b = -3, .expected = -1 }, // floor(-10/-3)=3 -> -10+9 + .{ .a = 7, .b = 7, .expected = 0 }, + .{ .a = 0, .b = 5, .expected = 0 }, + .{ .a = 7, .b = 5, .expected = 2 }, + .{ .a = -7, .b = 5, .expected = 3 }, + }; + for (cases) |c| { + var a = try Rational.initInt(alloc, c.a); + defer a.deinit(); + var b = try Rational.initInt(alloc, c.b); + defer b.deinit(); + var m = try Rational.mod(alloc, a, b); + defer m.deinit(); + + var expected = try Rational.initInt(alloc, c.expected); + defer expected.deinit(); + if (!try Rational.eql(alloc, m, expected)) { + const got = try m.toFractionString(alloc); + defer alloc.free(got); + std.debug.print("mod({d}, {d}): expected {d}, got {s}\n", .{ c.a, c.b, c.expected, got }); + return error.ModMismatch; + } + // The result must carry the sign of the divisor (or be zero). + if (c.expected != 0) { + try testing.expectEqual(c.b < 0, m.isNegative()); + } + } +} + +test "mod: fractional operands" { + var a = try Rational.parseDecimal(alloc, "7.5"); + defer a.deinit(); + var b = try Rational.parseDecimal(alloc, "2"); + defer b.deinit(); + var m = try Rational.mod(alloc, a, b); + defer m.deinit(); + try expectFrac("3/2", m); +} + +test "mod: by zero errors" { + var a = try Rational.initInt(alloc, 1); + defer a.deinit(); + var z = try Rational.initZero(alloc); + defer z.deinit(); + try testing.expectError(Error.DivisionByZero, Rational.mod(alloc, a, z)); +} + +test "factorial: small values" { + const cases = [_]struct { n: u64, out: []const u8 }{ + .{ .n = 0, .out = "1" }, + .{ .n = 1, .out = "1" }, + .{ .n = 5, .out = "120" }, + .{ .n = 10, .out = "3628800" }, + }; + for (cases) |c| { + var f = try Rational.factorial(alloc, c.n); + defer f.deinit(); + try expectFrac(c.out, f); + } +} + +test "factorial: unbounded past the f64 limit of 170" { + // 171! overflows f64 to infinity; exactly this is why the old evaluator + // rejected it. Exact arithmetic has no such wall. + var f = try Rational.factorial(alloc, 171); + defer f.deinit(); + try testing.expect(f.isInteger()); + const s = try f.toFractionString(alloc); + defer alloc.free(s); + try testing.expect(s.len > 300); // 171! has 310 digits + try testing.expect(std.math.isPositiveInf(f.toFloat(alloc))); +} + +test "factorial: 20 is exact where f64 is already lossy" { + var f = try Rational.factorial(alloc, 20); + defer f.deinit(); + try expectFrac("2432902008176640000", f); +} + +test "factorial: absurd inputs are rejected" { + try testing.expectError(Error.ExponentTooLarge, Rational.factorial(alloc, 20_001)); +}