Floats as written, never fused; fma and muladd; f32/f64 bits, f16 and bf16

- a * b + c is two roundings on every processor and --cpu: the IR has no fast-math flags, and clang
  compiles it with -ffp-contract=on (checked in the assembly: never fused, also with FMA hardware)
- fma(a, b, c) / v.fma(b, c): one rounding, always (llvm.fma)
- muladd(a, b, c) / v.muladd(b, c): the fastest way here (llvm.fmuladd), what v.fma was before
- f32_bits / f32_from_bits / f64_bits / f64_from_bits: bitcasts, also of vectors
- f16_to_f32 / f32_to_f16 (LLVM half), bf16_to_f32 / f32_to_bf16 (integer ops, ties to even, NaN kept)
- docs (en, pl): what is guaranteed about floats, the new built-ins

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Psc42moPLwP43Cmb9oZqSR

docs/en/cheatsheet.md

@@ -259,7 +259,8 @@ Types: `f32x4/8/16`, `f64x2/4/8`, `i8x16/32/64`, `i16x8/16/32`, `i32x4/8/16`, `i
 | `v[i]`, `v[i] = x`, `v.len` | a lane; the number of lanes |
 | `min(v, w)`, `max`, `abs(v)`, `clamp(v, lo, hi)` | lane by lane |
 | `v.sum()`, `v.min()`, `v.max()`, `v.dot(w)` | one number from all the lanes |
-| `v.sqrt()`, `v.floor()`, `v.ceil()`, `v.round()`, `v.fma(b, c)` | every lane (float lanes; `fma` = `v * b + c`) |
+| `v.sqrt()`, `v.floor()`, `v.ceil()`, `v.round()` | every lane (float lanes) |
+| `v.fma(b, c)`, `v.muladd(b, c)` | `v * b + c`: one rounding always / the fastest way here |
 | `v.store(list, i)` | lanes into `list[i]` … |
 
 ---
@@ -324,6 +325,9 @@ Automatic only when nothing is lost: smaller → bigger number, any int → floa
 | `shell(command)` | run a command, its exit code |
 | `exit(code)` | stop the program now (also ends an `or { }` block) |
 | `min(a, b, …)`, `max(a, b, …)`, `abs(x)`, `clamp(x, lo, hi)` | for any numbers (your own function of that name wins) |
+| `fma(a, b, c)`, `muladd(a, b, c)` | `a * b + c` with one rounding / fused where the processor can (`a * b + c` written out is never fused) |
+| `f32_bits(x)`, `f32_from_bits(u)`, `f64_bits`, `f64_from_bits` | a float's bits as `u32`/`u64` and back, unchanged |
+| `f16_to_f32(h)`, `f32_to_f16(x)`, `bf16_to_f32(b)`, `f32_to_bf16(x)` | half and bfloat16 numbers (in a `u16`) ↔ `f32`; also on vectors |
 | `error(message)` | `return error("…")` in a `-> T!` function |
 | `addr(x)` | pointer to x; `addr(fn)` = C callback |
 | `assert(cond)`, `assert(cond, "why")` | stop with a runtime error when cond is false (in a `test` block: that test fails) |

docs/en/for-engineers.md

@@ -308,7 +308,8 @@ type goes; a call `f32x8(…)` goes to `vecMake()`.
 | `min`, `max`, `abs` | `llvm.minnum`/`maxnum`, `llvm.smin`/`umin`…, `llvm.fabs`/`llvm.abs` |
 | `v.sum()`, `.min()`, `.max()` | `llvm.vector.reduce.fadd` (with `reassoc`), `.add`, `.fmin`, `.smin`… |
 | `v.sqrt()`, `.floor()`, `.ceil()`, `.round()` | the `llvm.*` intrinsic of that name |
-| `v.fma(b, c)` | `llvm.fmuladd` (fused where the target has FMA; `llvm.fma` would be a libm call per lane without it) |
+| `v.fma(b, c)`, `fma(a, b, c)` | `llvm.fma` (one rounding; a libm call per lane without FMA hardware) |
+| `v.muladd(b, c)`, `muladd(a, b, c)` | `llvm.fmuladd` (fused where the target has FMA, `fmul` + `fadd` where not) |
 | `v.store(list, i)` | bounds check, `store <8 x float> …, align 4` |
 
 Integer `/` and `%` are refused: neither x86 nor ARM has a lane-wise integer divide, and LLVM would scalarise it.
@@ -507,6 +508,19 @@ test.
 | Missing return | UB if used | compile error |
 | Order of evaluation | mostly unspecified | left to right; in an assignment the right side first |
 | Shift ≥ width, out-of-range float→int | UB | poison, as in LLVM (not checked) |
+| `a * b + c` on floats | may be fused into an FMA (`-ffp-contract`) | never fused: two roundings |
+
+### Floating point
+
+Float operations are emitted as plain `fadd`/`fsub`/`fmul`/`fdiv`/`fneg` with **no fast-math flags**, so LLVM may not
+reassociate, contract or assume no NaNs; `linker.fch` also compiles the IR with `-ffp-contract=on`, which fuses only the
+`llvm.fmuladd` calls Finch asks for. The only flagged instruction is `v.sum()`'s `llvm.vector.reduce.fadd` (`reassoc`).
+`fma` → `llvm.fma`, `muladd` → `llvm.fmuladd`, on scalars and vectors (`fmaCall`, `fusedOp`).
+
+`f32_bits`/`f32_from_bits`/`f64_…` are `bitcast` (also between vectors). `f16_to_f32` is `bitcast i16 → half` + `fpext`,
+`f32_to_f16` is `fptrunc` + `bitcast`: one instruction with F16C or on ARM64, `__extendhfsf2`/`__truncsfhf2` from the C
+compiler's runtime library elsewhere. `bf16_to_f32` is `zext` + `shl 16` + `bitcast`; `f32_to_bf16` rounds to nearest even
+in integer ops (`u + 0x7fff + ((u >> 16) & 1)`, then `>> 16`) and keeps NaNs quiet (`| 0x40`), so it needs no library.
 | Dangling pointer from `addr()` / `new` + `free` | UB | still UB |
 
 ### The chain of calls on runtime errors

docs/en/for-technicians.md

@@ -345,6 +345,40 @@ str(p)           // a C char pointer as text
 **Fixed-size numbers wrap around**, like real hardware: a `u8` holding 255 plus 1 becomes 0.
 Calculations happen in the type of the values involved, so `u8(250) + 10` is `4`, not `260`.
 
+### Decimal numbers: what is guaranteed
+
+`f32` and `f64` are IEEE 754 numbers, and Finch computes them **exactly as written**, the same on every processor and
+with every `--cpu`:
+
+- Each operation rounds once, in the order written. `a * b - c * d` is two multiplies and a subtract, never folded into a
+  fused multiply-add (C compilers often do that by default, which changes the last digits).
+- Nothing is reordered or simplified: `(a + b) + c` is not turned into `a + (b + c)`, and `x * 0.0` is not assumed to be 0
+  (it is NaN when `x` is). The one exception is `v.sum()` on a vector, which adds the lanes in the fastest order.
+
+When speed matters more than the last bit, ask for it by name:
+
+| | |
+|---|---|
+| `fma(a, b, c)` | `a * b + c` with a single rounding, always. One instruction where the processor has FMA (with `--cpu`, on ARM64); a slower library call elsewhere. |
+| `muladd(a, b, c)` | `a * b + c` the fastest way on this processor: fused where it can be, a multiply and an add where not. |
+
+Both work on numbers and on vectors (`v.fma(b, c)`, `v.muladd(b, c)`).
+
+**Numbers as bits.** File formats and model weights store numbers as raw bits. These look at the bits without
+converting the value (all of them also work lane by lane on vectors):
+
+| | |
+|---|---|
+| `f32_bits(x)` / `f32_from_bits(u)` | `f32` ↔ the `u32` with the same bits (`f32_bits(1.0)` is `0x3f800000`) |
+| `f64_bits(x)` / `f64_from_bits(u)` | `f64` ↔ `u64` |
+| `f16_to_f32(h)` / `f32_to_f16(x)` | half precision, kept in a `u16`, ↔ `f32` (rounded to the nearest) |
+| `bf16_to_f32(b)` / `f32_to_bf16(x)` | bfloat16, kept in a `u16`, ↔ `f32` (rounded to the nearest, ties to even) |
+
+```c
+[]u16 weights = read_weights()               // f16 numbers from a file
+v := f16_to_f32(u16x8(weights, 0))           // 8 of them as f32, at once
+```
+
 ---
 
 ## 7. Operators
@@ -697,6 +731,8 @@ Using `.value`, `.field` or `[i]` on a `null` pointer **stops the program** with
 | `exit(code)` | Ends the program right away. |
 | `addr(x)` | A pointer to `x`; `addr(function)` gives C a function pointer (callback). |
 | `new(value)`, `free(p)` | Manual heap memory (see [Memory](#13-memory-who-frees-what)). |
+| `fma(a, b, c)`, `muladd(a, b, c)` | `a * b + c` with one rounding / the fastest way (see [Decimal numbers](#decimal-numbers-what-is-guaranteed)). |
+| `f32_bits(x)`, `f32_from_bits(u)`, `f16_to_f32(h)`, … | A number's bits, and half-precision numbers (see [Decimal numbers](#decimal-numbers-what-is-guaranteed)). |
 | `int(x)`, `str(x)`, `u8(x)`, … | Conversions (see [Conversions](#conversions)). `str(x)` works on every type and gives the same text `print` shows. `int(text) or 0` handles text that isn't a number. |
 
 ---
@@ -1370,7 +1406,8 @@ arrays, maps and function parameters, and `print` and `str` show it.
 | `v.sum()`, `v.min()`, `v.max()` | all the lanes made into one number |
 | `v.dot(w)` | `(v * w).sum()` |
 | `v.sqrt()`, `v.floor()`, `v.ceil()`, `v.round()` | every lane (float lanes) |
-| `v.fma(b, c)` | `v * b + c`, as one instruction where the processor has one |
+| `v.fma(b, c)` | `v * b + c` with one rounding, always |
+| `v.muladd(b, c)` | `v * b + c` the fastest way here: one instruction where the processor can |
 | `v.store(list, i)` | writes the lanes to `list[i]` … (or through a `ptr`) |
 
 `<` and the other comparisons would give an answer per lane, so they aren't there; use `min`, `max` or `v[i] < w[i]`.
@@ -1382,7 +1419,7 @@ fn dot([]f32 a, []f32 b) -> f32 {
     acc := f32x8(0.0)
     i := 0
     while i + 8 <= a.len {
-        acc = f32x8(a, i).fma(f32x8(b, i), acc)
+        acc = f32x8(a, i).muladd(f32x8(b, i), acc)
         i += 8
     }
     total := acc.sum()

docs/pl/dla-inzynierow.md

@@ -314,7 +314,8 @@ wszędzie tam, gdzie może stać typ; wywołanie `f32x8(…)` trafia do `vecMake
 | `min`, `max`, `abs` | `llvm.minnum`/`maxnum`, `llvm.smin`/`umin`…, `llvm.fabs`/`llvm.abs` |
 | `v.sum()`, `.min()`, `.max()` | `llvm.vector.reduce.fadd` (z `reassoc`), `.add`, `.fmin`, `.smin`… |
 | `v.sqrt()`, `.floor()`, `.ceil()`, `.round()` | funkcja wewnętrzna `llvm.*` o tej nazwie |
