Velocity Reviews > sqrtl(), etc. for Cygwin

# sqrtl(), etc. for Cygwin

glen herrmannsfeldt
Guest
Posts: n/a

 11-08-2012
Philip Lantz <(E-Mail Removed)> wrote:

(snip, I wrote)

>> > return(x/r+(r-x/r)/2);

> Mathematically, ((r + x/r)/2) and (x/r + (r - x/r)/2) are equivalent. Is
> there a numerical reason to prefer one over the other? (Obviously, if
> there are no such concerns, one might as well use the simpler one.)

In base 2 floating point, either one works. The final divide just
changes the exponent. In any other base, though, the simpler one can
lose precision.

In two digit decimal, (32+34)/2 is 33, but (64+6/2 is 65.

64+68 is 13e1 and 13e1/2 is 65.

It is used for IBM base 16 floating point and, will be needed
for the now part of the IEEE standard, decimal float.

-- glen

Phil Carmody
Guest
Posts: n/a

 11-12-2012
Philip Lantz <(E-Mail Removed)> writes:
> Steven G. Kargl wrote:
> > glen herrmannsfeldt wrote:
> > > Steven G. Kargl wrote:
> > >> #include <math.h>
> > >> long double
> > >> foo_sqrtl(long double x)
> > >> {
> > >> long double r;
> > >> r = sqrt((double)x);
> > >> return ((r + (x / r))/2);
> > >> }

> >
> > > ... the last one should be
> > >
> > > return(x/r+(r-x/r)/2);

> >
> > I beg to differ.

>
> Mathematically, ((r + x/r)/2) and (x/r + (r - x/r)/2) are equivalent. Is
> there a numerical reason to prefer one over the other? (Obviously, if
> there are no such concerns, one might as well use the simpler one.)

Given that r and x/r are very close (on the assumption that there's a good
initial estimate, which holds in this case), the (r - x/r) expression is a
case of catastrophic loss of precision, which is normally something you
want to avoid like the plague. Gut feel would make me prefer the original.

Phil
--
Regarding TSA regulations:
How are four small bottles of liquid different from one large bottle?
Because four bottles can hold the components of a binary liquid explosive,
whereas one big bottle can't. -- camperdave responding to MacAndrew on /.

glen herrmannsfeldt
Guest
Posts: n/a

 11-12-2012
Phil Carmody <(E-Mail Removed)> wrote:
> Philip Lantz <(E-Mail Removed)> writes:
>> Steven G. Kargl wrote:

(snip)

>> > >> return ((r + (x / r))/2);

(snip, then I wrote)

>> > > return(x/r+(r-x/r)/2);

(snip)

>> Mathematically, ((r + x/r)/2) and (x/r + (r - x/r)/2) are equivalent. Is
>> there a numerical reason to prefer one over the other? (Obviously, if
>> there are no such concerns, one might as well use the simpler one.)

> Given that r and x/r are very close (on the assumption that there's a good
> initial estimate, which holds in this case), the (r - x/r) expression is a
> case of catastrophic loss of precision, which is normally something you
> want to avoid like the plague. Gut feel would make me prefer the original.

That is certainly true in a large number of algorithms, but you should
look more carefully at this one.

Yes the difference can have a large error relative to its (possibly
small) magnitude, but the absolute error is the same as the two
operands if the exponents are the same.

Note that the more cancellation of low digits, means that you are
closer to the correct square root, and so is better, not worse!

Unless the sum is kept with extra precision, in non-binary (base 2)
floating point, a low digit can be lost in the sum that can't be
restored in the divide by two.

It is the lower (relative) precision of the difference that allows for
no loss on the divide by two.

-- glen

Phil Carmody
Guest
Posts: n/a

 11-13-2012
glen herrmannsfeldt <(E-Mail Removed)> writes:
> Phil Carmody <(E-Mail Removed)> wrote:
> > Philip Lantz <(E-Mail Removed)> writes:
> >> Steven G. Kargl wrote:

>
> (snip)
>
> >> > >> return ((r + (x / r))/2);

>
> (snip, then I wrote)
>
> >> > > return(x/r+(r-x/r)/2);

>
> (snip)
>
> >> Mathematically, ((r + x/r)/2) and (x/r + (r - x/r)/2) are equivalent. Is
> >> there a numerical reason to prefer one over the other? (Obviously, if
> >> there are no such concerns, one might as well use the simpler one.)

>
> > Given that r and x/r are very close (on the assumption that there's a good
> > initial estimate, which holds in this case), the (r - x/r) expression is a
> > case of catastrophic loss of precision, which is normally something you
> > want to avoid like the plague. Gut feel would make me prefer the original.

>
> That is certainly true in a large number of algorithms, but you should
> look more carefully at this one.
>
> Yes the difference can have a large error relative to its (possibly
> small) magnitude, but the absolute error is the same as the two
> operands if the exponents are the same.
>
> Note that the more cancellation of low digits, means that you are
> closer to the correct square root, and so is better, not worse!

Agreed. Hence my "normally".

> Unless the sum is kept with extra precision, in non-binary (base 2)
> floating point, a low digit can be lost in the sum that can't be
> restored in the divide by two.

Anyone not using binary FP deserves all ills that come to them!

> It is the lower (relative) precision of the difference that allows for
> no loss on the divide by two.

The thing I like about the sum is its symmetry.
With the half-the-difference method, there are two ways of approaching the
middle.
(r + (x/r-r))/2
(x/r + (r-x/r))/2
I don't see why one should chose one rather the other. And in particular
why you should have chosen the one with the longer expression.

Phil
--
Regarding TSA regulations:
How are four small bottles of liquid different from one large bottle?
Because four bottles can hold the components of a binary liquid explosive,
whereas one big bottle can't. -- camperdave responding to MacAndrew on /.

glen herrmannsfeldt
Guest
Posts: n/a

 11-13-2012
Phil Carmody <(E-Mail Removed)> wrote:

(snip, I wrote)
>> >> > > return(x/r+(r-x/r)/2);

(snip)
>> >> Mathematically, ((r + x/r)/2) and (x/r + (r - x/r)/2) are equivalent. Is
>> >> there a numerical reason to prefer one over the other? (Obviously, if
>> >> there are no such concerns, one might as well use the simpler one.)

>> > Given that r and x/r are very close (on the assumption that there's a good
>> > initial estimate, which holds in this case), the (r - x/r) expression is a
>> > case of catastrophic loss of precision, which is normally something you
>> > want to avoid like the plague. Gut feel would make me prefer the original.

(snip)
>> Unless the sum is kept with extra precision, in non-binary (base 2)
>> floating point, a low digit can be lost in the sum that can't be
>> restored in the divide by two.

> Anyone not using binary FP deserves all ills that come to them!

>> It is the lower (relative) precision of the difference that allows for
>> no loss on the divide by two.

> The thing I like about the sum is its symmetry.
> With the half-the-difference method, there are two ways of approaching the
> middle.

> (r + (x/r-r))/2
> (x/r + (r-x/r))/2

> I don't see why one should chose one rather the other. And in particular
> why you should have chosen the one with the longer expression.

Once the values are in registers, it doesn't make so much difference,
but yes, I probably would do the former, to avoid the possibility of
doing the divide twice.

Assembler. IBM S/360 uses base 16 floating point. All the techniques
to avoid precision loss were well studied at the time, but maybe
now forgotten. Now, IEEE has standardized decimal float, so all
these will be brought back out again.

-- glen

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are Off Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post vc Firefox 10 02-28-2004 01:05 AM cygwin-help@cygwin.com Python 0 09-05-2003 04:42 PM cygwin-help@cygwin.com Python 1 09-05-2003 07:46 AM cygwin-announce-help@cygwin.com Python 0 09-05-2003 01:29 AM cygwin-help@cygwin.com Python 0 09-04-2003 06:34 PM