Contents

1 Introduction

Urbit is a deterministic functional operating system whose state is a pure function of its event log (~sorreg-namtyv et al., 2016). An event log replay must recover the state exactly, on any conforming runtime, indefinitely. This is not a performance nicety; it is the property that makes the system auditable, upgradeable, and peer-verifiable. Determinism is a load-bearing property all the way down.

Numerical computing is one of the most slippery places for determinism to leak through its abstractions. Standard ieee 754 floating-point arithmetic is not associative, as an example, but the situation is much more dire. The ieee 754 standard (2019) specifies the result of each individual operation under a given rounding mode, but a real numerical program is a composition of operations, and the standard says little about composition. A fused multiply-add computes \(a\cdot b + c\) with a single rounding where two separate operations would round twice; a vectorized dot product sums partial results in an order chosen by the hardware; and a threaded blas reassociates a reduction differently depending on how many cores happened to be available.¹ Each of these is standards-conformant while none is deterministic across hosts.

The treatment of floating-point representation in Hoon by ~lagrev-nocfep (2024) established the foundation for Urbit’s deterministic numerical computing stack. This article expounds both the state of the art of Urbit’s basic numerical packaging and the next layer: dense linear algebra, where the nondeterminism of conventional blas and ieee 754 floating-point mathematics are acute and consequential. We introduce SoftBLAS,² a software-defined blas package; SoftUnum,³ a software-defined universal number (unum) package; and other components of a C-based numerical stack that is deterministic and compliant with the Hoon specification. We further describe the “ray” data structure which facilitates the type dispatch for different numeric representations while preserving identity between Nock-specified atoms and C-recognized arrays.

The central conceit is architectural. In conventional high-performance computing, C (or Fortran) takes primacy: the C implementation is what runs, and correctness is in some sense whatever the C does.⁴ Urbit inverts this relationship. The Hoon/Nock isa implementation is the specification; the C jet is an optimization that is obligated to produce the specification’s exact output for every input, and should be discarded silently whenever it cannot. Determinism then is not merely a testable desideratum but a mandate whose violation the Nock-side architecture cannot express.

2 Background and Literature

Although the solution of partial differential equations is one of the oldest significant applications of computing (Charney, Fjortoft, and Neumann, 1950), no standard library emerged for decades—unsurprising given the churn of hardware architectures, the diversity of numerical methods, the lack of a common interface definition, the rapid evolution of programming languages, and the active development of linear algebraic algorithms.⁵ Given that situation, rigorous bitwise determinism was a concern of vanishingly small importance, whether cross-platform or even for serial runs on the same machine. However, over time the advent of high-performance computing, the emergence of blas as a pseudo-standard, and the increasing importance of reproducibility in numerical computing have made determinism both more pressing and more achievable.

2.1 BLAS as a pseudo-standard

The Basic Linear Algebra Subprograms originate with Lawson, Hanson, Kincaid, and Krogh’s 1979 specification of vector operations on strided arrays (Lawson et al., 1979), later extended to matrix-vector (Level 2) and matrix-matrix (Level 3) routines. The levels track both chronology and (naïve) asymptotic cost: Level 1 is \(O(n)\), Level 2 is \(O(n^2)\), Level 3 is \(O(n^3)\). blas is less a standard than a pseudo-standard—a widely honored set of names and calling conventions (?axpy, ?gemv, ?gemm, with type prefixes s, d, c, z) with many independent implementations. (The original blas predated floating-point standardization as ieee 754 by several years.) That plurality is exactly the problem for a consensus system: two conforming blas libraries routinely disagree in the low bits, and the same library disagrees with itself across machines and thread counts. Despite this fuzziness, blas quickly took its place as the bedrock layer of numerical and scientific computation.

2.2 Sources of floating-point nondeterminism

Cross-host divergence in a blas computation has a small number of recurring causes:

• Reduction order. Floating-point addition is not associative, so a dot product or a gemm inner loop summed in a different order yields a different rounded result; tiled and multithreaded implementations choose that order at run time.
• Fused multiply-add. Hardware that contracts \(a\cdot b + c\) to one rounding produces different bits than two roundings, and whether contraction happens is a compiler- and target-dependent decision.
• Vectorization and instruction selection. simd horizontal sums, x87 80-bit intermediates, and differing libm-style kernels for fast paths all perturb the last bits.
• Transcendental functions. The elementary functions are not specified to the bit by ieee 754 at all, and every platform libm differs in its implementation.

~lagrev-nocfep further enumerated other reasons that contemporary floating-point mathematics cannot be considered deterministic even when single-threaded on a particular hardware target and instruction set architecture.⁶

2.3 Reproducibility

Numerical determinism, or for our purposes synonymously reproducibility, requires that a computation yield the same result to the bit regardless of the platform on which it is carried out. This would not necessarily require the same algorithmic implementation in all cases (since all correct algorithms should converge on the true value given enough runway), but in practice floating-point quirks force this consideration in almost all cases.

Projects like ReproBLAS have yielded reproducible summation (independent of summation order) through error-free transformations and pre-rounding (Demmel and Nguyen, 2013). Intel’s Conditional Numerical Reproducibility offers run-to-run consistency within a fixed code path under documented constraints (Intel Corporation, 2023). The plasma and magma projects (Agullo et al., 2009) took the opposite design stance for performance: dynamic, dependency-driven scheduling that deliberately does not fix reduction order, which is excellent for throughput and fatal for bit-reproducibility. On the scalar side, correctly rounded elementary-function libraries—cr-libm (Daramy et al., 2003), rlibm (Lim and Nagarakatte, 2021), and the core-math project (core-math Project, 2022)—show that the transcendentals can be pinned down to the last bit, at a cost.

Urbit is distinguished by its requirement that determinism be not one option among many but a first-class feature of the OS and runtime. This requirement is structurally enforced by the Nock isa–jetted runtime relationship. Jets utilize relatively faster code paths and reductions but produce the identical result as the slow path: the specification is authoritative. Thus linear algebra techniques, even those following the blas pseudo-standard, must follow a reproducible fast path on all supported architectures.

3 Deterministic Architecture

A Nock computation is a pure function from noun to noun. A jet is a native routine registered against a named Hoon core; when the interpreter evaluates that core, it may run the jet instead of the Nock. The iron discipline that makes this sound is simple: the jet must produce the same noun as the Hoon, for every input, or it is a bug in the runtime. A divergence between jet and specification is termed a jet mismatch and is treated as a correctness failure, not a numerical tolerance.

The Hoon (or its compiled result as Nock) is therefore the specification, in the strong sense that it defines what the right answer is. The C jet is an accelerant with no semantic authority. In some cases, the jet may decline: if a native routine returns the sentinel u3_none, the interpreter falls through and evaluates the Hoon directly. The consequence for numerics is that an arithmetic kind with no jet, or with a jet that punts on some inputs, is not unsupported—it is merely slow. The specification still runs, and the answer is still defined and still deterministic, because Hoon over exact nouns is deterministic by construction.

This crystallizes, or perhaps even inverts, the usual relationship between correctness and performance. In conventional hpc, one has an optimized C or Fortran program, the optimized C arguably is the program and its source reference implementation is documentation (to say nothing of its documentation). There may of course be bugs in the compiled code, but the optimized program is the concrete artifact. By contrast in Urbit, the optimized C program must justify itself against an independent authoritative definition. Determinism follows from two facts: the specification is a pure function of its inputs, and the jet is contractually bit-identical to the specification. A host without the jet computes the same bits more slowly; a divergent jet is non-conforming and must be repaired.

For our linked implementations like SoftBLAS, this means the library is never permitted to be “approximately” the specification. A gemm jet that summed its inner products in a hardware-chosen order would be a mismatch even if every individual rounding were ieee 754-correct, because the Hoon specification sums them in one fixed order and the jet must reproduce that number.⁷ Our blas must satisfy the constraint of determinism to be admissible as a jet at all; thus, SoftBLAS was developed on top of SoftFloat (Hauser, 2018). This yields correctness and consistency rather than hardware-optimized performance.

3.1 The numerical stack