-| `v.fma(b, c)` | `llvm.fmuladd` (jedna instrukcja tam, gdzie procesor ma FMA; `llvm.fma` bez FMA wołałoby funkcję z libm dla każdego elementu) |
+| `v.fma(b, c)`, `fma(a, b, c)` | `llvm.fma` (jedno zaokrąglenie; bez sprzętowego FMA wywołanie libm dla każdego elementu) |
+| `v.muladd(b, c)`, `muladd(a, b, c)` | `llvm.fmuladd` (jedna instrukcja tam, gdzie procesor ma FMA, a gdzie nie, `fmul` + `fadd`) |
 | `v.store(lista, i)` | sprawdzenie zakresu, `store <8 x float> …, align 4` |
 
 Dzielenie całkowite `/` i `%` jest odrzucane: ani x86, ani ARM nie mają instrukcji dzielącej liczby całkowite element
@@ -515,6 +516,21 @@ Testem akceptacyjnym jest `tests/lsp_test.py`.
 | Brak `return` | UB, jeśli wynik jest użyty | błąd kompilacji |
 | Kolejność obliczania | w dużej mierze nieokreślona | od lewej do prawej; w przypisaniu najpierw prawa strona |
 | Przesunięcie o ≥ szerokość, liczba zmiennoprzecinkowa poza zakresem przy zamianie na całkowitą | UB | poison, jak w LLVM (niesprawdzane) |
