phypot - Pygmy Hippotause ?

Started by Paul Matthewsover 16 years ago8 messages

plm@netspace.net.au

over 16 years ago

1 attachment(s)

Another attempt at replacing the current HYPOT macro with a much better
behaved function. I've added comments addressing those sections of code
that tripped people up before. As well as explaining some of the design
decisions. Feedback appreciated.

Attachments:

patchfile_phypottext/plain; name=patchfile_phypotDownload

Index: src/backend/utils/adt/geo_ops.c
===================================================================
RCS file: /projects/cvsroot/pgsql/src/backend/utils/adt/geo_ops.c,v
retrieving revision 1.103
diff -c -r1.103 geo_ops.c
*** src/backend/utils/adt/geo_ops.c	28 Jul 2009 09:47:59 -0000	1.103
--- src/backend/utils/adt/geo_ops.c	28 Aug 2009 08:11:00 -0000
***************
*** 32,37 ****
--- 32,38 ----
   * Internal routines
   */
  
+ static double phypot(double x, double y);
  static int	point_inside(Point *p, int npts, Point *plist);
  static int	lseg_crossing(double x, double y, double px, double py);
  static BOX *box_construct(double x1, double x2, double y1, double y2);
***************
*** 825,831 ****
  	box_cn(&a, box1);
  	box_cn(&b, box2);
  
! 	PG_RETURN_FLOAT8(HYPOT(a.x - b.x, a.y - b.y));
  }
  
  
--- 826,832 ----
  	box_cn(&a, box1);
  	box_cn(&b, box2);
  
! 	PG_RETURN_FLOAT8(phypot(a.x-b.x, a.y-b.y));
  }
  
  
***************
*** 1971,1977 ****
  	Point	   *pt1 = PG_GETARG_POINT_P(0);
  	Point	   *pt2 = PG_GETARG_POINT_P(1);
  
! 	PG_RETURN_FLOAT8(HYPOT(pt1->x - pt2->x, pt1->y - pt2->y));
  }
  
  double
--- 1972,1978 ----
  	Point	   *pt1 = PG_GETARG_POINT_P(0);
  	Point	   *pt2 = PG_GETARG_POINT_P(1);
  
! 	PG_RETURN_FLOAT8(phypot(pt1->x - pt2->x, pt1->y - pt2->y));
  }
  
  double
***************
*** 1979,1987 ****
  {
  #ifdef GEODEBUG
  	printf("point_dt- segment (%f,%f),(%f,%f) length is %f\n",
! 	pt1->x, pt1->y, pt2->x, pt2->y, HYPOT(pt1->x - pt2->x, pt1->y - pt2->y));
  #endif
! 	return HYPOT(pt1->x - pt2->x, pt1->y - pt2->y);
  }
  
  Datum
--- 1980,1988 ----
  {
  #ifdef GEODEBUG
  	printf("point_dt- segment (%f,%f),(%f,%f) length is %f\n",
! 	pt1->x, pt1->y, pt2->x, pt2->y, phypot(pt1->x - pt2->x, pt1->y - pt2->y));
  #endif
! 	return phypot(pt1->x - pt2->x, pt1->y - pt2->y);
  }
  
  Datum
***************
*** 2444,2450 ****
  dist_pl_internal(Point *pt, LINE *line)
  {
  	return (line->A * pt->x + line->B * pt->y + line->C) /
! 		HYPOT(line->A, line->B);
  }
  
  Datum
--- 2445,2451 ----
  dist_pl_internal(Point *pt, LINE *line)
  {
  	return (line->A * pt->x + line->B * pt->y + line->C) /
! 		phypot(line->A, line->B);
  }
  
  Datum
***************
*** 4916,4922 ****
  										   PointPGetDatum(point)));
  	result->center.x = p->x;
  	result->center.y = p->y;
! 	result->radius *= HYPOT(point->x, point->y);
  
  	PG_RETURN_CIRCLE_P(result);
  }
--- 4917,4923 ----
  										   PointPGetDatum(point)));
  	result->center.x = p->x;
  	result->center.y = p->y;
! 	result->radius *= phypot(point->x, point->y);
  
  	PG_RETURN_CIRCLE_P(result);
  }
***************
*** 4936,4942 ****
  										   PointPGetDatum(point)));
  	result->center.x = p->x;
  	result->center.y = p->y;
! 	result->radius /= HYPOT(point->x, point->y);
  
  	PG_RETURN_CIRCLE_P(result);
  }
--- 4937,4943 ----
  										   PointPGetDatum(point)));
  	result->center.x = p->x;
  	result->center.y = p->y;
! 	result->radius /= phypot(point->x, point->y);
  
  	PG_RETURN_CIRCLE_P(result);
  }
***************
*** 5401,5403 ****
--- 5402,5468 ----
  
  	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. To obtain Single 
+  * UNIX Specification behaviour, swap the order of the isnan() and isinf() 
+  * test sections.
+  *
+  * The HYPOT macro this function is replacing did not check for, or 
+  * signal on overflow. In addition no part of geo_ops checks for overflow,
+  * underflow or NaN's. This function maintains the same behaviour.
+  *-----------------------------------------------------------------------*/
+ double phypot( double x, double y )
+ {
+     double yx;
+ 
+     /* As per IEEE Std 1003.1 */
+     if( isinf(x) || isinf(y) ) 
+       return get_float8_infinity();
+ 
+     /* As per IEEE Std 1003.1 */
+     if( isnan(x) || isnan(y) )
+       return get_float8_nan();
+ 
+     /* If required, swap x and y */
+     x = fabs(x);
+     y = fabs(y);    
+     if (x < y) {
+         double temp = x;
+         x = y;
+         y = temp;
+     }
+ 
+     /* As x is the larger value, this must be the correct answer. Also
+      * avoids division by zero. */
+     if( x == 0.0 )
+         return 0.0;
+ 
+     /* Trivial case. */
+     if( y == 0.0 )
+         return x;
+     
+     /* Determine the hypotenuse */
+     yx = y/x;
+     return x*sqrt(1.0+yx*yx);
+ }
Index: src/include/utils/geo_decls.h
===================================================================
RCS file: /projects/cvsroot/pgsql/src/include/utils/geo_decls.h,v
retrieving revision 1.55
diff -c -r1.55 geo_decls.h
*** src/include/utils/geo_decls.h	1 Jan 2009 17:24:02 -0000	1.55
--- src/include/utils/geo_decls.h	28 Aug 2009 08:11:05 -0000
***************
*** 50,56 ****
  #define FPge(A,B)				((A) >= (B))
  #endif
  
- #define HYPOT(A, B)				sqrt((A) * (A) + (B) * (B))
  
  /*---------------------------------------------------------------------
   * Point - (x,y)
--- 50,55 ----