On Mon, Aug 24, 2009 at 12:47:42PM -0500, Kevin Grittner wrote: > David Fetter <da...@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 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 <da...@fetter.org> http://fetter.org/ Phone: +1 415 235 3778 AIM: dfetter666 Yahoo!: dfetter Skype: davidfetter XMPP: david.fet...@gmail.com Remember to vote! Consider donating to Postgres: http://www.postgresql.org/about/donate -- Sent via pgsql-hackers mailing list (pgsql-hackers@postgresql.org) To make changes to your subscription: http://www.postgresql.org/mailpref/pgsql-hackers