Documentation
¶
Overview ¶
Package dd is a minimal double-double ("float128-lite") arithmetic used where a special function must carry more than 53 bits through an intermediate step: the argument of an exponential, a long product, a cancellation. A value is Hi + Lo with |Lo| <= ulp(Hi)/2.
The error-free transformations are the classical ones: TwoSum (Knuth, TAOCP vol. 2, 4.2.2), Fast2Sum (Dekker 1971) and TwoProd through a fused multiply-add (Ogita, Rump and Oishi, "Accurate sum and dot product", SIAM J. Sci. Comput. 26(6), 2005). Algorithms for the compound operations follow Joldes, Muller and Popescu, "Tight and rigorous error bounds for basic building blocks of double-word arithmetic", ACM TOMS 44(2), 2017 (DWPlusDW: Algorithm 6, DWTimesDW: Algorithm 12, DWDivDW: Algorithm 18).
This package is a copy of go-numerics/special's internal/dd (BSD-3-Clause), which an internal package cannot share across modules.
Index ¶
- Constants
- func Exp(a DD) float64
- func FastTwoSum(a, b float64) (s, e float64)
- func TwoProd(a, b float64) (p, e float64)
- func TwoSum(a, b float64) (s, e float64)
- type DD
- func Add(a, b DD) DD
- func AddF(a DD, b float64) DD
- func Div(a, b DD) DD
- func DivF(a DD, b float64) DD
- func ExpMant(a DD) (m DD, k int)
- func F(x float64) DD
- func Ldexp(a DD, k int) DD
- func Log(x float64) DD
- func LogAccurate(x float64) DD
- func Mul(a, b DD) DD
- func MulF(a DD, b float64) DD
- func MulFF(a, b float64) DD
- func Neg(a DD) DD
- func Sqrt(a DD) DD
- func Sub(a, b DD) DD
Constants ¶
const ( Ln2Hi = 6.93147180559945286227e-01 Ln2Lo = 2.31904681384629955842e-17 )
ln2 split as a double-double: Ln2Hi + Ln2Lo = ln 2 to 106 bits.
Variables ¶
This section is empty.
Functions ¶
func Exp ¶
Exp returns exp(a.Hi + a.Lo) rounded to float64: the double-double mantissa from ExpMant is accurate to about 2^-76, so the result is correctly rounded except in cases closer to a rounding boundary than that.
func FastTwoSum ¶
FastTwoSum requires |a| >= |b| (or a == 0).
func TwoProd ¶
TwoProd returns p, e with p = fl(a*b) and p + e = a*b exactly (barring underflow).
The explicit float64 conversion is load-bearing: the Go specification allows a*b to be fused with a later addition (an FMA on arm64, ppc64le, s390x, riscv64, loong64), which would break every error-free transformation built on p. A conversion forces the rounding.
Types ¶
type DD ¶
type DD struct{ Hi, Lo float64 }
DD is an unevaluated sum Hi + Lo.
func Add ¶
Add returns a + b (DWPlusDW, relative error <= 3u^2).
It is written out (rather than through TwoSum/FastTwoSum) to stay within the inliner's budget: it is the most frequent operation in the package.
func ExpMant ¶
ExpMant returns m, k with exp(a) = m * 2^k, m in about [0.99, 2.02] as a double-double accurate to about 2^-76 relative, so that a caller can multiply further factors in before the single final rounding. k is not limited to the float64 exponent range: a scaled result may be combined with others before math.Ldexp decides overflow or underflow. For |a.Hi| > 1e7 it returns m = +Inf or m = 0 with k = 0 (NaN for NaN).
Method (Tang, "Table-driven implementation of the exponential function in IEEE floating-point arithmetic", ACM TOMS 15(2), 1989, carried to double-double): a = (64k + j) ln2/64 + r, |r| <= ln2/128, the reduction exact (ln2/64 split in three parts, the first with 32 significant bits so that n times it is exact for |n| < 2^21, the second through TwoProd); exp(r) = 1 + r + r²/2 + r³ P(r) with 1 + r + r²/2 in double-double and P (to r^9/9!, truncation below 2^-86) in float64; 2^(j/64) from a 300-bit mpmath table (zz_exptable.go).
func Log ¶
Log returns log(x) for finite x > 0 as a double-double accurate to about 2^-72 relative (measured 2^-72.9; the float64 tail limits it).
Method (table-driven, after Tang, "Table-driven implementation of the logarithm function in IEEE floating-point arithmetic", ACM TOMS 16(4), 1990): x = 2^k m with m in [1/sqrt2, sqrt2); c = j/128 is the grid point nearest m, and log m = log c + 2 atanh(s), s = (m-c)/(m+c). m-c is exact (Sterbenz), s is formed in double-double, |s| < 1/360, and 2 atanh(s) = 2s + s³ (2/3 + 2s²/5 + 2s⁴/7) with the second part (below 2^-16 of the first) in float64. log c comes from a 300-bit mpmath table (zz_logtable.go); k ln2 is added exactly.
Log(0) = -Inf, Log(+Inf) = +Inf, NaN for x < 0 or NaN (as math.Log).
func LogAccurate ¶
LogAccurate is Log carried further, to about 2^-104 relative: the atanh series keeps its terms through s⁴/5 in double-double and the quotient s comes from a full double-double division. It costs about five times Log and is for the few places where the result is a difference of logarithms that cancel to below 2^-30 of their size (Digamma near its zeros).
Source Files