+| `a * b + c` na liczbach zmiennoprzecinkowych | może zostać złożone w FMA (`-ffp-contract`) | nigdy nie jest łączone: dwa zaokrąglenia |
+
+### Liczby zmiennoprzecinkowe
+
+Działania na liczbach zmiennoprzecinkowych to zwykłe `fadd`/`fsub`/`fmul`/`fdiv`/`fneg` **bez flag fast-math**, więc LLVM
+nie może zmieniać kolejności działań, łączyć ich ani zakładać braku NaN; `linker.fch` kompiluje też IR z
+`-ffp-contract=on`, które łączy tylko wywołania `llvm.fmuladd`, o które prosi Finch. Jedyna instrukcja z flagą to
+`llvm.vector.reduce.fadd` w `v.sum()` (`reassoc`). `fma` → `llvm.fma`, `muladd` → `llvm.fmuladd`, na liczbach i wektorach
+(`fmaCall`, `fusedOp`).
+
+`f32_bits`/`f32_from_bits`/`f64_…` to `bitcast` (także między wektorami). `f16_to_f32` to `bitcast i16 → half` + `fpext`,
+`f32_to_f16` to `fptrunc` + `bitcast`: jedna instrukcja z F16C albo na ARM64, a gdzie indziej `__extendhfsf2`/`__truncsfhf2`
+z biblioteki uruchomieniowej kompilatora C. `bf16_to_f32` to `zext` + `shl 16` + `bitcast`; `f32_to_bf16` zaokrągla do
+najbliższej parzystej na liczbach całkowitych (`u + 0x7fff + ((u >> 16) & 1)`, potem `>> 16`) i zostawia NaN cichym
+(`| 0x40`), więc nie potrzebuje żadnej biblioteki.
 | Wiszący wskaźnik z `addr()` / `new` + `free` | UB | nadal UB |
 
 ### Ścieżka wywołań przy błędach w trakcie działania

docs/pl/dla-technikow.md

@@ -348,6 +348,41 @@ str(p)           // wskaźnik na znaki z C jako tekst
 **Liczby o stałym rozmiarze się zawijają**, jak w prawdziwym sprzęcie: `u8` równe 255 plus 1 daje 0.
 Obliczenia odbywają się w typie użytych wartości, więc `u8(250) + 10` to `4`, a nie `260`.
 
