diff --git a/lagoon/README.md b/lagoon/README.md index 9bc59e8..18ce236 100644 --- a/lagoon/README.md +++ b/lagoon/README.md @@ -23,15 +23,18 @@ We envision four libraries and associated jet code living in this repository: ## Type System -- `%real` IEEE 754 float (currently supported) -- `%cplx` IEEE 754 float/BLAS packed complex (planned) -- `%uint` unsigned integers (currently supported) -- `%int2` signed twos-complement integers (planned) -- `%sint` signed ZigZag integers (planned) -- `%unum` → subdivide into one of: - - `%post` posits (32-bit or 16-bit) (planned) - - `%vald` valids (32-bit or 16-bit) (planned) -- `%fixp` fixed-precision (planned) +The element `kind` lives in `/sur/lagoon` (the old `%real` is now `%i754`): + +- `%i754` IEEE 754 float — `@rh`/`@rs`/`@rd`/`@rq` (supported) +- `%uint` unsigned integers (supported) +- `%int2` signed two's-complement integers, `/lib/twoc` (supported) +- `%unum` unum/posits — `@rpb`/`@rph`/`@rps`/`@rpd`, `/lib/unum` (supported) +- `%cplx` BLAS-packed complex — `@ch`/`@cs`/`@cd`/`@cq`, `/lib/complex` (supported) +- `%fixp` fixed-point Q a.b, `/lib/fixed`; precision `[a b]` in `meta.tail` (supported) + +(`%sint` ZigZag integers and `%vald` valids remain possible future additions.) +Arms added since the `%real`-only release below include `++dotc` (Hermitian dot) +and `++conj` (elementwise conjugate); Saloon adds eigendecomposition (`++eig`). ## Fixed-Point Library `/lib/fixed` diff --git a/lagoon/desk/lib/lagoon.hoon b/lagoon/desk/lib/lagoon.hoon index 2f76dd5..dc1ed7b 100644 --- a/lagoon/desk/lib/lagoon.hoon +++ b/lagoon/desk/lib/lagoon.hoon @@ -355,6 +355,10 @@ (^rip a b) -- :: + :: +check: validate a ray's data length against its shape. The data atom + :: carries a 0x1 most-significant "pin" above the elements (so leading-zero + :: elements aren't lost), hence (met bloq data) counts prod(shape)+1 blocks; + :: this asserts prod(shape) == (met ...) - 1. Every public arm calls it. ++ check |= =ray ^- ? @@ -426,8 +430,9 @@ :: %i754 ?+ kind !! - :: %i754 -> %uint - :: XXX will incorrectly convert negative values to %uint + :: %i754 -> %uint. KNOWN LIMITATION: a negative %i754 value has no + :: %uint representation; +toi yields a negative @s and the (div ... 2) + :: wraps it rather than clamping/erroring. Callers must pre-clamp. %uint %- en-ray :- [shape.meta.ray bloq %uint tail.meta.ray] diff --git a/lagoon/desk/sur/lagoon.hoon b/lagoon/desk/sur/lagoon.hoon index 18210e4..adc6cb0 100644 --- a/lagoon/desk/sur/lagoon.hoon +++ b/lagoon/desk/sur/lagoon.hoon @@ -12,7 +12,10 @@ $: shape=(list @) :: list of dimension lengths =bloq :: logarithm of bitwidth =kind :: name of data type - tail=* :: placeholder for future data (jet convenience) + tail=* :: per-kind specialization data (~ unless a kind needs + :: it). Deliberately last and untyped so a kind can add + :: data without rewriting the jet axes. %fixp uses it for + :: the [a=@ b=@] fixed-point precision; other kinds: ~. == :: +$ kind :: $kind: type of array scalars @@ -20,7 +23,7 @@ %uint :: unsigned integer %int2 :: 2s-complement integer (/lib/twoc) %unum :: unum/posit (/lib/unum) @rpb @rph @rps @rpd - %cplx :: complex (/lib/complex) @cs/@cd; bloq = total width + %cplx :: complex (/lib/complex) @ch/@cs/@cd/@cq; bloq=total width %fixp :: fixed-point Q a.b (/lib/fixed); prec [a b] in meta.tail == :: diff --git a/libmath/desk/lib/math.hoon b/libmath/desk/lib/math.hoon index 80be23d..5676f62 100644 --- a/libmath/desk/lib/math.hoon +++ b/libmath/desk/lib/math.hoon @@ -92,10 +92,11 @@ ++ huge `@rs`0x7f80.0000 :: 3.40282346638528859812e+38 :: +tiny: @rs :: - :: Returns the value of the smallest representable normal number. + :: Returns the smallest representable positive (subnormal) number, 2^-149. + :: (Not the smallest NORMAL, which is 2^-126 = .1.1754944e-38.) :: Examples :: > tiny - :: .1.1754944e-38 + :: .1e-45 :: Source ++ tiny `@rs`0x1 :: 1.40129846432481707092e-45 :: @@ -874,22 +875,20 @@ :: +round-bankers: @rs -> @rs :: :: Returns the floating-point atom rounded to the nearest integer, with - :: ties rounded to the nearest even integer. + :: ties rounded to the nearest even integer (banker's rounding). This is + :: exactly what +toi does under round-nearest (%n), which ties to even, so we + :: force %n regardless of the door's configured mode. :: Examples - :: > (round-bankers .1) - :: .1 :: > (round-bankers .1.5) :: .2 + :: > (round-bankers .2.5) + :: .2 :: > (round-bankers .1.49) :: .1 :: Source ++ round-bankers |= x=@rs ^- @rs - =/ int (san (need (toi x))) - =/ dcm (sub x int) - ?: (lth dcm .0.5) - int - (add int .1) + (san (need (~(toi ^rs %n) x))) -- :: double precision ++ rd @@ -1556,10 +1555,10 @@ |= z=@rd ^- @rd :: filter out non-finite arguments :: check infinities - ?: =(z 0x7ff0.0000.0000.0000) `@rd`0x7ff0.0000.0000.0000 :: exp(+inf) -> inf - ?: =(z 0xfff0.0000.0000.0000) .~0.0 :: exp(-inf) -> 0 + ?: =(z 0x7ff0.0000.0000.0000) `@rd`0x7ff0.0000.0000.0000 :: log(+inf) -> inf + ?: =(z 0xfff0.0000.0000.0000) `@rd`0x7ff8.0000.0000.0000 :: log(-inf) -> NaN :: check NaN - ?. (^gte (dis 0x7ff8.0000.0000.0000 z) 0) `@rd`0x7ff8.0000.0000.0000 :: exp(NaN) -> NaN + ?. (^gte (dis 0x7ff8.0000.0000.0000 z) 0) `@rd`0x7ff8.0000.0000.0000 :: log(NaN) -> NaN :: otherwise, use Taylor series =/ p .~0 =/ po .~-1 @@ -1742,23 +1741,21 @@ (div rnd-mantissa scaling) :: +round-bankers: @rs -> @rs :: - :: Returns the floating-point atom rounded to the nearest integer, with - :: ties rounded to the nearest even integer. + :: Returns the floating-point atom rounded to the nearest integer, with ties + :: rounded to the nearest even integer (banker's rounding). This is exactly + :: what +toi does under round-nearest (%n), which ties to even, so we force + :: %n regardless of the door's configured mode. :: Examples - :: > (round-bankers .1) - :: .1 - :: > (round-bankers .1.5) - :: .2 - :: > (round-bankers .1.49) - :: .1 + :: > (round-bankers .~1.5) + :: .~2 + :: > (round-bankers .~2.5) + :: .~2 + :: > (round-bankers .~1.49) + :: .~1 :: Source ++ round-bankers |= x=@rd ^- @rd - =/ int (san (need (toi x))) - =/ dcm (sub x int) - ?: (lth dcm .~0.5) - int - (add int .~1) + (san (need (~(toi ^rd %n) x))) -- :: half precision ++ rh @@ -1849,10 +1846,11 @@ ++ huge `@rh`0x7bff :: 6.55e+04 :: +tiny: @rh :: - :: Returns the value of the smallest representable normal number. + :: Returns the smallest representable positive (subnormal) number, 2^-24. + :: (Not the smallest NORMAL, which is 2^-14 = .~~6.10e-05.) :: Examples :: > tiny - :: .~~6.10e-05 + :: .~~6e-8 :: Source ++ tiny `@rh`0x1 :: 6e-08 :: diff --git a/libmath/desk/tests/lib/math-rounding.hoon b/libmath/desk/tests/lib/math-rounding.hoon new file mode 100644 index 0000000..467ab19 --- /dev/null +++ b/libmath/desk/tests/lib/math-rounding.hoon @@ -0,0 +1,23 @@ +/+ *test, math +:::: /tests/lib/math-rounding -- round-bankers (ties-to-even) + log(-inf) +:: +^| +|% +:: round-bankers ties to the nearest EVEN integer (banker's rounding). +++ test-round-bankers ^- tang + ;: weld + %+ expect-eq !>(`@rs`.2) !>((~(round-bankers rs:math [%n .1e-5]) .2.5)) :: -> 2 (even) + %+ expect-eq !>(`@rs`.2) !>((~(round-bankers rs:math [%n .1e-5]) .1.5)) :: -> 2 (even) + %+ expect-eq !>(`@rs`.4) !>((~(round-bankers rs:math [%n .1e-5]) .3.5)) :: -> 4 (even) + %+ expect-eq !>(`@rs`.1) !>((~(round-bankers rs:math [%n .1e-5]) .1.49)) :: -> 1 + %+ expect-eq !>(`@rd`.~2) !>((~(round-bankers rd:math [%n .~1e-10]) .~2.5)) + == +:: log(-inf) is NaN (was 0.0 in the rd door -- a copy-paste from exp). +++ test-log-ninf ^- tang + ;: weld + %+ expect-eq !>(`@rd`0x7ff8.0000.0000.0000) + !>((~(log rd:math [%n .~1e-10]) `@rd`0xfff0.0000.0000.0000)) + %+ expect-eq !>(`@rs`0x7fc0.0000) + !>((~(log rs:math [%n .1e-5]) `@rs`0xff80.0000)) + == +-- diff --git a/saloon/README.md b/saloon/README.md index bb13090..69a98d0 100644 --- a/saloon/README.md +++ b/saloon/README.md @@ -36,6 +36,17 @@ Logical functions: - `all-close` (pass-through from Lagoon `++all`) - `any-close` (pass-through from Lagoon `++any`) +Linear algebra (over Lagoon arrays): + +- `++eig`, eigendecomposition of a **symmetric** (`%i754`) or **Hermitian** + (`%cplx`) matrix via cyclic Jacobi → `[vals=ray vecs=ray]` (real eigenvalues, + orthonormal/unitary eigenvectors as columns). Dispatches on `kind`. +- `++eigvals`, eigenvalues only (1-D ray). +- `++eigvecs`, eigenvectors only. + +Set rounding mode and tolerance with `++sake` before calling `++eig` (the bare +`++sa` default `rtol` is unusable); `rtol`'s width must match the component. + ## References - Milton Abramowitz & Irene Stegun, _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_. 1964–2010. diff --git a/saloon/desk/lib/saloon.hoon b/saloon/desk/lib/saloon.hoon index 24462b6..f016977 100644 --- a/saloon/desk/lib/saloon.hoon +++ b/saloon/desk/lib/saloon.hoon @@ -80,7 +80,7 @@ :: Returns the BOOLEAN comparison of two floating-point rays, greater than :: Examples :: > (gth:sa:sa (ones:la:la [~[5 1] 5 %i754 ~]) (en-ray:la:la [~[5 1] 5 %i754 ~] ~[.1 .0 .1 .2 .0])) - :: [meta=[shape=~[5 1] bloq=5 kind=%i754 fxp=~] data=data=0x1.0000.0000.3f80.0000.0000.0000.0000.0000.3f80.0000] + :: [meta=[shape=~[5 1] bloq=5 kind=%i754 fxp=~] data=0x1.0000.0000.3f80.0000.0000.0000.0000.0000.3f80.0000] :: > ;;((list (list @rs)) data:(de-ray:la:la (gth:sa:sa (ones:la:la [~[5 1] 5 %i754 ~]) (en-ray:la:la [~[5 1] 5 %i754 ~] ~[.1 .0 .1 .2 .0])))) :: [i=~[.0] t=[i=~[.1] t=~[~[.0] ~[.0] ~[.1]]]] :: Source @@ -90,7 +90,7 @@ :: Returns the BOOLEAN comparison of two floating-point rays, greater than or equal to :: Examples :: > (gte:sa:sa (ones:la:la [~[5 1] 5 %i754 ~]) (en-ray:la:la [~[5 1] 5 %i754 ~] ~[.1 .0 .1 .2 .0])) - :: [meta=[shape=~[5 1] bloq=5 kind=%i754 fxp=~] data=data=0x1.3f80.0000.3f80.0000.3f80.0000.0000.0000.3f80.0000] + :: [meta=[shape=~[5 1] bloq=5 kind=%i754 fxp=~] data=0x1.3f80.0000.3f80.0000.3f80.0000.0000.0000.3f80.0000] :: > ;;((list (list @rs)) data:(de-ray:la:la (gte:sa:sa (ones:la:la [~[5 1] 5 %i754 ~]) (en-ray:la:la [~[5 1] 5 %i754 ~] ~[.1 .0 .1 .2 .0])))) :: [i=~[.1] t=[i=~[.1] t=~[~[.1] ~[.0] ~[.1]]]] :: Source @@ -101,14 +101,14 @@ :: Alias for +gte. :: Examples :: > (geq:sa:sa (ones:la:la [~[5 1] 5 %i754 ~]) (en-ray:la:la [~[5 1] 5 %i754 ~] ~[.1 .0 .1 .2 .0])) - :: [meta=[shape=~[5 1] bloq=5 kind=%i754 fxp=~] data=data=0x1.3f80.0000.3f80.0000.3f80.0000.0000.0000.3f80.0000] + :: [meta=[shape=~[5 1] bloq=5 kind=%i754 fxp=~] data=0x1.3f80.0000.3f80.0000.3f80.0000.0000.0000.3f80.0000] :: > ;;((list (list @rs)) data:(de-ray:la:la (geq:sa:sa (ones:la:la [~[5 1] 5 %i754 ~]) (en-ray:la:la [~[5 1] 5 %i754 ~] ~[.1 .0 .1 .2 .0])))) :: [i=~[.1] t=[i=~[.1] t=~[~[.1] ~[.0] ~[.1]]]] :: Source ++ geq gte:(lake rnd) :: +neq: [$ray $ray] -> $ray :: - :: Returns the BOOLEAN comparison of two floating-point rays, equal to + :: Returns the BOOLEAN comparison of two floating-point rays, not equal to :: Examples :: > (neq:la:la (ones:la:la [~[5 1] 5 %i754 ~]) (en-ray:la:la [~[5 1] 5 %i754 ~] ~[.1 .0 .1 .2 .0])) :: [meta=[shape=~[5 1] bloq=5 kind=%i754 fxp=~] data=0x1.0000.0000.3f80.0000.0000.0000.3f80.0000.3f80.0000]