Re: Slaying the HYPOTamus

On Mon, Aug 24, 2009 at 11:14:19PM +1000, Paul Matthews wrote:
> This is to go with the hypot() patch I posted the other day.
>
> As other code, such as found in adt/float.c and adt/numeric.c, simply
> assumes that isnan() exists, despite it being a C99 function, I have
> assumed the same.
>
> The below code should be placed into a file called src/port/hypot.c.
>
> Unfortunately I do not know how to drive configure and all the other
> associated build magics. Could some kind soul please implemented that
> portion. (Or shove me in the right direction)

Comments below :)

> #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;
> }

How about putting:

if (x >= 0.0 && y == 0.0)
return x;

here?

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);
> }
> }

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