8000 Replace the naive HYPOT() macro with a standards-conformant hypotenuse · postgrespro/postgres_cluster@db04f2b · GitHub
[go: up one dir, main page]
Skip to content

Commit db04f2b

Browse files
tglsfdctglsfdc
committed
Replace the naive HYPOT() macro with a standards-conformant hypotenuse
function. This avoids unnecessary overflows and probably gives a more accurate result as well. Paul Matthews, reviewed by Andrew Geery
1 parent 90a391c commit db04f2b

File tree

2 files changed

+64
-3
lines changed

2 files changed

+64
-3
lines changed

‎src/backend/utils/adt/geo_ops.c

Lines changed: 61 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@
88
*
99
*
1010
* IDENTIFICATION
11-
* $PostgreSQL: pgsql/src/backend/utils/adt/geo_ops.c,v 1.108 2010/02/26 02:01:08 momjian Exp $
11+
* $PostgreSQL: pgsql/src/backend/utils/adt/geo_ops.c,v 1.109 2010/08/03 21:21:03 tgl Exp $
1212
*
1313
*-------------------------------------------------------------------------
1414
*/
@@ -5410,3 +5410,63 @@ plist_same(int npts, Point *p1, Point *p2)
54105410

54115411
return FALSE;
54125412
}
5413+
5414+
5415+
/*-------------------------------------------------------------------------
5416+
* Determine the hypotenuse.
5417+
*
5418+
* If required, x and y are swapped to make x the larger number. The
5419+
* traditional formula of x^2+y^2 is rearranged to factor x outside the
5420+
* sqrt. This allows computation of the hypotenuse for significantly
5421+
* larger values, and with a higher precision than when using the naive
5422+
* formula. In particular, this cannot overflow unless the final result
5423+
* would be out-of-range.
5424+
*
5425+
* sqrt( x^2 + y^2 ) = sqrt( x^2( 1 + y^2/x^2) )
5426+
* = x * sqrt( 1 + y^2/x^2 )
5427+
* = x * sqrt( 1 + y/x * y/x )
5428+
*
5429+
* It is expected that this routine will eventually be replaced with the
5430+
* C99 hypot() function.
5431+
*
5432+
* This implementation conforms to IEEE Std 1003.1 and GLIBC, in that the
5433+
* case of hypot(inf,nan) results in INF, and not NAN.
5434+
*-----------------------------------------------------------------------
5435+
*/
5436+
double
5437+
pg_hypot(double x, double y)
5438+
{
5439+
double yx;
5440+
5441+
/* Handle INF and NaN properly */
5442+
if (isinf(x) || isinf(y))
5443+
return get_float8_infinity();
5444+
5445+
if (isnan(x) || isnan(y))
5446+
return get_float8_nan();
5447+
5448+
/* Else, drop any minus signs */
5449+
x = fabs(x);
5450+
y = fabs(y);
5451+
5452+
/* Swap x and y if needed to make x the larger one */
5453+
if (x < y)
5454+
{
5455+
double temp = x;
5456+
5457+
x = y;
5458+
y = temp;
5459+
}
5460+
5461+
/*
5462+
* If y is zero, the hypotenuse is x. This test saves a few cycles in
5463+
* such cases, but more importantly it also protects against
5464+
* divide-by-zero errors, since now x >= y.
5465+
*/
5466+
if (y == 0.0)
5467+
return x;
5468+
5469+
/* Determine the hypotenuse */
5470+
yx = y / x;
5471+
return x * sqrt(1.0 + (yx * yx));
5472+
}

‎src/include/utils/geo_decls.h

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@
66
* Portions Copyright (c) 1996-2010, PostgreSQL Global Development Group
77
* Portions Copyright (c) 1994, Regents of the University of California
88
*
9-
* $PostgreSQL: pgsql/src/include/utils/geo_decls.h,v 1.57 2010/01/14 16:31:09 teodor Exp $
9+
* $PostgreSQL: pgsql/src/include/utils/geo_decls.h,v 1.58 2010/08/03 21:21:03 tgl Exp $
1010
*
1111
* NOTE
1212
* These routines do *not* use the float types from adt/.
@@ -50,7 +50,7 @@
5050
#define FPge(A,B) ((A) >= (B))
5151
#endif
5252

53-
#define HYPOT(A, B) sqrt((A) * (A) + (B) * (B))
53+
#define HYPOT(A, B) pg_hypot(A, B)
5454

5555
/*---------------------------------------------------------------------
5656
* Point - (x,y)
@@ -211,6 +211,7 @@ extern Datum point_div(PG_FUNCTION_ARGS);
211211
/* private routines */
212212
extern double point_dt(Point *pt1, Point *pt2);
213213
extern double point_sl(Point *pt1, Point *pt2);
214+
extern double pg_hypot(double x, double y);
214215

215216
/* public lseg routines */
216217
extern Datum lseg_in(PG_FUNCTION_ARGS);

0 commit comments

Comments
 (0)
0