README

DoubleDouble

English | 日本語

Extended-precision floating-point arithmetic for Delphi, implemented as unevaluated sums of IEEE 754 Double values: a two-component double-double record (TDDouble, ~106-bit significand) and a four-component quad-double record (TQDouble), both with full operator overloading. A FireMonkey (FMX) test-bench application accompanies the two units.

利用ライブラリ

• LUX ：Foundational mathematics library for the LUXOPHIA projects.

1. Overview

• TDDouble — a double-double number $x = x_{hi} + x_{lo}$ built on error-free transformations: Knuth's TwoSum, Veltkamp splitting, and Dekker's TwoProd. Provides + - * /, unary + -, all six comparison operators, implicit cast from Double, explicit cast to Double, ToString / Parse, and elementary functions (TDD_Abs, TDD_Sqr, TDD_Sqrt, TDD_Pow, TDD_Exp, TDD_Ln, TDD_Sin, TDD_Cos).
• TQDouble — a quad-double number stored as four components hi + me1 + me2 + lo. Sums and products are accumulated as floating-point expansions and compressed back to four components by a sort-then-accumulate renormalization (Normalize4) in the style of Shewchuk's expansion arithmetic. Overloaded elementary functions (Abso, Pow2, Roo2, Power, Exp, Ln, Sin, Cos) are provided.
• Division in both types is performed by Newton–Raphson iteration on the reciprocal, seeded with a double-precision approximation.
• The FMX form (Main.pas / Main.fmx) is a minimal harness with tabs (D1 / D2 / D4) and image panes for visual experiments with the types.

2. Mathematical Background

2.1 Representation

A double-double value represents a real number as the unevaluated sum of two doubles,

x = x_{hi} + x_{lo}, \qquad |x_{lo}| \le \tfrac{1}{2}\,\mathrm{ulp}(x_{hi}),

which yields an effective significand of about $2 \times 53 = 106$ bits (roughly 31 decimal digits). A quad-double extends this to four components,

x = x_{hi} + x_{me1} + x_{me2} + x_{lo},

targeting an effective significand of about $212$ bits.

2.2 Error-free transformations

TwoSum (Knuth [2]). For any two doubles $a, b$, the floating-point sum and its exact rounding error are recovered branch-free in six operations:

s = a \oplus b,\quad v = s \ominus a,\quad e = \bigl(a \ominus (s \ominus v)\bigr) \oplus (b \ominus v), \qquad \text{(2.1)}

so that $a + b = s + e$ exactly ($\oplus, \ominus, \otimes$ denote correctly rounded floating-point operations). TDDouble.TwoSum implements (2.1) verbatim.

Fast2Sum (Dekker [1]). When $|a| \ge |b|$, the cheaper form

s = a \oplus b,\quad e = b \ominus (s \ominus a) \qquad \text{(2.2)}

is also exact. TQDouble.Normalize4 relies on (2.2), which is why it first sorts its inputs by descending magnitude (see §2.4).

Veltkamp splitting (Dekker [1]). With the constant $C = 2^{27} + 1 = 134217729$, a double $a$ is split into two non-overlapping halves of at most 26–27 significant bits:

t = C \otimes a,\quad a_h = t \ominus (t \ominus a),\quad a_l = a \ominus a_h, \qquad \text{(2.3)}

with $a = a_h + a_l$ exactly. This is TDDouble.Split.

TwoProd (Dekker [1]). Using (2.3) on both factors, the exact product error is recovered without FMA:

p = a \otimes b,\quad e = \bigl((a_h b_h - p) + a_h b_l + a_l b_h\bigr) + a_l b_l, \qquad \text{(2.4)}

so that $a \cdot b = p + e$ exactly. This is TDDouble.TwoProd; no FMA instruction is used anywhere in the code.

2.3 Double-double operations (TDDouble)

• Addition (AddDD): the high parts and low parts of the operands are combined with (2.1), the intermediate error terms are re-absorbed by a cascade of further TwoSum calls (six in total), and the result is returned as two components. Subtraction is $x + (-y)$.
• Multiplication (MulDD): the leading term is computed with (2.4) as $p_1 + p_2 = x_{hi} \cdot y_{hi}$; the cross terms $x_{hi} y_{lo} + x_{lo} y_{hi}$ are combined by a TwoSum-based renormalization (Normalize2), then $p_2$, the cross terms, and $x_{lo} y_{lo}$ are accumulated with the double-double addition above.
• Division (DivDD): the reciprocal $z \approx 1/y$ is seeded with the double-precision quotient and refined by two Newton–Raphson steps

z_{k+1} = z_k + z_k\,(1 - y\,z_k), \qquad \text{(2.5)}

each step performed entirely in double-double arithmetic (quadratic convergence), then $x / y = x \cdot z$. Division by zero raises an exception.

• Elementary functions: TDD_Sqrt refines the double-precision seed with two Newton steps $r \leftarrow r - (r^2 - x)/(2r)$; TDD_Pow uses binary (square-and-multiply) exponentiation for integer exponents; TDD_Exp, TDD_Sin, TDD_Cos use truncated Taylor series (6 correction terms); TDD_Ln inverts TDD_Exp by Newton iteration $y \leftarrow y + (x - e^y)/e^y$.

2.4 Quad-double operations (TQDouble)

