From aeeb6de59dd8a608f5f9fe259d1bb429a93a6601 Mon Sep 17 00:00:00 2001 From: Joseph Ferano Date: Thu, 17 Sep 2026 21:44:39 +0700 Subject: [PATCH] The maths is the whole family now, at both widths The prelude's declare surface was five f32 functions, and the five were there because somebody needed each one. Everything else a caller wanted was written as a declare at the top of their own file -- the identical libm call with none of the caveats written down. So the rest of libm is here: tan, the three inverses, the three logarithms, exp, fmod, hypot, cbrt, fabs, and an f64 face for every one of them including the five that already existed. A declare is a line, a symbol already on the link, and nothing in either backend, which is why this was cheap enough to do completely rather than one function at a time. The f64 half is not decoration. f32 is what a position is; f64 is what a measurement is -- the clock, parse-f64, format-f64, any sum over more than a few thousand terms -- and having only the f32 face forced a cast down and back at each of those boundaries, which is where the precision went. The paragraph the sqrt note draws for itself is now drawn once for the family: IEEE-754 specifies sqrt, fabs, floor, ceil, round and fmod as exact or correctly rounded, so those agree bit for bit across glibc, musl and wasi-libc; it requires nothing of the rest, so the sand-grid rule covers all of them unchanged. floor, ceil and round are Flan at f32 and libm at f64, and that is not an inconsistency: the f32 bodies work because every f32 with a fraction fits in an i32, and at f64 that trick is gone. abs-i32 and abs-i64 are Flan, one per width because min and max are builtins and no generic covers the numeric types. pi and tau at both widths, written out rather than derived so the compiler rounds each literal once. programs/math3.flan covers it at values that are exact in binary, so nothing pins one libm's last bit. The -O0 case is the one that matters: at -O2 LLVM folds a call over two literals and leaves no symbol to resolve, which is how a missing -lm hid the first time. --- lib/prelude.ml | 147 +++++++++++++++++++++++++++++++++++++++ test/programs/math3.flan | 92 ++++++++++++++++++++++++ test/test_acceptance.ml | 17 +++++ web/index.html | 34 +++++---- 4 files changed, 276 insertions(+), 14 deletions(-) create mode 100644 test/programs/math3.flan diff --git a/lib/prelude.ml b/lib/prelude.ml index 3ed65e6..353e69f 100644 --- a/lib/prelude.ml +++ b/lib/prelude.ml @@ -676,6 +676,153 @@ let source = {flan| (declare atan2-f32 [y f32 x f32] f32 "atan2f") (declare pow-f32 [x f32 y f32] f32 "powf") +;; ── The rest of libm, and both widths ───────────────────────────────── +;; +;; The five above were the whole of it for a long time, and the reason they +;; were is the reason the rest are here now: every one of these is a line, a +;; symbol that is already on the link, and nothing in the compiler. A program +;; that wanted a logarithm wrote the `declare` at the top of its own file — +;; which is the identical call with none of the caveats written down and +;; nobody's name on it. +;; +;; **The f64 half is not decoration.** f32 is what a position and a colour +;; are, and f64 is what a *measurement* is: the clock below is nanoseconds in +;; an i64 and seconds in an f64, parse-f64 and format-f64 are both f64, and a +;; sum over more than a few thousand f32 terms has already lost the low bits +;; the answer was about. Having only the f32 face forced a cast down and back +;; at every one of those boundaries, and a cast down is where the precision +;; went. +;; +;; The split below is the one the sqrt paragraph draws, applied to the whole +;; family, and it is the only thing here worth knowing before calling: +;; +;; **Exact on every target.** IEEE-754 specifies these as exact operations +;; or as correctly rounded, so the answer is the same bit pattern under +;; glibc, musl and wasi-libc, and a hash taken across targets may be routed +;; through them. sqrt, fabs, floor, ceil, round, fmod. +;; +;; **Not.** IEEE-754 requires nothing of these and the three libms do +;; differ in the last bit. The sand-grid rule from the sin/cos paragraph +;; above covers all of them without change: a value compared across targets +;; must not have been through one. Everything else here. + +(declare sqrt-f64 [x f64] f64 "sqrt") + +;; Magnitude, and the f32/f64 pair is libm's because fabs is a sign-bit clear +;; that the compiler folds into one instruction — cheaper than the branch a +;; Flan body would be, and right for -0.0 and for NaN, which a (< x 0.0) test +;; is not: -0.0 is not less than zero, so the branch returns it unchanged and +;; hands back a negative zero from a function named abs. +(declare abs-f32 [x f32] f32 "fabsf") +(declare abs-f64 [x f64] f64 "fabs") + +;; The f64 faces of the three rounding functions floor-f32, ceil-f32 and +;; round-f32 are libm's rather than Flan's, and that is not an inconsistency. +;; Those three are Flan because of a cast: (i32 x) is the whole of floor-f32's +;; body, and it works precisely because every f32 with a fractional part fits +;; in an i32. At f64 it does not — the exact range runs to 2^53 and i64's cast +;; would have to carry its own guard — so the trick that made them free is not +;; available and the libm call is both shorter and exact. +(declare floor-f64 [x f64] f64 "floor") +(declare ceil-f64 [x f64] f64 "ceil") +(declare round-f64 [x f64] f64 "round") + +;; Remainder, and it is C's fmod and not a modulo: the sign follows the +;; *dividend*, so (fmod-f32 -1.0 3.0) is -1.0 and not 2.0. An angle wrapped +;; into [0, tau) therefore needs the add-and-fmod-again that every wrap +;; function has, and this is the line where that is written down rather than +;; discovered. It is exact — the result is the true remainder, representable +;; by construction — so it belongs to the first group above. +(declare fmod-f32 [x f32 y f32] f32 "fmodf") +(declare fmod-f64 [x f64 y f64] f64 "fmod") + +;; The trigonometric family, in full and at both widths. tan is separate from +;; (/ (sin-f32 x) (cos-f32 x)) for the reason atan2 is separate from a +;; division: near pi/2 the quotient is a ratio of two small errors and tanf +;; is not. +(declare tan-f32 [x f32] f32 "tanf") +(declare sin-f64 [x f64] f64 "sin") +(declare cos-f64 [x f64] f64 "cos") +(declare tan-f64 [x f64] f64 "tan") + +;; The inverses. asin and acos answer NaN outside [-1, 1] rather than +;; clamping, which is what catches a dot product that drifted to 1.0000001 +;; through rounding — clamp it at the call site, on purpose, and the drift is +;; visible instead of silently becoming an angle of zero. +(declare asin-f32 [x f32] f32 "asinf") +(declare acos-f32 [x f32] f32 "acosf") +(declare atan-f32 [x f32] f32 "atanf") +(declare asin-f64 [x f64] f64 "asin") +(declare acos-f64 [x f64] f64 "acos") +(declare atan-f64 [x f64] f64 "atan") +(declare atan2-f64 [y f64 x f64] f64 "atan2") + +;; Logarithms and the exponential. log is the natural one, as in C and unlike +;; the spreadsheet convention — log2 and log10 are the other two and are named +;; for their bases, so nothing here is ambiguous. log2 is not (/ (log x) +;; (log 2.0)): it is exact at every power of two, which is the whole reason a +;; bit-width or an octave is computed with it. +;; +;; All four answer -inf at zero and NaN below it rather than signalling. A +;; condition per logarithm would cost a handler search on a path whose callers +;; are loops over samples, and NaN is the value that propagates to wherever +;; the caller does check. +(declare log-f32 [x f32] f32 "logf") +(declare log2-f32 [x f32] f32 "log2f") +(declare log10-f32 [x f32] f32 "log10f") +(declare exp-f32 [x f32] f32 "expf") +(declare log-f64 [x f64] f64 "log") +(declare log2-f64 [x f64] f64 "log2") +(declare log10-f64 [x f64] f64 "log10") +(declare exp-f64 [x f64] f64 "exp") +(declare pow-f64 [x f64 y f64] f64 "pow") + +;; hypot over (sqrt-f32 (+ (* x x) (* y y))) because the obvious form +;; overflows on inputs the answer does not: the square of an f32 above ~1.8e19 +;; is infinity, so a distance between two far-apart points comes back inf when +;; the distance itself is perfectly representable. hypotf scales first. It +;; costs more than the naive form and is worth it exactly when the naive form +;; is wrong. +(declare hypot-f32 [x f32 y f32] f32 "hypotf") +(declare hypot-f64 [x f64 y f64] f64 "hypot") + +;; Cube root, and it is here because (pow-f32 x 0.33333334) is not it: pow +;; goes through a logarithm, which is undefined for a negative base, so the +;; obvious spelling answers NaN for every negative number where cbrt answers +;; the negative root. +(declare cbrt-f32 [x f32] f32 "cbrtf") +(declare cbrt-f64 [x f64] f64 "cbrt") + +;; Integer magnitude, one per width because there are no generics over the +;; numeric types and min and max are builtins rather than functions, so a +;; single abs is not expressible today. +;; +;; The most negative value of each width has no positive counterpart, and this +;; does not special-case it: the subtraction is the same subtraction written +;; anywhere else and meets whatever the build's overflow rule is. Saturating +;; to the maximum would be a wrong answer returned quietly, which is the one +;; thing this file does not do. +(defn abs-i32 [x i32] i32 + (if (< x 0) (- 0 x) x)) + +(defn abs-i64 [x i64] i64 + (if (< x 0) (- 0 x) x)) + +;; pi and tau at both widths, because a defconst has a type and a cast between +;; them is where digits go missing. tau is 2pi and is written out rather than +;; multiplied, so the f32 one is the nearest f32 to tau and not twice the +;; nearest f32 to pi — which is the same number here and is not guaranteed to +;; be for the derived form in general. +;; +;; Both are given to more digits than either width holds. That is deliberate: +;; the literal is rounded once, by the compiler, to the nearest value of the +;; declared type, which is the best available answer and is the same answer on +;; both targets. +(defconst pi-f32 f32 3.14159265358979323846) +(defconst pi-f64 f64 3.14159265358979323846) +(defconst tau-f32 f32 6.28318530717958647692) +(defconst tau-f64 f64 6.28318530717958647692) + ;; ── Byte classes ────────────────────────────────────────────────────── ;; ;; ASCII only, and deliberately: a byte is a byte here, there is no code point diff --git a/test/programs/math3.flan b/test/programs/math3.flan new file mode 100644 index 0000000..7f8a381 --- /dev/null +++ b/test/programs/math3.flan @@ -0,0 +1,92 @@ +;;;; The rest of libm, at both widths — what math.flan and math2.flan left out. +;;;; +;;;; The rule those two set holds here unchanged and is the only reason this +;;;; file looks the way it does: none of these is correctly rounded under +;;;; IEEE-754 except sqrt, fabs, the rounding three and fmod, so every value +;;;; below is one whose answer is exact in binary — zero, one, a power of two, +;;;; a perfect square, a perfect cube. A case that pinned glibc's last bit +;;;; would pass here and fail on wasi-libc. +;;;; +;;;; The acceptance table builds this at -O0 as well, and that run is the one +;;;; that matters: at -O2 LLVM folds a libm call over two literals and leaves +;;;; no symbol to resolve, which is how a missing -lm hid the first time. + +(defn show [x f32] () + (print x) + (print " ")) + +(defn show64 [x f64] () + (print x) + (print " ")) + +(defn main [] i32 + ;; The f32 half. tan, the three inverses, the three logarithms and exp. + (show (tan-f32 0.0)) ; 0 + (show (asin-f32 0.0)) ; 0 + (show (acos-f32 1.0)) ; 0 + (show (atan-f32 0.0)) ; 0 + (show (log-f32 1.0)) ; 0 + (show (log2-f32 8.0)) ; 3 + (show (log10-f32 1000.0)) ; 3 + (show (exp-f32 0.0)) ; 1 + (println "") + + (show (fmod-f32 7.0 4.0)) ; 3 + ;; The sign follows the dividend and not the divisor, which is the line a + ;; caller reaching for a modulo gets wrong. Written down as a case. + (show (fmod-f32 -1.0 3.0)) ; -1 + (show (hypot-f32 3.0 4.0)) ; 5 + (show (cbrt-f32 27.0)) ; 3 + ;; The negative root, where (pow-f32 x 0.33333334) would be NaN: pow goes + ;; through a logarithm and cbrt does not. + (show (cbrt-f32 -8.0)) ; -2 + (show (abs-f32 -2.5)) ; 2.5 + (println "") + + ;; The f64 half, at the same exact values. This block is the whole point of + ;; the f64 face existing: before it, every one of these was a cast down to + ;; f32 and back, and the cast down is where the precision went. + (show64 (sqrt-f64 16.0)) ; 4 + (show64 (sin-f64 0.0)) ; 0 + (show64 (cos-f64 0.0)) ; 1 + (show64 (tan-f64 0.0)) ; 0 + (show64 (asin-f64 0.0)) ; 0 + (show64 (acos-f64 1.0)) ; 0 + (show64 (atan-f64 0.0)) ; 0 + (show64 (atan2-f64 0.0 1.0)) ; 0 + (println "") + + (show64 (log-f64 1.0)) ; 0 + (show64 (log2-f64 1024.0)) ; 10 + (show64 (log10-f64 100.0)) ; 2 + (show64 (exp-f64 0.0)) ; 1 + (show64 (pow-f64 2.0 10.0)) ; 1024 + (show64 (fmod-f64 7.0 4.0)) ; 3 + (show64 (hypot-f64 3.0 4.0)) ; 5 + (show64 (cbrt-f64 8.0)) ; 2 + (show64 (abs-f64 -1.5)) ; 1.5 + (println "") + + ;; The f64 rounding family, which is libm's where the f32 one is Flan's — + ;; and it agrees with the Flan one where they overlap: half away from zero, + ;; so -2.5 goes to -3 and not to -2. + (show64 (floor-f64 -2.5)) ; -3 + (show64 (ceil-f64 -2.5)) ; -2 + (show64 (round-f64 -2.5)) ; -3 + (show64 (round-f64 2.5)) ; 3 + (println "") + + ;; Integer magnitude, one per width. + (print (abs-i32 -7)) (print " ") ; 7 + (print (abs-i64 (i64 -7))) (print " ") ; 7 + (print (abs-i32 7)) (print " ") ; 7 + ;; tau is 2pi at both widths. Pinning the relation rather than the digits is + ;; what catches a constant written to too few of them. + (print (= tau-f32 (* 2.0 pi-f32))) (print " ") + (print (= tau-f64 (* 2.0 pi-f64))) + (println "") + + ;; pi is the one value here that can be pinned without pinning a libm: it is + ;; a literal the compiler rounds, so it is the same on every target. + (println (and (> pi-f64 3.14159265) (< pi-f64 3.14159266))) + 0) diff --git a/test/test_acceptance.ml b/test/test_acceptance.ml index a99ab31..4457ba1 100644 --- a/test/test_acceptance.ml +++ b/test/test_acceptance.ml @@ -268,6 +268,23 @@ let () = outputs "atan2, pow and clamp" "programs/math2.flan" math2_out; outputs ~opt:"-O0" "atan2, pow and clamp, -O0" "programs/math2.flan" math2_out; + (* The rest of libm, at both widths. Same rule as math2 above and for the + same reason — every value is exact in binary — and the -O0 pass is + doing the same job: at -O2 LLVM folds a libm call over two literals and + leaves no symbol to resolve, so that run is the one proving all thirty + new declares actually link. *) + let math3_out = + "0 0 0 0 0 3 3 1 \n\ + 3 -1 5 3 -2 2.5 \n\ + 4 0 1 0 0 0 0 0 \n\ + 0 10 2 1 1024 3 5 2 1.5 \n\ + -3 -2 -3 3 \n\ + 7 7 7 true true\n\ + true\n" + in + outputs "the rest of libm, both widths" "programs/math3.flan" math3_out; + outputs ~opt:"-O0" "the rest of libm, both widths, -O0" + "programs/math3.flan" math3_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 diff --git a/web/index.html b/web/index.html index a4ab6ec..5536c84 100644 --- a/web/index.html +++ b/web/index.html @@ -1047,7 +1047,7 @@ takes the value as it is and prints the number it holds.

The prelude

-

The prelude is written in Flan, all but five lines of it, and prepended to every +

The prelude is written in Flan, all but its declare lines, and prepended to every program, so nothing in it needs importing. It holds no printing of its own: print and println are the compiler's, and write-stdout — the one output primitive — is what they are written @@ -1063,7 +1063,7 @@ over.

textsplit-on-byte, split-next!, split, lower-ascii, upper-ascii, to-lower, to-upper building bytesappend!, append-i64!, append-f64!, concat, join, repeat-bytes, replace-bytes, slices-new, format-f64 UTF-8decode-rune, rune-at, rune-count, rune-size, rune-start?, valid-utf8?, encode-rune! -numberssign-f32, lerp, clamp, floor-f32, ceil-f32, round-f32, and the five declares: sqrt-f32, sin-f32, cos-f32, atan2-f32, pow-f32 +numberssign-f32, lerp, clamp, floor-f32, ceil-f32, round-f32, abs-i32, abs-i64, the constants pi-f32, pi-f64, tau-f32, tau-f64, and libm through a declare at both widths: sqrt, abs, floor, ceil, round, fmod, sin, cos, tan, asin, acos, atan, atan2, log, log2, log10, exp, pow, hypot, cbrt — each spelled -f32 or -f64 randomrand-seed, rand-u32, rand-f32, rand-i32-range, rand-f32-range forms, for macrosform-nil, form-cons, form-append, form-rest, form-items, form-pair, form-sym?, form-is-sym?, gensym, and unless and into, which are macros written here rather than special forms the restpause, which signals the Pause condition the break loop stops on, and embed-find @@ -1084,18 +1084,24 @@ native and on wasm32. The parsers are ours too: "abc" and 12 for "12x", which are three wrong answers a caller cannot tell from a real 12.

-

Five functions in the file are not Flan, and they are libm's: -(declare sqrt-f32 [x f32] f32 "sqrtf") and the same line for -sinf, cosf, atan2f and powf. Every -other number here is reachable from the four operations and a cast; a square root is -not, and the usual trick of seeding Newton's method from the exponent bits needs a -bit-cast between f32 and u32 that the language does not have. -IEEE-754 makes sqrt correctly rounded, so libm gives the same bit pattern -on both targets anyway. The other four are not: IEEE-754 requires -nothing of sinf, cosf, atan2f or -powf, and glibc, musl and wasi-libc do differ in the last bit — so the -byte-identical-hash property the RNG exists for does not survive a hash routed through -any of them. Every link carries -lm.

+

The maths in the file is not Flan, and it is libm's: +(declare sqrt-f32 [x f32] f32 "sqrtf") and the same line for thirty-odd +more. Every other number here is reachable from the four operations and a cast; a square +root is not, and the usual trick of seeding Newton's method from the exponent bits needs +a bit-cast between f32 and u32 that the language does not have. +A declare is also the cheapest thing in the language to add — a line, a +symbol already on the link, and nothing in either backend — which is why the surface is +now the whole family at both widths rather than the five it started as.

+ +

One split is worth knowing before calling any of them. IEEE-754 +specifies sqrt, fabs, floor, ceil, +round and fmod as exact or correctly rounded, so those give the +same bit pattern under glibc, musl and wasi-libc. It requires nothing of the +rest — sin, cos, tan, the inverses, the +logarithms, exp, pow, hypot, cbrt — +and the three libms do differ in the last bit, so the byte-identical-hash property the +RNG exists for does not survive a value routed through any of them. Every link carries +-lm.

The primitives underneath are few — a primitive is the only thing implemented twice per backend: argv,