approaching better vectorization

@chriselrod @ffevotte I could use some help providing better vectorization over basic arithmetic ops for DoubleFloats. I have tried to follow the guidelines you have given without seeing much any change. Underlying addition and multiplication, there are three core routines. Being able to vectorize well over them would go a long way.

Multiplication is faster than addition of Double64s. Both are needed to accelerate matrix ops. Please ask questions.

```
# Double64s are structs

struct Double64
hi::Float64
lo::Float64
end
```

```
# the underlying support for arithmetic

"""
two_prod(a, b)

Computes `hi = fl(a*b)` and `lo = err(a*b)`.
"""
@inline function two_prod(a::T, b::T) where {T<:Float64}
hi = a * b
lo = fma(a, b, -s)
return hi, lo
end

"""
two_sum(a, b)

Computes `hi = fl(a+b)` and `lo = err(a+b)`.
"""
@inline function two_sum(a::T, b::T) where {T<:Float64}
hi = a + b
a1 = hi - b
b1 = hi - a1
lo = (a - a1) + (b - b1)
return hi, lo
end

# there is an alternative implementation of two_sum.
@inline function two_sum(a::T, b::T) where {T<:Float64}
hi = a + b
v = hi - a
lo = (a - (hi - v)) + (b - v)
return hi, lo
end

# this is a special case of two_sum when the relative magnitudes are known
# it is used heavily to gain some speed where allowable
"""
two_hilo_sum(a, b)

*unchecked* requirement `|a| ≥ |b|`
Computes `hi = fl(a+b)` and `lo = err(a+b)`.
"""
@inline function two_hilo_sum(a::T, b::T) where {T<:Float64}
hi = a + b
lo = b - (s - a)
return hi, lo
end

# this is the support for division (called `two_dvi` in the source)
"""
two_div(a, b)

Computes `hi = fl(a/b)` and `lo = err(a/b)`.
"""
@inline function two_div(a::T, b::T) where {T<:Float64}
hi = a / b
lo = fma(-hi, b, a)
lo /= b
return hi, lo
end
```

```
# The basic arithmetic `add` and `mul` for Double64s

# relative error < 3u² (called `add_dddd_dd` in the source)
# each Tuple gives the Double64 struct fields `(hi, lo)`
#
# >>>> It is easy to pass the structs instead if that is easier to vectorize
#
@inline function add(x::Tuple{T,T}, y::Tuple{T,T}) where T<:Float64
xhi, xlo = x
yhi, ylo = y
hi, lo = two_sum(xhi, yhi)
thi, tlo = two_sum(xlo, ylo)
c = lo + thi
hi, lo = two_hilo_sum(hi, c)
c = tlo + lo
hi, lo = two_hilo_sum(hi, c)
return hi, lo
end

# relative error <= 5u² (called `mul_dddd_dd` in the source)
@inline function mul(x::Tuple{T,T}, y::Tuple{T,T}) where T<:Float64
xhi, xlo = x
yhi, ylo = y
hi, lo = two_prod(xhi, yhi)
t = xlo * ylo
t = fma(xhi, ylo, t)
t = fma(xlo, yhi, t)
t = lo + t
hi, lo = two_hilo_sum(hi, t)
return hi, lo
end

# called `dvi_dddd_dd` in the source
@inline function divide(x::Tuple{T,T}, y::Tuple{T,T}) where T<:Float64
xhi, xlo = x
yhi, ylo = y
hi = xhi / yhi
uh, ul = two_prod(hi, yhi)
lo = ((((xhi - uh) - ul) + xlo) - hi*ylo)/yhi
hi,lo = two_hilo_sum(hi, lo)
return hi, lo
end
```