b = 1
c = logspace(-20,-1,50) # 50 values from 10^-20 to 10^-1, equally spaced on a log scale
x = b - sqrt(b^2 - c)   # x computed in the default double precision (~15 digits)
set_bigfloat_precision(256) # 256 bits, about 80 decimal digits
xbig = b - sqrt(b^2 - big(c)) # x computed in "bigfloat" precision (~80 digits)
err = float64(abs(x - xbig) ./ abs(xbig)) # relative error, converted to double-precision to plot

using PyPlot
loglog(c, err, "bo-")
xlabel(L"c")
ylabel(L"relative error in $x$")
title(L"standard quadratic formula $x = b - \sqrt{b^2 - c}$")





# 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, logspace(1,7,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
loglog(N, err, "bo-")
title("simple div2sum")
xlabel("number of summands")
ylabel("relative error")
grid()

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