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

[Yitzchak Scott-Thoennes]

On Wed, 1 Oct 2003, Fred Morris wrote:

> At 4:05 PM 10/1/03, Itay Furman wrote:
> >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?

Whatever the underlying C math library uses.  You might have better luck
using sqrt instead of ** (which will call the C pow() or powl()).
But I suspect decent libraries will have .5 special-cased in pow()
to end up the same as sqrt().

> I've been painfully reminded that that is the expected behavior of Perl on
> several occasions. This can bite you even when trying to figure out what
> day of the week it is. There is some Math:: stuff. A long time ago I
> actually took Numerical Analysis, from which I learned:
>
> Several techniques for minimizing loss of precision, which contribute to
> unreadable code in most circumstances.
>
> That zero is only zero when you say it is, theoretically speaking.. unless
> you confine yourself to integer math.. which I don't believe OTS Perl is
> capable of (however, see Math::BigInt).

perldoc integer

(not for sqrt or **, though).

> OTOH, prime numbers are a lot of fun.
>
>
> You probably want some elucidation of exactly what standards-conformant
> floating point representation Perl uses; I can't help you there.

Whatever C uses for a double or long double.

> >## Smallest vector that yields the flow
> >$ test_innerprod.pl 1 1
> >                                sos=2   residue=2.22044604925031e-16
> >## 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
> >
>
> Interesting. Why is it always the same? Anybody got a calculator (with that
> kind of precision) handy?

$perl -wle'print log(2.22044604925031e-16)/log(2)'
-52

i.e. 2.2204...e-16 == 2**-52

So you probably have IEEE doubles with a 53 bit mantissa and the lowest
bit is wrong.  I suspect your workaround will sometimes fail, too.
But I'm not sure what you are trying to accomplish--  x/(x**.5 * x**.5)
is supposed to always be one.  Are you sure you don't mean to be
getting the ratio of the sum of squares to the square of the sum?

> It does seem to me that your example points the way to a technique for
> determining the precision (which was probably covered, I don't remember):
> take 1, and add a value to it which decreases by O(n**x) where n is 10 or 2
> and x is a decreasing series integer value until there is no effect. That
> will tell you the resolution on the platform on which you're running... or
> else you'll get an underflow error. :-\
>
> Post your (readable) code and I'll run it on a few machines.

The mantissa will be at least `perl -V:nv_preserves_uv_bits` bits.
(nv_preserves_uv_bits will be the lesser of the number of bits in
your unsigned integers and the mantissa of your floating point values.)

(Most of this is out of perl's hands and dependent on your libc/math lib.
There are problems with some versions of perl with converting a string
to a floating point number that have been fixed, though.)

(Anybody bored should check out Ilya's recent post to clpmoderated
for the metaphor of the year (though unsuitable for Time magazine's
cover)).