Urbit natively supports ieee 754 floating-point mathematics in its %base distribution, and through the supplied numerical stack additionally supports unums (posits and quires); complex ieee 754; two’s complement integers; and fixed-point mathematics.

The numerical libraries are layered as follows; each layer is Hoon-normative with optional, independently verifiable C jets.

1. Libraries. Imported from /lib at compile time.

   1. /lib/math. The basic ieee 754 arithmetic and transcendental functions library supports four precisions (16-bit @rh, 32-bit @rs, 64-bit @rd, and 128-bit @rq) across a battery of functions like exponentiation and logarithm, trigonometric and inverse-trigonometric functions, and rounding and comparison utilities. Transcendental functions are implemented as Chebyshev minimax polynomials with Cody–Waite range reduction (Cody and Waite, 1980). Jets run on SoftFloat (Hauser, 2018). Originally composed by ~lagrev-nocfep and ~mopfel-winrux.
   2. /lib/twoc. A two’s-complement integer arithmetic package supporting scalar operations. Jetted with an internal library against gmp. Originally composed by ~litlep-nibbyt.
   3. /lib/fixed. Fixed-point (Q-format) arithmetic at power-of-two total widths. Unjetted beyond standard integer arithmetic jets because de facto integer mathematics hold. Originally composed by ~lagrev-nocfep.
   4. /lib/unum. 2022-standard posit arithmetic (Posit Working Group, 2022) at five widths (8-bit @rpb; 16-bit @rph; 32-bit @rps; 64-bit @rpd; 128-bit @rpq, i.e. posit8 through posit128), \(es=\text {2}\) throughout, with quire accumulation. See Section 6.1 on oracle coverage. Jetted with SoftUnum. Originally composed by ~lagrev-nocfep.
   5. /lib/complex. Complex arithmetic over the four ieee 754 precisions, riding on the scalar float layer. Originally composed by ~lagrev-nocfep.
2. Lagoon (Linear AlGebra in hOON). A rank-\(N\) array library (a $ray) carrying a $kind field that dispatches to the appropriate arithmetic: %i754, %uint, %int2, %unum, %cplx, %fixp. Storage is row-major and identical in memory layout to underlying runtime arrays; the %i754 Level-3 path is jetted through SoftBLAS and the %unum path is jetted through SoftUnum. Composed by ~lagrev-nocfep and ~mopfel-winrux.
3. Saloon (Scientific ALgorithms in hOON). Scientific algorithms over Lagoon rays. Types correspond to the Lagoon field. Element-wise transcendentals are dispatched through /lib/math, and dense routines such as the Jacobi symmetric eigendecomposition.⁸ Composed by ~lagrev-nocfep and ~mopfel-winrux.
4. Maroon (MAchine leaRning in hOON). An inchoate machine learning library built atop Lagoon and Saloon.

The single most consequential design choice is Lagoon’s $kind dispatch with fall-through (Listing 1). One basic array type serves six arithmetics; adding a seventh requires only a Hoon implementation, and a jet for it is optional and separately auditable. This is architecturally distinct from NumPy (which supplies one dtype per array with no in-band kind tag); from Julia (which has compile-time specialization per concrete type); and from the apl/J/K tradition (which supply a uniform numeric tower with no explicit kind dispatch). Because an un-jetted kind still runs its Hoon, the specification is authoritative for every kind, jetted or not. Two other features of $ray’s design call attention: the msb pinned 1 bit, which preserves otherwise-stripped leading zeroes in atoms and thus makes the atom representation exactly match the expected C memory layout of an array; and the tail=* arbitrary noun specialization field, which preserves jet axis addresses in case of future extension.⁹

+$  ray               ::  $ray:  n-dimensional array 
  $:  =meta           ::  descriptor 
      data=@ux        ::  data, row-major order 
  == 
5:: 
+$  prec  [a=@ b=@]   ::  fixed-point precision 
+$  meta              ::  $meta:  metadata for a $ray 
  $:  shape=(list @)  ::  list of dimension lengths 
      =bloq           ::  logarithm of bitwidth 
10      =kind           ::  name of data type 
      tail=*          ::  per-kind specialization data 
  == 
