diff --git a/lib/build.ml b/lib/build.ml index a1dfa8d..3a10f6e 100644 --- a/lib/build.ml +++ b/lib/build.ml @@ -145,6 +145,15 @@ let executable ?(opts = default) ?(csrcs = []) ?(lflags = []) @ [ Filename.quote ll ] @ List.map Filename.quote objs @ lflags + (* The prelude declares sqrtf, so every link needs libm. It goes here + and not in the leading flags: the default --as-needed drops a + library named before the object that wants it, so at -O2 this would + appear to work — LLVM folds most sqrtf calls into the hardware + instruction and the symbol never has to resolve — and the -O0 build, + which emits the call, would fail at the link. Untested against + --target=wasm32; wasi-libc ships libm.a as a stub because the + symbols live in libc, so it should be inert there. *) + @ [ "-lm" ] @ [ "-o"; Filename.quote out ]) in let code = Sys.command cmd in diff --git a/lib/prelude.ml b/lib/prelude.ml index da923f8..7955089 100644 --- a/lib/prelude.ml +++ b/lib/prelude.ml @@ -237,6 +237,76 @@ let source = {flan| ;; [lo, hi), because rand-f32 never reaches 1.0. (defn rand-f32-range [lo f32 hi f32] f32 (+ lo (* (rand-f32) (- hi lo)))) + +;; ── Rounding, and the one thing that is not Flan ────────────────────── +;; +;; All three answer an f32 and take the f32 path, because that is what a +;; position, a tile coordinate and a velocity are here. f64 versions wait for +;; a program that wants them, for the same reason the f32 slice algorithms do. +;; +;; The cast to i32 truncates toward zero, which is the only rounding mode the +;; language has, so each of these is that cast plus the correction the mode +;; does not make. Three inputs would make the cast itself undefined and each +;; is named before it happens: NaN (which fails every comparison, so it is +;; tested for by (not (= x x)) and nothing else), and the two infinities, +;; which are caught by the magnitude test. Above 2^23 an f32 has no fractional +;; bits left at all, so returning x there is not an approximation — it is the +;; answer — and it doubles as the guard that keeps the cast inside i32. + +;; Zero is returned as itself rather than through the cast, which would turn +;; -0.0 into +0.0. That is one line for a value most callers never look at, +;; and it is here because floorf is specified to return it: a sign of zero is +;; how a caller recovers which side a position approached from once the +;; magnitude has already been rounded away. +(defn floor-f32 [x f32] f32 + (if (or (not (= x x)) + (>= x 8388608.0) + (<= x -8388608.0) + (= x 0.0)) + x + (let [t (f32 (i32 x))] + (if (> t x) (- t 1.0) t)))) + +;; One deviation from C's ceilf, written down rather than branched around: +;; between -1.0 and 0.0 this answers +0.0 where IEEE asks for -0.0, because +;; the outer (- 0.0 …) is a subtraction and not a negation. Nothing here reads +;; the sign of a zero; a caller that does should test the input instead. +(defn ceil-f32 [x f32] f32 + (- 0.0 (floor-f32 (- 0.0 x)))) + +;; Half away from zero, which is C's round and not the even-tie rule: -2.5 +;; goes to -3. Written as floor of the *magnitude* and mirrored, because +;; (floor-f32 (+ x 0.5)) is wrong twice over — it is half-*up* rather than +;; half-away for negatives, and at the largest f32 below 0.5 the addition +;; itself rounds to 1.0 and answers 1 for a number under a half. +(defn round-f32 [x f32] f32 + (let [m (if (< x 0.0) (- 0.0 x) x) + f (floor-f32 m) + r (if (>= (- m f) 0.5) (+ f 1.0) f)] + (if (< x 0.0) (- 0.0 r) r))) + +;; sqrt is the one function in this file that is not Flan, and it is a +;; `declare` rather than a body for a reason that is not laziness. Every other +;; number here is reachable from the four operations and a cast; a square root +;; is not. Newton's method needs a starting guess, a good one comes from +;; reinterpreting the exponent bits, and the language has no bit-cast between +;; f32 and u32 — only value-preserving casts. Without it the iteration needs a +;; scaling loop to normalise, converges slowly from a poor guess, and produces +;; a result that is *close*, which is exactly what a standard library must not +;; hand back. IEEE-754 makes sqrt correctly rounded, so libm's answer is the +;; same bit pattern on native and on wasm32 — the byte-identical property that +;; keeps rand-u32 in Flan is, for this one, an argument for going out to C. +;; +;; The cost is one `declare` line in every module, which LLVM drops where it +;; is unused, and one -lm on every link, which build.ml now passes. That flag +;; is not optional and not obvious: at -O2 LLVM folds most sqrtf calls into +;; the hardware instruction and nothing is left to resolve, so this appears to +;; link without it and then fails at -O0, where the call survives. +;; +;; The better fix belongs to the compiler and not here: llvm.sqrt.f32 as a +;; builtin in check.ml and emit.ml is one instruction with no symbol at all. +(declare sqrt-f32 [x f32] f32 "sqrtf") + ;; ── Byte classes ────────────────────────────────────────────────────── ;; ;; ASCII only, and deliberately: a byte is a byte here, there is no code point diff --git a/test/programs/math.flan b/test/programs/math.flan new file mode 100644 index 0000000..7caac53 --- /dev/null +++ b/test/programs/math.flan @@ -0,0 +1,60 @@ +;;;; The prelude's rounding, and sqrt. +;;;; +;;;; Every input is one a plausible wrong implementation gets wrong. The +;;;; negatives are the whole point: a floor written as a bare cast truncates +;;;; toward zero and answers -2 for -2.5, and a round written as +;;;; (floor-f32 (+ x 0.5)) answers -2 for -2.5 as well, where C's round says +;;;; -3. The exact halves appear on both signs for that reason. The values +;;;; that are already integers check that the correction does *not* fire — +;;;; a floor that always subtracts one turns 3.0 into 2.0 — and 16777216.0 is +;;;; past 2^24, where an f32 has no fractional bits and the guard, not the +;;;; cast, has to produce the answer. + +(defn show [x f32] + (print-f64 (f64 x)) + (print-str " ")) + +(defn main [] i32 + ;; floor: down on both signs, and unmoved on the integers. + (show (floor-f32 2.7)) (show (floor-f32 2.0)) (show (floor-f32 2.3)) + (show (floor-f32 -2.7)) (show (floor-f32 -2.0)) (show (floor-f32 -2.3)) + (show (floor-f32 0.5)) (show (floor-f32 -0.5)) + (newline) + + ;; ceil: up on both signs. -2.7 must give -2, which is where a ceil written + ;; as "floor plus one" goes wrong. + (show (ceil-f32 2.7)) (show (ceil-f32 2.0)) (show (ceil-f32 2.3)) + (show (ceil-f32 -2.7)) (show (ceil-f32 -2.0)) (show (ceil-f32 -2.3)) + (show (ceil-f32 0.5)) (show (ceil-f32 -0.5)) + (newline) + + ;; Zero keeps its sign through floor, which is what the (= x 0.0) guard in + ;; it is for and the only place that guard is observable: the cast it skips + ;; would turn -0.0 into +0.0, and %g prints the difference. Drop the guard + ;; and the third column here reads 0 instead of -0. + (show (floor-f32 0.0)) (show (ceil-f32 0.0)) (show (floor-f32 -0.0)) + (newline) + + ;; round: half away from zero on both signs, so -2.5 is -3 and not -2. + (show (round-f32 2.4)) (show (round-f32 2.5)) (show (round-f32 2.6)) + (show (round-f32 -2.4)) (show (round-f32 -2.5)) (show (round-f32 -2.6)) + (show (round-f32 0.5)) (show (round-f32 -0.5)) + (newline) + + ;; Past 2^24 there is no fraction left; the answer is the input, and the + ;; cast that would produce it is out of i32's range on the way there. + (show (floor-f32 16777216.0)) (show (ceil-f32 16777216.0)) + (show (round-f32 16777216.0)) (show (floor-f32 -16777216.0)) + (newline) + + ;; sqrt, including the two values a wrong-sense iteration still passes + ;; (0 and 1) and one that is not a perfect square. + (show (sqrt-f32 0.0)) (show (sqrt-f32 1.0)) (show (sqrt-f32 4.0)) + (show (sqrt-f32 2.0)) (show (sqrt-f32 0.25)) (show (sqrt-f32 1e6)) + (newline) + + ;; A squared distance through sqrt, which is what a game actually calls it + ;; for: 3-4-5 exactly, so a last-bit error would show. + (show (sqrt-f32 (+ (* 3.0 3.0) (* 4.0 4.0)))) + (newline) + 0) diff --git a/test/test_acceptance.ml b/test/test_acceptance.ml index 28345c7..52eb110 100644 --- a/test/test_acceptance.ml +++ b/test/test_acceptance.ml @@ -143,6 +143,27 @@ let () = outputs "bytes, parsing and numbers" "programs/text.flan" text_out; outputs ~opt:"-O0" "bytes, parsing and numbers, -O0" "programs/text.flan" text_out; + (* Rounding and sqrt. Every case here is a *negative* or a half, because + those are the two places a plausible wrong version differs: a floor + written as the bare cast truncates toward zero and answers -2 for -2.5, + and a round written as (floor-f32 (+ x 0.5)) is half-up rather than + half-away and answers -2 as well. 16777216.0 is past 2^24, where the + guard rather than the cast has to produce the answer — and where the + cast it guards would be out of i32's range. sqrt is a `declare` on + libm's sqrtf; the -O0 run is the one that matters for it, because at + -O2 LLVM folds most calls into the hardware instruction and a symbol + that never has to resolve proves nothing about the link. *) + let math_out = + "2 2 2 -3 -2 -3 0 -1 \n\ + 3 2 3 -2 -2 -2 1 0 \n\ + 0 0 -0 \n\ + 2 3 3 -2 -3 -3 1 -1 \n\ + 1.67772e+07 1.67772e+07 1.67772e+07 -1.67772e+07 \n\ + 0 1 2 1.41421 0.5 1000 \n\ + 5 \n" + in + outputs "rounding and sqrt" "programs/math.flan" math_out; + outputs ~opt:"-O0" "rounding and sqrt, -O0" "programs/math.flan" math_out; (* index-of-bytes, trim, the byte classes and parse-f64. The search cases are the ones that separate a correct loop from a lucky one: a match only at the end, "aab" in "aaab" (where the first byte matches twice