a more precise distance algorithm

[Steven D'Aprano]

On Tue, 26 May 2015 05:21 am, ravas wrote:

> I read an interesting comment:
> """
> The coolest thing I've ever discovered about Pythagorean's Theorem is an
> alternate way to calculate it. If you write a program that uses the
> distance form c = sqrt(a^2 + b^2) you will suffer from the lose of half of
> your available precision because the square root operation is last. A more
> accurate calculation is c = a * sqrt(1 + b^2 / a^2). If a is less than b,
> you should swap them and of course handle the special case of a = 0. """


Let's compare three methods. 


import math
import random

def naive(a, b):
    return math.sqrt(a**2 + b**2)

def alternate(a, b):
    a, b = min(a, b), max(a, b)
    if a == 0:  return b
    if b == 0:  return a
    return a * math.sqrt(1 + b**2 / a**2)

counter = 0
print("Type Ctrl-C to exit")
while True:
    counter += 1
    a = random.uniform(0, 1000)
    b = random.uniform(0, 1000)
    d1 = naive(a, b) 
    d2 = alternate(a, b)
    d3 = math.hypot(a, b)
    if not (d1 == d2 == d3):
        print("mismatch after %d trials" % counter)
        print("naive:", d1)
        print("alternate:", d2)
        print("hypot:", d3)
        break



When I run that, I get:

mismatch after 3 trials
naive: 1075.6464259886257
alternate: 1075.6464259886257
hypot: 1075.6464259886254

A second run gives:

mismatch after 3 trials
naive: 767.3916150255787
alternate: 767.3916150255789
hypot: 767.3916150255787


which shows that:

(1) It's not hard to find mismatches;
(2) It's not obvious which of the three methods is more accurate.



-- 
Steven