1.5e7 # a floating-point value 1.5 × 10⁷

x = 1/49 # division of two integers produces a floating-point value

x * 49

1 - x * 49

2.0^-52, eps(), eps(1.0), eps(Float64) # these are all the same thing

# 1.5 is exactly represented in binary floating point:
big(1.5), 1.5 == 3//2

# 0.1 is *not* exactly represented
big(0.1), 0.1 == 1//10

function my_cumsum(x)
    y = similar(x) # allocate an array of the same type and size as x
    y[1] = x[1]
    for i = 2:length(x)
        y[i] = y[i-1] + x[i]
    end
    return y
end

eps(Float32), eps(Float64)

x = rand(Float32, 10^7) # 10^7 single-precision values uniform in [0,1)
@time y = my_cumsum(x)
yexact = my_cumsum(Float64.(x)) # same thing in double precision
err = abs.(y .- yexact) ./ abs.(yexact) # relative error in y

using PyPlot
n = 1:10:length(err) # downsample by 10 for plotting speed
loglog(n, err[n])
ylabel("relative error")
xlabel("# summands")
# plot a √n line for comparison
loglog([1,length(err)], sqrt.([1,length(err)]) * 1e-7, "k--")
text(1e3,1e-5, L"\sim \sqrt{n}")
title("naive cumsum implementation")

VERSION < v"0.3.6" && warn("cumsum is less accurate for Julia < 0.3.6")
@time y2 = cumsum(x)
err2 = abs.(y2 .- yexact) ./ abs.(yexact)
loglog(n, err2[n])
ylabel("relative error")
xlabel("# summands")
title("built-in cumsum function")
loglog(n, sqrt.(log.(n)) * 1e-7, "k--")

rounding(Float32)

] add SetRounding

using SetRounding

errup = setrounding(Float32, RoundUp) do
    # error in single-precision (Float32) sum, rounding temporarily set to RoundUp
    abs.(my_cumsum(x) .- yexact) ./ abs.(yexact) # relative error in y
end

loglog(n, errup[n])
ylabel("relative error")
xlabel("# summands")
# plot an O(n) line for comparison
loglog([1,length(errup)], [1,length(errup)] * 1e-7, "k--")
text(1e3,4e-4, L"\sim n")
title("naive cumsum implementation, rounded up")

1e300 # okay: 10³⁰⁰ is in the representable range

(1e300)^2 # overflows

floatmax(Float64), floatmax(Float32)

1e-300 # okay

(1e-300)^2 # underflows to +0

floatmin(Float64), floatmin(Float32)

-1e-300 * 1e-300 # underflows to -0

+0.0 == -0.0

1 / +0.0, 1 / -0.0

signbit(+0.0), signbit(-0.0)

1 / (1 / -Inf)

0 * Inf, Inf / Inf, 0 / 0

sqrt(-1.0)

sqrt(-1.0 + 0im)