robust complex division benchmark

bench_cdiv.jl

function bench_smith(n,d)
    for i = 1:length(n)
        smith_cdiv(n[i], d[i])
    end
end

function bench_robust(n,d)
    for i = 1:length(n)
        robust_cdiv(n[i], d[i])
    end
end
function bench_naive_test(n,d)
    for i = 1:length(n)
        naive_test_cdiv(n[i], d[i])
    end
end

function naive_test_cdiv(a::Complex128, b::Complex128)
    c = naive_cdiv(a,b)
    cr,ci = reim(c)
    if !isfinite(cr) || !isfinite(ci) || (cr == 0.0 && real(a) != 0.0) || (ci == 0.0 && imag(a) != 0.0)
        return robust_cdiv(a,b)
    end
    c
end

naive_cdiv(a::Complex, b::Complex) = (a * conj(b)) / abs2(b)
function smith_cdiv(a::Complex, b::Complex)
    are = real(a); aim = imag(a); bre = real(b); bim = imag(b)
    if abs(bre) <= abs(bim)
        if isinf(bre) && isinf(bim)
            r = sign(bre)/sign(bim)
        else
            r = bre / bim
        end
        den = bim + r*bre
        complex((are*r + aim)/den, (aim*r - are)/den)
    else
        if isinf(bre) && isinf(bim)
            r = sign(bim)/sign(bre)
        else
            r = bim / bre
        end
        den = bre + r*bim
        complex((are + aim*r)/den, (aim - are*r)/den)
    end
end

# robust complex division for double precision
# the first step is to scale variables if appropriate ,then do calculations
# in a way that avoids over/underflow (subfuncs 1 and 2), then undo the scaling.
# scaling variable s and other techniques
# based on arxiv.1210.4539
#             a + i*b
#  p + i*q = ---------
#             c + i*d
function robust_cdiv(z::Complex128, w::Complex128)
    a, b = reim(z); c, d = reim(w)
    half = 0.5
    two = 2.0
    ab = max(abs(a), abs(b))
    cd = max(abs(c), abs(d))
    ov = realmax(a)
    un = realmin(a)
    ϵ = eps(Float64)
    bs = two/(ϵ*ϵ)
    s = 1.0
    ab >= half*ov  && (a=half*a; b=half*b; s=two*s ) # scale down a,b
    cd >= half*ov  && (c=half*c; d=half*d; s=s*half) # scale down c,d
    ab <= un*two/ϵ && (a=a*bs; b=b*bs; s=s/bs      ) # scale up a,b
    cd <= un*two/ϵ && (c=c*bs; d=d*bs; s=s*bs      ) # scale up c,d
    abs(d)<=abs(c) ? ((p,q)=robust_cdiv1(a,b,c,d)  ) : ((p,q)=robust_cdiv1(b,a,d,c); q=-q)
    return Complex128(p*s,q*s) # undo scaling
end
function robust_cdiv1(a::Float64, b::Float64, c::Float64, d::Float64)
    r = d/c
    t = 1.0/(c+d*r)
    p = robust_cdiv2(a,b,c,d,r,t)
    q = robust_cdiv2(b,-a,c,d,r,t)
    return p,q
end
function robust_cdiv2(a::Float64, b::Float64, c::Float64, d::Float64, r::Float64, t::Float64)
    if r != 0
        br = b*r
        return (br != 0 ? (a+br)*t : a*t + (b*t)*r)
    else
        return (a + d*(b/c)) * t
    end
end



z7 = Complex{Float64}(3.898125604559113300e289, 8.174961907852353577e295)
z9 = Complex{Float64}(0.001953125, -0.001953125)
z10 = Complex{Float64}( 1.02951151789360578e-84, 6.97145987515076231e-220)
harddivs = ((1.0+im*1.0, 1.0+im*2^1023.0, 2^-1023.0-im*2^-1023.0), #1
      (1.0+im*1.0, 2^-1023.0+im*2^-1023.0, 2^1023.0+im*0.0), #2
      (2^1023.0+im*2^-1023.0, 2^677.0+im*2^-677.0, 2^346.0-im*2^-1008.0), #3
      (2^1023.0+im*2^1023.0, 1.0+im*1.0, 2^1023.0+im*0.0), #4
      (2^1020.0+im*2^-844., 2^656.0+im*2^-780.0, 2^364.0-im*2^-1072.0), #5
      (2^-71.0+im*2^1021., 2^1001.0+im*2^-323.0, 2^-1072.0+im*2^20.0), #6
      (2^-347.0+im*2^-54., 2^-1037.0+im*2^-1058.0, z7), #7
      (2^-1074.0+im*2^-1074., 2^-1073.0+im*2^-1074., 0.6+im*0.2), #8
      (2^1015.0+im*2^-989., 2^1023.0+im*2^1023.0, z9), #9
      (2^-622.0+im*2^-1071., 2^-343.0+im*2^-798.0, z10)#10
      )

# calculate "accurate bits" in range 0:53 by algorithm given in arxiv.1210.4539
function sb_accuracy(x,expected)
  min(logacc(real(x),real(expected)),
    logacc(imag(x),imag(expected)) )
end
relacc(x,expected) = abs(x-expected)/abs(expected)
function logacc(x::Float64,expected::Float64)
  x == expected && (return 53)
  expected == 0 && (return 0)
  (x == Inf || x == -Inf) && (return 0)
  isnan(x) && (return 0)
  ra = relacc(BigFloat(x),BigFloat(expected))
  max(ifloor(-log2(ra)),0)
end
# the robust division algorithm should have 53 or 52
# bits accuracy for each of the hard divisions
[sb_accuracy(h[1]/h[2],h[3]) for h in harddivs]

n = rand(Complex{Float64},1000000)
d = rand(Complex{Float64},1000000)
bench_smith(n,d)
bench_robust(n,d)
print("robust/current: ")
@time bench_robust(n,d)
println("")
println([sb_accuracy(robust_cdiv(h[1],h[2]),h[3]) for h in harddivs]')
print("smith/old: ")
@time bench_smith(n,d)
println("")
println([sb_accuracy(smith_cdiv(h[1],h[2]),h[3]) for h in harddivs]')
print("naive_test: ")
[sb_accuracy(naive_test_cdiv(h[1],h[2]),h[3]) for h in harddivs]'
@time bench_naive_test(n,d)
println("")
println([sb_accuracy(naive_test_cdiv(h[1],h[2]),h[3]) for h in harddivs]')

output.txt

robust/current: elapsed time: 0.098058714 seconds (64006848 bytes allocated)

53	53	53	53	53	53	53	52	53	53

smith/old: elapsed time: 0.019170227 seconds (48 bytes allocated)

53	53	0	0	0	53	41	0	0	5

naive_test: elapsed time: 0.01587455 seconds (88712 bytes allocated)

53	53	53	53	53	53	53	52	53	53