Marc:
I've started implementing the codes in Kahan's paper "Branch Cuts for Complex Elementary Functions, or Much Ado About Nothing's Sign Bit". The first nontrivial code is for complex magnitude, which he calls cabs.
I've met Kahan a number of times, mainly when I was a graduate student, and he's a very bright guy. A bit scary, actually.
Anyway, I've now compared his code to our current magnitude code for 1,000,000 pairs of floating-point numbers exponentially distributed (the logarithms are distributed uniformly) between the smallest positive flonum and the largest finite flonum.
A correctly-rounded result has an error of less than 1/2 unit in the last place (ulp). I recorded the number of times the answers were incorrect (but still within one ulp of the true answer) and the number of times the error was greater than 1 ulp. The correct answers were computed using my computable reals code.
The results were as follows:
Current Gambit magnitude:
Incorrect, but within 1 ulp: 7909 Incorrect, error > 1 ulp: 2132
Kahan's cabs:
Incorrect, but within 1 ulp: 150 Incorrect, error > 1 ulp: 0
Like I said, a bit scary.
Brad
Afficher les réponses par date
On 11/18/2013 06:21 PM, Bradley Lucier wrote:
Anyway, I've now compared his code to our current magnitude code for 1,000,000 pairs of floating-point numbers exponentially distributed (the logarithms are distributed uniformly) between the smallest positive flonum and the largest finite flonum.
A correctly-rounded result has an error of less than 1/2 unit in the last place (ulp). I recorded the number of times the answers were incorrect (but still within one ulp of the true answer) and the number of times the error was greater than 1 ulp. The correct answers were computed using my computable reals code.
The results were as follows:
Current Gambit magnitude:
Incorrect, but within 1 ulp: 7909 Incorrect, error > 1 ulp: 2132
Kahan's cabs:
Incorrect, but within 1 ulp: 150 Incorrect, error > 1 ulp: 0
Here's a more rigorous test. I set epsilon equal to the smallest positive number such that (= (- 1.0 epsilon) 1.0) is true, and ran 1,000,000 random tests with each real and imaginary part exponentially distributed between (sqrt epsilon) and 1:
Current Gambit Magnitude:
Incorrect, but within 1 ulp: 315879 Incorrect, error > 1 ulp: 78337
Kahan's cabs:
Incorrect, but within 1 ulp: 10959 Incorrect, error > 1 ulp: 0
Put like this, where the current code returns a (very slightly) incorrect answer about 40% of the time, it should be replaced.
Brad
-
Bradley Lucier