:: 
+$  kind              ::  $kind:  type of array scalars 
15  $?  %i754           ::  IEEE 754 float 
      %uint           ::  unsigned integer 
      %int2           ::  2s-complement integer (/lib/twoc) 
      %unum           ::  unum/posit (/lib/unum) 
      %cplx           ::  complex (/lib/complex) 
20      %fixp           ::  fixed-point Q a.b (/lib/fixed) 
  == 
:: 
+$  baum              ::  $baum:  ndray with metadata 
  $:  =meta           :: 
25      data=ndray      :: 
  == 
:: 
+$  ndray             ::  $ndray:  n-dim array as list 
    $@  @             ::  single item 
30    (list ndray)      ::  nonempty list of children

3.2 SoftBLAS and SoftUnum

At the current time, two external numerical libraries facilitate the runtime jets: SoftBLAS and SoftUnum. Both are software-defined, deterministic, and cross-platform.¹⁰

3.2.1 SoftBLAS

SoftBLAS is built on Berkeley SoftFloat 3e (Hauser, 2018). Following SoftFloat, it requires a C99-compliant host and passes all floating-point quantities by value. Every arithmetic primitive used by every routine is a SoftFloat call: there is no path to a hardware floating-point instruction, and therefore no path to hardware-dependent rounding, contraction, or intermediate precision. For reasons enumerated by ~lagrev-nocfep (2024), this is one of the only current methods which guarantees that the output of a floating-point computation is a pure function of its inputs, and therefore deterministic across hosts.¹¹

SoftBLAS covers the four ieee 754 binary precisions used across the stack: half (h), single (s), double (d), and quadruple (q).¹² Routines follow the blas prefixing convention extended to the half and quad types. Real-valued routines across all three levels supply a subset sufficient for Lagoon requirements as of the cited release; complex-valued routines are in progress but sufficient to support Lagoon’s requirements.

The rounding mode is a global, softblas_roundingMode (a char typedef declared in softblas.h) interpreted by the underlying SoftFloat layer.¹³ The admissible values are the SoftFloat 3e rounding modes (round to nearest even, round toward zero, round toward \(+\infty \), round toward \(-\infty \), and round to nearest-away).¹⁴ The default is round to nearest-even, %n. The Hoon specification names the rounding mode in the %base distribution, and the jet is required to honor it.

SoftBLAS satisfies three conditions that conventional blas libraries do not:

1. Fixed reduction order. Every reduction (dot products, gemm inner loops, norms) sums in one canonical order, identical to the order the Hoon specification uses. No tiling or threading is permitted to reassociate it. (Urbit runtimes are, in any case, single-threaded at the time of writing.)
2. No contraction. Multiply-add is two roundings, never a fused one, because the specification rounds twice. There is no fma path, which as ~lagrev-nocfep (2024) exhibited is prone to nondeterminism.
3. Software rounding. All rounding is performed by SoftFloat in software, so it is identical on every host regardless of native floating-point behavior.

3.2.2 SoftUnum

SoftUnum is a software-defined universal number (unum) library that implements the 2022 posit standard (Posit Working Group, 2022). It supports five widths (8, 16, 32, 64, and 128 bits) with \(es=\text {2}\) throughout. The quire accumulator is implemented in software and is bit-identical across hosts. The library is jetted into Lagoon’s %unum kind.

4 Types

The mathematical base types are enumerated in Listing 2. These have occasioned the introduction of several new atom auras, given in Table 2. We have not introduced an aura for two’s-complement integers at the current time, but these will likely nest under @s (signed integer) in the future. The @rqb through @rqq auras are quire accumulators for the corresponding posit widths, and are intended for internal algorithms not general use.¹⁵

+$  kind              ::  $kind:  type of array scalars 
  $?  %i754           ::  IEEE 754 float 
      %uint           ::  unsigned integer 
      %int2           ::  2s-complement integer (/lib/twoc) 
5      %unum           ::  unum/posit (/lib/unum) 
      %cplx           ::  complex (/lib/complex) 
      %fixp           ::  fixed-point Q a.b (/lib/fixed) 
  ==

