Crate zenith_float

Expand description

zenith-float implements arbitrary-precision software floating-point numbers (ExactNum) and software IEEE-754 binary32/binary64 (Ieee32 / Ieee64). All arithmetic uses integer limbs. The library does not use hardware floating-point for calculations.

Repository guides: doc/GETTING_STARTED.md (short path) and doc/HELP.md (longer tutorial).

§Introduction

Numbers

The number is defined by the data type ExactNum. Each finite number consists of an array of words representing the mantissa, exponent, and sign. ExactNum can also be Inf (positive infinity), -Inf (negative infinity) or NaN (not-a-number).

ExactNum creation operations take bit precision as an argument. Precision is always rounded up to the nearest word. For example, if you specify a precision of 1 bit, then it will be converted to 64 bits when one word has a size of 64 bits. If you specify a precision of 65 bits, the resulting precision will be 128 bits (2 words), and so on.

Most operations take the rounding mode as an argument. The operation will typically internally result in a number with more precision than necessary. Before the result is returned to the user, the result is rounded according to the rounding mode and reduced to the expected precision.

The result of an operation is marked as inexact if some of the bits were rounded when producing the result, or if any of the operation’s arguments were marked as inexact. The information about exactness is used to achieve correct rounding.

ExactNum can be parsed from a string and formatted into a string using binary, octal, decimal, or hexadecimal representation.

Numbers can be subnormal. Usually any number is normalized: the most significant bit of the mantissa is set to 1. If the result of the operation has the smallest possible exponent, then normalization cannot be performed, and some significant bits of the mantissa may become 0. This allows for a more gradual transition to zero.

Error handling

In case of an error, such as memory allocation error, ExactNum takes the value NaN. ExactNum::err() can be used to get the associated error in this situation.

Constants

Constants such as pi or the Euler number have arbitrary precision and are evaluated lazily and then cached in the constants cache. Some functions expect constants cache as parameter.

Rounding

ExactNum methods that take a rounding mode other than RoundingMode::None round to the requested precision. RoundingMode::None skips that step and may keep extra bits. expr! raises working precision to compensate cancellation; it does not itself perform correct rounding.

§Examples

The example below computes Pi with precision 1024, rounding to even, using expr!.

use zenith_float::Consts;
use zenith_float::RoundingMode;
use zenith_float::ctx::Context;
use zenith_float::expr;

// Create a context with precision 1024, rounding to the nearest even,
// and exponent range from -100000 to 100000.
let mut ctx = Context::new(1024, RoundingMode::ToEven,
    Consts::new().expect("Constants cache initialized"),
    -100000, 100000);

// Compute pi: pi = 6*arctan(1/sqrt(3))
let pi = expr!(6 * atan(1 / sqrt(3)), &mut ctx);

// Use library's constant value for verifying the result.
let pi_lib = ctx.const_pi();

// Compare computed constant with library's constant
assert_eq!(pi.cmp(&pi_lib), Some(0));

// Print using decimal radix.
#[cfg(feature="std")]
println!("{}", pi);

// output: 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458699748e+0

The example below computes value of Pi with precision 1024 rounded to the nearest even number using ExactNum directly. We will take care of the error in this case.

use zenith_float::ExactNum;
use zenith_float::Consts;
use zenith_float::RoundingMode;

// Precision with some space for error.
let p = 1024 + 8;

// The results of computations will not be rounded.
// That will be more performant, even though it may give an incorrectly rounded result.
let rm = RoundingMode::None;

// Initialize mathematical constants cache
let mut cc = Consts::new().expect("An error occured when initializing constants");

// Compute pi: pi = 6*arctan(1/sqrt(3))
let six = ExactNum::from_word(6, 1);
let three = ExactNum::from_word(3, p);

let n = three.sqrt(p, rm);
let n = n.reciprocal(p, rm);
let n = n.atan(p, rm, &mut cc);
let mut pi = six.mul(&n, p, rm);

// Reduce precision to 1024 and round to the nearest even number.
pi.set_precision(1024, RoundingMode::ToEven).expect("Precision updated");

// Use library's constant for verifying the result
let pi_lib = cc.pi(1024, RoundingMode::ToEven);

// Compare computed constant with library's constant
assert_eq!(pi.cmp(&pi_lib), Some(0));

// Print using decimal radix.
#[cfg(feature="std")]
println!("{}", pi);

// output: 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458699748e+0

§Performance recommendations

When small error is acceptable because of rounding it is recommended to do all computations with RoundingMode::None, and use ExactNum::set_precision or ExactNum::round with a specific rounding mode just once for the final result.

§no_std

The library can work without the standard library provided there is a memory allocator. The standard library dependency is activated by the feature std. The feature std is active by default and must be excluded when specifying dependency, e.g.:

[dependencies]
zenith-float = { version = "1.0.1", default-features = false }

Modules§

ctx: Context is used in expressions returning ExactNum.

