crnumerics 0.10.0

crdist adds double-precision uniform and exponential draws from its
Philox stream. crdist_uniform gives (2m + 1)·2⁻⁵³ from two words, exact,
and crdist_exponential the correctly rounded −log of it. Like crdist's
other samplers they are specified bit for bit, so a double-precision
simulation, such as the Gillespie algorithm, gets the same draws on every
machine; crdist-check computes them again with MPFR and compares every
bit. In Python, Stream.uniform and Stream.exponential take
dtype=np.float64, and pip install crdistr brings them.

A new example, examples/ssa-xplat, runs four classic reaction networks
through libRoadRunner, COPASI and GillesPy2 on eleven platforms, beside a
reference built on crdist. The stochastic trajectories were the same
everywhere. The same models solved as ODEs differed between x86 and Arm,
and COPASI's ODE results differed from one run to the next on macOS and
Windows.

ival, crlin, crfix, crnn and lowp are unchanged but for the version
strings, now 0.10.0.

crnumerics 0.9.0

crfix adds SiLU and the exact GELU (v·Φ(v), not the tanh approximation),
correctly rounded in fixed point and in TFLite's affine formats like its
other functions, and checked against MPFR on every 16-bit input. It adds
layer normalization computed exactly: each result is the exact value
rounded once, settled by comparing integer products, with no
approximation anywhere. And it adds unsigned containers, TFLite's uint8
among them, for results and for lookup tables.

pip install crfix brings all of it to Python. Its tests check SiLU, GELU
and layer normalization against a second reference, independent of MPFR.
The Cortex-M0 build still has no floating point and gives the host's
bits.

The Python packages carry keywords, classifiers, and links to the project
page and changelog on PyPI: crfix from this release, the others from
their next.

ival, crdist, crlin, crnn and lowp are unchanged but for the version
strings, now 0.9.0.

crnumerics 0.8.0

crfix is on PyPI: pip install crfix gives its functions, lookup tables,
TFLite's affine formats and softmax on NumPy integer arrays. Its tests
use a second reference, independent of MPFR: Python's decimal and exact
fractions. Array functions in C (crfix_eval_n, crfix_eval_affine_n,
crfix_div_n) back it, and python -m crfix table prints a lookup table
as C source.

The affine formats and softmax have fast paths: an int8 lookup table in
13 to 26 µs instead of 0.9 to 2.4 ms, and softmax at 0.1 µs an element,
with the same results. Whatever the 64-bit enclosures leave undecided
still goes to the slow path.

crfix runs on a Cortex-M0. The build has no floating point and no
instruction beyond ARMv6-M, and a digest of 23 million results matches
the host's on 32-bit soft-float ARM and on the M0 itself, under qemu.

Two fixes to 0.7.0's affine formats and softmax:

• With an output zero point, CRFIX_ZERO rounded y + zero point toward 0
  rather than y, as crfix.h specifies. With TFLite's logistic output,
  most inexact entries were one too large. The other modes were right.
• exp and exp2 values from 2^62 up were saturated whatever the output
  scale, which was wrong for output scales above 2^30.

ival, crdist, crlin, crnn and lowp are unchanged but for the version
strings, now 0.8.0.

crnumerics 0.7.0

crfix is new: exp, exp2, log, log2, sigmoid, tanh, sqrt, rsqrt,
reciprocal and division on fixed-point numbers, correctly rounded in four
modes and saturated, in integer arithmetic only, so it runs where there
is no FPU. Every 16-bit input, every pair of formats and every mode is
checked against MPFR, on x86-64 and aarch64. The few inputs a fast 64-bit
enclosure cannot decide go to a slow path of interval arithmetic that is
correct by construction, so 32-bit inputs are correctly rounded too.
It also takes TFLite's affine formats, any scale and zero point:
correctly rounded int8 tables for logistic, tanh, exp and log in about
2 ms each, and softmax.

crlin solves banded and sparse systems: crlin_solve_band and
crlin_solve_sparse, and solve_banded and solve_sparse in Python in
scipy's layouts. A proven lower bound on the smallest eigenvalue of A,
or of A^T A, takes the place of an inverse. pip install crlin brings
them.

crsum rounds short sums three times faster, with the same results.

ival, crdist, crnn and lowp are unchanged but for ival's and crdist's
version strings, now 0.7.0; ival's soname is still libival.so.0.

crnumerics 0.6.0

crlin is new: verified linear algebra in binary64. crlin_solve returns,
for each component of the solution of A x = b, an interval proved to
hold the exact solution of the system as given, typically an ulp wide,
or says that the proof did not go through (A singular, or too
ill-conditioned). It is Rump's method with every bound taken from exact
sums and products or from an a priori error bound in round-to-nearest,
so it never switches the rounding mode, and a threaded BLAS can't break
it: any BLAS's dgemm, or a GPU's, can do the heavy products. In the
default mode a solve costs about 10 times an unverified one.

Beyond one system: several right-hand sides, interval data, least
squares and minimum-norm solutions, very ill-conditioned systems (to
condition numbers past 1e30), proofs of positive definiteness, and an
enclosure of every eigenvalue of a symmetric matrix. Each is checked
against exact rational answers.

pip install crlin brings the Python binding over NumPy, with wheels for
Linux (x86-64, aarch64) and macOS (Apple Silicon).

ival, crsum, crdist, crnn and lowp are unchanged but for ival's and
crdist's version strings, now 0.6.0; ival's soname is still libival.so.0.

