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,