Re: phypot - Pygmy Hippotause ? - Mailing list pgsql-hackers

Kevin Grittner wrote:
>
> The first test seems unnecessary if you have the second.
> x >= 0, so x can't be zero unless y is, too.
> Returning x on y == 0.0 will return 0.0 whenever x == 0.0.
>
> -Kevin
>
Wish granted. :-)

--
--
Fools ignore complexity. Pragmatists suffer it.
Some can avoid it. Geniuses remove it.

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    29 Aug 2009 02:47:14 -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,5465 ----

      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;
+     }
+
+     /* If x is also 0.0, then 0.0 is returned. This is the correct result
+      * and avoids division by zero. When x > 0.0 we can trivially avoid
+      * the calculation of the hypotenuse. */
+     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    29 Aug 2009 02:47:19 -0000
***************
*** 50,56 ****
  #define FPge(A,B)                ((A) >= (B))
  #endif

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

  /*---------------------------------------------------------------------
   * Point - (x,y)
--- 50,55 ----