0003-geo-float-v2.patch
application/octet-stream
Filename: 0003-geo-float-v2.patch
Type: application/octet-stream
Part: 2
Message:
Re: [PATCH] Improve geometric types
Patch
Format: format-patch
Series: patch v2-0003
Subject: geo-float-v2
| File | + | − |
|---|---|---|
| src/backend/access/gist/gistproc.c | 79 | 77 |
| src/backend/utils/adt/geo_ops.c | 284 | 262 |
| src/backend/utils/adt/geo_spgist.c | 18 | 18 |
| src/include/utils/float.h | 13 | 0 |
| src/include/utils/geo_decls.h | 18 | 22 |
From e4a1939dc75fb607245643fb928db508746f7bc2 Mon Sep 17 00:00:00 2001
From: Emre Hasegeli <emre@hasegeli.com>
Date: Sat, 27 May 2017 16:17:42 -0400
Subject: [PATCH 3/4] geo-float-v2
Use the built-in float datatype to implement geometric types
This will provide:
* Check for underflow and overflow
* Check for division by zero
* Handle NaNs consistently
The patch also replaces all occurrences of "double" as "float8". They
are the same, but were randomly spread on the same file.
---
src/backend/access/gist/gistproc.c | 156 +++++------
src/backend/utils/adt/geo_ops.c | 546 +++++++++++++++++++------------------
src/backend/utils/adt/geo_spgist.c | 36 +--
src/include/utils/float.h | 13 +
src/include/utils/geo_decls.h | 40 ++-
5 files changed, 412 insertions(+), 379 deletions(-)
diff --git a/src/backend/access/gist/gistproc.c b/src/backend/access/gist/gistproc.c
index c1ea88f90d..3dbd4aaa36 100644
--- a/src/backend/access/gist/gistproc.c
+++ b/src/backend/access/gist/gistproc.c
@@ -16,20 +16,21 @@
*-------------------------------------------------------------------------
*/
#include "postgres.h"
#include <float.h>
#include <math.h>
#include "access/gist.h"
#include "access/stratnum.h"
#include "utils/builtins.h"
+#include "utils/float.h"
#include "utils/geo_decls.h"
static bool gist_box_leaf_consistent(BOX *key, BOX *query,
StrategyNumber strategy);
static bool rtree_internal_consistent(BOX *key, BOX *query,
StrategyNumber strategy);
/* Minimum accepted ratio of split */
#define LIMIT_RATIO 0.3
@@ -48,55 +49,56 @@ rt_box_union(BOX *n, const BOX *a, const BOX *b)
n->high.x = float8_max(a->high.x, b->high.x);
n->high.y = float8_max(a->high.y, b->high.y);
n->low.x = float8_min(a->low.x, b->low.x);
n->low.y = float8_min(a->low.y, b->low.y);
}
/*
* Size of a BOX for penalty-calculation purposes.
* The result can be +Infinity, but not NaN.
*/
-static double
+static float8
size_box(const BOX *box)
{
/*
* Check for zero-width cases. Note that we define the size of a zero-
* by-infinity box as zero. It's important to special-case this somehow,
* as naively multiplying infinity by zero will produce NaN.
*
* The less-than cases should not happen, but if they do, say "zero".
*/
if (float8_le(box->high.x, box->low.x) ||
float8_le(box->high.y, box->low.y))
return 0.0;
/*
* We treat NaN as larger than +Infinity, so any distance involving a NaN
* and a non-NaN is infinite. Note the previous check eliminated the
* possibility that the low fields are NaNs.
*/
if (isnan(box->high.x) || isnan(box->high.y))
return get_float8_infinity();
- return (box->high.x - box->low.x) * (box->high.y - box->low.y);
+ return float8_mul(float8_mi(box->high.x, box->low.x),
+ float8_mi(box->high.y, box->low.y));
}
/*
* Return amount by which the union of the two boxes is larger than
* the original BOX's area. The result can be +Infinity, but not NaN.
*/
-static double
+static float8
box_penalty(const BOX *original, const BOX *new)
{
BOX unionbox;
rt_box_union(&unionbox, original, new);
- return size_box(&unionbox) - size_box(original);
+ return float8_mi(size_box(&unionbox), size_box(original));
}
/*
* The GiST Consistent method for boxes
*
* Should return false if for all data items x below entry,
* the predicate x op query must be FALSE, where op is the oper
* corresponding to strategy in the pg_amop table.
*/
Datum
@@ -284,74 +286,74 @@ fallbackSplit(GistEntryVector *entryvec, GIST_SPLITVEC *v)
/*
* Represents information about an entry that can be placed to either group
* without affecting overlap over selected axis ("common entry").
*/
typedef struct
{
/* Index of entry in the initial array */
int index;
/* Delta between penalties of entry insertion into different groups */
- double delta;
+ float8 delta;
} CommonEntry;
/*
* Context for g_box_consider_split. Contains information about currently
* selected split and some general information.
*/
typedef struct
{
int entriesCount; /* total number of entries being split */
BOX boundingBox; /* minimum bounding box across all entries */
/* Information about currently selected split follows */
bool first; /* true if no split was selected yet */
- double leftUpper; /* upper bound of left interval */
- double rightLower; /* lower bound of right interval */
+ float8 leftUpper; /* upper bound of left interval */
+ float8 rightLower; /* lower bound of right interval */
float4 ratio;
float4 overlap;
int dim; /* axis of this split */
- double range; /* width of general MBR projection to the
+ float8 range; /* width of general MBR projection to the
* selected axis */
} ConsiderSplitContext;
/*
* Interval represents projection of box to axis.
*/
typedef struct
{
- double lower,
+ float8 lower,
upper;
} SplitInterval;
/*
* Interval comparison function by lower bound of the interval;
*/
static int
interval_cmp_lower(const void *i1, const void *i2)
{
- double lower1 = ((const SplitInterval *) i1)->lower,
+ float8 lower1 = ((const SplitInterval *) i1)->lower,
lower2 = ((const SplitInterval *) i2)->lower;
return float8_cmp_internal(lower1, lower2);
}
/*
* Interval comparison function by upper bound of the interval;
*/
static int
interval_cmp_upper(const void *i1, const void *i2)
{
- double upper1 = ((const SplitInterval *) i1)->upper,
+ float8 upper1 = ((const SplitInterval *) i1)->upper,
upper2 = ((const SplitInterval *) i2)->upper;
return float8_cmp_internal(upper1, upper2);
}
/*
* Replace negative (or NaN) value with zero.
*/
static inline float
non_negative(float val)
@@ -360,28 +362,28 @@ non_negative(float val)
return val;
else
return 0.0f;
}
/*
* Consider replacement of currently selected split with the better one.
*/
static inline void
g_box_consider_split(ConsiderSplitContext *context, int dimNum,
- double rightLower, int minLeftCount,
- double leftUpper, int maxLeftCount)
+ float8 rightLower, int minLeftCount,
+ float8 leftUpper, int maxLeftCount)
{
int leftCount,
rightCount;
float4 ratio,
overlap;
- double range;
+ float8 range;
/*
* Calculate entries distribution ratio assuming most uniform distribution
* of common entries.
*/
if (minLeftCount >= (context->entriesCount + 1) / 2)
{
leftCount = minLeftCount;
}
else
@@ -390,52 +392,54 @@ g_box_consider_split(ConsiderSplitContext *context, int dimNum,
leftCount = maxLeftCount;
else
leftCount = context->entriesCount / 2;
}
rightCount = context->entriesCount - leftCount;
/*
* Ratio of split - quotient between size of lesser group and total
* entries count.
*/
- ratio = ((float4) Min(leftCount, rightCount)) /
- ((float4) context->entriesCount);
+ ratio = float4_div(Min(leftCount, rightCount), context->entriesCount);
- if (ratio > LIMIT_RATIO)
+ if (float4_gt(ratio, LIMIT_RATIO))
{
bool selectthis = false;
/*
* The ratio is acceptable, so compare current split with previously
* selected one. Between splits of one dimension we search for minimal
* overlap (allowing negative values) and minimal ration (between same
* overlaps. We switch dimension if find less overlap (non-negative)
* or less range with same overlap.
*/
if (dimNum == 0)
- range = context->boundingBox.high.x - context->boundingBox.low.x;
+ range = float8_mi(context->boundingBox.high.x,
+ context->boundingBox.low.x);
else
- range = context->boundingBox.high.y - context->boundingBox.low.y;
+ range = float8_mi(context->boundingBox.high.y,
+ context->boundingBox.low.y);
- overlap = (leftUpper - rightLower) / range;
+ overlap = float8_div(float8_mi(leftUpper, rightLower), range);
/* If there is no previous selection, select this */
if (context->first)
selectthis = true;
else if (context->dim == dimNum)
{
/*
* Within the same dimension, choose the new split if it has a
* smaller overlap, or same overlap but better ratio.
*/
- if (overlap < context->overlap ||
- (overlap == context->overlap && ratio > context->ratio))
+ if (float4_lt(overlap, context->overlap) ||
+ (float4_eq(overlap, context->overlap) &&
+ float4_gt(ratio, context->ratio)))
selectthis = true;
}
else
{
/*
* Across dimensions, choose the new split if it has a smaller
* *non-negative* overlap, or same *non-negative* overlap but
* bigger range. This condition differs from the one described in
* the article. On the datasets where leaf MBRs don't overlap
* themselves, non-overlapping splits (i.e. splits which have zero
@@ -444,55 +448,49 @@ g_box_consider_split(ConsiderSplitContext *context, int dimNum,
* non-overlapping splits (i.e. having lowest negative overlap)
* appears to be in the same dimension as in the previous split.
* Therefore MBRs appear to be very prolonged along another
* dimension, which leads to bad search performance. Using range
* as the second split criteria makes MBRs more quadratic. Using
* *non-negative* overlap instead of overlap as the first split
* criteria gives to range criteria a chance to matter, because
* non-overlapping splits are equivalent in this criteria.
*/
if (non_negative(overlap) < non_negative(context->overlap) ||
- (range > context->range &&
+ (float4_gt(range, context->range) &&
non_negative(overlap) <= non_negative(context->overlap)))
selectthis = true;
}
if (selectthis)
{
/* save information about selected split */
context->first = false;
context->ratio = ratio;
context->range = range;
context->overlap = overlap;
context->rightLower = rightLower;
context->leftUpper = leftUpper;
context->dim = dimNum;
}
}
}
/*
* Compare common entries by their deltas.
- * (We assume the deltas can't be NaN.)
*/
static int
common_entry_cmp(const void *i1, const void *i2)
{
- double delta1 = ((const CommonEntry *) i1)->delta,
+ float8 delta1 = ((const CommonEntry *) i1)->delta,
delta2 = ((const CommonEntry *) i2)->delta;
- if (delta1 < delta2)
- return -1;
- else if (delta1 > delta2)
- return 1;
- else
- return 0;
+ return float8_cmp_internal(delta1, delta2);
}
/*
* --------------------------------------------------------------------------
* Double sorting split algorithm. This is used for both boxes and points.
*
* The algorithm finds split of boxes by considering splits along each axis.
* Each entry is first projected as an interval on the X-axis, and different
* ways to split the intervals into two groups are considered, trying to
* minimize the overlap of the groups. Then the same is repeated for the
@@ -552,21 +550,21 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
else
adjustBox(&context.boundingBox, box);
}
/*
* Iterate over axes for optimal split searching.
*/
context.first = true; /* nothing selected yet */
for (dim = 0; dim < 2; dim++)
{
- double leftUpper,
+ float8 leftUpper,
rightLower;
int i1,
i2;
/* Project each entry as an interval on the selected axis. */
for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
{
box = DatumGetBoxP(entryvec->vector[i].key);
if (dim == 0)
{
@@ -749,21 +747,21 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
*rightBox = *(box); \
v->spl_right[v->spl_nright++] = off; \
} while(0)
/*
* Distribute entries which can be distributed unambiguously, and collect
* common entries.
*/
for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
{
- double lower,
+ float8 lower,
upper;
/*
* Get upper and lower bounds along selected axis.
*/
box = DatumGetBoxP(entryvec->vector[i].key);
if (context.dim == 0)
{
lower = box->low.x;
upper = box->high.x;
@@ -804,31 +802,31 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
/*
* Distribute "common entries", if any.
*/
if (commonEntriesCount > 0)
{
/*
* Calculate minimum number of entries that must be placed in both
* groups, to reach LIMIT_RATIO.
*/
- int m = ceil(LIMIT_RATIO * (double) nentries);
+ int m = ceil(LIMIT_RATIO * (float8) nentries);
/*
* Calculate delta between penalties of join "common entries" to
* different groups.
*/
for (i = 0; i < commonEntriesCount; i++)
{
box = DatumGetBoxP(entryvec->vector[commonEntries[i].index].key);
- commonEntries[i].delta = Abs(box_penalty(leftBox, box) -
- box_penalty(rightBox, box));
+ commonEntries[i].delta = Abs(float8_mi(box_penalty(leftBox, box),
+ box_penalty(rightBox, box)));
}
/*
* Sort "common entries" by calculated deltas in order to distribute
* the most ambiguous entries first.
*/
qsort(commonEntries, commonEntriesCount, sizeof(CommonEntry), common_entry_cmp);
/*
* Distribute "common entries" between groups.
@@ -1128,24 +1126,24 @@ gist_circle_compress(PG_FUNCTION_ARGS)
{
GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
GISTENTRY *retval;
if (entry->leafkey)
{
CIRCLE *in = DatumGetCircleP(entry->key);
BOX *r;
r = (BOX *) palloc(sizeof(BOX));
- r->high.x = in->center.x + in->radius;
- r->low.x = in->center.x - in->radius;
- r->high.y = in->center.y + in->radius;
- r->low.y = in->center.y - in->radius;
+ r->high.x = float8_pl(in->center.x, in->radius);
+ r->low.x = float8_mi(in->center.x, in->radius);
+ r->high.y = float8_pl(in->center.y, in->radius);
+ r->low.y = float8_mi(in->center.y, in->radius);
retval = (GISTENTRY *) palloc(sizeof(GISTENTRY));
gistentryinit(*retval, PointerGetDatum(r),
entry->rel, entry->page,
entry->offset, FALSE);
}
else
retval = entry;
PG_RETURN_POINTER(retval);
}
@@ -1169,24 +1167,24 @@ gist_circle_consistent(PG_FUNCTION_ARGS)
*recheck = true;
if (DatumGetBoxP(entry->key) == NULL || query == NULL)
PG_RETURN_BOOL(FALSE);
/*
* Since the operators require recheck anyway, we can just use
* rtree_internal_consistent even at leaf nodes. (This works in part
* because the index entries are bounding boxes not circles.)
*/
- bbox.high.x = query->center.x + query->radius;
- bbox.low.x = query->center.x - query->radius;
- bbox.high.y = query->center.y + query->radius;
- bbox.low.y = query->center.y - query->radius;
+ bbox.high.x = float8_pl(query->center.x, query->radius);
+ bbox.low.x = float8_mi(query->center.x, query->radius);
+ bbox.high.y = float8_pl(query->center.y, query->radius);
+ bbox.low.y = float8_mi(query->center.y, query->radius);
result = rtree_internal_consistent(DatumGetBoxP(entry->key),
&bbox, strategy);
PG_RETURN_BOOL(result);
}
/**************************************************
* Point ops
**************************************************/
@@ -1237,80 +1235,84 @@ gist_point_fetch(PG_FUNCTION_ARGS)
entry->offset, FALSE);
PG_RETURN_POINTER(retval);
}
#define point_point_distance(p1,p2) \
DatumGetFloat8(DirectFunctionCall2(point_distance, \
PointPGetDatum(p1), PointPGetDatum(p2)))
-static double
+static float8
computeDistance(bool isLeaf, BOX *box, Point *point)
{
- double result = 0.0;
+ float8 result = 0.0;
if (isLeaf)
{
/* simple point to point distance */
result = point_point_distance(point, &box->low);
}
- else if (point->x <= box->high.x && point->x >= box->low.x &&
- point->y <= box->high.y && point->y >= box->low.y)
+ else if (float8_le(point->x, box->high.x) &&
+ float8_ge(point->x, box->low.x) &&
+ float8_le(point->y, box->high.y) &&
+ float8_ge(point->y, box->low.y))
{
/* point inside the box */
result = 0.0;
}
- else if (point->x <= box->high.x && point->x >= box->low.x)
+ else if (float8_le(point->x, box->high.x) &&
+ float8_ge(point->x, box->low.x))
{
/* point is over or below box */
- Assert(box->low.y <= box->high.y);
- if (point->y > box->high.y)
- result = point->y - box->high.y;
- else if (point->y < box->low.y)
- result = box->low.y - point->y;
+ Assert(float8_le(box->low.y, box->high.y));
+ if (float8_gt(point->y, box->high.y))
+ result = float8_mi(point->y, box->high.y);
+ else if (float8_lt(point->y, box->low.y))
+ result = float8_mi(box->low.y, point->y);
else
elog(ERROR, "inconsistent point values");
}
- else if (point->y <= box->high.y && point->y >= box->low.y)
+ else if (float8_le(point->y, box->high.y) &&
+ float8_ge(point->y, box->low.y))
{
/* point is to left or right of box */
- Assert(box->low.x <= box->high.x);
- if (point->x > box->high.x)
- result = point->x - box->high.x;
- else if (point->x < box->low.x)
- result = box->low.x - point->x;
+ Assert(float8_lt(box->low.x, box->high.x));
+ if (float8_gt(point->x, box->high.x))
+ result = float8_mi(point->x, box->high.x);
+ else if (float8_lt(point->x, box->low.x))
+ result = float8_mi(box->low.x, point->x);
else
elog(ERROR, "inconsistent point values");
}
else
{
/* closest point will be a vertex */
Point p;
- double subresult;
+ float8 subresult;
result = point_point_distance(point, &box->low);
subresult = point_point_distance(point, &box->high);
- if (result > subresult)
+ if (float8_gt(result, subresult))
result = subresult;
p.x = box->low.x;
p.y = box->high.y;
subresult = point_point_distance(point, &p);
- if (result > subresult)
+ if (float8_gt(result, subresult))
result = subresult;
p.x = box->high.x;
p.y = box->low.y;
subresult = point_point_distance(point, &p);
- if (result > subresult)
+ if (float8_gt(result, subresult))
result = subresult;
}
return result;
}
static bool
gist_point_consistent_internal(StrategyNumber strategy,
bool isLeaf, BOX *key, Point *query)
{
@@ -1391,24 +1393,24 @@ gist_point_consistent(PG_FUNCTION_ARGS)
* Instead we write a non-fuzzy overlap test. The same code
* will also serve for leaf-page tests, since leaf keys have
* high == low.
*/
BOX *query,
*key;
query = PG_GETARG_BOX_P(1);
key = DatumGetBoxP(entry->key);
- result = (key->high.x >= query->low.x &&
- key->low.x <= query->high.x &&
- key->high.y >= query->low.y &&
- key->low.y <= query->high.y);
+ result = (float8_ge(key->high.x, query->low.x) &&
+ float8_le(key->low.x, query->high.x) &&
+ float8_ge(key->high.y, query->low.y) &&
+ float8_le(key->low.y, query->high.y));
*recheck = false;
}
break;
case PolygonStrategyNumberGroup:
{
POLYGON *query = PG_GETARG_POLYGON_P(1);
result = DatumGetBool(DirectFunctionCall5(
gist_poly_consistent,
PointerGetDatum(entry),
@@ -1417,22 +1419,22 @@ gist_point_consistent(PG_FUNCTION_ARGS)
0, PointerGetDatum(recheck)));
if (GIST_LEAF(entry) && result)
{
/*
* We are on leaf page and quick check shows overlapping
* of polygon's bounding box and point
*/
BOX *box = DatumGetBoxP(entry->key);
- Assert(box->high.x == box->low.x
- && box->high.y == box->low.y);
+ Assert(float8_eq(box->high.x, box->low.x) &&
+ float8_eq(box->high.y, box->low.y));
result = DatumGetBool(DirectFunctionCall2(
poly_contain_pt,
PolygonPGetDatum(query),
PointPGetDatum(&box->high)));
*recheck = false;
}
}
break;
case CircleStrategyNumberGroup:
{
@@ -1446,22 +1448,22 @@ gist_point_consistent(PG_FUNCTION_ARGS)
0, PointerGetDatum(recheck)));
if (GIST_LEAF(entry) && result)
{
/*
* We are on leaf page and quick check shows overlapping
* of polygon's bounding box and point
*/
BOX *box = DatumGetBoxP(entry->key);
- Assert(box->high.x == box->low.x
- && box->high.y == box->low.y);
+ Assert(float8_eq(box->high.x, box->low.x) &&
+ float8_eq(box->high.y, box->low.y));
result = DatumGetBool(DirectFunctionCall2(
circle_contain_pt,
CirclePGetDatum(query),
PointPGetDatum(&box->high)));
*recheck = false;
}
}
break;
default:
elog(ERROR, "unrecognized strategy number: %d", strategy);
@@ -1470,21 +1472,21 @@ gist_point_consistent(PG_FUNCTION_ARGS)
}
PG_RETURN_BOOL(result);
}
Datum
gist_point_distance(PG_FUNCTION_ARGS)
{
GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
- double distance;
+ float8 distance;
StrategyNumber strategyGroup = strategy / GeoStrategyNumberOffset;
switch (strategyGroup)
{
case PointStrategyNumberGroup:
distance = computeDistance(GIST_LEAF(entry),
DatumGetBoxP(entry->key),
PG_GETARG_POINT_P(1));
break;
default:
@@ -1499,25 +1501,25 @@ gist_point_distance(PG_FUNCTION_ARGS)
/*
* The inexact GiST distance method for geometric types that store bounding
* boxes.
*
* Compute lossy distance from point to index entries. The result is inexact
* because index entries are bounding boxes, not the exact shapes of the
* indexed geometric types. We use distance from point to MBR of index entry.
* This is a lower bound estimate of distance from point to indexed geometric
* type.
*/
-static double
+static float8
gist_bbox_distance(GISTENTRY *entry, Datum query,
StrategyNumber strategy, bool *recheck)
{
- double distance;
+ float8 distance;
StrategyNumber strategyGroup = strategy / GeoStrategyNumberOffset;
/* Bounding box distance is always inexact. */
*recheck = true;
switch (strategyGroup)
{
case PointStrategyNumberGroup:
distance = computeDistance(false,
DatumGetBoxP(entry->key),
@@ -1533,32 +1535,32 @@ gist_bbox_distance(GISTENTRY *entry, Datum query,
Datum
gist_circle_distance(PG_FUNCTION_ARGS)
{
GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
Datum query = PG_GETARG_DATUM(1);
StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
/* Oid subtype = PG_GETARG_OID(3); */
bool *recheck = (bool *) PG_GETARG_POINTER(4);
- double distance;
+ float8 distance;
distance = gist_bbox_distance(entry, query, strategy, recheck);
PG_RETURN_FLOAT8(distance);
}
Datum
gist_poly_distance(PG_FUNCTION_ARGS)
{
GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
Datum query = PG_GETARG_DATUM(1);
StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
/* Oid subtype = PG_GETARG_OID(3); */
bool *recheck = (bool *) PG_GETARG_POINTER(4);
- double distance;
+ float8 distance;
distance = gist_bbox_distance(entry, query, strategy, recheck);
PG_RETURN_FLOAT8(distance);
}
diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index 865698738e..dbff7abf64 100644
--- a/src/backend/utils/adt/geo_ops.c
+++ b/src/backend/utils/adt/geo_ops.c
@@ -41,55 +41,55 @@ static void point_sub_internal(Point *result, Point *pt1, Point *pt2);
static void point_mul_internal(Point *result, Point *pt1, Point *pt2);
static void point_div_internal(Point *result, Point *pt1, Point *pt2);
static bool point_eq_internal(Point *pt1, Point *pt2);
static float8 point_dt(Point *pt1, Point *pt2);
static float8 point_sl(Point *pt1, Point *pt2);
static int point_inside(Point *p, int npts, Point *plist);
/* Methods on box */
static void box_construct(BOX *result, float8 x1, float8 x2, float8 y1, float8 y2);
static void box_cn(Point *center, BOX *box);
static bool box_ov(BOX *box1, BOX *box2);
-static double box_ar(BOX *box);
-static double box_ht(BOX *box);
-static double box_wd(BOX *box);
+static float8 box_ar(BOX *box);
+static float8 box_ht(BOX *box);
+static float8 box_wd(BOX *box);
/* Methods on circle */
-static double circle_ar(CIRCLE *circle);
+static float8 circle_ar(CIRCLE *circle);
/* Methods on line */
static void line_construct_pm(LINE *result, Point *pt, float8 m);
static void line_construct_pts(LINE *line, Point *pt1, Point *pt2);
static bool line_interpt_internal(Point *result, LINE *l1, LINE *l2);
static float8 line_calculate_point(LINE *line, Point *pt);
/* Methods on line segment */
static bool lseg_interpt_line_internal(Point *result, LSEG *lseg, LINE *line);
static bool lseg_interpt_internal(Point *result, LSEG *l1, LSEG *l2);
static bool lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start);
static float8 lseg_dt(LSEG *l1, LSEG *l2);
-static int lseg_crossing(double x, double y, double px, double py);
+static int lseg_crossing(float8 x, float8 y, float8 px, float8 py);
/* Others */
static bool on_ps_internal(Point *pt, LSEG *lseg);
static void make_bound_box(POLYGON *poly);
static bool plist_same(int npts, Point *p1, Point *p2);
-static double single_decode(char *num, char **endptr_p,
+static float8 single_decode(char *num, char **endptr_p,
const char *type_name, const char *orig_string);
static void single_encode(float8 x, StringInfo str);
-static void pair_decode(char *str, double *x, double *y, char **endptr_p,
+static void pair_decode(char *str, float8 *x, float8 *y, char **endptr_p,
const char *type_name, const char *orig_string);
static void pair_encode(float8 x, float8 y, StringInfo str);
static int pair_count(char *s, char delim);
static void path_decode(char *str, bool opentype, int npts, Point *p,
bool *isopen, char **endptr_p,
const char *type_name, const char *orig_string);
static char *path_encode(enum path_delim path_delim, int npts, Point *pt);
static void statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2);
-static double dist_pl_internal(Point *pt, LINE *line);
-static double dist_ps_internal(Point *pt, LSEG *lseg);
-static double dist_ppoly_internal(Point *pt, POLYGON *poly);
+static float8 dist_pl_internal(Point *pt, LINE *line);
+static float8 dist_ps_internal(Point *pt, LSEG *lseg);
+static float8 dist_ppoly_internal(Point *pt, POLYGON *poly);
/*
* Delimiters for input and output strings.
* LDELIM, RDELIM, and DELIM are left, right, and separator delimiters, respectively.
* LDELIM_EP, RDELIM_EP are left and right delimiters for paths with endpoints.
*/
#define LDELIM '('
#define RDELIM ')'
@@ -118,38 +118,38 @@ static double dist_ppoly_internal(Point *pt, POLYGON *poly);
*
* For boxes, the points are opposite corners with the first point at the top right.
* For closed paths and polygons, the points should be reordered to allow
* fast and correct equality comparisons.
*
* XXX perhaps points in complex shapes should be reordered internally
* to allow faster internal operations, but should keep track of input order
* and restore that order for text output - tgl 97/01/16
*/
-static double
+static float8
single_decode(char *num, char **endptr_p,
const char *type_name, const char *orig_string)
{
return float8in_internal(num, endptr_p, type_name, orig_string);
} /* single_decode() */
static void
single_encode(float8 x, StringInfo str)
{
char *xstr = float8out_internal(x);
appendStringInfoString(str, xstr);
pfree(xstr);
} /* single_encode() */
static void
-pair_decode(char *str, double *x, double *y, char **endptr_p,
+pair_decode(char *str, float8 *x, float8 *y, char **endptr_p,
const char *type_name, const char *orig_string)
{
bool has_delim;
while (isspace((unsigned char) *str))
str++;
if ((has_delim = (*str == LDELIM)))
str++;
*x = float8in_internal(str, &str, type_name, orig_string);
@@ -351,33 +351,33 @@ pair_count(char *s, char delim)
* External format: (two corners of box)
* "(f8, f8), (f8, f8)"
* also supports the older style "(f8, f8, f8, f8)"
*/
Datum
box_in(PG_FUNCTION_ARGS)
{
char *str = PG_GETARG_CSTRING(0);
BOX *box = (BOX *) palloc(sizeof(BOX));
bool isopen;
- double x,
+ float8 x,
y;
path_decode(str, false, 2, &(box->high), &isopen, NULL, "box", str);
/* reorder corners if necessary... */
- if (box->high.x < box->low.x)
+ if (float8_lt(box->high.x, box->low.x))
{
x = box->high.x;
box->high.x = box->low.x;
box->low.x = x;
}
- if (box->high.y < box->low.y)
+ if (float8_lt(box->high.y, box->low.y))
{
y = box->high.y;
box->high.y = box->low.y;
box->low.y = y;
}
PG_RETURN_BOX_P(box);
}
/* box_out - convert a box to external form.
@@ -391,38 +391,38 @@ box_out(PG_FUNCTION_ARGS)
}
/*
* box_recv - converts external binary format to box
*/
Datum
box_recv(PG_FUNCTION_ARGS)
{
StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
BOX *box;
- double x,
+ float8 x,
y;
box = (BOX *) palloc(sizeof(BOX));
box->high.x = pq_getmsgfloat8(buf);
box->high.y = pq_getmsgfloat8(buf);
box->low.x = pq_getmsgfloat8(buf);
box->low.y = pq_getmsgfloat8(buf);
/* reorder corners if necessary... */
- if (box->high.x < box->low.x)
+ if (float8_lt(box->high.x, box->low.x))
{
x = box->high.x;
box->high.x = box->low.x;
box->low.x = x;
}
- if (box->high.y < box->low.y)
+ if (float8_lt(box->high.y, box->low.y))
{
y = box->high.y;
box->high.y = box->low.y;
box->low.y = y;
}
PG_RETURN_BOX_P(box);
}
/*
@@ -441,31 +441,31 @@ box_send(PG_FUNCTION_ARGS)
pq_sendfloat8(&buf, box->low.y);
PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
}
/* box_construct - fill in a new box.
*/
static void
box_construct(BOX *result, float8 x1, float8 x2, float8 y1, float8 y2)
{
- if (x1 > x2)
+ if (float8_gt(x1, x2))
{
result->high.x = x1;
result->low.x = x2;
}
else
{
result->high.x = x2;
result->low.x = x1;
}
- if (y1 > y2)
+ if (float8_gt(y1, y2))
{
result->high.y = y1;
result->low.y = y2;
}
else
{
result->high.y = y2;
result->low.y = y1;
}
}
@@ -777,54 +777,54 @@ box_center(PG_FUNCTION_ARGS)
Point *result = (Point *) palloc(sizeof(Point));
box_cn(result, box);
PG_RETURN_POINT_P(result);
}
/* box_ar - returns the area of the box.
*/
-static double
+static float8
box_ar(BOX *box)
{
- return box_wd(box) * box_ht(box);
+ return float8_mul(box_wd(box), box_ht(box));
}
/* box_cn - stores the centerpoint of the box into *center.
*/
static void
box_cn(Point *center, BOX *box)
{
- center->x = (box->high.x + box->low.x) / 2.0;
- center->y = (box->high.y + box->low.y) / 2.0;
+ center->x = float8_div(float8_pl(box->high.x, box->low.x), 2.0);
+ center->y = float8_div(float8_pl(box->high.y, box->low.y), 2.0);
}
/* box_wd - returns the width (length) of the box
* (horizontal magnitude).
*/
-static double
+static float8
box_wd(BOX *box)
{
- return box->high.x - box->low.x;
+ return float8_mi(box->high.x, box->low.x);
}
/* box_ht - returns the height of the box
* (vertical magnitude).
*/
-static double
+static float8
box_ht(BOX *box)
{
- return box->high.y - box->low.y;
+ return float8_mi(box->high.y, box->low.y);
}
/*----------------------------------------------------------
* Funky operations.
*---------------------------------------------------------*/
/* box_intersect -
* returns the overlapping portion of two boxes,
* or NULL if they do not intersect.
@@ -834,24 +834,24 @@ box_intersect(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
BOX *result;
if (!box_ov(box1, box2))
PG_RETURN_NULL();
result = (BOX *) palloc(sizeof(BOX));
- result->high.x = Min(box1->high.x, box2->high.x);
- result->low.x = Max(box1->low.x, box2->low.x);
- result->high.y = Min(box1->high.y, box2->high.y);
- result->low.y = Max(box1->low.y, box2->low.y);
+ result->high.x = float8_min(box1->high.x, box2->high.x);
+ result->low.x = float8_max(box1->low.x, box2->low.x);
+ result->high.y = float8_min(box1->high.y, box2->high.y);
+ result->low.y = float8_max(box1->low.y, box2->low.y);
PG_RETURN_BOX_P(result);
}
/* box_diagonal -
* returns a line segment which happens to be the
* positive-slope diagonal of "box".
*/
Datum
@@ -982,65 +982,65 @@ line_send(PG_FUNCTION_ARGS)
/* line_construct_pm()
* point-slope
*/
static void
line_construct_pm(LINE *result, Point *pt, float8 m)
{
if (m == DBL_MAX)
{
/* vertical - use "x = C" */
- result->A = -1;
- result->B = 0;
+ result->A = -1.0;
+ result->B = 0.0;
result->C = pt->x;
}
else
{
/* use "mx - y + yinter = 0" */
result->A = m;
result->B = -1.0;
- result->C = pt->y - m * pt->x;
+ result->C = float8_mi(pt->y, float8_mul(m, pt->x));
}
}
/*
* Fill already-allocated LINE struct from two points on the line
*/
static void
line_construct_pts(LINE *line, Point *pt1, Point *pt2)
{
if (FPeq(pt1->x, pt2->x))
{ /* vertical */
/* use "x = C" */
- line->A = -1;
- line->B = 0;
+ line->A = -1.0;
+ line->B = 0.0;
line->C = pt1->x;
#ifdef GEODEBUG
printf("line_construct_pts- line is vertical\n");
#endif
}
else if (FPeq(pt1->y, pt2->y))
{ /* horizontal */
/* use "y = C" */
- line->A = 0;
- line->B = -1;
+ line->A = 0.0;
+ line->B = -1.0;
line->C = pt1->y;
#ifdef GEODEBUG
printf("line_construct_pts- line is horizontal\n");
#endif
}
else
{
/* use "mx - y + yinter = 0" */
line->A = point_sl(pt2, pt1);
line->B = -1.0;
- line->C = pt1->y - line->A * pt1->x;
+ line->C = float8_mi(pt1->y, float8_mul(line->A, pt1->x));
/* on some platforms, the preceding expression tends to produce -0 */
if (line->C == 0.0)
line->C = 0.0;
#ifdef GEODEBUG
printf("line_construct_pts- line is neither vertical nor horizontal (diffs x=%.*g, y=%.*g\n",
DBL_DIG, (pt2->x - pt1->x), DBL_DIG, (pt2->y - pt1->y));
#endif
}
}
@@ -1060,21 +1060,22 @@ line_construct_pp(PG_FUNCTION_ARGS)
/*
* Calculate the line equation for a point
*
* This returns the result of the line equation Ax + By + C. The result
* needs to be 0 for the point to be on the line.
*/
static float8
line_calculate_point(LINE *line, Point *pt)
{
- return line->A * pt->x + line->B * pt->y + line->C;
+ return float8_pl(float8_pl(float8_mul(line->A, pt->x),
+ float8_mul(line->B, pt->y)), line->C);
}
/*----------------------------------------------------------
* Relative position routines.
*---------------------------------------------------------*/
Datum
line_intersect(PG_FUNCTION_ARGS)
{
@@ -1097,21 +1098,22 @@ Datum
line_perp(PG_FUNCTION_ARGS)
{
LINE *l1 = PG_GETARG_LINE_P(0);
LINE *l2 = PG_GETARG_LINE_P(1);
if (FPzero(l1->A))
PG_RETURN_BOOL(FPzero(l2->B));
else if (FPzero(l1->B))
PG_RETURN_BOOL(FPzero(l2->A));
- PG_RETURN_BOOL(FPeq(((l1->A * l2->B) / (l1->B * l2->A)), -1.0));
+ PG_RETURN_BOOL(FPeq(float8_div(float8_mul(l1->A, l2->B),
+ float8_mul(l1->B, l2->A)), -1.0));
}
Datum
line_vertical(PG_FUNCTION_ARGS)
{
LINE *line = PG_GETARG_LINE_P(0);
PG_RETURN_BOOL(FPzero(line->B));
}
@@ -1121,34 +1123,34 @@ line_horizontal(PG_FUNCTION_ARGS)
LINE *line = PG_GETARG_LINE_P(0);
PG_RETURN_BOOL(FPzero(line->A));
}
Datum
line_eq(PG_FUNCTION_ARGS)
{
LINE *l1 = PG_GETARG_LINE_P(0);
LINE *l2 = PG_GETARG_LINE_P(1);
- double k;
+ float8 ratio;
if (!FPzero(l2->A))
- k = l1->A / l2->A;
+ ratio = float8_div(l1->A, l2->A);
else if (!FPzero(l2->B))
- k = l1->B / l2->B;
+ ratio = float8_div(l1->B, l2->B);
else if (!FPzero(l2->C))
- k = l1->C / l2->C;
+ ratio = float8_div(l1->C, l2->C);
else
- k = 1.0;
+ ratio = 1.0;
- PG_RETURN_BOOL(FPeq(l1->A, k * l2->A) &&
- FPeq(l1->B, k * l2->B) &&
- FPeq(l1->C, k * l2->C));
+ PG_RETURN_BOOL(FPeq(l1->A, float8_mul(ratio, l2->A)) &&
+ FPeq(l1->B, float8_mul(ratio, l2->B)) &&
+ FPeq(l1->C, float8_mul(ratio, l2->C)));
}
/*----------------------------------------------------------
* Line arithmetic routines.
*---------------------------------------------------------*/
/* line_distance()
* Distance between two lines.
*/
@@ -1156,21 +1158,21 @@ Datum
line_distance(PG_FUNCTION_ARGS)
{
LINE *l1 = PG_GETARG_LINE_P(0);
LINE *l2 = PG_GETARG_LINE_P(1);
float8 result;
Point tmp;
if (line_interpt_internal(NULL, l1, l2))
PG_RETURN_FLOAT8(0.0);
if (FPzero(l1->B)) /* vertical? */
- PG_RETURN_FLOAT8(fabs(l1->C - l2->C));
+ PG_RETURN_FLOAT8(fabs(float8_mi(l1->C, l2->C)));
point_construct(&tmp, 0.0, l1->C);
result = dist_pl_internal(&tmp, l2);
PG_RETURN_FLOAT8(result);
}
/* line_interpt()
* Point where two lines l1, l2 intersect (if any)
*/
Datum
line_interpt(PG_FUNCTION_ARGS)
@@ -1193,41 +1195,41 @@ line_interpt(PG_FUNCTION_ARGS)
* parallel), false if they do not. This also sets the intersection point
* to *result, if it is not NULL.
*
* NOTE: If the lines are identical then we will find they are parallel
* and report "no intersection". This is a little weird, but since
* there's no *unique* intersection, maybe it's appropriate behavior.
*/
static bool
line_interpt_internal(Point *result, LINE *l1, LINE *l2)
{
- double x,
+ float8 x,
y;
if (FPzero(l1->B)) /* l1 vertical? */
{
if (FPzero(l2->B)) /* l2 vertical? */
return false;
x = l1->C;
- y = (l2->A * x + l2->C);
+ y = float8_pl(float8_mul(l2->A, x), l2->C);
}
else
{
- if (FPeq(l2->A, l1->A * (l2->B / l1->B)))
+ if (FPeq(l2->A, float8_mul(l1->A, float8_div(l2->B, l1->B))))
return false;
if (FPzero(l2->B))
x = l2->C;
else
- x = (l1->C - l2->C) / (l2->A - l1->A);
- y = (l1->A * x + l1->C);
+ x = float8_div(float8_mi(l1->C, l2->C), float8_mi(l2->A, l1->A));
+ y = float8_pl(float8_mul(l1->A, x), l1->C);
}
if (result != NULL)
point_construct(result, x, y);
#ifdef GEODEBUG
printf("line_interpt- lines are A=%.*g, B=%.*g, C=%.*g, A=%.*g, B=%.*g, C=%.*g\n",
DBL_DIG, l1->A, DBL_DIG, l1->B, DBL_DIG, l1->C, DBL_DIG, l2->A, DBL_DIG, l2->B, DBL_DIG, l2->C);
printf("line_interpt- lines intersect at (%.*g,%.*g)\n", DBL_DIG, x, DBL_DIG, y);
#endif
@@ -1254,36 +1256,35 @@ line_interpt_internal(Point *result, LINE *l1, LINE *l2)
* "(xcoord, ycoord),... "
* "[xcoord, ycoord,... ]"
* Also support older format:
* "(closed, npts, xcoord, ycoord,... )"
*---------------------------------------------------------*/
Datum
path_area(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P(0);
- double area = 0.0;
+ float8 area = 0.0;
int i,
j;
if (!path->closed)
PG_RETURN_NULL();
for (i = 0; i < path->npts; i++)
{
j = (i + 1) % path->npts;
- area += path->p[i].x * path->p[j].y;
- area -= path->p[i].y * path->p[j].x;
+ area = float8_pl(area, float8_mul(path->p[i].x, path->p[j].y));
+ area = float8_mi(area, float8_mul(path->p[i].y, path->p[j].x));
}
- area *= 0.5;
- PG_RETURN_FLOAT8(area < 0.0 ? -area : area);
+ PG_RETURN_FLOAT8(float8_div(fabs(area), 2.0));
}
Datum
path_in(PG_FUNCTION_ARGS)
{
char *str = PG_GETARG_CSTRING(0);
PATH *path;
bool isopen;
char *s;
@@ -1500,31 +1501,31 @@ path_npoints(PG_FUNCTION_ARGS)
PG_RETURN_INT32(path->npts);
}
Datum
path_close(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P_COPY(0);
- path->closed = TRUE;
+ path->closed = true;
PG_RETURN_PATH_P(path);
}
Datum
path_open(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P_COPY(0);
- path->closed = FALSE;
+ path->closed = false;
PG_RETURN_PATH_P(path);
}
/* path_inter -
* Does p1 intersect p2 at any point?
* Use bounding boxes for a quick (O(n)) check, then do a
* O(n^2) iterative edge check.
*/
@@ -1540,33 +1541,33 @@ path_inter(PG_FUNCTION_ARGS)
LSEG seg1,
seg2;
if (p1->npts <= 0 || p2->npts <= 0)
PG_RETURN_BOOL(false);
b1.high.x = b1.low.x = p1->p[0].x;
b1.high.y = b1.low.y = p1->p[0].y;
for (i = 1; i < p1->npts; i++)
{
- b1.high.x = Max(p1->p[i].x, b1.high.x);
- b1.high.y = Max(p1->p[i].y, b1.high.y);
- b1.low.x = Min(p1->p[i].x, b1.low.x);
- b1.low.y = Min(p1->p[i].y, b1.low.y);
+ b1.high.x = float8_max(p1->p[i].x, b1.high.x);
+ b1.high.y = float8_max(p1->p[i].y, b1.high.y);
+ b1.low.x = float8_min(p1->p[i].x, b1.low.x);
+ b1.low.y = float8_min(p1->p[i].y, b1.low.y);
}
b2.high.x = b2.low.x = p2->p[0].x;
b2.high.y = b2.low.y = p2->p[0].y;
for (i = 1; i < p2->npts; i++)
{
- b2.high.x = Max(p2->p[i].x, b2.high.x);
- b2.high.y = Max(p2->p[i].y, b2.high.y);
- b2.low.x = Min(p2->p[i].x, b2.low.x);
- b2.low.y = Min(p2->p[i].y, b2.low.y);
+ b2.high.x = float8_max(p2->p[i].x, b2.high.x);
+ b2.high.y = float8_max(p2->p[i].y, b2.high.y);
+ b2.low.x = float8_min(p2->p[i].x, b2.low.x);
+ b2.low.y = float8_min(p2->p[i].y, b2.low.y);
}
if (!box_ov(&b1, &b2))
PG_RETURN_BOOL(false);
/* pairwise check lseg intersections */
for (i = 0; i < p1->npts; i++)
{
int iprev;
if (i > 0)
@@ -1642,21 +1643,21 @@ path_distance(PG_FUNCTION_ARGS)
{
if (!p2->closed)
continue;
jprev = p2->npts - 1; /* include the closure segment */
}
statlseg_construct(&seg1, &p1->p[iprev], &p1->p[i]);
statlseg_construct(&seg2, &p2->p[jprev], &p2->p[j]);
tmp = lseg_dt(&seg1, &seg2);
- if (!have_min || tmp < min)
+ if (!have_min || float8_lt(tmp, min))
{
min = tmp;
have_min = true;
}
}
}
if (!have_min)
PG_RETURN_NULL();
@@ -1681,21 +1682,21 @@ path_length(PG_FUNCTION_ARGS)
if (i > 0)
iprev = i - 1;
else
{
if (!path->closed)
continue;
iprev = path->npts - 1; /* include the closure segment */
}
- result += point_dt(&path->p[iprev], &path->p[i]);
+ result = float8_pl(result, point_dt(&path->p[iprev], &path->p[i]));
}
PG_RETURN_FLOAT8(result);
}
/***********************************************************************
**
** Routines for 2D points.
**
***********************************************************************/
@@ -1867,43 +1868,48 @@ point_distance(PG_FUNCTION_ARGS)
{
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
PG_RETURN_FLOAT8(point_dt(pt1, pt2));
}
static float8
point_dt(Point *pt1, Point *pt2)
{
+ float8 result;
+
+ result = float8_hypot(float8_mi(pt1->x, pt2->x), float8_mi(pt1->y, pt2->y));
+
#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));
+ pt1->x, pt1->y, pt2->x, pt2->y, result);
#endif
- return HYPOT(pt1->x - pt2->x, pt1->y - pt2->y);
+
+ return result;
}
Datum
point_slope(PG_FUNCTION_ARGS)
{
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
PG_RETURN_FLOAT8(point_sl(pt1, pt2));
}
static float8
point_sl(Point *pt1, Point *pt2)
{
- return (FPeq(pt1->x, pt2->x)
- ? (double) DBL_MAX
- : (pt1->y - pt2->y) / (pt1->x - pt2->x));
+ if (FPeq(pt1->x, pt2->x))
+ return DBL_MAX;
+ return float8_div(float8_mi(pt1->y, pt2->y), float8_mi(pt1->x, pt2->x));
}
/***********************************************************************
**
** Routines for 2D line segments.
**
***********************************************************************/
/*----------------------------------------------------------
@@ -2043,35 +2049,35 @@ lseg_parallel(PG_FUNCTION_ARGS)
* '((0,0),(0,1))'::lseg ?-| '((0,0),(1,0))'::lseg
* So, modified it to check explicitly for slope of vertical line
* returned by point_sl() and the results seem better.
* - thomas 1998-01-31
*/
Datum
lseg_perp(PG_FUNCTION_ARGS)
{
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
- double m1,
+ float8 m1,
m2;
m1 = point_sl(&l1->p[0], &l1->p[1]);
m2 = point_sl(&l2->p[0], &l2->p[1]);
#ifdef GEODEBUG
printf("lseg_perp- slopes are %g and %g\n", m1, m2);
#endif
if (FPzero(m1))
PG_RETURN_BOOL(FPeq(m2, DBL_MAX));
else if (FPzero(m2))
PG_RETURN_BOOL(FPeq(m1, DBL_MAX));
- PG_RETURN_BOOL(FPeq(m1 / m2, -1.0));
+ PG_RETURN_BOOL(FPeq(float8_div(m1, m2), -1.0));
}
Datum
lseg_vertical(PG_FUNCTION_ARGS)
{
LSEG *lseg = PG_GETARG_LSEG_P(0);
PG_RETURN_BOOL(FPeq(lseg->p[0].x, lseg->p[1].x));
}
@@ -2161,52 +2167,47 @@ lseg_distance(PG_FUNCTION_ARGS)
LSEG *l2 = PG_GETARG_LSEG_P(1);
PG_RETURN_FLOAT8(lseg_dt(l1, l2));
}
/* lseg_dt()
* Distance between two line segments.
* Must check both sets of endpoints to ensure minimum distance is found.
* - thomas 1998-02-01
*/
-static double
+static float8
lseg_dt(LSEG *l1, LSEG *l2)
{
- double result,
- d;
+ float8 result;
if (lseg_interpt_internal(NULL, l1, l2))
return 0.0;
- d = dist_ps_internal(&l1->p[0], l2);
- result = d;
- d = dist_ps_internal(&l1->p[1], l2);
- result = Min(result, d);
- d = dist_ps_internal(&l2->p[0], l1);
- result = Min(result, d);
- d = dist_ps_internal(&l2->p[1], l1);
- result = Min(result, d);
+ result = dist_ps_internal(&l1->p[0], l2);
+ result = float8_min(result, dist_ps_internal(&l1->p[1], l2));
+ result = float8_min(result, dist_ps_internal(&l2->p[0], l1));
+ result = float8_min(result, dist_ps_internal(&l2->p[1], l1));
return result;
}
Datum
lseg_center(PG_FUNCTION_ARGS)
{
LSEG *lseg = PG_GETARG_LSEG_P(0);
Point *result;
result = (Point *) palloc(sizeof(Point));
- result->x = (lseg->p[0].x + lseg->p[1].x) / 2.0;
- result->y = (lseg->p[0].y + lseg->p[1].y) / 2.0;
+ result->x = float8_div(float8_pl(lseg->p[0].x, lseg->p[1].x), 2.0);
+ result->y = float8_div(float8_pl(lseg->p[0].y, lseg->p[1].y), 2.0);
PG_RETURN_POINT_P(result);
}
static bool
lseg_interpt_internal(Point *result, LSEG *l1, LSEG *l2)
{
Point interpt;
LINE tmp1,
tmp2;
@@ -2290,58 +2291,58 @@ lseg_interpt(PG_FUNCTION_ARGS)
*/
Datum
dist_pl(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
LINE *line = PG_GETARG_LINE_P(1);
PG_RETURN_FLOAT8(dist_pl_internal(pt, line));
}
-static double
+static float8
dist_pl_internal(Point *pt, LINE *line)
{
- return fabs(line_calculate_point(line, pt) /
- HYPOT(line->A, line->B));
+ return float8_div(fabs(line_calculate_point(line, pt)),
+ float8_hypot(line->A, line->B));
}
/*
* Distance from a point to a lseg
*/
Datum
dist_ps(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
LSEG *lseg = PG_GETARG_LSEG_P(1);
PG_RETURN_FLOAT8(dist_ps_internal(pt, lseg));
}
-static double
+static float8
dist_ps_internal(Point *pt, LSEG *lseg)
{
- double m; /* slope of perp. */
- double result,
+ float8 m; /* slope of perp. */
+ float8 result,
tmpdist;
Point interpt;
LINE ln;
/*
* Construct a line perpendicular to the input segment and through the
* input point
*/
- if (lseg->p[1].x == lseg->p[0].x)
- m = 0;
- else if (lseg->p[1].y == lseg->p[0].y)
- m = (double) DBL_MAX; /* slope is infinite */
+ if (float8_eq(lseg->p[0].x, lseg->p[1].x))
+ m = 0.0;
+ else if (float8_eq(lseg->p[0].y, lseg->p[1].y))
+ m = DBL_MAX; /* slope is infinite */
else
- m = -1.0 / point_sl(&lseg->p[0], &lseg->p[1]);
+ m = float8_div(-1.0, point_sl(&lseg->p[0], &lseg->p[1]));
line_construct_pm(&ln, pt, m);
#ifdef GEODEBUG
printf("dist_ps- line is A=%g B=%g C=%g from (point) slope (%f,%f) %g\n",
ln->A, ln->B, ln->C, pt->x, pt->y, m);
#endif
/*
* Calculate distance to the line segment or to the nearest endpoint of
* the segment.
@@ -2353,24 +2354,22 @@ dist_ps_internal(Point *pt, LSEG *lseg)
/* yes, so use distance to the intersection point */
result = point_dt(pt, &interpt);
#ifdef GEODEBUG
printf("dist_ps- distance is %f to intersection point is (%f,%f)\n",
result, tmp.x, tmp.y);
#endif
}
else
{
/* no, so use distance to the nearer endpoint */
- result = point_dt(pt, &lseg->p[0]);
- tmpdist = point_dt(pt, &lseg->p[1]);
- if (tmpdist < result)
- result = tmpdist;
+ result = float8_min(point_dt(pt, &lseg->p[0]),
+ point_dt(pt, &lseg->p[1]));
}
return result;
}
/*
* Distance from a point to a path
*/
Datum
dist_ppath(PG_FUNCTION_ARGS)
@@ -2408,21 +2407,21 @@ dist_ppath(PG_FUNCTION_ARGS)
iprev = i - 1;
else
{
if (!path->closed)
continue;
iprev = path->npts - 1; /* include the closure segment */
}
statlseg_construct(&lseg, &path->p[iprev], &path->p[i]);
tmp = dist_ps_internal(pt, &lseg);
- if (!have_min || tmp < result)
+ if (!have_min || float8_lt(tmp, result))
{
result = tmp;
have_min = true;
}
}
break;
}
PG_RETURN_FLOAT8(result);
}
@@ -2446,33 +2445,28 @@ dist_pb(PG_FUNCTION_ARGS)
}
/*
* Distance from a lseg to a line
*/
Datum
dist_sl(PG_FUNCTION_ARGS)
{
LSEG *lseg = PG_GETARG_LSEG_P(0);
LINE *line = PG_GETARG_LINE_P(1);
- float8 result,
- d2;
+ float8 result;
if (lseg_interpt_line_internal(NULL, lseg, line))
result = 0.0;
else
- {
- result = dist_pl_internal(&lseg->p[0], line);
- d2 = dist_pl_internal(&lseg->p[1], line);
/* XXX shouldn't we take the min not max? */
- if (d2 > result)
- result = d2;
- }
+ result = float8_max(dist_pl_internal(&lseg->p[0], line),
+ dist_pl_internal(&lseg->p[1], line));
PG_RETURN_FLOAT8(result);
}
/*
* Distance from a lseg to a box
*/
Datum
dist_sb(PG_FUNCTION_ARGS)
{
@@ -2514,25 +2508,22 @@ dist_lb(PG_FUNCTION_ARGS)
* Distance from a circle to a polygon
*/
Datum
dist_cpoly(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
POLYGON *poly = PG_GETARG_POLYGON_P(1);
float8 result;
/* calculate distance to center, and subtract radius */
- result = dist_ppoly_internal(&circle->center, poly);
-
- result -= circle->radius;
- if (result < 0)
- result = 0;
+ result = float8_min(float8_mi(dist_ppoly_internal(&circle->center, poly),
+ circle->radius), 0.0);
PG_RETURN_FLOAT8(result);
}
/*
* Distance from a point to a polygon
*/
Datum
dist_ppoly(PG_FUNCTION_ARGS)
{
@@ -2550,21 +2541,21 @@ dist_polyp(PG_FUNCTION_ARGS)
{
POLYGON *poly = PG_GETARG_POLYGON_P(0);
Point *point = PG_GETARG_POINT_P(1);
float8 result;
result = dist_ppoly_internal(point, poly);
PG_RETURN_FLOAT8(result);
}
-static double
+static float8
dist_ppoly_internal(Point *pt, POLYGON *poly)
{
float8 result;
float8 d;
int i;
LSEG seg;
if (point_inside(pt, poly->npts, poly->p) != 0)
{
#ifdef GEODEBUG
@@ -2587,21 +2578,21 @@ dist_ppoly_internal(Point *pt, POLYGON *poly)
for (i = 0; (i < poly->npts - 1); i++)
{
seg.p[0].x = poly->p[i].x;
seg.p[0].y = poly->p[i].y;
seg.p[1].x = poly->p[i + 1].x;
seg.p[1].y = poly->p[i + 1].y;
d = dist_ps_internal(pt, &seg);
#ifdef GEODEBUG
printf("dist_ppoly_internal- segment %d distance is %f\n", (i + 1), d);
#endif
- if (d < result)
+ if (float8_lt(d, result))
result = d;
}
return result;
}
/*---------------------------------------------------------------------
* interpt_
* Intersection point of objects.
@@ -2655,41 +2646,41 @@ lseg_interpt_line_internal(Point *result, LSEG *lseg, LINE *line)
/* close_pl -
* The intersection point of a perpendicular of the line
* through the point.
*/
Datum
close_pl(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
LINE *line = PG_GETARG_LINE_P(1);
Point *result;
- double invm;
+ float8 invm;
LINE tmp;
result = (Point *) palloc(sizeof(Point));
if (FPzero(line->B)) /* vertical? */
{
result->x = line->C;
result->y = pt->y;
PG_RETURN_POINT_P(result);
}
if (FPzero(line->A)) /* horizontal? */
{
result->x = pt->x;
result->y = line->C;
PG_RETURN_POINT_P(result);
}
/* drop a perpendicular and find the intersection point */
/* invert and flip the sign on the slope to get a perpendicular */
- invm = line->B / line->A;
+ invm = float8_div(line->B, line->A);
line_construct_pm(&tmp, pt, invm);
line_interpt_internal(result, &tmp, line);
PG_RETURN_POINT_P(result);
}
/* close_ps()
* Closest point on line segment to specified point.
* Take the closest endpoint if the point is left, right,
* above, or below the segment, otherwise find the intersection
@@ -2698,21 +2689,21 @@ close_pl(PG_FUNCTION_ARGS)
* Some tricky code here, relying on boolean expressions
* evaluating to only zero or one to use as an array index.
* bug fixes by gthaker@atl.lmco.com; May 1, 1998
*/
Datum
close_ps(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
LSEG *lseg = PG_GETARG_LSEG_P(1);
Point *result;
- double invm;
+ float8 invm;
int xh,
yh;
LINE tmp;
result = (Point *) palloc(sizeof(Point));
#ifdef GEODEBUG
printf("close_sp:pt->x %f pt->y %f\nlseg(0).x %f lseg(0).y %f lseg(1).x %f lseg(1).y %f\n",
pt->x, pt->y, lseg->p[0].x, lseg->p[0].y,
lseg->p[1].x, lseg->p[1].y);
@@ -2754,35 +2745,35 @@ close_ps(PG_FUNCTION_ARGS)
point_construct(result, pt->x, lseg->p[0].y);
PG_RETURN_POINT_P(result);
}
/*
* vert. and horiz. cases are down, now check if the closest point is one
* of the end points or someplace on the lseg.
*/
- invm = -1.0 / point_sl(&(lseg->p[0]), &(lseg->p[1]));
+ invm = float8_div(-1.0, point_sl(&lseg->p[0], &lseg->p[1]));
line_construct_pm(&tmp, &lseg->p[!yh], invm); /* lower edge of the
* "band" */
- if (pt->y < (tmp.A * pt->x + tmp.C))
+ if (pt->y < float8_pl(float8_mul(tmp.A, pt->x), tmp.C))
{ /* we are below the lower edge */
*result = lseg->p[!yh]; /* below the lseg, take lower end pt */
#ifdef GEODEBUG
printf("close_ps below: tmp A %f B %f C %f\n",
tmp->A, tmp->B, tmp->C);
#endif
PG_RETURN_POINT_P(result);
}
line_construct_pm(&tmp, &lseg->p[yh], invm); /* upper edge of the
* "band" */
- if (pt->y > (tmp.A * pt->x + tmp.C))
+ if (pt->y > float8_pl(float8_mul(tmp.A, pt->x), tmp.C))
{ /* we are below the lower edge */
*result = lseg->p[yh]; /* above the lseg, take higher end pt */
#ifdef GEODEBUG
printf("close_ps above: tmp A %f B %f C %f\n",
tmp->A, tmp->B, tmp->C);
#endif
PG_RETURN_POINT_P(result);
}
/*
@@ -2818,32 +2809,32 @@ close_lseg(PG_FUNCTION_ARGS)
{
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
Point *result;
double dist,
dist0,
dist1;
dist0 = dist_ps_internal(&l1->p[0], l2);
dist1 = dist_ps_internal(&l1->p[1], l2);
- dist = Min(dist0, dist1);
+ dist = float8_min(dist0, dist1);
- if (dist_ps_internal(&l2->p[0], l1) < dist)
+ if (float8_lt(dist_ps_internal(&l2->p[0], l1), dist))
{
result = DatumGetPointP(DirectFunctionCall2(close_ps,
PointPGetDatum(&l2->p[0]),
LsegPGetDatum(l1)));
result = DatumGetPointP(DirectFunctionCall2(close_ps,
PointPGetDatum(result),
LsegPGetDatum(l2)));
}
- else if (dist_ps_internal(&l2->p[1], l1) < dist)
+ else if (float8_lt(dist_ps_internal(&l2->p[1], l1), dist))
{
result = DatumGetPointP(DirectFunctionCall2(close_ps,
PointPGetDatum(&l2->p[1]),
LsegPGetDatum(l1)));
result = DatumGetPointP(DirectFunctionCall2(close_ps,
PointPGetDatum(result),
LsegPGetDatum(l2)));
}
else
{
@@ -2858,52 +2849,55 @@ close_lseg(PG_FUNCTION_ARGS)
* Closest point on or in box to specified point.
*/
Datum
close_pb(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
BOX *box = PG_GETARG_BOX_P(1);
LSEG lseg,
seg;
Point point;
- double dist,
+ float8 dist,
d;
if (DatumGetBool(DirectFunctionCall2(on_pb,
PointPGetDatum(pt),
BoxPGetDatum(box))))
PG_RETURN_POINT_P(pt);
/* pairwise check lseg distances */
point.x = box->low.x;
point.y = box->high.y;
statlseg_construct(&lseg, &box->low, &point);
dist = dist_ps_internal(pt, &lseg);
statlseg_construct(&seg, &box->high, &point);
- if ((d = dist_ps_internal(pt, &seg)) < dist)
+ d = dist_ps_internal(pt, &seg);
+ if (float8_lt(d, dist))
{
dist = d;
memcpy(&lseg, &seg, sizeof(lseg));
}
point.x = box->high.x;
point.y = box->low.y;
statlseg_construct(&seg, &box->low, &point);
- if ((d = dist_ps_internal(pt, &seg)) < dist)
+ d = dist_ps_internal(pt, &seg);
+ if (float8_lt(d, dist))
{
dist = d;
memcpy(&lseg, &seg, sizeof(lseg));
}
statlseg_construct(&seg, &box->high, &point);
- if ((d = dist_ps_internal(pt, &seg)) < dist)
+ d = dist_ps_internal(pt, &seg);
+ if (float8_lt(d, dist))
{
dist = d;
memcpy(&lseg, &seg, sizeof(lseg));
}
PG_RETURN_DATUM(DirectFunctionCall2(close_ps,
PointPGetDatum(pt),
LsegPGetDatum(&lseg)));
}
@@ -2926,21 +2920,21 @@ close_sl(PG_FUNCTION_ARGS)
float8 d1,
d2;
result = (Point *) palloc(sizeof(Point));
if (lseg_interpt_line_internal(result, lseg, line))
PG_RETURN_POINT_P(result);
d1 = dist_pl_internal(&lseg->p[0], line);
d2 = dist_pl_internal(&lseg->p[1], line);
- if (d1 < d2)
+ if (float8_lt(d1, d2))
*result = lseg->p[0];
else
*result = lseg->p[1];
PG_RETURN_POINT_P(result);
#endif
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("function \"close_sl\" not implemented")));
@@ -2960,77 +2954,80 @@ close_ls(PG_FUNCTION_ARGS)
float8 d1,
d2;
result = (Point *) palloc(sizeof(Point));
if (lseg_interpt_line_internal(result, lseg, line))
PG_RETURN_POINT_P(result);
d1 = dist_pl_internal(&lseg->p[0], line);
d2 = dist_pl_internal(&lseg->p[1], line);
- if (d1 < d2)
+ if (float8_lt(d1, d2))
*result = lseg->p[0];
else
*result = lseg->p[1];
PG_RETURN_POINT_P(result);
}
/* close_sb()
* Closest point on or in box to line segment.
*/
Datum
close_sb(PG_FUNCTION_ARGS)
{
LSEG *lseg = PG_GETARG_LSEG_P(0);
BOX *box = PG_GETARG_BOX_P(1);
Point point;
LSEG bseg,
seg;
- double dist,
+ float8 dist,
d;
/* segment intersects box? then just return closest point to center */
if (DatumGetBool(DirectFunctionCall2(inter_sb,
LsegPGetDatum(lseg),
BoxPGetDatum(box))))
{
box_cn(&point, box);
PG_RETURN_DATUM(DirectFunctionCall2(close_ps,
PointPGetDatum(&point),
LsegPGetDatum(lseg)));
}
/* pairwise check lseg distances */
point.x = box->low.x;
point.y = box->high.y;
statlseg_construct(&bseg, &box->low, &point);
dist = lseg_dt(lseg, &bseg);
statlseg_construct(&seg, &box->high, &point);
- if ((d = lseg_dt(lseg, &seg)) < dist)
+ d = lseg_dt(lseg, &seg);
+ if (float8_lt(d, dist))
{
dist = d;
memcpy(&bseg, &seg, sizeof(bseg));
}
point.x = box->high.x;
point.y = box->low.y;
statlseg_construct(&seg, &box->low, &point);
- if ((d = lseg_dt(lseg, &seg)) < dist)
+ d = lseg_dt(lseg, &seg);
+ if (float8_lt(d, dist))
{
dist = d;
memcpy(&bseg, &seg, sizeof(bseg));
}
statlseg_construct(&seg, &box->high, &point);
- if ((d = lseg_dt(lseg, &seg)) < dist)
+ d = lseg_dt(lseg, &seg);
+ if (float8_lt(d, dist))
{
dist = d;
memcpy(&bseg, &seg, sizeof(bseg));
}
/* OK, we now have the closest line segment on the box boundary */
PG_RETURN_DATUM(DirectFunctionCall2(close_lseg,
LsegPGetDatum(lseg),
LsegPGetDatum(&bseg)));
}
@@ -3078,75 +3075,80 @@ on_ps(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
LSEG *lseg = PG_GETARG_LSEG_P(1);
PG_RETURN_BOOL(on_ps_internal(pt, lseg));
}
static bool
on_ps_internal(Point *pt, LSEG *lseg)
{
- return FPeq(point_dt(pt, &lseg->p[0]) + point_dt(pt, &lseg->p[1]),
+ return FPeq(float8_pl(point_dt(pt, &lseg->p[0]),
+ point_dt(pt, &lseg->p[1])),
point_dt(&lseg->p[0], &lseg->p[1]));
}
+
Datum
on_pb(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
BOX *box = PG_GETARG_BOX_P(1);
- PG_RETURN_BOOL(pt->x <= box->high.x && pt->x >= box->low.x &&
- pt->y <= box->high.y && pt->y >= box->low.y);
+ PG_RETURN_BOOL(float8_le(pt->x, box->high.x) &&
+ float8_ge(pt->x, box->low.x) &&
+ float8_le(pt->y, box->high.y) &&
+ float8_ge(pt->y, box->low.y));
}
Datum
box_contain_pt(PG_FUNCTION_ARGS)
{
BOX *box = PG_GETARG_BOX_P(0);
Point *pt = PG_GETARG_POINT_P(1);
- PG_RETURN_BOOL(pt->x <= box->high.x && pt->x >= box->low.x &&
- pt->y <= box->high.y && pt->y >= box->low.y);
+ PG_RETURN_BOOL(float8_le(pt->x, box->high.x) &&
+ float8_ge(pt->x, box->low.x) &&
+ float8_le(pt->y, box->high.y) &&
+ float8_ge(pt->y, box->low.y));
}
/* on_ppath -
* Whether a point lies within (on) a polyline.
* If open, we have to (groan) check each segment.
* (uses same algorithm as for point intersecting segment - tgl 1997-07-09)
* If closed, we use the old O(n) ray method for point-in-polygon.
* The ray is horizontal, from pt out to the right.
* Each segment that crosses the ray counts as an
* intersection; note that an endpoint or edge may touch
* but not cross.
* (we can do p-in-p in lg(n), but it takes preprocessing)
*/
Datum
on_ppath(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
PATH *path = PG_GETARG_PATH_P(1);
int i,
n;
- double a,
+ float8 a,
b;
/*-- OPEN --*/
if (!path->closed)
{
n = path->npts - 1;
a = point_dt(pt, &path->p[0]);
for (i = 0; i < n; i++)
{
b = point_dt(pt, &path->p[i + 1]);
- if (FPeq(a + b,
- point_dt(&path->p[i], &path->p[i + 1])))
+ if (FPeq(float8_pl(a, b), point_dt(&path->p[i], &path->p[i + 1])))
PG_RETURN_BOOL(true);
a = b;
}
PG_RETURN_BOOL(false);
}
/*-- CLOSED --*/
PG_RETURN_BOOL(point_inside(pt, path->npts, path->p) != 0);
}
@@ -3204,24 +3206,24 @@ inter_sl(PG_FUNCTION_ARGS)
*/
Datum
inter_sb(PG_FUNCTION_ARGS)
{
LSEG *lseg = PG_GETARG_LSEG_P(0);
BOX *box = PG_GETARG_BOX_P(1);
BOX lbox;
LSEG bseg;
Point point;
- lbox.low.x = Min(lseg->p[0].x, lseg->p[1].x);
- lbox.low.y = Min(lseg->p[0].y, lseg->p[1].y);
- lbox.high.x = Max(lseg->p[0].x, lseg->p[1].x);
- lbox.high.y = Max(lseg->p[0].y, lseg->p[1].y);
+ lbox.low.x = float8_min(lseg->p[0].x, lseg->p[1].x);
+ lbox.low.y = float8_min(lseg->p[0].y, lseg->p[1].y);
+ lbox.high.x = float8_max(lseg->p[0].x, lseg->p[1].x);
+ lbox.high.y = float8_max(lseg->p[0].y, lseg->p[1].y);
/* nothing close to overlap? then not going to intersect */
if (!box_ov(&lbox, box))
PG_RETURN_BOOL(false);
/* an endpoint of segment is inside box? then clearly intersects */
if (DatumGetBool(DirectFunctionCall2(on_pb,
PointPGetDatum(&lseg->p[0]),
BoxPGetDatum(box))) ||
DatumGetBool(DirectFunctionCall2(on_pb,
@@ -3302,38 +3304,38 @@ inter_lb(PG_FUNCTION_ARGS)
* make_bound_box - create the bounding box for the input polygon
*------------------------------------------------------------------*/
/*---------------------------------------------------------------------
* Make the smallest bounding box for the given polygon.
*---------------------------------------------------------------------*/
static void
make_bound_box(POLYGON *poly)
{
int i;
- double x1,
+ float8 x1,
y1,
x2,
y2;
if (poly->npts > 0)
{
x2 = x1 = poly->p[0].x;
y2 = y1 = poly->p[0].y;
for (i = 1; i < poly->npts; i++)
{
- if (poly->p[i].x < x1)
+ if (float8_lt(poly->p[i].x, x1))
x1 = poly->p[i].x;
- if (poly->p[i].x > x2)
+ if (float8_gt(poly->p[i].x, x2))
x2 = poly->p[i].x;
- if (poly->p[i].y < y1)
+ if (float8_lt(poly->p[i].y, y1))
y1 = poly->p[i].y;
- if (poly->p[i].y > y2)
+ if (float8_gt(poly->p[i].y, y2))
y2 = poly->p[i].y;
}
box_construct(&poly->boundbox, x1, x2, y1, y2);
}
else
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("cannot create bounding box for empty polygon")));
}
@@ -3461,21 +3463,21 @@ poly_send(PG_FUNCTION_ARGS)
* the right most point of A left of the left most point
* of B?
*-------------------------------------------------------*/
Datum
poly_left(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
- result = polya->boundbox.high.x < polyb->boundbox.low.x;
+ result = float8_lt(polya->boundbox.high.x, polyb->boundbox.low.x);
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
@@ -3484,21 +3486,21 @@ poly_left(PG_FUNCTION_ARGS)
* the right most point of A at or left of the right most point
* of B?
*-------------------------------------------------------*/
Datum
poly_overleft(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
- result = polya->boundbox.high.x <= polyb->boundbox.high.x;
+ result = float8_le(polya->boundbox.high.x, polyb->boundbox.high.x);
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
@@ -3507,21 +3509,21 @@ poly_overleft(PG_FUNCTION_ARGS)
* the left most point of A right of the right most point
* of B?
*-------------------------------------------------------*/
Datum
poly_right(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
- result = polya->boundbox.low.x > polyb->boundbox.high.x;
+ result = float8_gt(polya->boundbox.low.x, polyb->boundbox.high.x);
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
@@ -3530,21 +3532,21 @@ poly_right(PG_FUNCTION_ARGS)
* the left most point of A at or right of the left most point
* of B?
*-------------------------------------------------------*/
Datum
poly_overright(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
- result = polya->boundbox.low.x >= polyb->boundbox.low.x;
+ result = float8_ge(polya->boundbox.low.x, polyb->boundbox.low.x);
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
@@ -3553,21 +3555,21 @@ poly_overright(PG_FUNCTION_ARGS)
* the upper most point of A below the lower most point
* of B?
*-------------------------------------------------------*/
Datum
poly_below(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
- result = polya->boundbox.high.y < polyb->boundbox.low.y;
+ result = float8_lt(polya->boundbox.high.y, polyb->boundbox.low.y);
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
@@ -3576,21 +3578,21 @@ poly_below(PG_FUNCTION_ARGS)
* the upper most point of A at or below the upper most point
* of B?
*-------------------------------------------------------*/
Datum
poly_overbelow(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
- result = polya->boundbox.high.y <= polyb->boundbox.high.y;
+ result = float8_le(polya->boundbox.high.y, polyb->boundbox.high.y);
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
@@ -3599,21 +3601,21 @@ poly_overbelow(PG_FUNCTION_ARGS)
* the lower most point of A above the upper most point
* of B?
*-------------------------------------------------------*/
Datum
poly_above(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
- result = polya->boundbox.low.y > polyb->boundbox.high.y;
+ result = float8_gt(polya->boundbox.low.y, polyb->boundbox.high.y);
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
@@ -3622,21 +3624,21 @@ poly_above(PG_FUNCTION_ARGS)
* the lower most point of A at or above the lower most point
* of B?
*-------------------------------------------------------*/
Datum
poly_overabove(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
- result = polya->boundbox.low.y >= polyb->boundbox.low.y;
+ result = float8_ge(polya->boundbox.low.y, polyb->boundbox.low.y);
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
@@ -3833,22 +3835,22 @@ lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start)
}
if (res && !intersection)
{
Point p;
/*
* if X-intersection wasn't found then check central point of tested
* segment. In opposite case we already check all subsegments
*/
- p.x = (t.p[0].x + t.p[1].x) / 2.0;
- p.y = (t.p[0].y + t.p[1].y) / 2.0;
+ p.x = float8_div(float8_pl(t.p[0].x, t.p[1].x), 2.0);
+ p.y = float8_div(float8_pl(t.p[0].y, t.p[1].y), 2.0);
res = point_inside(&p, poly->npts, poly->p);
}
return res;
}
/*-----------------------------------------------------------------
* Determine if polygon A contains polygon B.
*-----------------------------------------------------------------*/
@@ -3975,66 +3977,68 @@ point_add(PG_FUNCTION_ARGS)
point_add_internal(result, p1, p2);
PG_RETURN_POINT_P(result);
}
static void
point_add_internal(Point *result, Point *pt1, Point *pt2)
{
point_construct(result,
- pt1->x + pt2->x,
- pt1->y + pt2->y);
+ float8_pl(pt1->x, pt2->x),
+ float8_pl(pt1->y, pt2->y));
}
Datum
point_sub(PG_FUNCTION_ARGS)
{
Point *p1 = PG_GETARG_POINT_P(0);
Point *p2 = PG_GETARG_POINT_P(1);
Point *result;
result = (Point *) palloc(sizeof(Point));
point_sub_internal(result, p1, p2);
PG_RETURN_POINT_P(result);
}
static void
point_sub_internal(Point *result, Point *pt1, Point *pt2)
{
point_construct(result,
- pt1->x - pt2->x,
- pt1->y - pt2->y);
+ float8_mi(pt1->x, pt2->x),
+ float8_mi(pt1->y, pt2->y));
}
Datum
point_mul(PG_FUNCTION_ARGS)
{
Point *p1 = PG_GETARG_POINT_P(0);
Point *p2 = PG_GETARG_POINT_P(1);
Point *result;
result = (Point *) palloc(sizeof(Point));
point_mul_internal(result, p1, p2);
PG_RETURN_POINT_P(result);
}
static void
point_mul_internal(Point *result, Point *pt1, Point *pt2)
{
point_construct(result,
- (pt1->x * pt2->x) - (pt1->y * pt2->y),
- (pt1->x * pt2->y) + (pt1->y * pt2->x));
+ float8_mi(float8_mul(pt1->x, pt2->x),
+ float8_mul(pt1->y, pt2->y)),
+ float8_pl(float8_mul(pt1->x, pt2->y),
+ float8_mul(pt1->y, pt2->x)));
}
Datum
point_div(PG_FUNCTION_ARGS)
{
Point *p1 = PG_GETARG_POINT_P(0);
Point *p2 = PG_GETARG_POINT_P(1);
Point *result;
result = (Point *) palloc(sizeof(Point));
@@ -4042,24 +4046,26 @@ point_div(PG_FUNCTION_ARGS)
point_div_internal(result, p1, p2);
PG_RETURN_POINT_P(result);
}
static void
point_div_internal(Point *result, Point *pt1, Point *pt2)
{
float8 div;
- div = (pt2->x * pt2->x) + (pt2->y * pt2->y);
+ div = float8_pl(float8_mul(pt2->x, pt2->x), float8_mul(pt2->y, pt2->y));
point_construct(result,
- ((pt1->x * pt2->x) + (pt1->y * pt2->y)) / div,
- ((pt2->x * pt1->y) - (pt2->y * pt1->x)) / div);
+ float8_div(float8_pl(float8_mul(pt1->x, pt2->x),
+ float8_mul(pt1->y, pt2->y)), div),
+ float8_div(float8_mi(float8_mul(pt1->y, pt2->x),
+ float8_mul(pt1->x, pt2->y)), div));
}
/***********************************************************************
**
** Routines for 2D boxes.
**
***********************************************************************/
Datum
@@ -4168,24 +4174,24 @@ point_box(PG_FUNCTION_ARGS)
*/
Datum
boxes_bound_box(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0),
*box2 = PG_GETARG_BOX_P(1),
*container;
container = (BOX *) palloc(sizeof(BOX));
- container->high.x = Max(box1->high.x, box2->high.x);
- container->low.x = Min(box1->low.x, box2->low.x);
- container->high.y = Max(box1->high.y, box2->high.y);
- container->low.y = Min(box1->low.y, box2->low.y);
+ container->high.x = float8_max(box1->high.x, box2->high.x);
+ container->low.x = float8_min(box1->low.x, box2->low.x);
+ container->high.y = float8_max(box1->high.y, box2->high.y);
+ container->low.y = float8_min(box1->low.y, box2->low.y);
PG_RETURN_BOX_P(container);
}
/***********************************************************************
**
** Routines for 2D paths.
**
***********************************************************************/
@@ -4495,21 +4501,21 @@ circle_in(PG_FUNCTION_ARGS)
if (*cp == LDELIM)
s = cp;
}
pair_decode(s, &circle->center.x, &circle->center.y, &s, "circle", str);
if (*s == DELIM)
s++;
circle->radius = single_decode(s, &s, "circle", str);
- if (circle->radius < 0)
+ if (float8_lt(circle->radius, 0.0))
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type %s: \"%s\"",
"circle", str)));
while (depth > 0)
{
if ((*s == RDELIM) || ((*s == RDELIM_C) && (depth == 1)))
{
depth--;
@@ -4562,21 +4568,21 @@ circle_recv(PG_FUNCTION_ARGS)
{
StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
CIRCLE *circle;
circle = (CIRCLE *) palloc(sizeof(CIRCLE));
circle->center.x = pq_getmsgfloat8(buf);
circle->center.y = pq_getmsgfloat8(buf);
circle->radius = pq_getmsgfloat8(buf);
- if (circle->radius < 0)
+ if (float8_lt(circle->radius, 0.0))
ereport(ERROR,
(errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
errmsg("invalid radius in external \"circle\" value")));
PG_RETURN_CIRCLE_P(circle);
}
/*
* circle_send - converts circle to binary format
*/
@@ -4613,144 +4619,146 @@ circle_same(PG_FUNCTION_ARGS)
/* circle_overlap - does circle1 overlap circle2?
*/
Datum
circle_overlap(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
- circle1->radius + circle2->radius));
+ float8_pl(circle1->radius, circle2->radius)));
}
/* circle_overleft - is the right edge of circle1 at or left of
* the right edge of circle2?
*/
Datum
circle_overleft(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
- PG_RETURN_BOOL(FPle((circle1->center.x + circle1->radius),
- (circle2->center.x + circle2->radius)));
+ PG_RETURN_BOOL(FPle(float8_pl(circle1->center.x, circle1->radius),
+ float8_pl(circle2->center.x, circle2->radius)));
}
/* circle_left - is circle1 strictly left of circle2?
*/
Datum
circle_left(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
- PG_RETURN_BOOL(FPlt((circle1->center.x + circle1->radius),
- (circle2->center.x - circle2->radius)));
+ PG_RETURN_BOOL(FPlt(float8_pl(circle1->center.x, circle1->radius),
+ float8_mi(circle2->center.x, circle2->radius)));
}
/* circle_right - is circle1 strictly right of circle2?
*/
Datum
circle_right(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
- PG_RETURN_BOOL(FPgt((circle1->center.x - circle1->radius),
- (circle2->center.x + circle2->radius)));
+ PG_RETURN_BOOL(FPgt(float8_mi(circle1->center.x, circle1->radius),
+ float8_pl(circle2->center.x, circle2->radius)));
}
/* circle_overright - is the left edge of circle1 at or right of
* the left edge of circle2?
*/
Datum
circle_overright(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
- PG_RETURN_BOOL(FPge((circle1->center.x - circle1->radius),
- (circle2->center.x - circle2->radius)));
+ PG_RETURN_BOOL(FPge(float8_mi(circle1->center.x, circle1->radius),
+ float8_mi(circle2->center.x, circle2->radius)));
}
/* circle_contained - is circle1 contained by circle2?
*/
Datum
circle_contained(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
- PG_RETURN_BOOL(FPle((point_dt(&circle1->center, &circle2->center) + circle1->radius), circle2->radius));
+ PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
+ float8_mi(circle2->radius, circle1->radius)));
}
/* circle_contain - does circle1 contain circle2?
*/
Datum
circle_contain(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
- PG_RETURN_BOOL(FPle((point_dt(&circle1->center, &circle2->center) + circle2->radius), circle1->radius));
+ PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
+ float8_mi(circle1->radius, circle2->radius)));
}
/* circle_below - is circle1 strictly below circle2?
*/
Datum
circle_below(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
- PG_RETURN_BOOL(FPlt((circle1->center.y + circle1->radius),
- (circle2->center.y - circle2->radius)));
+ PG_RETURN_BOOL(FPlt(float8_pl(circle1->center.y, circle1->radius),
+ float8_mi(circle2->center.y, circle2->radius)));
}
/* circle_above - is circle1 strictly above circle2?
*/
Datum
circle_above(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
- PG_RETURN_BOOL(FPgt((circle1->center.y - circle1->radius),
- (circle2->center.y + circle2->radius)));
+ PG_RETURN_BOOL(FPgt(float8_mi(circle1->center.y, circle1->radius),
+ float8_pl(circle2->center.y, circle2->radius)));
}
/* circle_overbelow - is the upper edge of circle1 at or below
* the upper edge of circle2?
*/
Datum
circle_overbelow(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
- PG_RETURN_BOOL(FPle((circle1->center.y + circle1->radius),
- (circle2->center.y + circle2->radius)));
+ PG_RETURN_BOOL(FPle(float8_pl(circle1->center.y, circle1->radius),
+ float8_pl(circle2->center.y, circle2->radius)));
}
/* circle_overabove - is the lower edge of circle1 at or above
* the lower edge of circle2?
*/
Datum
circle_overabove(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
- PG_RETURN_BOOL(FPge((circle1->center.y - circle1->radius),
- (circle2->center.y - circle2->radius)));
+ PG_RETURN_BOOL(FPge(float8_mi(circle1->center.y, circle1->radius),
+ float8_mi(circle2->center.y, circle2->radius)));
}
/* circle_relop - is area(circle1) relop area(circle2), within
* our accuracy constraint?
*/
Datum
circle_eq(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
@@ -4849,36 +4857,38 @@ circle_sub_pt(PG_FUNCTION_ARGS)
Datum
circle_mul_pt(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
Point *point = PG_GETARG_POINT_P(1);
CIRCLE *result;
result = (CIRCLE *) palloc(sizeof(CIRCLE));
point_mul_internal(&result->center, &circle->center, point);
- result->radius *= HYPOT(point->x, point->y);
+ result->radius = float8_mul(circle->radius,
+ float8_hypot(point->x, point->y));
PG_RETURN_CIRCLE_P(result);
}
Datum
circle_div_pt(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
Point *point = PG_GETARG_POINT_P(1);
CIRCLE *result;
result = (CIRCLE *) palloc(sizeof(CIRCLE));
point_div_internal(&result->center, &circle->center, point);
- result->radius /= HYPOT(point->x, point->y);
+ result->radius = float8_div(circle->radius,
+ float8_hypot(point->x, point->y));
PG_RETURN_CIRCLE_P(result);
}
/* circle_area - returns the area of the circle.
*/
Datum
circle_area(PG_FUNCTION_ARGS)
{
@@ -4888,21 +4898,21 @@ circle_area(PG_FUNCTION_ARGS)
}
/* circle_diameter - returns the diameter of the circle.
*/
Datum
circle_diameter(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
- PG_RETURN_FLOAT8(2 * circle->radius);
+ PG_RETURN_FLOAT8(float8_mul(circle->radius, 2.0));
}
/* circle_radius - returns the radius of the circle.
*/
Datum
circle_radius(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
@@ -4913,81 +4923,84 @@ circle_radius(PG_FUNCTION_ARGS)
/* circle_distance - returns the distance between
* two circles.
*/
Datum
circle_distance(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
float8 result;
- result = point_dt(&circle1->center, &circle2->center)
- - (circle1->radius + circle2->radius);
- if (result < 0)
- result = 0;
+ result = float8_mi(point_dt(&circle1->center, &circle2->center),
+ float8_pl(circle1->radius, circle2->radius));
+ if (float8_lt(result, 0.0))
+ result = 0.0;
+
PG_RETURN_FLOAT8(result);
}
Datum
circle_contain_pt(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
Point *point = PG_GETARG_POINT_P(1);
- double d;
+ float8 d;
d = point_dt(&circle->center, point);
- PG_RETURN_BOOL(d <= circle->radius);
+ PG_RETURN_BOOL(float8_le(d, circle->radius));
}
Datum
pt_contained_circle(PG_FUNCTION_ARGS)
{
Point *point = PG_GETARG_POINT_P(0);
CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
- double d;
+ float8 d;
d = point_dt(&circle->center, point);
- PG_RETURN_BOOL(d <= circle->radius);
+ PG_RETURN_BOOL(float8_le(d, circle->radius));
}
/* dist_pc - returns the distance between
* a point and a circle.
*/
Datum
dist_pc(PG_FUNCTION_ARGS)
{
Point *point = PG_GETARG_POINT_P(0);
CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
float8 result;
- result = point_dt(point, &circle->center) - circle->radius;
- if (result < 0)
- result = 0;
+ result = float8_mi(point_dt(point, &circle->center), circle->radius);
+ if (float8_lt(result, 0.0))
+ result = 0.0;
+
PG_RETURN_FLOAT8(result);
}
/*
* Distance from a circle to a point
*/
Datum
dist_cpoint(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
Point *point = PG_GETARG_POINT_P(1);
float8 result;
- result = point_dt(point, &circle->center) - circle->radius;
- if (result < 0)
- result = 0;
+ result = float8_mi(point_dt(point, &circle->center), circle->radius);
+ if (float8_lt(result, 0.0))
+ result = 0.0;
+
PG_RETURN_FLOAT8(result);
}
/* circle_center - returns the center point of the circle.
*/
Datum
circle_center(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
Point *result;
@@ -4995,24 +5008,24 @@ circle_center(PG_FUNCTION_ARGS)
result = (Point *) palloc(sizeof(Point));
result->x = circle->center.x;
result->y = circle->center.y;
PG_RETURN_POINT_P(result);
}
/* circle_ar - returns the area of the circle.
*/
-static double
+static float8
circle_ar(CIRCLE *circle)
{
- return M_PI * (circle->radius * circle->radius);
+ return float8_mul(float8_mul(circle->radius, circle->radius), M_PI);
}
/*----------------------------------------------------------
* Conversion operators.
*---------------------------------------------------------*/
Datum
cr_circle(PG_FUNCTION_ARGS)
{
@@ -5027,65 +5040,65 @@ cr_circle(PG_FUNCTION_ARGS)
result->radius = radius;
PG_RETURN_CIRCLE_P(result);
}
Datum
circle_box(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
BOX *box;
- double delta;
+ float8 delta;
box = (BOX *) palloc(sizeof(BOX));
- delta = circle->radius / sqrt(2.0);
+ delta = float8_div(circle->radius, sqrt(2.0));
- box->high.x = circle->center.x + delta;
- box->low.x = circle->center.x - delta;
- box->high.y = circle->center.y + delta;
- box->low.y = circle->center.y - delta;
+ box->high.x = float8_pl(circle->center.x, delta);
+ box->low.x = float8_mi(circle->center.x, delta);
+ box->high.y = float8_pl(circle->center.y, delta);
+ box->low.y = float8_mi(circle->center.y, delta);
PG_RETURN_BOX_P(box);
}
/* box_circle()
* Convert a box to a circle.
*/
Datum
box_circle(PG_FUNCTION_ARGS)
{
BOX *box = PG_GETARG_BOX_P(0);
CIRCLE *circle;
circle = (CIRCLE *) palloc(sizeof(CIRCLE));
- circle->center.x = (box->high.x + box->low.x) / 2;
- circle->center.y = (box->high.y + box->low.y) / 2;
+ circle->center.x = float8_div(float8_pl(box->high.x, box->low.x), 2.0);
+ circle->center.y = float8_div(float8_pl(box->high.y, box->low.y), 2.0);
circle->radius = point_dt(&circle->center, &box->high);
PG_RETURN_CIRCLE_P(circle);
}
Datum
circle_poly(PG_FUNCTION_ARGS)
{
int32 npts = PG_GETARG_INT32(0);
CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
POLYGON *poly;
int base_size,
size;
int i;
- double angle;
- double anglestep;
+ float8 angle;
+ float8 anglestep;
if (FPzero(circle->radius))
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("cannot convert circle with radius zero to polygon")));
if (npts < 2)
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("must request at least 2 points")));
@@ -5096,27 +5109,30 @@ circle_poly(PG_FUNCTION_ARGS)
/* Check for integer overflow */
if (base_size / npts != sizeof(poly->p[0]) || size <= base_size)
ereport(ERROR,
(errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
errmsg("too many points requested")));
poly = (POLYGON *) palloc0(size); /* zero any holes */
SET_VARSIZE(poly, size);
poly->npts = npts;
- anglestep = (2.0 * M_PI) / npts;
+ anglestep = float8_div(2.0 * M_PI, npts);
for (i = 0; i < npts; i++)
{
- angle = i * anglestep;
- poly->p[i].x = circle->center.x - (circle->radius * cos(angle));
- poly->p[i].y = circle->center.y + (circle->radius * sin(angle));
+ angle = float8_mul(anglestep, i);
+
+ poly->p[i].x = float8_mi(circle->center.x,
+ float8_mul(circle->radius, cos(angle)));
+ poly->p[i].y = float8_pl(circle->center.y,
+ float8_mul(circle->radius, sin(angle)));
}
make_bound_box(poly);
PG_RETURN_POLYGON_P(poly);
}
/* poly_circle - convert polygon to circle
*
* XXX This algorithm should use weighted means of line segments
@@ -5135,29 +5151,30 @@ poly_circle(PG_FUNCTION_ARGS)
errmsg("cannot convert empty polygon to circle")));
circle = (CIRCLE *) palloc(sizeof(CIRCLE));
circle->center.x = 0;
circle->center.y = 0;
circle->radius = 0;
for (i = 0; i < poly->npts; i++)
{
- circle->center.x += poly->p[i].x;
- circle->center.y += poly->p[i].y;
+ circle->center.x = float8_pl(circle->center.x, poly->p[i].x);
+ circle->center.y = float8_pl(circle->center.y, poly->p[i].y);
}
- circle->center.x /= poly->npts;
- circle->center.y /= poly->npts;
+ circle->center.x = float8_div(circle->center.x, poly->npts);
+ circle->center.y = float8_div(circle->center.y, poly->npts);
for (i = 0; i < poly->npts; i++)
- circle->radius += point_dt(&poly->p[i], &circle->center);
- circle->radius /= poly->npts;
+ circle->radius = float8_pl(circle->radius, point_dt(&poly->p[i],
+ &circle->center));
+ circle->radius = float8_div(circle->radius, poly->npts);
PG_RETURN_CIRCLE_P(circle);
}
/***********************************************************************
**
** Private routines for multiple types.
**
***********************************************************************/
@@ -5171,45 +5188,45 @@ poly_circle(PG_FUNCTION_ARGS)
* http://hopf.math.northwestern.edu/index.html
* Description of algorithm: http://www.linuxjournal.com/article/2197
* http://www.linuxjournal.com/article/2029
*/
#define POINT_ON_POLYGON INT_MAX
static int
point_inside(Point *p, int npts, Point *plist)
{
- double x0,
+ float8 x0,
y0;
- double prev_x,
+ float8 prev_x,
prev_y;
int i = 0;
- double x,
+ float8 x,
y;
int cross,
total_cross = 0;
if (npts <= 0)
return 0;
/* compute first polygon point relative to single point */
- x0 = plist[0].x - p->x;
- y0 = plist[0].y - p->y;
+ x0 = float8_mi(plist[0].x, p->x);
+ y0 = float8_mi(plist[0].y, p->y);
prev_x = x0;
prev_y = y0;
/* loop over polygon points and aggregate total_cross */
for (i = 1; i < npts; i++)
{
/* compute next polygon point relative to single point */
- x = plist[i].x - p->x;
- y = plist[i].y - p->y;
+ x = float8_mi(plist[i].x, p->x);
+ y = float8_mi(plist[i].y, p->y);
/* compute previous to current point crossing */
if ((cross = lseg_crossing(x, y, prev_x, prev_y)) == POINT_ON_POLYGON)
return 2;
total_cross += cross;
prev_x = x;
prev_y = y;
}
@@ -5227,69 +5244,74 @@ point_inside(Point *p, int npts, Point *plist)
/* lseg_crossing()
* Returns +/-2 if line segment crosses the positive X-axis in a +/- direction.
* Returns +/-1 if one point is on the positive X-axis.
* Returns 0 if both points are on the positive X-axis, or there is no crossing.
* Returns POINT_ON_POLYGON if the segment contains (0,0).
* Wow, that is one confusing API, but it is used above, and when summed,
* can tell is if a point is in a polygon.
*/
static int
-lseg_crossing(double x, double y, double prev_x, double prev_y)
+lseg_crossing(float8 x, float8 y, float8 prev_x, float8 prev_y)
{
- double z;
+ float8 z;
int y_sign;
if (FPzero(y))
{ /* y == 0, on X axis */
if (FPzero(x)) /* (x,y) is (0,0)? */
return POINT_ON_POLYGON;
else if (FPgt(x, 0))
{ /* x > 0 */
if (FPzero(prev_y)) /* y and prev_y are zero */
/* prev_x > 0? */
- return FPgt(prev_x, 0) ? 0 : POINT_ON_POLYGON;
- return FPlt(prev_y, 0) ? 1 : -1;
+ return FPgt(prev_x, 0.0) ? 0 : POINT_ON_POLYGON;
+ return FPlt(prev_y, 0.0) ? 1 : -1;
}
else
{ /* x < 0, x not on positive X axis */
if (FPzero(prev_y))
/* prev_x < 0? */
- return FPlt(prev_x, 0) ? 0 : POINT_ON_POLYGON;
+ return FPlt(prev_x, 0.0) ? 0 : POINT_ON_POLYGON;
return 0;
}
}
else
{ /* y != 0 */
/* compute y crossing direction from previous point */
- y_sign = FPgt(y, 0) ? 1 : -1;
+ y_sign = FPgt(y, 0.0) ? 1 : -1;
if (FPzero(prev_y))
/* previous point was on X axis, so new point is either off or on */
- return FPlt(prev_x, 0) ? 0 : y_sign;
- else if (FPgt(y_sign * prev_y, 0))
+ return FPlt(prev_x, 0.0) ? 0 : y_sign;
+ else if ((y_sign < 0 && FPlt(prev_y, 0.0)) ||
+ (y_sign > 0 && FPgt(prev_y, 0.0)))
/* both above or below X axis */
return 0; /* same sign */
else
{ /* y and prev_y cross X-axis */
- if (FPge(x, 0) && FPgt(prev_x, 0))
+ if (FPge(x, 0.0) && FPgt(prev_x, 0.0))
/* both non-negative so cross positive X-axis */
return 2 * y_sign;
- if (FPlt(x, 0) && FPle(prev_x, 0))
+ if (FPlt(x, 0.0) && FPle(prev_x, 0.0))
/* both non-positive so do not cross positive X-axis */
return 0;
/* x and y cross axises, see URL above point_inside() */
- z = (x - prev_x) * y - (y - prev_y) * x;
+ z = float8_mi(float8_mul(float8_mi(x, prev_x), y),
+ float8_mul(float8_mi(y, prev_y), x));
if (FPzero(z))
return POINT_ON_POLYGON;
- return FPgt((y_sign * z), 0) ? 0 : 2 * y_sign;
+ if ((y_sign < 0 && FPlt(z, 0.0)) ||
+ (y_sign > 0 && FPgt(z, 0.0)))
+ return 0;
+ return 2 * y_sign;
}
}
}
static bool
plist_same(int npts, Point *p1, Point *p2)
{
int i,
ii,
diff --git a/src/backend/utils/adt/geo_spgist.c b/src/backend/utils/adt/geo_spgist.c
index c800bb1338..f62be1061e 100644
--- a/src/backend/utils/adt/geo_spgist.c
+++ b/src/backend/utils/adt/geo_spgist.c
@@ -76,38 +76,38 @@
#include "access/spgist.h"
#include "access/stratnum.h"
#include "catalog/pg_type.h"
#include "utils/fmgrprotos.h"
#include "utils/geo_decls.h"
/*
* Comparator for qsort
*
* We don't need to use the floating point macros in here, because this
- * is going only going to be used in a place to effect the performance
+ * is only going to be used in a place to effect the performance
* of the index, not the correctness.
*/
static int
compareDoubles(const void *a, const void *b)
{
- double x = *(double *) a;
- double y = *(double *) b;
+ float8 x = *(float8 *) a;
+ float8 y = *(float8 *) b;
if (x == y)
return 0;
return (x > y) ? 1 : -1;
}
typedef struct
{
- double low;
- double high;
+ float8 low;
+ float8 high;
} Range;
typedef struct
{
Range left;
Range right;
} RangeBox;
typedef struct
{
@@ -121,30 +121,30 @@ typedef struct
* The quadrant is 8 bit unsigned integer with 4 least bits in use.
* This function accepts BOXes as input. They are not casted to
* RangeBoxes, yet. All 4 bits are set by comparing a corner of the box.
* This makes 16 quadrants in total.
*/
static uint8
getQuadrant(BOX *centroid, BOX *inBox)
{
uint8 quadrant = 0;
- if (inBox->low.x > centroid->low.x)
+ if (float8_gt(inBox->low.x, centroid->low.x))
quadrant |= 0x8;
- if (inBox->high.x > centroid->high.x)
+ if (float8_gt(inBox->high.x, centroid->high.x))
quadrant |= 0x4;
- if (inBox->low.y > centroid->low.y)
+ if (float8_gt(inBox->low.y, centroid->low.y))
quadrant |= 0x2;
- if (inBox->high.y > centroid->high.y)
+ if (float8_gt(inBox->high.y, centroid->high.y))
quadrant |= 0x1;
return quadrant;
}
/*
* Get RangeBox using BOX
*
* We are turning the BOX to our structures to emphasize their function
* of representing points in 4D space. It also is more convenient to
@@ -167,21 +167,21 @@ getRangeBox(BOX *box)
/*
* Initialize the traversal value
*
* In the beginning, we don't have any restrictions. We have to
* initialize the struct to cover the whole 4D space.
*/
static RectBox *
initRectBox(void)
{
RectBox *rect_box = (RectBox *) palloc(sizeof(RectBox));
- double infinity = get_float8_infinity();
+ float8 infinity = get_float8_infinity();
rect_box->range_box_x.left.low = -infinity;
rect_box->range_box_x.left.high = infinity;
rect_box->range_box_x.right.low = -infinity;
rect_box->range_box_x.right.high = infinity;
rect_box->range_box_y.left.low = -infinity;
rect_box->range_box_y.left.high = infinity;
@@ -410,40 +410,40 @@ spg_box_quad_choose(PG_FUNCTION_ARGS)
* point as the median of the coordinates of the boxes.
*/
Datum
spg_box_quad_picksplit(PG_FUNCTION_ARGS)
{
spgPickSplitIn *in = (spgPickSplitIn *) PG_GETARG_POINTER(0);
spgPickSplitOut *out = (spgPickSplitOut *) PG_GETARG_POINTER(1);
BOX *centroid;
int median,
i;
- double *lowXs = palloc(sizeof(double) * in->nTuples);
- double *highXs = palloc(sizeof(double) * in->nTuples);
- double *lowYs = palloc(sizeof(double) * in->nTuples);
- double *highYs = palloc(sizeof(double) * in->nTuples);
+ float8 *lowXs = palloc(sizeof(float8) * in->nTuples);
+ float8 *highXs = palloc(sizeof(float8) * in->nTuples);
+ float8 *lowYs = palloc(sizeof(float8) * in->nTuples);
+ float8 *highYs = palloc(sizeof(float8) * in->nTuples);
/* Calculate median of all 4D coordinates */
for (i = 0; i < in->nTuples; i++)
{
BOX *box = DatumGetBoxP(in->datums[i]);
lowXs[i] = box->low.x;
highXs[i] = box->high.x;
lowYs[i] = box->low.y;
highYs[i] = box->high.y;
}
- qsort(lowXs, in->nTuples, sizeof(double), compareDoubles);
- qsort(highXs, in->nTuples, sizeof(double), compareDoubles);
- qsort(lowYs, in->nTuples, sizeof(double), compareDoubles);
- qsort(highYs, in->nTuples, sizeof(double), compareDoubles);
+ qsort(lowXs, in->nTuples, sizeof(float8), compareDoubles);
+ qsort(highXs, in->nTuples, sizeof(float8), compareDoubles);
+ qsort(lowYs, in->nTuples, sizeof(float8), compareDoubles);
+ qsort(highYs, in->nTuples, sizeof(float8), compareDoubles);
median = in->nTuples / 2;
centroid = palloc(sizeof(BOX));
centroid->low.x = lowXs[median];
centroid->high.x = highXs[median];
centroid->low.y = lowYs[median];
centroid->high.y = highYs[median];
diff --git a/src/include/utils/float.h b/src/include/utils/float.h
index c26de1cca4..c893f79a4d 100644
--- a/src/include/utils/float.h
+++ b/src/include/utils/float.h
@@ -269,20 +269,33 @@ float8_div(float8 val1, float8 val2)
ereport(ERROR,
(errcode(ERRCODE_DIVISION_BY_ZERO),
errmsg("division by zero")));
result = val1 / val2;
check_float8_val(result, isinf(val1) || isinf(val2), val1 == 0.0);
return result;
}
+static inline float8
+float8_hypot(float8 val1, float8 val2)
+{
+ float8 result;
+
+ result = hypot(val1, val2);
+ check_float8_val(result, isinf(val1) || isinf(val2),
+ val1 == 0.0 && val2 == 0.0);
+ Assert(isinf(result) || result >= 0.0);
+
+ return result;
+}
+
/*
* Routines for NaN-aware comparisons
*
* We consider all NANs to be equal and larger than any non-NAN. This is
* somewhat arbitrary; the important thing is to have a consistent sort
* order.
*/
static inline bool
float4_eq(float4 val1, float4 val2)
diff --git a/src/include/utils/geo_decls.h b/src/include/utils/geo_decls.h
index f845dffcb7..527036f1fc 100644
--- a/src/include/utils/geo_decls.h
+++ b/src/include/utils/geo_decls.h
@@ -1,23 +1,20 @@
/*-------------------------------------------------------------------------
*
* geo_decls.h - Declarations for various 2D constructs.
*
*
* Portions Copyright (c) 1996-2017, PostgreSQL Global Development Group
* Portions Copyright (c) 1994, Regents of the University of California
*
* src/include/utils/geo_decls.h
*
- * NOTE
- * These routines do *not* use the float types from adt/.
- *
* XXX These routines were not written by a numerical analyst.
*
* XXX I have made some attempt to flesh out the operators
* and data types. There are still some more to do. - tgl 97/04/19
*
*-------------------------------------------------------------------------
*/
#ifndef GEO_DECLS_H
#define GEO_DECLS_H
@@ -28,44 +25,43 @@
/*--------------------------------------------------------------------
* Useful floating point utilities and constants.
*-------------------------------------------------------------------*/
#define EPSILON 1.0E-06
#ifdef EPSILON
#define FPzero(A) (fabs(A) <= EPSILON)
-#define FPeq(A,B) (fabs((A) - (B)) <= EPSILON)
-#define FPne(A,B) (fabs((A) - (B)) > EPSILON)
-#define FPlt(A,B) ((B) - (A) > EPSILON)
-#define FPle(A,B) ((A) - (B) <= EPSILON)
-#define FPgt(A,B) ((A) - (B) > EPSILON)
-#define FPge(A,B) ((B) - (A) <= EPSILON)
+#define FPeq(A, B) (float8_le(fabs(float8_mi(A, B)), EPSILON))
+#define FPne(A, B) (float8_gt(fabs(float8_mi(A, B)), EPSILON))
+#define FPlt(A, B) (float8_gt(float8_mi(B, A), EPSILON))
+#define FPle(A, B) (float8_le(float8_mi(A, B), EPSILON))
+#define FPgt(A, B) (float8_gt(float8_mi(A, B), EPSILON))
+#define FPge(A, B) (float8_le(float8_mi(B, A), EPSILON))
#else
-#define FPzero(A) ((A) == 0)
-#define FPeq(A,B) ((A) == (B))
-#define FPne(A,B) ((A) != (B))
-#define FPlt(A,B) ((A) < (B))
-#define FPle(A,B) ((A) <= (B))
-#define FPgt(A,B) ((A) > (B))
-#define FPge(A,B) ((A) >= (B))
+#define FPzero(A) ((A) == 0.0)
+#define FPeq(A, B) (float8_eq(A, B))
+#define FPne(A, B) (float8_ne(A, B))
+#define FPlt(A, B) (float8_lt(A, B))
+#define FPle(A, B) (float8_le(A, B))
+#define FPgt(A, B) (float8_gt(A, B))
+#define FPge(A, B) (float8_ge(A, B))
#endif
-#define HYPOT(A, B) hypot(A, B)
/*---------------------------------------------------------------------
* Point - (x,y)
*-------------------------------------------------------------------*/
typedef struct
{
- double x,
+ float8 x,
y;
} Point;
/*---------------------------------------------------------------------
* LSEG - A straight line, specified by endpoints.
*-------------------------------------------------------------------*/
typedef struct
{
Point p[2];
@@ -73,66 +69,66 @@ typedef struct
/*---------------------------------------------------------------------
* PATH - Specified by vertex points.
*-------------------------------------------------------------------*/
typedef struct
{
int32 vl_len_; /* varlena header (do not touch directly!) */
int32 npts;
int32 closed; /* is this a closed polygon? */
- int32 dummy; /* padding to make it double align */
+ int32 dummy; /* padding to make it float8 align */
Point p[FLEXIBLE_ARRAY_MEMBER];
} PATH;
/*---------------------------------------------------------------------
* LINE - Specified by its general equation (Ax+By+C=0).
*-------------------------------------------------------------------*/
typedef struct
{
- double A,
+ float8 A,
B,
C;
} LINE;
/*---------------------------------------------------------------------
* BOX - Specified by two corner points, which are
* sorted to save calculation time later.
*-------------------------------------------------------------------*/
typedef struct
{
Point high,
low; /* corner POINTs */
} BOX;
/*---------------------------------------------------------------------
- * POLYGON - Specified by an array of doubles defining the points,
+ * POLYGON - Specified by an array of float8s defining the points,
* keeping the number of points and the bounding box for
* speed purposes.
*-------------------------------------------------------------------*/
typedef struct
{
int32 vl_len_; /* varlena header (do not touch directly!) */
int32 npts;
BOX boundbox;
Point p[FLEXIBLE_ARRAY_MEMBER];
} POLYGON;
/*---------------------------------------------------------------------
* CIRCLE - Specified by a center point and radius.
*-------------------------------------------------------------------*/
typedef struct
{
Point center;
- double radius;
+ float8 radius;
} CIRCLE;
/*
* fmgr interface macros
*
* Path and Polygon are toastable varlena types, the others are just
* fixed-size pass-by-reference types.
*/
#define DatumGetPointP(X) ((Point *) DatumGetPointer(X))
--
2.11.0 (Apple Git-81)