Work in progress towards robust complex division algorithm in Julia

output.txt

The archive paper provides 10 example hard complex divions with answers.  It also provides an algorith for measuring the accuracy of the result in "bits", which I have implemented.  The robust cdiv algorithm is said to get 53 bits for all but #8 on the hard problems, where it gets 52 bits. I reproduce all of the results (also the default julia / reproduces results for smith's algorithm), except I get 0 bits on #5.  I've tracked this down to the division of b/c returning zero inside robust_cdiv2, after r is also zero. But it's not clear how to fix it.


cdiv #1 ,53 bits accurate, 1.1125369292536007e-308 - 1.1125369292536007e-308im
cdiv #2 ,53 bits accurate, 8.98846567431158e307 + 0.0im
cdiv #3 ,53 bits accurate, 1.4334366349937947e104 - 3.645561009778199e-304im
cdiv #4 ,53 bits accurate, 8.98846567431158e307 + 0.0im
cdiv #5 ,0 bits accurate, 3.757668132438133e109 - 2.0e-323im
cdiv #6 ,53 bits accurate, 2.0e-323 + 1.048576e6im
cdiv #7 ,53 bits accurate, 3.8981256045591133e289 + 8.174961907852354e295im
cdiv #8 ,52 bits accurate, 0.6000000000000001 + 0.2im
cdiv #9 ,53 bits accurate, 0.001953125 - 0.001953125im
cdiv #10 ,53 bits accurate, 1.0295115178936058e-84 + 6.971459875150762e-220im
/ #1 ,53 bits accurate, 1.1125369292536007e-308 - 1.1125369292536007e-308im
/ #2 ,53 bits accurate, 8.98846567431158e307 + 0.0im
/ #3 ,0 bits accurate, 1.4334366349937947e104 + 0.0im
/ #4 ,0 bits accurate, Inf + 0.0im
/ #5 ,0 bits accurate, 3.757668132438133e109 + 0.0im
/ #6 ,53 bits accurate, 2.0e-323 + 1.048576e6im
/ #7 ,41 bits accurate, 3.89812560456e289 + 8.174961907854212e295im
/ #8 ,0 bits accurate, 0.5 + 0.5im
/ #9 ,0 bits accurate, 0.0 - 0.0im
/ #10 ,5 bits accurate, 1.0295115178936058e-84 + 7.082117968407124e-220im

robust_cdiv.jl

function cdiv(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
# see arxiv.1210.4539 and a c version at link
# http://forge.scilab.org/index.php/p/compdiv/source/tree/HEAD/src/c/compdiv.c
#             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 && println("#1 $a $b")
	ab >= half*ov && (a=half*a; b=half*b; s=two*s )
	#cd >= half*ov && println("#2 $a $b")
	cd >= half*ov && (c=half*c; d=half*d; s=s*half)
	#ab <= un*two/ϵ && println("#3 $a $b")
	ab <= un*two/ϵ && (a=a*bs; b=b*bs; s=s/bs      )
	#cd <= un*two/ϵ && println("#4 $a $b")
	cd <= un*two/ϵ && (c=c*bs; d=d*bs; s=s*bs      )
	#abs(d)<=abs(c)  && (c=c*bs; d=d*bs; s=s*bs     ) # in LAPACK but not arxiv
	#abs(d)<=abs(c) ? println("path 1 #1") : println("path1 #2")
	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
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)
	#println("r $r t $t")
	p::T=robust_cdiv2(a,b,c,d,r,t)
	q::T=robust_cdiv2(b,-a,c,d,r,t)
	#println("cdiv1 p,q $p, $q")
	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
		#println("path 2 #1")
		#br != zero ? println("path 3 #1") : println("path 3 #2")
		return (br != zero ? (a+br)*t : a*t + (b*t)*r)
	else
		#println("path 2 #2 a $a b $b c $c d $d")
		return (a + d*(b/c)) * t
	end
end

# hard divisions from Fig 6 of arxiv paper
z7 = (3.898125604559113300e289, 8.174961907852353577e295)
z9 = (0.001953125, -0.001953125)
z10 = ( 1.02951151789360578e-84, 6.97145987515076231e-220)
harddivs = ((1.0, 1.0, 1.0, 2^1023.0, 2^-1023.0, -2^-1023.0), #1
			(1.0, 1.0, 2^-1023.0, 2^-1023.0, 2^1023.0, 0.0), #2
			(2^1023.0, 2^-1023.0, 2^677.0, 2^-677.0, 2^346.0, -2^-1008.0), #3
			(2^1023.0, 2^1023.0, 1.0, 1.0, 2^1023.0, 0.0), #4
			(2^1020.0, 2^-844., 2^656., 2^-780., 2^364., 2^-1072.), #5
			(2^-71., 2^1021., 2^1001., 2^-323., 2^-1072., 2^20.), #6
			(2^-347.,2^-54., 2^-1037., 2^-1058., z7[1], z7[2]), #7
			(2^-1074., 2^-1074., 2^-1073., 2^-1074., 0.6, 0.2), #8
			(2^1015., 2^-989., 2^1023., 2^1023., z9[1], z9[2]), #9
			(2^-622., 2^-1071., 2^-343., 2^-798., z10[1], z10[2])#10
			)

function 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

function ccdiv(f,harddivs)
	for i = 1:length(harddivs)
		h = harddivs[i]
		x=Complex{Float64}(h[1],h[2])
		y=Complex{Float64}(h[3],h[4])
		z=Complex{Float64}(h[5],h[6])
		out = f(x,y)
		acc = accuracy(out,z)
		println("$f #$i ,$acc bits accurate, $out")
	end	
end

ccdiv(cdiv,harddivs)
ccdiv(/,harddivs)

It seems like the comparison value for case 5 is incorrect. If I plug this into wolfram alpha, the imaginary part comes back negative. In your test code, the expected imaginary part is positive. In other words, with that fix, it seems to work for me!