SPUG: S.O.S: Sum-Of-Squares (nickname:Help)

[Fred Morris]

I suppose I should be nice and answer this question, even though it makes
my head hurt. However, there's a fundamental lesson here that languages
cannot solve thinking problems, and leaving that unanswered made my head
continue to hurt. So, I sincerely hope that this explanation is useful.


In a nutshell, there is no rational representation for the square root of 2.

>The problem: for some cases I get a small, non-zero deviation.
>My question: does any one have some insight or rules-of-a-thumb as to
>the expected numerical precision in basic arithmetic operations under Perl?
>
>I should point out that I have found a simple workaround (see further below),
>but would like some insight for future reference.
>
>        Thanks,
>        Itay
>
>[...]
>#!/usr/local/bin/perl -w

BTW, Perl lives at /usr/bin/perl on all three of my test boxes.

>use strict;
>my @elem = @ARGV;               # Elements of vector.
>my $sos  = 0;                   # Sum-Of-Squares
>$sos    += $_**2 foreach (@elem);
>print "\t"x4, "sos=$sos\tresidue=",
>  1-$sos/($sos**0.5 * $sos**0.5), "\n";
>exit;

I get divide by zero on all platforms without any arguments. :-p

>#### Some examples ####
>
>## OK
>$ test_innerprod.pl 1 1 1 1
>                                sos=4   residue=0
>$ test_innerprod.pl 1 0 0 0
>                                sos=1   residue=0
>## Wrong !
>$ test_innerprod.pl 1 1 1 0
>                                sos=3   residue=-2.22044604925031e-16
>## Number of 1's doesn't matter, necessarily
>$ test_innerprod.pl 1 1 0 0
>                                sos=2   residue=2.22044604925031e-16
>## Order doesn't matter
>$ test_innerprod.pl 1 0 1 0
>                                sos=2   residue=2.22044604925031e-16
>## Length doesn't matter
>$ test_innerprod.pl 1 1 0
>                                sos=2   residue=2.22044604925031e-16
>## OK
>$ test_innerprod.pl 1 0 0
>                                sos=1   residue=0
>## Smallest vector that yields the flow
>$ test_innerprod.pl 1 1
>                                sos=2   residue=2.22044604925031e-16

This should have been the clue: recast this statement as "simple scalars do
not leave a residue."

>## These don't have to be 1's
>$ test_innerprod.pl 0.5 0.5 0 0
>                                sos=0.5 residue=2.22044604925031e-16
>$ test_innerprod.pl 0.5 0.5
>                                sos=0.5 residue=2.22044604925031e-16
>

Without going exhaustively through all of these cases, I get the same
values on SuSE 8.2/Perl 5.8, SuSE 6.4/Perl 5.5.3, Solaris (Sparc)/Perl
5.5.3.

Earlier I observed that the residue was always the same. [3,3,3,1] (sos 28)
is a case where it's not.

Here's what's going on: Things which yield a rational (and representable)
square root are reversible, and things which don't aren't. This is not
surprising. (In fact, I'm more surprised that [.1,.1,.1,.1] doesn't leave a
residue, since 1/10 is a repeating decimal in base 2. I'm sure if I thought
about it long enough, I'd divine the reason.)

>#### My workaround ####
>
>$ cat test_innerprod.pl
>#!/usr/local/bin/perl -w
>use strict;
>my @elem = @ARGV;               # Elements of vector.
>my $sos  = 0;                   # Sum-Of-Squares
>$sos    += $_**2 foreach (@elem);
>print "\t"x4, "sos=$sos\tresidue=",
>  #################################################
>  # Take the square-root AFTER making the product #
>  # Compare with previous version                 #
>  #   1-$sos/($sos**0.5 * $sos**0.5)              #
>  #################################################
>  1-$sos/($sos*$sos)**0.5, "\n";
>exit;

That's perfectly reasonable.


As a general philosophy, it is best to assume that when you're trying to do
"high-precision" math with floating points... you can't! The language is
irrelevant. Languages (and implementations of languages) do vary in their
handling of representational errors, but that's a different.


A general practice is to round or floor at important points based on
analysis of the algorithm (although refactoring is another valuable
technique). Never assume that the result of a floating point calc is
accurate to the limits of the mantissa; never represent as much, either.

Analysis (and appropriate rounding/flooring) should allow you to assert
that the result of the calc is accurate to some number of places. Flooring
is often preferred to rounding because it can allow the direction of the
error to be predicted better. (Just make sure you know the difference
between "flooring" towards 0, and flooring towards -1*infinity.)


HTH.

--

Fred Morris
m3047 at inwa.net