• Addition (AddQD): the eight components of both operands are concatenated and compressed by Normalize4.
• Renormalization (Normalize4): the input terms are sorted by descending absolute value, then accumulated left-to-right with (2.2) in the manner of Shewchuk's GrowExpansion [4], keeping at most four components hi, me1, me2, lo. The magnitude sort establishes the $|a| \ge |b|$ precondition of Fast2Sum.
• Multiplication (MulQD): all $4 \times 4 = 16$ pairwise component products $x_i , y_j$ are formed as ordinary double products and compressed by Normalize4. (The individual products are rounded — no TwoProd is applied at this stage — so the result is an approximate quad-double, accurate in its leading components.)
• Division (DivQD): the same Newton–Raphson reciprocal iteration (2.5) as in the double-double case, with two steps in quad-double arithmetic.
• Elementary functions: Roo2 (square root) by Newton iteration, Power by binary exponentiation, Exp / Sin / Cos by truncated Taylor series, Ln by Newton inversion of Exp — structurally identical to their TDDouble counterparts.

3. Architecture

Program structure (unit dependencies)

???DoubleDouble.dpr       ????????? (FMX application)
  ??????Main.pas / Main.fmx
     ??????TForm1 :TabControl (D1 / D2 / D4) + Images
        ??????uses LUX.D1.DDouble.pas
        ??????uses LUX.D1.QDouble.pas

Unit contents

???LUX.D1.DDouble.pas
  ??????TDDouble = record
  ???  ??????hi, lo :Double
  ???  ??????EFT primitives
  ???  ???  ??????TwoSum     ????????? (2.1)
  ???  ???  ??????Split      ????????? (2.3)
  ???  ???  ??????TwoProd    ????????? (2.4)
  ???  ???  ??????Normalize2
  ???  ??????core ops
  ???  ???  ??????AddDD SubDD MulDD DivDD
  ???  ??????operators
  ???     ??????+ - * /     ????????? (unary + -)
  ???     ??????= <> < > <= >=
  ???     ??????Implicit(Double)
  ???     ??????Explicit ??? Double
  ???     ??????ToString / Parse
  ??????elementary functions
     ??????TDD_Abs TDD_Sqr TDD_Sqrt TDD_Pow
     ??????TDD_Exp TDD_Ln TDD_Sin TDD_Cos

???LUX.D1.QDouble.pas
  ??????TQDouble = record
  ???  ??????hi, me1, me2, lo :Double
  ???  ??????renormalization
  ???  ???  ??????Normalize4 ????????? (sort + Fast2Sum (2.2))
  ???  ??????core ops
  ???  ???  ??????AddQD SubQD MulQD DivQD
  ???  ??????operators
  ???     ??????+ - * /     ????????? (unary + -)
  ???     ??????Implicit(Double)
  ???     ??????Explicit ??? Double
  ???     ??????ToString
  ??????elementary functions (overloads)
     ??????Abso Pow2 Roo2 Power
     ??????Exp Ln Sin Cos

File layout:

???DoubleDouble/
  ??????DoubleDouble.dpr      ????????? FMX application entry point
  ??????DoubleDouble.dproj    ????????? RAD Studio project (FMX, ProjectVersion 20.2)
  ??????Main.pas / Main.fmx   ????????? test-bench form (tabs D1 / D2 / D4)
  ??????_LIBRARY/
     ??????LUX.D1.DDouble.pas ????????? TDDouble: double-double arithmetic
     ??????LUX.D1.QDouble.pas ????????? TQDouble: quad-double arithmetic
     ??????LUXOPHIA/
        ??????LUX/            ????????? git subtree: https://github.com/LUXOPHIA/LUX

4. Usage

uses LUX.D1.DDouble, LUX.D1.QDouble;

var
   A, B, C :TDDouble;
   Q       :TQDouble;
begin
     // Double → TDDouble is implicit
     A := 1.0;
     B := 3.0;

     C := A / B;                       // Newton–Raphson division
     C := C * B - A;                   // residual, ≈ 0 to ~106 bits

     ShowMessage( C.ToString );        // formatted via ToString

     C := TDD_Sqrt( TDDouble( 2.0 ) ); // √2 in double-double

     // Quad-double
     Q := 1.0;
     Q := Q / TQDouble( 7.0 );
     ShowMessage( Q.ToString );

     // Explicit cast back to Double (rounds)
     ShowMessage( FloatToStr( Double( C ) ) );
end;

5. Building

• Open DoubleDouble.dproj in RAD Studio 12 (Athens) or later (ProjectVersion 20.2, FireMonkey framework).
• Enabled target platforms in the project file: Win32, Win64, macOS 64-bit, macOS ARM 64-bit, iOS Device 64-bit, iOS Simulator ARM 64-bit, Android 32/64-bit.
• The two arithmetic units have no dependencies beyond System.SysUtils and System.Math, so they can be dropped into any Delphi project (VCL, FMX, or console) as-is.

6. References

1. T. J. Dekker, "A floating-point technique for extending the available precision", Numerische Mathematik 18, pp. 224–242, 1971.
2. D. E. Knuth, The Art of Computer Programming, Vol. 2: Seminumerical Algorithms, 3rd ed., Addison-Wesley, 1997. (TwoSum, §4.2.2)
3. Y. Hida, X. S. Li, D. H. Bailey, "Algorithms for Quad-Double Precision Floating Point Arithmetic", Proc. 15th IEEE Symposium on Computer Arithmetic (ARITH-15), 2001. QD library
4. J. R. Shewchuk, "Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates", Discrete & Computational Geometry 18, pp. 305–363, 1997.
5. Wikipedia: Double-double arithmetic.

💖 Embarcadero Delphi

Integrated Development Environment (IDE) for Creating Native Cross-Platform Apps.

Free Download: Delphi Community Edition