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