L4(x,y) = ...

L4(1e-100,0.0)

L4(1e+100,0.0)

cotdiff(x,y) = ...

cotdiff(1.0, 1e-20)

cotdiff(big(1.0), big(1e-20)) # compute in BigFloat precision



# Sum x[first:last].  This function works, but is a little slower than we would like.
function div2sum(x, first=1, last=length(x))
    n = last - first + 1;
    if n < 2
        s = zero(eltype(x))
        for i = first:last
            s += x[i]
        end
        return s
    else
        mid = div(first + last, 2) # find middle as (first+last)/2, rounding down
        return div2sum(x, first, mid) + div2sum(x, mid+1, last)
    end
end

# check its accuracy for a set logarithmically spaced n's.  Since div2sum is slow,
# we won't go to very large n or use too many points
N = round.(Int, 10 .^ range(1,7,length=50)) # 50 points from 10¹ to 10⁷
err = Float64[]
for n in N
    x = rand(Float32, n)
    xdouble = Float64.(x)
    push!(err, abs(div2sum(x) - sum(xdouble)) / abs(sum(xdouble)))
end

using PyPlot
loglog(N, err, "bo-")
title("simple div2sum")
xlabel("number of summands")
ylabel("relative error")
grid()

x = rand(Float32, 10^7)
@time div2sum(x)
@time div2sum(x)
@time div2sum(x)
@time sum(x)
@time sum(x)
@time sum(x)