mirrors/postgres
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Replace the naive HYPOT() macro with a standards-conformant hypotenuse

Browse Source 
function.  This avoids unnecessary overflows and probably gives a more
accurate result as well.

Paul Matthews, reviewed by Andrew Geery
This commit is contained in:
Tom Lane 2010-08-03 21:21:03 +00:00
parent 90a391c645
commit db04f2b322
2 changed files with 64 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

62
src/backend/utils/adt/geo_ops.c
View File

@ -8,7 +8,7 @@
*
*
* IDENTIFICATION
* $PostgreSQL: pgsql/src/backend/utils/adt/geo_ops.c,v 1.108 2010/02/26 02:01:08 momjian Exp $
* $PostgreSQL: pgsql/src/backend/utils/adt/geo_ops.c,v 1.109 2010/08/03 21:21:03 tgl Exp $
*
*-------------------------------------------------------------------------
*/
@ -5410,3 +5410,63 @@ plist_same(int npts, Point *p1, Point *p2)
return FALSE;
}
/*-------------------------------------------------------------------------
* Determine the hypotenuse.
*
* If required, x and y are swapped to make x the larger number. The
* traditional formula of x^2+y^2 is rearranged to factor x outside the
* sqrt. This allows computation of the hypotenuse for significantly
* larger values, and with a higher precision than when using the naive
* formula. In particular, this cannot overflow unless the final result
* would be out-of-range.
*
* 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 )
*
* It is expected that this routine will eventually be replaced with the
* C99 hypot() function.
*
* This implementation conforms to IEEE Std 1003.1 and GLIBC, in that the
* case of hypot(inf,nan) results in INF, and not NAN.
*-----------------------------------------------------------------------
*/
double
pg_hypot(double x, double y)
{
double yx;
/* Handle INF and NaN properly */
if (isinf(x) || isinf(y))
return get_float8_infinity();
if (isnan(x) || isnan(y))
return get_float8_nan();
/* Else, drop any minus signs */
x = fabs(x);
y = fabs(y);
/* Swap x and y if needed to make x the larger one */
if (x < y)
{
double temp = x;
x = y;
y = temp;
}
/*
* If y is zero, the hypotenuse is x. This test saves a few cycles in
* such cases, but more importantly it also protects against
* divide-by-zero errors, since now x >= y.
*/
if (y == 0.0)
return x;
/* Determine the hypotenuse */
yx = y / x;
return x * sqrt(1.0 + (yx * yx));
}

5
src/include/utils/geo_decls.h
View File

@ -6,7 +6,7 @@
* Portions Copyright (c) 1996-2010, PostgreSQL Global Development Group
* Portions Copyright (c) 1994, Regents of the University of California
*
* $PostgreSQL: pgsql/src/include/utils/geo_decls.h,v 1.57 2010/01/14 16:31:09 teodor Exp $
* $PostgreSQL: pgsql/src/include/utils/geo_decls.h,v 1.58 2010/08/03 21:21:03 tgl Exp $
*
* NOTE
* These routines do *not* use the float types from adt/.
@ -50,7 +50,7 @@
#define FPge(A,B) ((A) >= (B))
#endif
#define HYPOT(A, B) sqrt((A) * (A) + (B) * (B))
#define HYPOT(A, B) pg_hypot(A, B)
/*---------------------------------------------------------------------
* Point - (x,y)
@ -211,6 +211,7 @@ extern Datum point_div(PG_FUNCTION_ARGS);
/* private routines */
extern double point_dt(Point *pt1, Point *pt2);
extern double point_sl(Point *pt1, Point *pt2);
extern double pg_hypot(double x, double y);
/* public lseg routines */
extern Datum lseg_in(PG_FUNCTION_ARGS);