Macros§

cexpr: Computes an expression with the specified precision and rounding mode.
exact: Computes an expression with the specified precision and rounding mode.
expr: Computes an expression with the specified precision and rounding mode.
fbig: Computes an expression with the specified precision and rounding mode.

Structs§

Ball: Enclosure mid ± rad used to certify that a rounded midpoint is unique.
CachedFBig: A float stored at extra working precision so later requests at lower (or equal) precision reuse the cache instead of recomputing.
ComplexBall: Disk enclosure in (\mathbb{C}): center mid, radius rad.
ConstCacheInfo: Snapshot of how many mantissa bits of each constant are currently cached.
Consts: Constants cache contains arbitrary-precision mathematical constants.
ExactComplex: Complex value re + i·im with software-limb real and imaginary parts.
ExactInt: Arbitrary-precision signed integer. Distinct from ExactNum, which is floating-point.
ExactNum: A floating point number of arbitrary precision.
ExactNumArray: Batch of ExactNum values at a shared precision p, row-major.
ExactNumPoly: Dense univariate polynomial c₀ + c₁ x + ⋯ + cₙ xⁿ.
ExactRational: Exact rational num/den. Distinct from ExactNum, which is floating-point.
Ieee32: Software IEEE-754 binary32 (24-bit significand, 8-bit exponent). Stored as u32 bits.
Ieee64: Software IEEE-754 binary64 (53-bit significand, 11-bit exponent). Stored as u64 bits.
Ieee32Array: Contiguous binary32 lanes (u32 bits), row-major.
Ieee64Array: Contiguous binary64 lanes (u64 bits), row-major.
InlineBinaryBuffer: Fixed 16-byte stack buffer for an inlined ExactNum.
Radix: Radix for parse/format (bases 2 through 36).
RadixFloat: A floating-point value with an associated parse/format radix.
SharedConsts: Thread-safe Consts for batch evaluation across threads (std only). Callers still take &mut Consts inside SharedConsts::with; the mutex serializes cache fills.

Enums§

Error: Possible errors.
RoundingMode: Rounding modes.
Sign: Sign.

Constants§

BINARY_FLAG_ARRAY: ExactNumArray record.
BINARY_FLAG_HEAP_NEG: Heap finite, negative.
BINARY_FLAG_HEAP_POS: Heap finite, positive.
BINARY_FLAG_INF_NEG: INF_NEG.
BINARY_FLAG_INF_POS: INF_POS.
BINARY_FLAG_NAN_BARE: NaN with no associated Error.
BINARY_FLAG_NAN_DIV0: NaN: Error::DivisionByZero.
BINARY_FORMAT_VERSION: First-byte version stored in every record.
BINARY_HEADER_LEN: Header bytes of a heap or array record (same size as the inline record).
BINARY_INLINE_LEN: Bytes in InlineBinaryBuffer / ExactNum::to_inline_bytes.
BINARY_INLINE_MANT_BITS: Mantissa bits that fit in the inline record (2 × 32).
BINARY_INLINE_U32_WORDS: u32 limbs in the inline mantissa field.
BINARY_MAX_ELEMS: Maximum array elements accepted by ExactNumArray::from_bytes.
BINARY_MAX_U32: Maximum u32 limbs accepted by ExactNum::from_bytes.
CHEBYSHEV_MAX_DEGREE: Maximum number of Chebyshev coefficients chebyshev_coeffs will compute.
CSV_MAX_COLS: Maximum columns accepted in one CSV row.
CSV_MAX_ROWS: Maximum data rows accepted by Ieee64Array::from_csv_str.
DSP_MAX_POINTS: Maximum real length for dct / idct / dst / idst and the window generators. Transforms use a 2N-point FFT, so this is half of FFT_MAX_POINTS.
EXPONENT_BIT_SIZE: The size of exponent type in bits.
EXPONENT_MAX: Maximum exponent value.
EXPONENT_MIN: Minimum exponent value.
IEEE_SIMD_LANE_WIDTH: Number of u32 lanes in one integer SIMD vector (128-bit register). Binary64 uses IEEE_SIMD_LANE_WIDTH / 2 u64 lanes. Scalar fallback uses the same width so wrappers can size buffers without cfg.
INF_NEG: Negative infinity.
INF_POS: Positive infinity.
INLINE_WORDS: Number of words kept on the stack before allocating.
JACOBI_AGM_MAX: AGM / descending Landen steps for am, sn, cn, dn.
MAX_PREC_RETRY: Maximum extra correct-rounding retries.
NAN: Not a number.
ODE_MAX_STEPS: Maximum accepted steps for any solver in this module.
ODE_MIN_STEP: Default minimum step exponent: h_min = 2^{ODE_MIN_STEP}.
ORTHOPOLY_N_MAX: Maximum degree for every family in this module. Larger n is NaN.
POLLARD_RHO_ITER_MAX: Maximum f evaluations in one pollard_rho run before None.
POLY_COMPANION_CLOSED_DEG: Highest degree whose companion eigenvalues have a closed form on this type.
PROPTEST_CASES: Cases per proptest property under cargo test (no mpfr-tests required).
QUADRATURE_MAX_NODES: Maximum number of Gauss nodes. Larger n is None.
ROOT_DEFAULT_TOL: Default absolute tolerance exponent: tol = 2^{ROOT_DEFAULT_TOL}.
ROOT_MAX_ITER: Maximum iterations for every method in this module.
TANH_SINH_LEVELS_MAX: Maximum tanh–sinh step halvings after the precision-based h.
WORD_BASE: Base of words.
WORD_BIT_SIZE: Size of a word in bits.
WORD_MAX: Maximum value of a word.
WORD_SIGNIFICANT_BIT: Word with the most significant bit set.

