Re: Slaying the HYPOTamus - Mailing list pgsql-hackers

From David Fetter
Subject Re: Slaying the HYPOTamus
Date
Msg-id 20090824140713.GC5896@fetter.org
Whole thread Raw
In response to Re: Slaying the HYPOTamus  (Paul Matthews <plm@netspace.net.au>)
Responses Re: Slaying the HYPOTamus  ("Kevin Grittner" <Kevin.Grittner@wicourts.gov>)
Re: Slaying the HYPOTamus  (Sam Mason <sam@samason.me.uk>)
List pgsql-hackers
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@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


pgsql-hackers by date:

Previous
From: Robert Haas
Date:
Subject: Re: 8.5 release timetable, again
Next
From: Andrew Dunstan
Date:
Subject: Re: hba load error and silent mode