RFC: robust division for Complex{Float64}

base/complex.jl

@@ -197,7 +197,54 @@ function /{T<:Real}(a::T, b::Complex{T})
         complex(a/den, -a*r/den)
     end
 end
-
+function /(a::Complex{Float64}, b::Complex{Float64})
+     r,i = robust_cdiv(real(a), imag(a), real(b), imag(b))
+     Complex{Float64}(r,i)
+end
+# robust_cdiv performs complex division in real arithmetic
+# the first step is to scale variables if appropriate
+# then do calculations in 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{T<:Float64}(a::T,b::T,c::T,d::T)
+    const half::T = 0.5
+    const two::T = 2.0
+    const ab::T = max(abs(a), abs(b))
+    const cd::T = max(abs(c), abs(d))
+    const ov::T = realmax(T)
+    const un::T = realmin(T)
+    const ϵ::T = eps(T)
+    const bs=two/(ϵ*ϵ)
+    s::T = 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 p*s,q*s # undo scaling
+end
+function robust_cdiv1{T<:Float64}(a::T,b::T,c::T,d::T)
+    const one::T = 1.0
+    const r::T=d/c
+    const t::T=one/(c+d*r)
+    p::T=robust_cdiv2(a,b,c,d,r,t)
+    q::T=robust_cdiv2(b,-a,c,d,r,t)
+    return p,q
+end
+function robust_cdiv2{T<:Float64}(a::T,b::T,c::T,d::T,r::T,t::T)
+    const zero::T = 0.0
+    if r != zero
+        const br::T = b*r
+        return (br != zero ? (a+br)*t : a*t + (b*t)*r)
+    else
+        return (a + d*(b/c)) * t
+    end
+end
+
 function ssqs{T<:FloatingPoint}(x::T, y::T)
     k::Int = 0
     ρ = x*x + y*y

test/complex.jl

@@ -598,3 +598,46 @@
 @test complex(2,2)^2 === complex(0,8)
 @test_throws complex(2,2)^(-2)
 @test complex(2.0,2.0)^(-2) === complex(0.0, -0.125)
+
+# robust division of Float64
+# hard complex divisions from Fig 6 of arxiv.1210.4539
+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
+@test 52 <= minimum([sb_accuracy(h[1]/h[2],h[3]) for h in harddivs])
+
+# division of non-Float64
+function cdiv_test(a,b)
+  c=convert(Complex{Float64},a)/convert(Complex{Float64},b)
+  50 <= sb_accuracy(c,convert(Complex{Float64},a/b))
+end
+@test cdiv_test(complex(1//2, 3//4), complex(17//13, 4//5))
+@test cdiv_test(complex(1,2), complex(8997,2432))