Re: Slaying the HYPOTamus

On Mon, Aug 24, 2009 at 12:47:42PM -0500, Kevin Grittner wrote:
> David Fetter <david@fetter.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
formulaeof x^2+y^2 is rearranged* to bring x outside the sqrt. This allows computation of the hypotenuse* for much
largermagnitudes 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;
returnx*sqrt(1.0+yx*yx);
 
}

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

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