a more precise distance algorithm

[Oscar Benjamin]

On 27 May 2015 at 19:00, Brian Blais <bblais at gmail.com> wrote:
> On Mon, May 25, 2015 at 11:11 PM, Steven D'Aprano <steve at pearwood.info> wrote:
>>
>> Let's compare three methods.
>>
>> 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)
>
>
>>     d1 = naive(a, b)
>>     d2 = alternate(a, b)
>>     d3 = math.hypot(a, b)
>>
>
>> which shows that:
>>
>> (1) It's not hard to find mismatches;
>> (2) It's not obvious which of the three methods is more accurate.
>>
>
> Bottom line: they all suck.  :)
>
> I ran the program you posted, and, like you, got the following two examples:
>
> for fun in [naive, alternate, math.hypot]:
>     print '%.20f' % fun(222.44802484683657,680.255801504161)
>
> 715.70320611153294976248
> 715.70320611153283607564
> 715.70320611153283607564
>
> and
>
> for fun in [naive, alternate, math.hypot]:
>     print '%.20f' % fun(376.47153302262484,943.1877995550265)
>
> 1015.54617837194291496417
> 1015.54617837194280127733
> 1015.54617837194291496417
>
> but when comparing to Wolfram Alpha, which calculates these out many
> more decimal places, we have for the two cases:
>
> 715.7032061115328768204988784125331443593766145937358347357252...
> 715.70320611153294976248
> 715.70320611153283607564
> 715.70320611153283607564
>
> 1015.546178371942943007625196455666280385821355370154991424749...
> 1015.54617837194291496417
> 1015.54617837194280127733
> 1015.54617837194291496417
>
> where all of the methods deviate at the 13/14 decimal place.

So they have 12/13 correct digits after the decimal point. Including
the digits before the decimal point they all have 15/16 correct
decimal digits. This is exactly what you should expect when using
double precision floating point. Feel free to print as many digits as
you like but the extra ones won't add any more accuracy.

The difference between the answers you showed are just 1 ULP:
>>> a = 715.70320611153283607564
>>> b = 715.70320611153294976248
>>> a
715.7032061115328
>>> a + 5e-14
715.7032061115328
>>> a + 6e-14
715.703206111533
>>> a + 6e-14 == b
True

Since you don't show the source numbers we can't know whether the true
result lies between these two or not.


Oscar