Traits§

FromExt: A trait for conversion with additional arguments.

Functions§

bisect: Bisection on [a, b]. None if f(a) and f(b) do not have opposite signs.
blackman_window: Symmetric Blackman: 0.42 − 0.5 cos(θ) + 0.08 cos(2θ), θ = 2πk/(n−1).
brent: Brent’s method (bisection + secant + inverse quadratic) on [a, b].
chebyshev_coeffs: Interpolation coefficients of f on [a, b] at the n Chebyshev–Gauss nodes.
chebyshev_error_bound: ℓ¹ tail Σ_{k≥1} |c_k| at precision p.
chebyshev_eval: Evaluate the Chebyshev expansion of coeffs at x ∈ [a, b].
clenshaw: Clenshaw recurrence for Σ_{k=0}^{n-1} c_k T_k(x) at precision p.
constant_time_eq: Timing-safe equality: always scans both slices; length mismatch is false.
dct: Discrete cosine transform, type II, via a 2N-point FFT.
dst: Discrete sine transform, type II, via a 2N-point FFT.
euler: Explicit Euler for y' = f(t, y) on [t0, t1] with n_steps equal steps.
fft_real: Unnormalized DFT of a real row or column. Output is (2, n) (row 0 real, row 1 imaginary). Same power-of-two rule as ExactNumArray::fft.
gauss_hermite: Gauss–Hermite quadrature: ∫_{-∞}^{∞} f(x) e^{-x²} dx with n nodes.
gauss_laguerre: Gauss–Laguerre quadrature: ∫₀^∞ f(x) e^{-x} dx with n nodes.
gauss_legendre: Gauss–Legendre quadrature of f on [a, b] with n nodes.
hamming_window: Symmetric Hamming: 0.54 − 0.46 cos(2πk/(n−1)). Endpoints are 0.08.
hann_window: Symmetric Hann: ½(1 − cos(2πk/(n−1))). n = 1 is [1].
hmac_sha256: HMAC-SHA-256 (key, data).
idct: Inverse of dct: type-III DCT scaled by 2/N so idct(dct(x)) = x.
idst: Inverse of dst: type-III DST scaled by 2/N so idst(dst(x)) = x.
ifft_real: Inverse of fft_real: ExactNumArray::ifft then the real row.
illinois: Illinois (modified regula falsi) on [a, b].
kaiser_window: Kaiser–Bessel: I_0(β √(1−t_k²)) / I_0(β) with t_k = (k−(n−1)/2)/((n−1)/2).
miller_rabin: Miller–Rabin on |n| with the given bases. n < 2 is false.
mod_inv: Modular inverse of a modulo modulus. None if not invertible.
mod_pow: base^exp mod modulus. exp < 0 or a zero modulus is None. |modulus| = 1 is 0.
newton: Newton–Raphson from x0 with exact derivative df.
ode_min_step: Minimum step 2^{ODE_MIN_STEP} at precision p.
pollard_rho: Brent Pollard ρ. A proper factor of |n|, or None if prime / cap hit.
rectangular_window: Rectangular window: n ones.
rk4: Classical RK4 for y' = f(t, y) on [t0, t1] with n_steps equal steps.
rk45_adaptive: Dormand–Prince RK5(4) with absolute / relative step control.
root_default_tol: Suggested absolute tolerance 2^{ROOT_DEFAULT_TOL} at precision p.
sha256: SHA-256 of data.
sha512: SHA-512 of data.
tanh_sinh: Tanh–sinh (double-exponential) quadrature of f on [a, b].
ziv_round: Evaluate compute at increasing working precision until the result rounds uniquely to p bits (try_set_precision). Same retry budget as the transcendental kernel (crate::MAX_PREC_RETRY).
ziv_round_vec: Same Ziv loop as ziv_round, with every input lifted to p_wrk before compute.

Type Aliases§

ConstCache: Alias for Consts: a progressive cache of π, e, ln 2, ln 10, √2, φ, and γ.
Exponent: An exponent.
Word: A word.