+### Liczby zmiennoprzecinkowe: co jest gwarantowane
+
+`f32` i `f64` to liczby w standardzie IEEE 754, a Finch liczy je **dokładnie tak, jak są zapisane**, tak samo na każdym
+procesorze i z każdym `--cpu`:
+
+- Każde działanie zaokrągla raz, w zapisanej kolejności. `a * b - c * d` to dwa mnożenia i odejmowanie, nigdy nie
+  złożone w jedno „mnożenie z dodawaniem” (FMA); kompilatory C często robią to domyślnie, co zmienia ostatnie cyfry.
+- Nic nie jest przestawiane ani upraszczane: `(a + b) + c` nie zamienia się w `a + (b + c)`, a `x * 0.0` nie jest
+  z góry uznawane za 0 (gdy `x` to NaN, wynik też jest NaN). Jedyny wyjątek to `v.sum()` na wektorze, które dodaje
+  elementy w najszybszej kolejności.
+
+Gdy szybkość jest ważniejsza niż ostatni bit, poproś o nią wprost:
+
+| | |
+|---|---|
+| `fma(a, b, c)` | `a * b + c` z jednym zaokrągleniem, zawsze. Jedna instrukcja tam, gdzie procesor ma FMA (z `--cpu`, na ARM64); gdzie indziej wolniejsza funkcja biblioteczna. |
+| `muladd(a, b, c)` | `a * b + c` najszybciej, jak się da na tym procesorze: jedną instrukcją, gdzie można, a gdzie nie, mnożeniem i dodawaniem. |
+
+Obie działają na liczbach i na wektorach (`v.fma(b, c)`, `v.muladd(b, c)`).
+
+**Liczby jako bity.** Formaty plików i wagi modeli zapisują liczby jako surowe bity. Te funkcje patrzą na bity bez
+przeliczania wartości (wszystkie działają też na wektorach, element po elemencie):
+
+| | |
+|---|---|
+| `f32_bits(x)` / `f32_from_bits(u)` | `f32` ↔ `u32` z tymi samymi bitami (`f32_bits(1.0)` to `0x3f800000`) |
+| `f64_bits(x)` / `f64_from_bits(u)` | `f64` ↔ `u64` |
+| `f16_to_f32(h)` / `f32_to_f16(x)` | liczba połówkowej precyzji, trzymana w `u16`, ↔ `f32` (zaokrąglenie do najbliższej) |
+| `bf16_to_f32(b)` / `f32_to_bf16(x)` | bfloat16, trzymany w `u16`, ↔ `f32` (zaokrąglenie do najbliższej, remisy do parzystej) |
+
+```c
+[]u16 wagi = wczytaj_wagi()                  // liczby f16 z pliku
+v := f16_to_f32(u16x8(wagi, 0))              // osiem z nich jako f32, naraz
+```
+
 ---
 
 ## 7. Operatory