Aura types have been named after the target width in the cases of complex numbers and quires.

5 Contributions

This deterministic numerical stack contributes several properties to the broader discussion of numerical determinism and reproducibility. Determinism is a structural property of Urbit: while reproducible floating point is a standing goal of several parties (Demmel and Nguyen (2013); Intel Corporation (2023)), Urbit contributes a first-class philosophical requirement.¹⁶

Although not novel, the library covers a broad range of bit widths and precisions, including half-precision and quadruple-precision, which are not widely supported in other numerical computing libraries. Likewise we provided faithfully-rounded elementary functions at most precisions and quantify their accuracy.

Lagoon’s heterogeneous-kind array design with a single $ray type is a novel approach to array libraries, allowing for in-band kind dispatch and fall-through to specification. As mentioned previously, this permits more flexibility and extensibility in the base representation than NumPy, Julia, or the apl/J/K tradition. The kind field allows for easy addition of new numeric types without requiring changes to the underlying array structure, and the fall-through mechanism ensures that un-jetted kinds can still be computed correctly using the Hoon specification.

6 Benchmarks

While a software-defined floating-point library will never match the performance of a hardware-accelerated one, we have nevertheless seen impressive speedups by refining the algorithms and the runtime jets. The benchmarks reported here are for a single-threaded, single-core execution on a modern 64-bit compiler cpu (AArch64, Apple M4 Max). They should serve as order-of-magnitude indications on contemporary hardware.

The first benchmark results tabulate what the jetted system buys over interpreting the Hoon specification directly. Table 3 reports jetted versus interpreted time for ieee 754 /lib/math. For most trials, 100,000 runs were evaluated and averaged.¹⁷ For the very slowest functions (e.g. half-precision Taylor series inverse sine), this number was modulated downwards substantially.¹⁸ The original /lib/math and Saloon implementations of transcendentals were naïve Taylor series expansions and serve as a worst-case demonstration. Accuracies are comparable for the number of Taylor series terms utilized to achieve an appropriate relative tolerance, but the Chebyshev basis functions are uniformly faster and more correct. All values reported are measured within a running Urbit ship (rather than from an independent test harness).

The /lib/unum transcendentals are currently implemented as Chebyshev basis expansions like the /lib/math operations. Results are included in Table 4.

6.1 Verification

The determinism architecture reduces correctness to one question per routine: does the jet equal the specification on every input? We answer it insofar as possible with oracles external to both the Hoon and C implementations. SoftFloat by and large served as the oracle for the ieee 754 floating-point operations (real and complex), while the posit mathematics were verified against two independent references: SoftPosit (Leong, 2018) (posit8, posit16, posit32) and the Universal library (Omtzigt et al., 2020) (posit64, posit128).

The quire’s defining property—exact accumulation until a single final rounding—is checked directly against an exact rational sum.

Elementary functions are checked against mpfr (Fousse et al., 2007) at high precision. At half precision the input space is \(\text {2}^\text {16}\) values and the test is exhaustive; at single precision the \(\text {2}^\text {32}\) space is likewise tested in full. At double and quadruple precision, where exhaustive testing is infeasible, we report maximum observed unit in the last place (ulp) error over a sampled and edge-weighted input set rather than a guarantee. We claim a faithful-rounding (\(\leq \text {1}\) ulp) result where the testing supports it and a measured bound where it does not.

We report the errors of the transcendental functions in terms of ulp error, which is a standard measure of floating-point accuracy.¹⁹ The maximum observed ulp error for each function is reported against an independent oracle, which is typically a high-precision library such as mpmath. The results are summarized in Tables 5 and 6.

The libraries instrumented with the Chebyshev basis are designed to produce “faithfully rounded functions” rather than “correctly rounded functions” for composed operations, while dedicated (non-composed) kernels frequently do better still. The distinction is that a faithfully rounded result is within 1 ulp of the exact result (Rump, Ogita, and Oishi, 2008), while a correctly rounded result is the exact representable value closest to the exact result. Faithful rounding is often sufficient for numerical applications and is generally faster to compute than correct rounding. Truly correct rounding is often prohibitively expensive, especially for transcendental functions, and can require much higher intermediate precision to attain.²⁰ Tables 5 and 6 report the maximum observed ulp error of every transcendental arm against an independent mpmath-checked oracle.²¹