crnumerics 0.5.0

crdist's samplers gain Poisson and binomial counts, and crdist now builds and runs on macOS (Apple Silicon), where pip install crdistr gets a ready-built wheel.

Poisson and binomial: inversion for small means, then Hormann's transformed rejection (PTRS for Poisson, as NumPy uses it; BTRS for the binomial, as TensorFlow and JAX use it), each specified in crdist.h as a fixed sequence of correctly rounded operations, so every machine gives the same counts for the same stream. The check recomputes the specification with MPFR and compares every draw, and tests the distributions against the exact probabilities.

crdist-check hash-all hashes every output of every function on all 2^32 binary32 inputs and compares with the digests of the builds that were checked against MPFR on every input. Equal digests carry that proof to another compiler or machine in minutes: they matched with gcc 13 and 14 on x86-64 and aarch64 Linux, and with Apple clang on an Apple M1.

ival, crsum, crnn and lowp are unchanged; ival's soname is still libival.so.0.

crnumerics 0.4.1

The Python package is now crdistr both on PyPI and as a module: pip install crdistr, then import crdistr (it was import crdist). "crdist" on PyPI is an unrelated project. The C library, crdist, is unchanged, as is everything else in this release but the version strings.

crnumerics 0.4.0

crdist grows from the inverse of the normal distribution into the normal distribution: four more functions, correctly rounded, the 16-bit formats, and samplers whose draws are the same bits on every machine.

New functions, correctly rounded to nearest: ndtr (Phi, R's pnorm), log_ndtr (log Phi, also where Phi underflows: log_ndtr(-40) is -804.6084), normpdf and logit. Every one of the 2^32 binary32 inputs of each was checked against MPFR, with none differing.

All seven functions now also take binary16 and bfloat16 (crdist_ndtr_f16, crdist_ndtr_bf16 and so on, on bit patterns), each correctly rounded in its own format from the same binary64 values; every 16-bit input was checked against MPFR. Rounding the binary32 result again to 16 bits would have been wrong on up to 150 inputs of a function.

Samplers: uniform, normal, exponential, gamma, beta, chi-squared and Student t, from a Philox4x32-10 stream (seed, substream, index). A variate has no exact value to round, so crdist.h specifies each one as a fixed sequence of correctly rounded operations, and every machine gives the same bits. The check recomputes that specification with MPFR, tests the distributions, and records fingerprints of fixed streams, which another implementation of the specification can check itself against.

make install now installs crdist too: libcrdist.so.0, libcrdist.a, crdist.h and crdist.pc for pkg-config (Linux).

From Python, the crdist module (pip install crdistr) gives the functions over NumPy arrays, the format following the dtype, and crdist.Stream with the samplers.

ival, crsum, crnn and lowp are unchanged; ival's soname is still libival.so.0.

crnumerics 0.3.0

New: crdist (dist/), the inverse error functions and the standard normal quantile in binary32, correctly rounded to nearest on every input: crdist_erfinvf, crdist_erfcinvf and crdist_probitf (probit is R's qnorm and SciPy's ndtri). We know of no other correctly rounded implementation: CORE-MATH and RLIBM have no inverse of erf. With a counter-based random number generator, inverse-CDF sampling through crdist gives the same normal draws on every machine.

Each is computed in binary64 from polynomials and CORE-MATH's correctly rounded log, within 2^-50 of the exact value, and rounded once to binary32 when that cannot go two ways. The 118 inputs where it can come from a table, made with MPFR from all 2^32 inputs of each function. On x86-64 and aarch64, every input of each function was checked against MPFR, with none differing, and the 2^-50 bound was checked on every input as well. A call takes 36 to 54 ns an element on a laptop's Ryzen 5. Every call runs in the default floating-point environment and restores the caller's, so neither the rounding mode nor flush-to-zero can change a result.

make builds build/libcrdist.a and libcrdist.so (the header is dist/crdist.h); like crsum and crnn, crdist is not part of make install yet. make check runs a sampled check; make crdist-check-all checks every input (2.5 hours on 64 threads of an x86-64 server).

ival, crsum, crnn and lowp are unchanged: ival's results and exported functions are 0.2.2's, and its soname is still libival.so.0.

crnumerics 0.2.2

ival now carries CORE-MATH at ea7a421 (it carried 040ee48): upstream changed ten of the functions it uses (cosh, cospi, exp, exp10, exp2, log10, log2, sinh, sinpi and tanpi). Old against new, in all four rounding modes, the results and exception flags are the same on 686,469,600 calls on x86-64 and aarch64, every worst case CORE-MATH lists for those functions included, so ival's intervals are the same as 0.2.1's.

The self-test's bfloat16 references are right when built by GCC 13 on AArch64 with the fp16 extension (-march=native on most current Arm cores): GCC 13 loads __bf16 constants as binary16 (GCC bug 111867, fixed in GCC 14), and 17 of the 42 reference functions were wrong, which the self-test reported. Those files are now compiled with fp16 off under GCC before 14. ival, crsum, crnn and lowp do not use the __bf16 type.

New checks: make ival-julia-ccall and make ival-julia-threads call ival from Julia, from one thread and from several at once, and compare every bound with IntervalArithmetic.jl's correctly rounded mode; Julia stays in round-to-nearest around every call.

Compatibility: the soname is still libival.so.0 and the exported functions are the same as 0.2.1's.