Tom Lane <[email protected]> wrote:
> I think the patch is good in principle
Since everyone seems to agree this is a good patch which needed
minor tweaks, and we haven't heard from the author, I just went
ahead and made the changes.
Andrew, could you take another look and see if you agree I've
covered it all before it gets marked ready for a committer?
-Kevin
*** a/src/backend/utils/adt/geo_ops.c
--- b/src/backend/utils/adt/geo_ops.c
***************
*** 5410,5412 **** plist_same(int npts, Point *p1, Point *p2)
--- 5410,5470 ----
return FALSE;
}
+
+
+ /*-------------------------------------------------------------------------
+ * Determine the hypotenuse.
+ *
+ * If required, x and y are swapped to make x the larger number. The
+ * traditional formulae 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 otherwise normally
+ * possible.
+ *
+ * Only arguments of > 1.27e308 are at risk of causing overflow. Whereas
+ * the naive approach limits arguments to < 9.5e153.
+ *
+ * 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
+ phypot(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 required, swap x and y */
+ if (x < y)
+ {
+ double temp = x;
+
+ x = y;
+ y = temp;
+ }
+
+ /*
+ * If y is zero, the hypotenuse is x, whether or not x is also zero.
This
+ * test protects against divide-by-zero errors.
+ */
+ if (y == 0.0)
+ return x;
+
+ /* Determine the hypotenuse */
+ yx = y / x;
+ return x * sqrt(1.0 + (yx * yx));
+ }
*** a/src/include/utils/geo_decls.h
--- b/src/include/utils/geo_decls.h
***************
*** 50,56 ****
#define FPge(A,B) ((A) >= (B))
#endif
! #define HYPOT(A, B) sqrt((A) * (A) + (B) * (B))
/*---------------------------------------------------------------------
* Point - (x,y)
--- 50,56 ----
#define FPge(A,B) ((A) >= (B))
#endif
! #define HYPOT(A, B) phypot((A), (B))
/*---------------------------------------------------------------------
* Point - (x,y)
***************
*** 211,216 **** extern Datum point_div(PG_FUNCTION_ARGS);
--- 211,217 ----
/* private routines */
extern double point_dt(Point *pt1, Point *pt2);
extern double point_sl(Point *pt1, Point *pt2);
+ extern double phypot(double x, double y);
/* public lseg routines */
extern Datum lseg_in(PG_FUNCTION_ARGS);
--
Sent via pgsql-hackers mailing list ([email protected])
To make changes to your subscription:
http://www.postgresql.org/mailpref/pgsql-hackers
