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

I assume this is not something we are supposed to apply.

---------------------------------------------------------------------------

Paul Matthews wrote:
> 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 ----

> 
> -- 
> Sent via pgsql-hackers mailing list (pgsql-hackers@postgresql.org)
> To make changes to your subscription:
> http://www.postgresql.org/mailpref/pgsql-hackers

--  Bruce Momjian  <bruce@momjian.us>        http://momjian.us EnterpriseDB
http://enterprisedb.comPG East:  http://www.enterprisedb.com/community/nav-pg-east-2010.do + If your life is a hard
drive,Christ can be your backup. +