@@ -700,6 +735,8 @@ Użycie `.value`, `.pole` albo `[i]` na wskaźniku `null` **zatrzymuje program**
 | `exit(kod)` | Natychmiast kończy program. |
 | `addr(x)` | Wskaźnik na `x`; `addr(funkcja)` daje C wskaźnik na funkcję (funkcję zwrotną). |
 | `new(wartosc)`, `free(p)` | Ręczna pamięć na stercie (zob. [Pamięć](#13-pamięć-kto-co-zwalnia)). |
+| `fma(a, b, c)`, `muladd(a, b, c)` | `a * b + c` z jednym zaokrągleniem / najszybciej (zob. [Liczby zmiennoprzecinkowe](#liczby-zmiennoprzecinkowe-co-jest-gwarantowane)). |
+| `f32_bits(x)`, `f32_from_bits(u)`, `f16_to_f32(h)`, … | Bity liczby i liczby połówkowej precyzji (zob. [Liczby zmiennoprzecinkowe](#liczby-zmiennoprzecinkowe-co-jest-gwarantowane)). |
 | `int(x)`, `str(x)`, `u8(x)`, … | Zamiany typów (zob. [Zamiana typów](#zamiana-typów)). `str(x)` działa na każdym typie i daje ten sam tekst, który pokazuje `print`. `int(tekst) or 0` obsługuje tekst, który nie jest liczbą. |
 
 ---
@@ -1377,7 +1414,8 @@ struktury, elementem tablicy, wartością w mapie i parametrem funkcji, a `print
 | `v.sum()`, `v.min()`, `v.max()` | wszystkie elementy złożone w jedną liczbę |
 | `v.dot(w)` | `(v * w).sum()` |
 | `v.sqrt()`, `v.floor()`, `v.ceil()`, `v.round()` | każdy element (liczby zmiennoprzecinkowe) |
-| `v.fma(b, c)` | `v * b + c`, jedną instrukcją tam, gdzie procesor ją ma |
+| `v.fma(b, c)` | `v * b + c` z jednym zaokrągleniem, zawsze |
+| `v.muladd(b, c)` | `v * b + c` najszybciej, jak się da: jedną instrukcją, gdzie procesor to potrafi |
 | `v.store(lista, i)` | zapisuje elementy do `lista[i]` … (albo przez `ptr`) |
 
 `<` i inne porównania dałyby osobną odpowiedź dla każdego elementu, więc ich nie ma; użyj `min`, `max` albo `v[i] < w[i]`.
@@ -1389,7 +1427,7 @@ fn iloczyn([]f32 a, []f32 b) -> f32 {
     suma := f32x8(0.0)
     i := 0
     while i + 8 <= a.len {
-        suma = f32x8(a, i).fma(f32x8(b, i), suma)
+        suma = f32x8(a, i).muladd(f32x8(b, i), suma)
         i += 8
     }
     wynik := suma.sum()

docs/pl/sciaga.md

@@ -259,7 +259,8 @@ Typy: `f32x4/8/16`, `f64x2/4/8`, `i8x16/32/64`, `i16x8/16/32`, `i32x4/8/16`, `i6
 | `v[i]`, `v[i] = x`, `v.len` | jeden element; liczba elementów |
 | `min(v, w)`, `max`, `abs(v)`, `clamp(v, lo, hi)` | element po elemencie |
 | `v.sum()`, `v.min()`, `v.max()`, `v.dot(w)` | jedna liczba ze wszystkich elementów |
-| `v.sqrt()`, `v.floor()`, `v.ceil()`, `v.round()`, `v.fma(b, c)` | każdy element (liczby zmiennoprzecinkowe; `fma` = `v * b + c`) |
+| `v.sqrt()`, `v.floor()`, `v.ceil()`, `v.round()` | każdy element (liczby zmiennoprzecinkowe) |
+| `v.fma(b, c)`, `v.muladd(b, c)` | `v * b + c`: zawsze z jednym zaokrągleniem / najszybciej na tym procesorze |
 | `v.store(lista, i)` | zapisuje elementy do `lista[i]` … |
 
 ---
@@ -324,6 +325,9 @@ Automatycznie tylko wtedy, gdy nic nie ginie: mniejsza liczba → większa, licz
 | `shell(polecenie)` | uruchamia polecenie, zwraca jego kod wyjścia |
 | `exit(kod)` | od razu kończy program (kończy też blok `or { }`) |
 | `min(a, b, …)`, `max(a, b, …)`, `abs(x)`, `clamp(x, lo, hi)` | dla dowolnych liczb (twoja funkcja o tej samej nazwie ma pierwszeństwo) |
+| `fma(a, b, c)`, `muladd(a, b, c)` | `a * b + c` z jednym zaokrągleniem / jedną instrukcją, gdzie procesor potrafi (zapisane `a * b + c` nigdy nie jest łączone) |
+| `f32_bits(x)`, `f32_from_bits(u)`, `f64_bits`, `f64_from_bits` | bity liczby zmiennoprzecinkowej jako `u32`/`u64` i z powrotem, bez zmian |
+| `f16_to_f32(h)`, `f32_to_f16(x)`, `bf16_to_f32(b)`, `f32_to_bf16(x)` | liczby połówkowej precyzji i bfloat16 (w `u16`) ↔ `f32`; także na wektorach |
 | `error(komunikat)` | `return error("…")` w funkcji `-> T!` |
 | `addr(x)` | wskaźnik na x; `addr(fn)` = funkcja zwrotna dla C |
 | `assert(warunek)`, `assert(warunek, "dlaczego")` | zatrzymuje program, gdy warunek jest fałszywy (w bloku `test`: ten test nie przechodzi) |

examples/simd.fch

@@ -13,12 +13,12 @@ fn dot_plain([]f32 a, []f32 b) -> f32 {
     return total
 }
 
-// eight at a time: f32x8(list, i) takes list[i] ... list[i+7], fma multiplies and adds in one step
+// eight at a time: f32x8(list, i) takes list[i] ... list[i+7]; muladd multiplies and adds, in one step where it can
 fn dot_simd([]f32 a, []f32 b) -> f32 {
     acc := f32x8(0.0)
     i := 0
     while i + 8 <= a.len {
-        acc = f32x8(a, i).fma(f32x8(b, i), acc)
+        acc = f32x8(a, i).muladd(f32x8(b, i), acc)
         i += 8
     }
     total := acc.sum()

src/docs.fch

@@ -26,6 +26,16 @@ fn builtins() -> []Doc {
         doc("max", "max(a, b, ...) -> number", "The biggest of the numbers, in the type they all fit."),
         doc("abs", "abs(x) -> number", "x without its minus sign (any number type)."),
         doc("clamp", "clamp(x, lo, hi) -> number", "x, but at least lo and at most hi."),
+        doc("fma", "fma(a, b, c) -> float", "a * b + c with one rounding, always (numbers or vectors). a * b + c written out is never fused: two roundings."),
+        doc("muladd", "muladd(a, b, c) -> float", "a * b + c the fastest way here: one fused instruction where the processor has one, else a multiply and an add."),
+        doc("f32_bits", "f32_bits(f32 x) -> u32", "The bits of x, unchanged (also f32xN -> u32xN). f64_bits gives a u64."),
+        doc("f32_from_bits", "f32_from_bits(u32 bits) -> f32", "The f32 with these bits (also u32xN -> f32xN). f64_from_bits takes a u64."),
+        doc("f64_bits", "f64_bits(float x) -> u64", "The bits of x, unchanged (also f64xN -> u64xN)."),
+        doc("f64_from_bits", "f64_from_bits(u64 bits) -> float", "The float with these bits (also u64xN -> f64xN)."),
+        doc("f16_to_f32", "f16_to_f32(u16 half) -> f32", "A half-precision number, kept in a u16, as an f32 (also u16xN -> f32xN)."),
+        doc("f32_to_f16", "f32_to_f16(f32 x) -> u16", "x as a half-precision number in a u16, rounded to the nearest (also f32xN -> u16xN)."),
+        doc("bf16_to_f32", "bf16_to_f32(u16 b) -> f32", "A bfloat16, kept in a u16, as an f32 (also u16xN -> f32xN)."),
+        doc("f32_to_bf16", "f32_to_bf16(f32 x) -> u16", "x as a bfloat16 in a u16, rounded to the nearest, ties to even (also f32xN -> u16xN)."),
         doc("atomic_add", "atomic_add(ptr[int] p, n) -> int", "Adds n to the number at p in one step other threads can't split; gives back what was there before."),
         doc("atomic_sub", "atomic_sub(ptr[int] p, n) -> int", "Takes n from the number at p in one step; gives back what was there before."),
         doc("atomic_load", "atomic_load(ptr[int] p) -> int", "Reads the number at p, as another thread last wrote it."),
@@ -107,7 +117,8 @@ fn vecMethods() -> []Doc {
         doc("ceil", "ceil() -> V", "Every lane rounded up (f32/f64 lanes)."),
         doc("round", "round() -> V", "Every lane rounded to the nearest whole number, halves away from 0 (f32/f64 lanes)."),
         doc("abs", "abs() -> V", "Every lane without its minus sign."),
-        doc("fma", "fma(V b, V c) -> V", "v * b + c, as one instruction where the processor has one (f32/f64 lanes)."),
+        doc("fma", "fma(V b, V c) -> V", "v * b + c with one rounding, always (f32/f64 lanes). One instruction with FMA hardware, slower without."),
+        doc("muladd", "muladd(V b, V c) -> V", "v * b + c the fastest way here: fused into one instruction where the processor can (f32/f64 lanes)."),
         doc("store", "store([]T list, int i)", "Writes the lanes into list[i], list[i+1], ... (or through a ptr[T]).")
     ]
 }

src/gen.fch

@@ -5572,6 +5572,12 @@ fn call(ptr[Gen] g, ast.Node c) -> Val {
     if ei >= 0 {
         return enumFrom(g, c, ei)
     }
+    if ["fma", "muladd"].contains(n) {
+        return fmaCall(g, c)        // Finch's, also for f32 and vectors (C's fma is the same for one double)
+    }
+    if ["f32_bits", "f32_from_bits", "f64_bits", "f64_from_bits", "f16_to_f32", "f32_to_f16", "bf16_to_f32", "f32_to_bf16"].contains(n) && !g.ci.fns.has(n) {
+        return bitsCall(g, c)
+    }
     if g.ci.fns.has(n) && !["min", "max", "abs", "clamp"].contains(n) {
         return callC(g, c, g.ci.fns[n])        // (C's abs only takes an int; Finch's takes any number)
     }
@@ -6647,8 +6653,14 @@ fn laneBytes(ast.Type vt) -> str {
     return str(bitsOf(vt.elem[0]) / 8)
 }
 
-// v8f32, v32i8: how LLVM's vector intrinsics are named for this type
+// v8f32, v32i8 (or f32, i64 for a plain number): how LLVM's intrinsics are named for this type
 fn vecSuffix(ast.Type vt) -> str {
+    if vt.kind != "vec" {
+        if isFloat(vt) {
+            return "f" + str(bitsOf(vt))
+        }
+        return "i" + str(bitsOf(vt))
+    }
     if isFloat(vt.elem[0]) {
         return "v" + str(vt.count) + "f" + str(bitsOf(vt.elem[0]))
     }
@@ -6984,17 +6996,15 @@ fn vecMethod(ptr[Gen] g, ast.Node m, Val v) -> Val {
         need(g, m, 0)
         return vecAbs(g, v)
     }
-    if n == "sqrt" || n == "floor" || n == "ceil" || n == "round" || n == "fma" {
+    if n == "sqrt" || n == "floor" || n == "ceil" || n == "round" || n == "fma" || n == "muladd" {
         if !isFloat(lane) {
             fail(g, m, n + "() is for f32 and f64 lanes, but this is " + ast.show(vt))
         }
-        if n == "fma" {
+        if n == "fma" || n == "muladd" {
             need(g, m, 2)
-            b := asVec(g, expr(g, m.kids[1]), vt, "fma", m.kids[1])
+            b := asVec(g, expr(g, m.kids[1]), vt, n, m.kids[1])
-            c := asVec(g, expr(g, m.kids[2]), vt, "fma", m.kids[2])
+            c := asVec(g, expr(g, m.kids[2]), vt, n, m.kids[2])
-            // fmuladd: one fused instruction where the processor has one (FMA, NEON), a multiply and an add where
+            return Val(v: fusedOp(g, n, vt, v.v, b.v, c.v), ty: vt)
-            // it doesn't (llvm.fma would call a slow library function there, to round only once)
-            return Val(v: vecIntrinsic(g, lt, "fmuladd", vt, [lt, lt, lt], [v.v, b.v, c.v], ""), ty: vt)
         }
         need(g, m, 0)
         return Val(v: vecIntrinsic(g, lt, n, vt, [lt], [v.v], ""), ty: vt)
@@ -7028,6 +7038,178 @@ fn vecMethod(ptr[Gen] g, ast.Node m, Val v) -> Val {
         emit(g, "store " + lt + " " + v.v + ", ptr " + ep + ", align " + laneBytes(vt))
         return Val(ty: T("void"))
     }
-    fail(g, m, "a " + ast.show(vt) + " has no method '" + n + "' (it has: sum, min, max, dot, abs, sqrt, floor, ceil, round, fma, store)")
+    fail(g, m, "a " + ast.show(vt) + " has no method '" + n + "' (it has: sum, min, max, dot, abs, sqrt, floor, ceil, round, fma, muladd, store)")
     return v
 }
+
+// ---------- fused multiply-add, and numbers as bits ----------
+//
+// Finch never fuses on its own: a * b + c is a multiply and an add, two roundings, on every processor and with every
+// --cpu (the IR has no fast-math flags, and clang compiles it with -ffp-contract=on, which fuses only what asks to be).
+// Fusing is asked for by name:
+//   fma(a, b, c)      one rounding, always (one instruction with FMA hardware, a library call without it)
+//   muladd(a, b, c)   whichever is faster here: fused where the processor can, a multiply and an add where it can't
+
+// the llvm.fma / llvm.fmuladd call for a number or a vector type
+fn fusedOp(ptr[Gen] g, str name, ast.Type t, str a, str b, str c) -> str {
+    lt := llty(g, t)
+    which := "fma"
+    if name == "muladd" {
+        which = "fmuladd"
+    }
+    return vecIntrinsic(g, lt, which, t, [lt, lt, lt], [a, b, c], "")
+}
+
+fn fmaCall(ptr[Gen] g, ast.Node c) -> Val {
+    n := c.name
+    if g.index {
+        noteRef(g, c, c.col, n.len, docs.hover(docs.find(docs.builtins(), n)), "", 0, 0)
+    }
+    checkArgs(g, c, 0, 3, n)
+    []Val vals
+    for a in c.kids {
+        v := expr(g, a)
+        if !isNumber(v.ty) && v.ty.kind != "vec" {
+            fail(g, a, n + "(...) is for numbers, but this is " + ast.show(v.ty))
+        }
+        vals.push(v)
+    }
+    // the type all three take: a vector if one is, else the float type they fit
+    t := T("void")
+    for v in vals {
+        if v.ty.kind == "vec" {
+            t = v.ty
+        }
+    }
+    if t.kind == "vec" {
+        if !isFloat(t.elem[0]) {
+            fail(g, c, n + "(...) is for f32 and f64 lanes, but this is " + ast.show(t))
+        }
+        []str vs
+        for i in 0..3 {
+            vs.push(asVec(g, vals[i], t, n, c.kids[i]).v)
+        }
+        return Val(v: fusedOp(g, n, t, vs[0], vs[1], vs[2]), ty: t)
+    }
+    // f32 when an f32 takes part and nothing wider does (numbers written in the code adapt), else float
+    narrow := false
+    wide := false
+    for v in vals {
+        if v.ty.kind == "f32" {
+            narrow = true
+        } else if !v.lit {
+            wide = true
+        }
+    }
+    t = T("float")
+    if narrow && !wide {
+        t = T("f32")
+    }
+    []str vs
+    for i in 0..3 {
+        vs.push(coerce(g, vals[i], t, c.kids[i], "each number of " + n + "(...)").v)
+    }
+    return Val(v: fusedOp(g, n, t, vs[0], vs[1], vs[2]), ty: t)
+}
+
+// the same shape as t (a number, or a vector of as many lanes) with lanes of another kind
+fn shaped(ast.Type t, str lane) -> ast.Type {
+    if t.kind != "vec" {
+        return T(lane)
+    }
+    out := ast.Type(kind: "vec", count: t.count)
+    out.elem.push(T(lane))
+    return out
+}
+
+// a constant in every lane of t (or just the number)
+fn konst(ptr[Gen] g, ast.Type t, str n) -> str {
+    if t.kind != "vec" {
+        return n
+    }
+    return splat(g, n, t)
+}
+
+// f32_bits(x) -> u32 and back, f64 the same; f16 and bf16 (as the u16 that holds them) to f32 and back.
+// Each works on a number or on a vector, lane by lane: f16_to_f32(u16x8) is an f32x8.
+fn bitsCall(ptr[Gen] g, ast.Node c) -> Val {
+    n := c.name
+    if g.index {
+        noteRef(g, c, c.col, n.len, docs.hover(docs.find(docs.builtins(), n)), "", 0, 0)
+    }
+    checkArgs(g, c, 0, 1, n)
+    v := expr(g, c.kids[0])
+    from := ""          // the lane type it takes
+    to := ""
+    if n == "f32_bits" || n == "f32_to_f16" || n == "f32_to_bf16" {
+        from = "f32"
+        to = "u32"
+        if n != "f32_bits" {
+            to = "u16"
+        }
+    } else if n == "f64_bits" {
+        from = "float"
+        to = "u64"
+    } else if n == "f32_from_bits" {
+        from = "u32"
+        to = "f32"
+    } else if n == "f64_from_bits" {
+        from = "u64"
+        to = "float"
+    } else {
+        from = "u16"
+        to = "f32"
+    }
+    lane := v.ty
+    if v.ty.kind == "vec" {
+        lane = v.ty.elem[0]
+    }
+    // the bits may come as a signed number of the same size (i32 for u32), and a number written in the code fits
+    if isInt(T(from)) && isInt(lane) && bitsOf(lane) == bitsOf(T(from)) {
+        v = Val(v: v.v, ty: shaped(v.ty, from))
+    } else if v.ty.kind != "vec" && v.lit {
+        v = coerce(g, v, T(from), c.kids[0], "the value of " + n + "(...)")
+    } else if lane.kind != from {
+        fail(g, c.kids[0], n + "(...) takes " + ast.show(T(from)) + " (or a vector of them), but this is " + ast.show(v.ty))
+    }
+    src := shaped(v.ty, from)
+    out := shaped(v.ty, to)
+    il := llty(g, src)
+    ol := llty(g, out)
+    if n == "f32_bits" || n == "f64_bits" || n == "f32_from_bits" || n == "f64_from_bits" {
+        return Val(v: op(g, "bitcast " + il + " " + v.v + " to " + ol), ty: out)
+    }
+    if n == "bf16_to_f32" {
+        // bf16 is the top half of an f32
+        wide := shaped(v.ty, "u32")
+        wl := llty(g, wide)
+        z := op(g, "zext " + il + " " + v.v + " to " + wl)
+        sh := op(g, "shl " + wl + " " + z + ", " + konst(g, wide, "16"))
+        return Val(v: op(g, "bitcast " + wl + " " + sh + " to " + ol), ty: out)
+    }
+    if n == "f32_to_bf16" {
+        // the top half, rounded to the nearest (ties to even); a NaN stays a NaN
+        wide := shaped(v.ty, "u32")
+        wl := llty(g, wide)
+        u := op(g, "bitcast " + il + " " + v.v + " to " + wl)
+        hi := op(g, "lshr " + wl + " " + u + ", " + konst(g, wide, "16"))
+        odd := op(g, "and " + wl + " " + hi + ", " + konst(g, wide, "1"))
+        bias := op(g, "add " + wl + " " + odd + ", " + konst(g, wide, "32767"))
+        rounded := op(g, "lshr " + wl + " " + op(g, "add " + wl + " " + u + ", " + bias) + ", " + konst(g, wide, "16"))
+        quiet := op(g, "or " + wl + " " + hi + ", " + konst(g, wide, "64"))
+        isNan := op(g, "fcmp uno " + il + " " + v.v + ", " + v.v)
+        pickd := op(g, "select " + llty(g, shaped(v.ty, "bool")) + " " + isNan + ", " + wl + " " + quiet + ", " + wl + " " + rounded)
+        return Val(v: op(g, "trunc " + wl + " " + pickd + " to " + ol), ty: out)
+    }
+    // f16: LLVM's half type (one instruction with F16C or on ARM64, a small library function elsewhere)
+    hl := "half"
+    if v.ty.kind == "vec" {
+        hl = "<" + str(v.ty.count) + " x half>"
+    }
+    if n == "f16_to_f32" {
+        h := op(g, "bitcast " + il + " " + v.v + " to " + hl)
+        return Val(v: op(g, "fpext " + hl + " " + h + " to " + ol), ty: out)
+    }
+    h := op(g, "fptrunc " + il + " " + v.v + " to " + hl)
+    return Val(v: op(g, "bitcast " + hl + " " + h + " to " + ol), ty: out)
+}

src/linker.fch

@@ -391,7 +391,8 @@ fn build(Build b, platform.Target tg, cimport.CImports ci) -> int {
         target = " --target=" + tg.triple
     }
     log := b.out + ".log"
-    obj := irObject(clang, target + " " + b.opt + b.extra, b.ir, b.name + "\n" + tg.triple, log) or {
+    // -ffp-contract=on: only what asks to be fused is (muladd); a * b + c stays two roundings, whatever clang's default
+    obj := irObject(clang, target + " " + b.opt + " -ffp-contract=on" + b.extra, b.ir, b.name + "\n" + tg.triple, log) or {
         print(err)
         host.remove(log)
         return 1

tests/run/floatbits.fch

@@ -0,0 +1,31 @@
+// a * b + c is never fused; fma rounds once; numbers as their bits; f16 and bf16
+
+fn plain(float a, float b, float c) -> float {
+    return a * b + c
+}
+
+fn main() {
+    // 0.1 * 10 is 1.0000000000000000555..., which rounds to 1 before the - 1 (two roundings), but not inside fma
+    print(plain(0.1, 10.0, -1.0), fma(0.1, 10.0, -1.0))
+    print(fma(f32(0.1), 10.0, -1.0), fma(2, 3, 4))
+    m := muladd(1.5, 2.0, 0.25)          // fused or not, exact here
+    print(m, f32x4(1.5).muladd(2.0, 0.25), fma(f32x4(0.1), 10.0, f32x4(-1.0))[0])
+    det := f32x4(0.1) * f32x4(10.0) - f32x4(1.0)
+    print(det[0])
+
+    // bits
+    print(f32_bits(1.0), f64_bits(1.0), f32_bits(-0.0))
+    print(f32_from_bits(0x40490fdb), f64_from_bits(f64_bits(2.5)), f32_from_bits(i32(-1082130432)))
+    u32 inf = 0x7f800000
+    print(f32_from_bits(inf), f32_bits(f32x4(1.0, 2.0, -2.0, 0.5)))
+
+    // f16
+    print(f16_to_f32(0x3c00), f16_to_f32(0xc000), f16_to_f32(0x7bff), f16_to_f32(0x0001))
+    print(f32_to_f16(0.1), f32_to_f16(65504.0), f32_to_f16(1.0e6), f16_to_f32(f32_to_f16(0.333)))
+    print(f16_to_f32(u16x8(0x3c00, 0x4000, 0x4200, 0x4400, 0, 0x8000, 0x3800, 0x3400)))
+
+    // bf16
+    print(bf16_to_f32(0x3f80), bf16_to_f32(0x4049), f32_to_bf16(3.14159))
+    print(f32_to_bf16(f32_from_bits(0x3f808000)), f32_to_bf16(f32_from_bits(0x3f818000)), f32_to_bf16(f32_from_bits(0x7fc00001)))
+    print(f32_to_bf16(f32x4(1.0, -2.0, 0.5, 3.14159)), bf16_to_f32(u16x4(0x3f80, 0xc000, 0x3f00, 0x4049)))
+}

tests/run/floatbits.out

@@ -0,0 +1,13 @@
+0 5.551115123125783e-17
+1.4901161e-08 10
+3.25 f32x4(3.25, 3.25, 3.25, 3.25) 1.4901161e-08
+0
+1065353216 4607182418800017408 2147483648
+3.1415927 2.5 -1
+inf u32x4(1065353216, 1073741824, 3221225472, 1056964608)
+1 -2 65504 5.9604645e-08
+11878 31743 31744 0.3330078
+f32x8(1, 2, 3, 4, 0, -0, 0.5, 0.25)
+1 3.140625 16457
+16256 16258 32704
+u16x4(16256, 49152, 16128, 16457) f32x4(1, -2, 0.5, 3.140625)

tests/run/vectors.fch

@@ -14,7 +14,7 @@ fn dot([]f32 a, []f32 b) -> f32 {
     acc := f32x8(0.0)
     i := 0
     while i + 8 <= a.len {
-        acc = f32x8(a, i).fma(f32x8(b, i), acc)
+        acc = f32x8(a, i).muladd(f32x8(b, i), acc)
         i += 8
     }
     total := acc.sum()