Re: Slaying the HYPOTamus

On Mon, Aug 24, 2009 at 12:47:42PM -0500, Kevin Grittner wrote:
> David Fetter <david(at)fetter(dot)org> wrote:
> > On Mon, Aug 24, 2009 at 11:14:19PM +1000, Paul Matthews wrote:
>
> > These next two lines are a teensy bit baroque. Is there some
> > significant speed increase that would justify them?
> >
> >> if (x == 0.0)
> >> return 0.0;
> >> else {
> >> yx = y/x;
> >> return x*sqrt(1.0+yx*yx);
> >> }
> >> }
>
> I think the reason is overflow. From the function comment:
>
> >> * The traditional formulae of x^2+y^2 is rearranged
> >> * to bring x outside the sqrt. This allows computation of the
> hypotenuse
> >> * for much larger magnitudes than otherwise normally possible.
>
> Although I don't see why the first part isn't:
>
> if (y == 0.0)
> return x;

D'oh!

Good point :)

So the code should read as follows?

#include <math.h>
#include "c.h"
#include "utils/builtins.h"

/*
* Find the hypotenuse. Firstly x and y are swapped, if required, to make
* x the larger number. The traditional formulae of x^2+y^2 is rearranged
* to bring x outside the sqrt. This allows computation of the hypotenuse
* for much larger magnitudes than otherwise normally possible.
*
* sqrt( x^2 + y^2 ) = sqrt( x^2( 1 + y^2/x^2) )
* = x * sqrt( 1 + y^2/x^2 )
* = x * sqrt( 1 + y/x * y/x )
*/
double hypot( double x, double y )
{
double yx;

if( isinf(x) || isinf(y) )
return get_float8_infinity();

if( isnan(x) || isnan(y) )
return get_float8_nan();

x = fabs(x);
y = fabs(y);

if (x < y) {
double temp = x;
x = y;
y = temp;
}
if (y == 0.0)
return x;
yx = y/x;
return x*sqrt(1.0+yx*yx);
}

Cheers,
David.
--
David Fetter <david(at)fetter(dot)org> http://fetter.org/
Phone: +1 415 235 3778 AIM: dfetter666 Yahoo!: dfetter
Skype: davidfetter XMPP: david(dot)fetter(at)gmail(dot)com

Remember to vote!
Consider donating to Postgres: http://www.postgresql.org/about/donate