*** 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);