repository·main·Indexed 41 days ago
https://github.com/juliamath/doublefloats.jlExtended precision floating-point and complex types for Julia, providing Double64, Double32, and Double16. It offers 85+ accurate bits and is designed to be more performant than Float128 or BigFloat for specific high-precision requirements. Includes support for complex types like ComplexDF64 and ComplexDF32, as well as the df macro to prevent precision loss when parsing literals.
Accuracy for Double64 is measured against a 512-bit BigFloat reference. Errors are reported in units of the Double64 working precision, where 1 ulp = 2⁻¹⁰⁴ (approximately 4.9e-32) relative.
A worst-case error of approximately 1 ulp indicates that the function is performing with the maximum accuracy allowed by the format.
A DoubleFloat{T} uses an unevaluated sum representation consisting of two floats, hi and lo:
value = hi + lo, with |lo| ≤ ulp(hi)/2
hi is the value rounded to type T.lo is the leading part of the rounding error.Double64, this provides a 106-bit significand (~32 decimal digits).T. Overflow and underflow occur at the same points as Float64.|value| < ~4e-292, the lo word becomes subnormal and precision tapers. Below ~1e-308, only Float64 precision remains.Double32 when you need approximately 14 significant decimal digits of precision. Like Double64, it is an immutable bitstype that offers much better performance than BigFloat for high-precision needs, though it is slower than native Float32.DoubleFloats offer high precision for elementary functions. For Double64 arguments within the range 0.0..2.0, you can expect:
abserr) of elementary functions to be 1e-30 or better.relerr) of elementary functions to be 1e-28 or better.When using reasonably sized values, successive DoubleFloat operations add no more than 10⋅𝘂² to the cumulative relative error (where 𝘂 is the relative rounding unit, typically 𝘂 = eps(x)/2). While relative error can accrue steadily over many operations (e.g., after 100,000 ops, relerr could approach 100,000 * 10⋅𝘂²), the algorithms are designed to be highly resilient.
DoubleFloats provide extended precision by using nonoverlapping pairs of standard floating-point types. They are magnitude ordered.
Double64 is a pair of Float64 with a 106-bit significand and 11-bit exponent.Double32 is a pair of Float32 with a 48-bit significand and 8-bit exponent.Double16 is a pair of Float16 with a 22-bit significand and 5-bit exponent.For complex numbers, the package provides:
ComplexDF64 (pair of Double64)ComplexDF32 (pair of Double32)ComplexDF16 (pair of Double16)The eigen function for DoubleFloat matrices follows LAPACK-route conventions:
(real, imag).Accuracy for Double32 is measured against a 512-bit BigFloat reference. Errors are reported in units of the Double32 working precision, where 1 ulp = 2⁻⁴⁶ (approximately 1.4e-14) relative.
A worst-case error of approximately 1 ulp indicates that the function is performing with the maximum accuracy allowed by the format.
Double64 is designed to limit the accrual of relative error during arithmetic operations. For reasonably sized values, a single operation should limit relative error to approximately $10 \cdot \mathbf{u}^2$, where $\mathbf{u}$ is the relative rounding unit (often eps(x)/2).
Double64 are considered reliable when mapped back to Float64 or Float32.DoubleFloats.jl provides high accuracy for transcendental functions (e.g., exp, log, sin, cos) with absolute and relative errors typically around $10^{-31}$.
When using isapprox, the package uses a default relative tolerance of rtol = eps(1.0)^(37/64).
The library uses specific acronyms to denote different representations and operations involving DoubleFloats. Understanding these prefixes helps in identifying the data format or operation being used:
op: Refers to operators such as add, sub, mul, and dvi.fp: Refers to standard floating point types (Float64, Float32, Float16).dd: Represents a DoubleFloat stored as a 2-tuple of (hi, lo) components.db: Represents a DoubleFloat stored as a struct DoubleFloat{::IEEEFloat}(hi, lo).acronym | meaning
--------|--------
op | operator (`add`, `sub`, `mul`, `dvi`)
fp | floating point (`Float64`, `Float32`, `Float16`)
dd | DoubleFloat as a 2-tuple (hi, lo)
db | DoubleFloat as a struct DoubleFloat{::IEEEFloat}(hi, lo)Arithmetic in DoubleFloats.jl is built upon exact building blocks to ensure high precision:
two_sum(a, b) / two_hilo_sum(a, b): Computes hi + lo == a + b exactly.two_prod(a, b): Computes hi + lo == a * b exactly using fma.Note on Accuracy: A zero lo word does not guarantee a value is exact. Algorithms must always include correction steps even if lo is zero to avoid losing precision (e.g., in inv operations).
DoubleFloats supports matrix functions with different requirements based on the algorithm used:
Dense Schur-based and scaling-and-squaring algorithms: These work on every square matrix, including defective (non-diagonalizable) ones. This applies to:
explogsqrtEigen-diagonalization: These require the argument to be diagonalizable. This applies to:
cbrtsinpi/cospiNote on Result Types: Real matrices return real results whenever the mathematical result is real. If a negative real eigenvalue makes the result genuinely complex, the return type will reflect that (e.g., a complex matrix).
using DoubleFloats
# Works even on defective matrices
exp(Double64[1 1; 0 1])