In Tables 5 and 6, the superscripts carry the following meanings:

• ^c marks a composed arm,²² whose error is inherited from its constituents’ chained roundings rather than directly controlled.
• ^e marks an exhaustive sweep over the entire input space. (Unmarked cells describe sampled ranges.)
• ^f marks where an imprecision in a previous version of the Chebyshev fit led to massive error in the reduced argument, in which the range reduction quotient for @rh was not representable in half precision (Table 7). This was replaced by widening the intermediate result to @rs, producing a faithful rounding for @rh as well. As this bug results from a naïve application of the Chebyshev fit to half-precision range reduction, we report it as a cautionary aside.
• ^k marks dedicated kernels for +tan (replacing the composed sin/cos ratio used at @rh/@rq). A dedicated kernel was investigated but found to regress @rh (the 11-bit mantissa leaves insufficient headroom) and to need a different technique entirely at @rq (the Chebyshev series coefficients grow rather than shrink at high degree causing numerical instability).
• ⁿ marks /lib/unum’s +cbt, which returns NaR for negative inputs and so is swept over positive inputs only.

The original Chebyshev-based implementation of the @rh transcendental functions had a subtle but catastrophic flaw in its range reduction. The @rh Chebyshev fits were correct, but the range reduction was performed entirely in half precision, which is insufficient for large arguments. @rh’s trigonometric arms did quarter-turn range reduction (\(q=\mathrm {round}(|x|\cdot \frac {\text {2}}{\pi })=|x|\cdot \frac {\text {4}}{\tau }\), \(r=|x|-q\cdot (\frac {\pi }{\text {2}})=|x|-q\cdot (\frac {\tau }{\text {4}})\)) entirely in native @rh arithmetic, giving every reduction constant and every intermediate product only an 11-bit-mantissa half-precision value. Half precision guarantees exact integer representation only up to 2,048; once \(q\) exceeds that value (\(|x|\) past \(\sim \)500), \(q\)’s own rounding error was amplified by the \(\frac {\pi }{\text {2}}=\frac {\tau }{\text {4}}\) leading term, injecting a multi-radian error into the reduced remainder. Table 7 reports the resulting damage, measured before the fix: catastrophic, and well within @rh’s ordinary input range (max magnitude \(\sim \)65,504), not an edge case. The fix widens to @rs (24-bit mantissa, exact \(q\) up to 2²⁴, vastly beyond the dynamic range of @rh) via +widen-hs (exact) and +narrow-sh (correctly rounded), primitives. These were verified exhaustively.

Every /lib/unum transcendental with a dedicated, non-composed kernel is exactly correctly rounded (0 ulp) at all three tested widths, a stronger guarantee than the ieee 754 /lib/math achieves at any precision, where the best case is faithful rounding around 0.5–1 ulp rather than 0 ulp. The composed arms (+pow, +cbt, and +pow-n beyond small positive exponents) lose this property in both libraries alike, landing in the 1–134 ulp range depending on precision and the number of chained roundings the composition involves—an inherent cost of composition rather than a library-specific weakness.

6.2 The price of determinism

The third field of benchmarking measures the overhead of the determinism guarantee relative to a conventional, nondeterministic, hardware blas. We report SoftBLAS against Apple’s Accelerate framework—a stock blas on macOS which dispatches to Apple Silicon’s amx matrix coprocessor—for sgemm and dgemm (single-threaded to remove Accelerate’s automatic multithreading as a confound). We expect software-defined SoftBLAS to be slower; here we quantify how that cost scales. The results are summarized in Table 8 and Figure 1.

^*At \(n\leq \text {20}\) Accelerate’s own run-to-run standard deviation is comparable to its mean, so the reported slowdown there is noise-dominated; the trend is reliable and monotonic above \(n=\text {30}\) (SoftBLAS’s timings are stable throughout since a fixed triple loop has no comparable startup-dominated regime).

