CheckKahanAngle.jl

module CheckKahanAngle

using PyPlot

function kahanAngle(a, b, c)
  b >= c >= 0 || c > b >= 0 || error("not a triangle")
  μ = ifelse(b >= c >= 0, c - (a - b), b - (a - c))
  2*atan(sqrt(((a-b)+c)*μ/((a+(b+c))*((a-c)+b))))
end

function f(v1, v2)
  a, b, c = norm(v1), norm(v2), norm(v1 - v2)
  kahanAngle(a, b, c)
end

function f2(v1, v2)
  v1, v2 = big.(v1), big.(v2)
  acos(dot(v1./norm(v1), v2./norm(v2)))
end

function fe{T<:AbstractFloat}(θ::T)
  v1, v2 = [1.0, 0.0], T.([cos(big(θ)), sin(big(θ))])
  T(abs(big(f(v1, v2)) / f2(big.(v1), big.(v2)) - 1))
end

function main()
  x = 10.^linspace(-16, -4, 1024)
  figure()
  loglog(x, fe.(pi/2-x))
  loglog(x, fe.(x))
end

end