Rounding from the one mode the language has, and sqrt from libm

floor, ceil and round over f32, which is what a position and a tile coordinate
are here. The only rounding mode available is the cast's truncation toward
zero, so each of these is that cast plus the correction the mode does not
make, and the content is which inputs make the cast itself undefined. NaN
fails every comparison, so it needs its own (not (= x x)) and nothing else
finds it; the infinities fall out of the magnitude test; and above 2^23 an f32
has no fractional bits left, which makes returning the input there the exact
answer and also the guard that keeps the cast inside i32.

round is half away from zero, written as floor of the magnitude and mirrored.
The obvious (floor-f32 (+ x 0.5)) is wrong twice: half-up rather than
half-away, so -2.5 comes out -2, and at the largest f32 below 0.5 the addition
alone rounds to 1.0 and answers 1 for a number under a half. Both are in the
table, which is why every case there is a negative or a half.

sqrt is the decision in this commit and it goes out to libm, which is a change
to the release link and so is said out loud. Every other number in the prelude
is reachable from the four operations and a cast; a square root is not.
Newton's method needs a starting guess, the good guess comes from
reinterpreting the exponent bits, and the language has only value-preserving
casts - no bit-cast between f32 and u32. Without one the iteration needs a
scaling loop to normalise and still produces a result that is merely close,
which is the one thing 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; for this function the byte-identical argument points at C rather
than away from it.

The cost is -lm on every link, and its placement matters. It goes after the
objects, not in the leading flags, because --as-needed drops a library named
before the object that wants it. Worse, at -O2 LLVM folds most sqrtf calls
into the hardware instruction and the symbol never has to resolve - so this
looked linked before the flag existed and failed only at -O0, which is exactly
why the table runs both. Untested against --target=wasm32: wasi-libc ships
libm.a as a stub because the symbols live in libc, so it should be inert
there, but nothing here exercises it.

The better fix is not in this lane. llvm.sqrt.f32 as a builtin in check.ml and
emit.ml is one instruction, no symbol and no flag, and it belongs to whoever
owns the compiler.

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

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

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)

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