Determinism does not levy a fixed tax, but its cost grows with problem size, from single digits at \(n=\text {10}\) to nearly four orders of magnitude by \(n=\text {1,000}\) (sgemm) or over three orders of magnitude (dgemm). SoftBLAS’s triple-nested loop over SoftFloat calls is a clean \(O(n^3)\) throughout the range tested (doubling \(n\) consistently multiplies time by \(\sim \!8\times \)), with no cache blocking or vectorization—each scalar multiply-add pays a full software floating-point call.²³ Accelerate is sub-cubic through the moderate range (cache blocking and simd/amx dispatch absorb work that would otherwise scale as \(n^3\)) before trending back toward cubic by \(n=\text {1,000}\). In throughput terms, SoftBLAS delivers roughly 0.12 gflop/s at \(n=\text {1,000}\) regardless of precision—consistent with SoftFloat call overhead, not memory bandwidth, being the bottleneck—against Accelerate’s approximately 1,477 gflop/s (sgemm) and 401 gflop/s (dgemm), the latter gap itself illustrating that Accelerate’s single-precision path better exploits the amx unit’s native format.

This sharply demonstrates a quantitative trade-off between determinism and performance introduced in Section 3.1: SoftBLAS is not merely “slower,” but slower by a margin that widens with the problem sizes—large dense linear algebra—where blas exists to be fast in the first place. A deterministic numerical stack should budget accordingly. A software blas that never touches a hardware floating-point instruction cannot reach a coprocessor-backed blas’s throughput, only close the distance to it, which argues in favor of new paradigms for deterministic numerical computing that can leverage hardware acceleration without sacrificing reproducibility.

7 Future Work

The complex-valued SoftBLAS routines are in progress and will complete the Level 1–3 coverage. Quire operations in /lib/unum are implemented within the generic posit core but are not yet surfaced as callable top-level arms; exposing them is a prerequisite both for user access and for routing them through the verification harness. The unum specification also defines valids, which are arithmetic intervals which have not yet seen widespread adoption.

Saloon will grow in methods adjacent to scientific computing: conjugate gradient minimization and other optimization algorithms; random number generators; and the like. The fixed-point library is relatively underdeveloped compared to the other arithmetics, and will be expanded to support more general Q-format widths and operations. More components of the various libraries need to be upgraded to a state-of-the-art Chebyshev basis and benchmarked for performance.²⁴

Finally, Maroon, the machine learning library, is in its infancy and will be developed to support a range of machine learning algorithms and techniques. Part of this field of consideration encompasses the utilization of gpu programming for llm-based computing within or adjacent to Urbit. gpu numerical nondeterminism and host portability presents a challenge of the same order of magnitude as the baseline reproducibility we have tackled over the past several years in this project. One can envision a ReproCL as a deterministic counterpart to OpenCL, but it is unclear whether sufficient warrant exists to justify the effort as of the mid-2020s. Heterogeneous compute would be benefitted by the adoption of alternatives to ieee 754 with better guarantees of portability: not just unums but bfloat16 and other newer types may answer the need.

8 Conclusion

Urbit as a deterministic state machine cannot treat numerical determinism or bitwise reproducibility as merely a tolerance. The architectural stack we have described here constitutes a partial answer to the challenge underlying ~lagrev-nocfep (2024): how can a numerical computation built on a particular stack be reliably reproducible? We maintain this determinism by painstaking construction. The computational cost of determinism is currently (and regrettably) non-negligible, while one can hope for future developments, protocols, and architectures that make it more natural and tenable.

References

Agullo, Emmanuel et al. (2009). “Numerical Linear Algebra on Emerging Architectures: The plasma and magma Projects.” In: Journal of Physics: Conference Series 180, p. 012037. doi: 10.1088/1742-6596/180/1/012037.

Charney, J., R. Fjortoft, and J. von Neumann (1950). “Numerical Integration of the Barotropic Vorticity Equation.” In: Tellus 2.4, pp. 237–254. doi: 10.1111/j.2153-3490.1950.tb00336.x. url: https://doi.org/10.1111/j.2153-3490.1950.tb00336.x.

