@doc """
  Splits a floating-point number `a` into two parts `a = h + l` such that the
  halves `h` and `l` are non-overlapping and `|x| ≥ |y|`. That is, both halves
  have at most `s-1` nonzero bits, where `s = ceil(p/2)` and `p` is the
  precision (`p = 24` for `Float32` and `p = 53` for `Float64`).

  References:

    * T.J. Dekker, A floating-point technique for extending the available
      precision, Numerische Mathematik 18, pp. 224–242, 1971.
  """ ->
# TODO: Extend the definition to floating-point arrays
for T in (Float32, Float64)
  @eval begin
    function fsplit(a::$T)
      # Float32: factor = 2^12 + 1 = 4097
      # Float64: factor = 2^27 + 1 = 134217729
      f = $(convert(T, 2^ceil(precision(T)/2)+1));
      c = a * f
      h = c - (c-a)
      l = a - h
      h, l
    end
  end
end


@doc """
  Computes the error-free multiplication `a * b = x + y`, where `x = fl(a * b)`,
  for floating-point numbers `a` and `b`. This algorithm uses all cross products
  between the halves generated by `fsplit` from `a` and `b`. It requires 17
  flops.

  References:

    * T.J. Dekker, A floating-point technique for extending the available
      precision, Numerische Mathematik 18, pp. 224–242, 1971.
  """ ->
# TODO: Extend the definition to floating-point arrays
function err_mul{T<:IEEE754}(a::T, b::T)
  x = a * b
  ah, al = fsplit(a)
  bh, bl = fsplit(b)
  y = al*bl - (((x - ah*bh) - al*bh) - ah*bl)
  x, y
end


@doc """
  Computes the error-free transformation `a * b = x + y`, where `x = fl(a * b)`,
  for floating-point numbers `a` and `b`. By using the fused multiply-add (FMA)
  operation, it only requires 2 flops instead of the 17 flops required by
  `err_mul` but the improvement in execution time typically a less than that.

  References:

    * T. Ogita, S.M. Rump, and S. Oishi, Accurate sum and dot product, SIAM
      Journal on Scientific Computing 26, pp. 1955–1988, 2005.
  """ ->
# TODO: Extend the definition to floating-point arrays
function err_fast_mul{T<:IEEE754}(a::T, b::T)
  x = a * b
  y = fma(a, b, -x)
  x, y
end
