Thread by @raymondh

Another building block for a #Python floating point ninja toolset:

def veltkamp_split(x):
'Exact split into two 26-bit precision components'
t = x * 134217729.0
hi = t - (t - x)
lo = x - hi
return hi, lo

csclub.uwaterloo.ca/~pbarfuss/dekk…

Input: one signed 53-bit precision float

Output: two signed 26-bit precision floats

Invariant: x == hi + lo

Constant: 134217729.0 == 2.0 ** 27 + 1.0

Example:

>>> hi, lo = veltkamp_split(pi)
>>> hi + lo == pi
True
>>> hi.hex()
'0x1.921fb58000000p+1'
>>> lo.hex()
'-0x1.dde9740000000p-26'

Note all the trailing zeros and the difference between the two exponents. Also both the lo and hi values are signed.

Use case: The 26-bit precision components can be multiplied losslessly (without rounding):

# Four exact components of e * pi:
>>> pi_hi*e_hi
8.539734226211163
>>> pi_hi*e_lo
7.207993525551209e-08
>>> pi_lo*e_hi
-7.561753118743836e-08
>>> pi_lo*e_lo
-6.382525038592985e-16

The payoff for Veltkamp-Dekker splitting and piecewise multiplication is that you can build quad precision arithmetic out of everyday double precision arithmetic.

The part that is pure magic: How do you split a 53 bit number into two 26 bit numbers? Where did the extra bit go?

Answer: The extra bit is stored in the sign bit of the "lo" component.

53 bits + 1 sign = 26 bits + 1 sign + 26 bits + 1 sign

All bits accounted for.🧐

Cool!😏