Cody, William J. and William M. Waite (1980). Software Manual for the Elementary Functions. Prentice-Hall.

core-math Project (2022) “core-math: Correctly Rounded Mathematical Functions”. url: https://core-math.gitlabpages.inria.fr/.

Daramy, Catherine et al. (2003). “cr-libm: A Correctly Rounded Elementary Function Library.” In: Proc. spie 5205, Advanced Signal Processing Algorithms, Architectures, and Implementations xiii. doi: 10.1117/12.505591.

Demmel, James and Hong Diep Nguyen (2013). “Fast Reproducible Floating-Point Summation.” In: 2013 ieee 21st Symposium on Computer Arithmetic (arith). ReproBLAS project, pp. 163–172. doi: 10.1109/ARITH.2013.9.

Fousse, Laurent et al. (2007). “mpfr: A Multiple-Precision Binary Floating-Point Library with Correct Rounding.” In: acm Transactions on Mathematical Software 33.2, 13:1–13:15. doi: 10.1145/1236463.1236468.

Free Software Foundation (2024). The gnu C Library Reference Manual. Section 19.7, Known Maximum Errors in Math Functions. url: https://www.gnu.org/software/libc/manual/html_node/Errors-in-Math-Functions.html.

Hauser, John R. (2018) “Berkeley SoftFloat, Release 3e”.

ieee (2019). ieee Standard for Floating-Point Arithmetic. ieee Std 754-2019. ieee. doi: 10.1109/IEEESTD.2019.8766229.

Intel Corporation (2023) “Conditional Numerical Reproducibility (cnr) in Intel oneAPI Math Kernel Library”.

~lagrev-nocfep, N. E. Davis (2024). “The Desert of the Reals: Floating-Point Arithmetic on Deterministic Systems.” In: Urbit Systems Technical Journal 1.1, pp. 93–131.

— (2026a). “Exact Uniform Random Variates for Tapered-Precision Posit Arithmetic.” In: The Seventh Conference on Next Generation Arithmetic 2026 (CoNGA’26). accepted for publication.

— (2026b) “SoftUnum: Software-Defined Universal Number Library for Deterministic Jetting”. url: https://github.com/sigilante/SoftUnum.

Lawson, C. L. et al. (1979). “Basic Linear Algebra Subprograms for Fortran Usage.” In: acm Transactions on Mathematical Software 5.3, pp. 308–323. doi: 10.1145/355841.355847.

Leong, S. H. (Cerlane) (2018) “SoftPosit: A Software Implementation of Posit Arithmetic”.

Lim, Jay P. and Santosh Nagarakatte (2021). “High Performance Correctly Rounded Math Libraries for 32-bit Floating Point Representations.” In: Proc. 42nd acm sigplan Conf. on Programming Language Design and Implementation (pldi), pp. 359–374. doi: 10.1145/3453483.3454049.

Muller, Jean-Michel et al. (2018). Handbook of Floating-Point Arithmetic, 2nd edition. Birkhäuser Boston, p. 632.

Omtzigt, E. Theodore L. et al. (2020). “Universal Numbers Library: Design and Implementation of a High- Performance Reproducible Number Systems Library.” In: arXiv preprint arXiv:2012.11011. url: https://arxiv.org/abs/2012.11011.

Posit Working Group (2022). Standard for Posit Arithmetic (2022). Posit Working Group. url: https://posithub.org/docs/posit_standard-2.pdf.

Rump, Siegfried M., Takeshi Ogita, and Shin’ichi Oishi (2008). “Accurate Floating-Point Summation Part I: Faithful Rounding.” In: siam Journal on Scientific Computing 31.1, pp. 189–224. doi: 10.1137/050645671. url: https://doi.org/10.1137/050645671.

~sorreg-namtyv, Curtis Yarvin et al. (2016). Urbit: A Solid-State Interpreter. Whitepaper. Tlon Corporation. url: https://media.urbit.org/whitepaper.pdf (visited on ~2024.1.25).

Urbit Foundation (2024) “SoftBLAS: Software-Defined blas for Deterministic Jetting”. url: https://github.com/urbit/SoftBLAS.

Footnotes