Crate multidouble 

Double-double and quad-double floating-point arithmetic, scalar and on 8-lane SIMD vectors.

A double-double is an unevaluated sum hi + lo of two f64 (about 32 significant digits); a quad-double is a sum of four (about 64 digits). The lane types run the same algorithms on 8 values at once, with bitwise the same results per lane as the scalar types, on every backend:

• v8::V8: 8 f64 in one AVX-512 register, two AVX2 registers, four SSE2 or four NEON registers, or a plain array (the cargo feature portable forces it), with fused multiply-add throughout.
• simd: the backends behind V8. simd::dispatch runs code written once over the backend on the best one the CPU has, detected at run time; native has the lane types on the backend the build’s target features select at compile time.
• dd::DD scalar and dv::DV 8-lane double-double; dv::DotAcc a sum-of-products accumulator.
• qd::QD scalar and qv::QV 8-lane quad-double (the algorithms of Hida, Li and Bailey’s QD library, with two repairs found by this crate’s tests), exact decimal parsing and printing; qd::QdAcc and qv::QDotAcc accumulators.
• exact: exact dyadic and decimal arithmetic, the reference the tests compare against.
• dual: forward-mode dual numbers over these lanes (dual::Du double-double value and tangent, dual::Dtn a double-double value with N f64 tangents, and f64 variants), behind the dual::Lane trait, so one generic code gives values and exact derivatives; dual::Counted counts the operations of such code.
• field: 8 B values of a lane type as B lane vectors (field::Fl), with Cauchy products, Horner evaluation and a constant double-double matrix applied to them (field::DMat): the arithmetic of power-series methods.
• linalg: dense LU with partial pivoting in double-double (and blocked in f64 for mixed precision), solves, transpose solves, solves refined to quad-double and Gram matrices; linalg::lu_rect_f64, the f64 LU of the first columns of a rectangular block (the step of a block LU), and linalg::StaircaseLu, the LU of a staircase system with one global parameter (the Newton system of multiple shooting) built from those steps.

All operations are IEEE round-to-nearest with no reassociation, which the error-free transformations need.

The cargo feature f128 adds conversions between f128 and dd::DD and qd::QD. It needs a nightly compiler, as f128 is unstable.

use multidouble::dd::DD;
use multidouble::qd::QD;

let third = DD::from_f64(1.0) / DD::from_f64(3.0);
assert!(((third * DD::from_f64(3.0)) - DD::from_f64(1.0)).to_f64().abs() < 1e-31);

let x: QD = "1.41421356237309504880168872420969807856967187537694807317667973799".parse().unwrap();
assert!((x * x - QD::from_f64(2.0)).to_f64().abs() < 1e-62);

Modules

dd

Scalar double-double (hi + lo, |lo| <= ulp(hi)/2, ~32 significant digits). Used for constants and set-up; the hot loops use the 8-lane dv::DV with the same algorithms.

dual

Lane types for code that is written once and run in several arithmetics, with or without derivatives: the trait Lane is 8 values of a computation (8 grid points, 8 independent problems, …), and its implementations are

dv

8 lanes of double-double: (hi, lo) as two V8. The algorithms are the QD library’s accurate (“IEEE”) add, the FMA two-product multiply and the three-quotient divide, so every operation is accurate to a few units of 2^-104 relative.

exact

Exact arithmetic on dyadic rationals ±m·2^e (m an unsigned big integer) and on m / 10^k, with no dependency.

field

Fields: N = 8 B values of a lane type as B lane vectors (Fl<T, B>; a function on a grid of 8 B points, or 8 B independent problems), with the power-series arithmetic that series methods for ODEs need on them: Cauchy products of coefficient sequences (cs, cs_range, cs_w), Horner evaluation (horner), and a constant double-double matrix applied to a field (DMat, for instance a spectral derivative). Every operation is the lane type’s, so code written over Fl<T, B> runs with values, dual numbers or quad-double lanes alike (dual::Lane).

linalg

Dense double-double LU with partial pivoting (LAPACK getf2 order), column-major with the leading dimension padded to a multiple of 8 so the trailing update runs on 8-row lane vectors, columns updated in parallel (rayon) for large matrices; the same factorization blocked in plain f64 (LAPACK’s right-looking getrf) for mixed-precision use, solves and transpose solves in double-double with either set of factors, solves refined to quad-double (DDLu::solve_refine_qd), and Gram matrices.

native

The lane types on the compile-time backend simd::Native, for code that does not dispatch.

qd

Scalar quad-double: an unevaluated sum x0 + x1 + x2 + x3 of four f64, renormalized so that the parts do not overlap (|x_{i+1}| <= ulp(x_i)/2), about 64 significant digits (2^-209 relative, the QD library’s epsilon). The 8-lane qv::QV runs the same algorithms with bitwise the same results.

qv

8 lanes of quad-double: four V8 registers (structure of arrays, part k of the 8 lanes in c[k]), with the algorithms of the scalar qd::QD and bitwise the same results lane by lane, on every v8 backend (AVX-512, 2 x AVX2, portable).

simd

Runtime choice of the 8-lane backend.

v8

8-lane f64 vector V8<S> on a backend S (crate::simd), chosen at run time by crate::simd::dispatch or at compile time as crate::simd::Native: