Remove the local hacks for float rounding and go

with system level functions, that are hopefully
more optimized. Also, this seems to provide more
correct results than the old code!
(Closes #3852)

git-svn-id: http://svn.osgeo.org/postgis/trunk@15772 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
index 1207d6e5add1aea8dae5a0567a0c8e21abb4fcca..702283d6464238b49adee4fabb883a0341d6fedd 100644 (file)
--- a/liblwgeom/liblwgeom.h.in
+++ b/liblwgeom/liblwgeom.h.in
@@ -651,7 +651,7 @@ extern uint32_t gserialized_get_type(const GSERIALIZED *g);
 extern uint32_t gserialized_max_header_size(void);
 
 /**
-* Returns the size in bytes of the header, from the start of the 
+* Returns the size in bytes of the header, from the start of the
 * object up to the type number.
 */
 extern uint32_t gserialized_header_size(const GSERIALIZED *gser);
@@ -710,10 +710,10 @@ extern int gserialized_ndims(const GSERIALIZED *gser);
 /**
 * Return -1 if g1 is "less than" g2, 1 if g1 is "greater than"
 * g2 and 0 if g1 and g2 are the "same". Equality is evaluated
-* with a memcmp and size check. So it is possible that two 
-* identical objects where one lacks a bounding box could be 
+* with a memcmp and size check. So it is possible that two
+* identical objects where one lacks a bounding box could be
 * evaluated as non-equal initially. Greater and less than
-* are evaluated by calculating a sortable key from the center 
+* are evaluated by calculating a sortable key from the center
 * point of the object bounds.
 */
 extern int gserialized_cmp(const GSERIALIZED *g1, const GSERIALIZED *g2);
@@ -1165,8 +1165,6 @@ extern void printLWTIN(LWTIN *tin);
 
 extern float  next_float_down(double d);
 extern float  next_float_up(double d);
-extern double next_double_down(float d);
-extern double next_double_up(float d);
 
 /* general utilities 2D */
 extern double  distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2);

diff --git a/liblwgeom/lwgeom_api.c b/liblwgeom/lwgeom_api.c
index 399b169f9d96baff449c92745cf34dae6681b964..b1b888fad8b8bf8ea7a3d71ef9054d266f013b16 100644 (file)
--- a/liblwgeom/lwgeom_api.c
+++ b/liblwgeom/lwgeom_api.c
@@ -47,115 +47,15 @@ lwgeom_version()
 }
 
 
-/**********************************************************************
- * BOX routines
- *
- * returns the float thats very close to the input, but <=
- *  handles the funny differences in float4 and float8 reps.
- **********************************************************************/
-
-typedef union
-{
-       float value;
-       uint32_t word;
-} ieee_float_shape_type;
-
-#define GET_FLOAT_WORD(i,d)                    \
-       do {                                    \
-               ieee_float_shape_type gf_u;     \
-               gf_u.value = (d);               \
-               (i) = gf_u.word;                \
-       } while (0)
-
-
-#define SET_FLOAT_WORD(d,i)                    \
-       do {                                    \
-               ieee_float_shape_type sf_u;     \
-               sf_u.word = (i);                \
-               (d) = sf_u.value;               \
-       } while (0)
-
-
-/*
- * Returns the next smaller or next larger float
- * from x (in direction of y).
- */
-static float
-nextafterf_custom(float x, float y)
-{
-       int hx,hy,ix,iy;
-
-       GET_FLOAT_WORD(hx,x);
-       GET_FLOAT_WORD(hy,y);
-       ix = hx&0x7fffffff;             /* |x| */
-       iy = hy&0x7fffffff;             /* |y| */
-
-       if ((ix>0x7f800000) ||   /* x is nan */
-               (iy>0x7f800000))     /* y is nan */
-               return x+y;
-       if (x==y) return y;              /* x=y, return y */
-       if (ix==0)
-       {
-               /* x == 0 */
-               SET_FLOAT_WORD(x,(hy&0x80000000)|1);/* return +-minsubnormal */
-               y = x*x;
-               if (y==x) return y;
-               else return x;   /* raise underflow flag */
-       }
-       if (hx>=0)
-       {
-               /* x > 0 */
-               if (hx>hy)
-               {
-                       /* x > y, x -= ulp */
-                       hx -= 1;
-               }
-               else
-               {
-                       /* x < y, x += ulp */
-                       hx += 1;
-               }
-       }
-       else
-       {
-               /* x < 0 */
-               if (hy>=0||hx>hy)
-               {
-                       /* x < y, x -= ulp */
-                       hx -= 1;
-               }
-               else
-               {
-                       /* x > y, x += ulp */
-                       hx += 1;
-               }
-       }
-       hy = hx&0x7f800000;
-       if (hy>=0x7f800000) return x+x;  /* overflow  */
-       if (hy<0x00800000)
-       {
-               /* underflow */
-               y = x*x;
-               if (y!=x)
-               {
-                       /* raise underflow flag */
-                       SET_FLOAT_WORD(y,hx);
-                       return y;
-               }
-       }
-       SET_FLOAT_WORD(x,hx);
-       return x;
-}
-
-
-float next_float_down(double d)
+inline float
+next_float_down(double d)
 {
        float result  = d;
 
-       if ( ((double) result) <=d)
+       if ( ((double)result) <=d )
                return result;
 
-       return nextafterf_custom(result, result - 1000000);
+       return nextafterf(result, -1*FLT_MAX);
 
 }
 
@@ -163,47 +63,18 @@ float next_float_down(double d)
  * Returns the float thats very close to the input, but >=.
  * handles the funny differences in float4 and float8 reps.
  */
-float
+inline float
 next_float_up(double d)
 {
        float result  = d;
 
-       if ( ((double) result) >=d)
-               return result;
-
-       return nextafterf_custom(result, result + 1000000);
-}
-
-
-/*
- * Returns the double thats very close to the input, but <.
- * handles the funny differences in float4 and float8 reps.
- */
-double
-next_double_down(float d)
-{
-       double result  = d;
-
-       if ( result < d)
+       if ( ((double)result) >=d )
                return result;
 
-       return nextafterf_custom(result, result - 1000000);
+       return nextafterf(result, FLT_MAX);
 }
 
-/*
- * Returns the double thats very close to the input, but >
- * handles the funny differences in float4 and float8 reps.
- */
-double
-next_double_up(float d)
-{
-       double result  = d;
 
-       if ( result > d)
-               return result;
-
-       return nextafterf_custom(result, result + 1000000);
-}
 
 
 /************************************************************************

diff --git a/regress/lwgeom_regress_expected b/regress/lwgeom_regress_expected
index 6bc193ed0b995b4c316287b165379d2ef1bd6580..5e20f9c6984e170ad0b28e27982f48737ee7ad67 100644 (file)
--- a/regress/lwgeom_regress_expected
+++ b/regress/lwgeom_regress_expected
@@ -15,9 +15,9 @@ BOX3D(0 0.1 -55,11 12 12)
 #3069|BOX(0 0,1 1)
 #3069|BOX(1 1,1 1)
 #3069|BOX(0 0,1 1)
-BoundingDiagonal1|SRID=4326;LINESTRING(999999986991104 999999986991104,999999986991104 999999986991104)
+BoundingDiagonal1|SRID=4326;LINESTRING(999999986991104 999999986991104,1.00000005409997e+15 1.00000005409997e+15)
 BoundingDiagonal2|SRID=4326;LINESTRING(1e+15 1e+15,1e+15 1e+15)
-BoundingDiagonal3|SRID=4326;LINESTRING(999999986991104 999999986991104,999999986991104 999999986991104)
+BoundingDiagonal3|SRID=4326;LINESTRING(999999986991104 999999986991104,1.00000005409997e+15 1.00000005409997e+15)
 BoundingDiagonal4|SRID=3857;LINESTRING(-1 -2 -8 2,1 2 3 9)
 BoundingDiagonal5|SRID=3857;LINESTRINGM(4 4 0,5 4 1)
 BoundingDiagonal6|SRID=3857;LINESTRINGM EMPTY