JosephFerano/flan
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

The maths is the whole family now, at both widths

Browse Source 
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.
This commit is contained in:
Joseph Ferano 2026-09-17 21:44:39 +07:00
parent 8a430bb50c
commit aeeb6de59d
4 changed files with 276 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

147
lib/prelude.ml
View File

@ -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

92
test/programs/math3.flan Normal file
View File

@ -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)

17
test/test_acceptance.ml
View File

@ -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

34
web/index.html
View File

@ -1047,7 +1047,7 @@ takes the value as it is and prints the number it holds.</p>
<h2 id="prelude">The prelude</h2>
<p>The prelude is written in Flan, all but five lines of it, and prepended to every
<p>The prelude is written in Flan, all but its <code>declare</code> lines, and prepended to every
program, so nothing in it needs importing. It holds no printing of its own:
<code>print</code> and <code>println</code> are the compiler's, and
<code>write-stdout</code> — the one output primitive — is what they are written
@ -1063,7 +1063,7 @@ over.</p>
<tr><td>text</td><td><code>split-on-byte</code>, <code>split-next!</code>, <code>split</code>, <code>lower-ascii</code>, <code>upper-ascii</code>, <code>to-lower</code>, <code>to-upper</code></td></tr>
<tr><td>building bytes</td><td><code>append!</code>, <code>append-i64!</code>, <code>append-f64!</code>, <code>concat</code>, <code>join</code>, <code>repeat-bytes</code>, <code>replace-bytes</code>, <code>slices-new</code>, <code>format-f64</code></td></tr>
<tr><td>UTF-8</td><td><code>decode-rune</code>, <code>rune-at</code>, <code>rune-count</code>, <code>rune-size</code>, <code>rune-start?</code>, <code>valid-utf8?</code>, <code>encode-rune!</code></td></tr>
<tr><td>numbers</td><td><code>sign-f32</code>, <code>lerp</code>, <code>clamp</code>, <code>floor-f32</code>, <code>ceil-f32</code>, <code>round-f32</code>, and the five <code>declare</code>s: <code>sqrt-f32</code>, <code>sin-f32</code>, <code>cos-f32</code>, <code>atan2-f32</code>, <code>pow-f32</code></td></tr>
<tr><td>numbers</td><td><code>sign-f32</code>, <code>lerp</code>, <code>clamp</code>, <code>floor-f32</code>, <code>ceil-f32</code>, <code>round-f32</code>, <code>abs-i32</code>, <code>abs-i64</code>, the constants <code>pi-f32</code>, <code>pi-f64</code>, <code>tau-f32</code>, <code>tau-f64</code>, and libm through a <code>declare</code> at both widths: <code>sqrt</code>, <code>abs</code>, <code>floor</code>, <code>ceil</code>, <code>round</code>, <code>fmod</code>, <code>sin</code>, <code>cos</code>, <code>tan</code>, <code>asin</code>, <code>acos</code>, <code>atan</code>, <code>atan2</code>, <code>log</code>, <code>log2</code>, <code>log10</code>, <code>exp</code>, <code>pow</code>, <code>hypot</code>, <code>cbrt</code> — each spelled <code>-f32</code> or <code>-f64</code></td></tr>
<tr><td>random</td><td><code>rand-seed</code>, <code>rand-u32</code>, <code>rand-f32</code>, <code>rand-i32-range</code>, <code>rand-f32-range</code></td></tr>
<tr><td>forms, for macros</td><td><code>form-nil</code>, <code>form-cons</code>, <code>form-append</code>, <code>form-rest</code>, <code>form-items</code>, <code>form-pair</code>, <code>form-sym?</code>, <code>form-is-sym?</code>, <code>gensym</code>, and <code>unless</code> and <code>into</code>, which are macros written here rather than special forms</td></tr>
<tr><td>the rest</td><td><code>pause</code>, which signals the <code>Pause</code> condition the break loop stops on, and <code>embed-find</code></td></tr>
@ -1084,18 +1084,24 @@ native and on wasm32. <strong>The parsers are ours too</strong>:
<code>"abc"</code> and 12 for <code>"12x"</code>, which are three wrong answers a caller
cannot tell from a real 12.</p>
<p>Five functions in the file are not Flan, and they are libm's:
<code>(declare sqrt-f32 [x f32] f32 "sqrtf")</code> and the same line for
<code>sinf</code>, <code>cosf</code>, <code>atan2f</code> and <code>powf</code>. 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 <code>f32</code> and <code>u32</code> that the language does not have.
IEEE-754 makes <code>sqrt</code> correctly rounded, so libm gives the same bit pattern
on both targets anyway. <strong>The other four are not</strong>: IEEE-754 requires
nothing of <code>sinf</code>, <code>cosf</code>, <code>atan2f</code> or
<code>powf</code>, 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 <code>-lm</code>.</p>
<p>The maths in the file is not Flan, and it is libm's:
<code>(declare sqrt-f32 [x f32] f32 "sqrtf")</code> 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 <code>f32</code> and <code>u32</code> that the language does not have.
A <code>declare</code> 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.</p>
<p><strong>One split is worth knowing before calling any of them.</strong> IEEE-754
specifies <code>sqrt</code>, <code>fabs</code>, <code>floor</code>, <code>ceil</code>,
<code>round</code> and <code>fmod</code> as exact or correctly rounded, so those give the
same bit pattern under glibc, musl and wasi-libc. <strong>It requires nothing of the
rest</strong> — <code>sin</code>, <code>cos</code>, <code>tan</code>, the inverses, the
logarithms, <code>exp</code>, <code>pow</code>, <code>hypot</code>, <code>cbrt</code> —
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
<code>-lm</code>.</p>
<p>The primitives underneath are few — a primitive is the only thing implemented
twice per backend: <code>argv</code>,