[SCM] PostGIS branch master updated. 3.7.0alpha1-14-gd7a16bb43
git at osgeo.org
git at osgeo.org
Mon Jul 6 02:53:10 PDT 2026
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "PostGIS".
The branch, master has been updated
via d7a16bb43295c653e86c69ae5eeace3661587ca3 (commit)
via 408572b681512693cadc9ce9e00d3f090bbb66d9 (commit)
via ff20f053438d1c082f05bc5834a44eca083aa7bb (commit)
via c0521bcc74cd4168f495909c4fa31a1cc4cd956e (commit)
via 0e16d8991f1356c0b735d6ab3046e036a7e08740 (commit)
via 5336ccf1fe35d63672357a11f0082f4b70c0bbf8 (commit)
via 5124ad4a50a4a914e4b188240997ee3f06e258d7 (commit)
via 6fff4b58bbf6b671fdb34a0bbaa70d1f8ab11b2c (commit)
via 6ea8421662ed55c8b895e1a6e4d4ef3edefde4c6 (commit)
via 8a03193ebc452f335c0f51c26dbdbb53b1a15904 (commit)
via d9525f032bba08eefdce60e55f333c101dace23d (commit)
via 9000ef67ed38616b10e06d8da7bd012f9d866aa0 (commit)
from 9b4d142f5b771557dcc19f10473f2553a3b16b3a (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit d7a16bb43295c653e86c69ae5eeace3661587ca3
Author: Darafei Praliaskouski <me at komzpa.net>
Date: Mon Jul 6 13:46:12 2026 +0400
NEWS for GH-897
References https://github.com/postgis/postgis/pull/1143
References https://github.com/postgis/postgis/pull/1147
References https://github.com/postgis/postgis/pull/1148
References https://github.com/postgis/postgis/pull/1149
References https://github.com/postgis/postgis/pull/1153
diff --git a/NEWS b/NEWS
index c4a2ebec8..6fda3ea30 100644
--- a/NEWS
+++ b/NEWS
@@ -171,6 +171,9 @@ To take advantage of all postgis_sfcgal extension features SFCGAL 2.3+ is needed
(Darafei Praliaskouski)
- GH-895, Qualify pg_class lookup during SRID/PROJ validation
(Darafei Praliaskouski)
+ - GH-897, [topology] Harden topology helper functions and
+ maintenance SQL against identifier, tolerance, bigint, and cleanup
+ regressions (Darafei Praliaskouski)
- GH-898, Fix ST_DFullyWithin indexed plans and large-coordinate
ring-orientation precision regressions (Darafei Praliaskouski)
- #4828, geometry_columns handles NOT VALID SRID checks without errors (Darafei Praliaskouski)
commit 408572b681512693cadc9ce9e00d3f090bbb66d9
Author: Darafei Praliaskouski <me at komzpa.net>
Date: Fri Jun 19 14:52:46 2026 +0400
liblwgeom: preserve scaled area subtraction tails
(cherry picked from commit 04ea138387f71044c700ceb28ec0719a7b3ef71c)
Closes https://github.com/postgis/postgis/pull/898
diff --git a/liblwgeom/cunit/cu_ptarray.c b/liblwgeom/cunit/cu_ptarray.c
index 933a5ca4f..f616a5928 100644
--- a/liblwgeom/cunit/cu_ptarray.c
+++ b/liblwgeom/cunit/cu_ptarray.c
@@ -479,6 +479,29 @@ static void test_ptarray_signed_area(void)
CU_ASSERT(area < 6e99);
CU_ASSERT_EQUAL(ptarray_isccw(line->points), 0);
lwline_free(line);
+
+ /* Detect residual tails after scaled determinant cancellation. */
+ line =
+ lwgeom_as_lwline(lwgeom_from_text("LINESTRING(3e120 1.0000000000000002e220,"
+ "3e155 -1e-50,3e155 -1.0000000000000003e140,"
+ "3e120 1.0000000000000002e220)"));
+ area = ptarray_signed_area(line->points);
+ CU_ASSERT(isfinite(area));
+ CU_ASSERT(area > 1e295);
+ CU_ASSERT(area < 2e295);
+ CU_ASSERT_EQUAL(ptarray_isccw(line->points), 0);
+ lwline_free(line);
+
+ line =
+ lwgeom_as_lwline(lwgeom_from_text("LINESTRING(3e120 1.0000000000000002e220,"
+ "3e155 -1.0000000000000003e140,3e155 -1e-50,"
+ "3e120 1.0000000000000002e220)"));
+ area = ptarray_signed_area(line->points);
+ CU_ASSERT(isfinite(area));
+ CU_ASSERT(area < -1e295);
+ CU_ASSERT(area > -2e295);
+ CU_ASSERT_EQUAL(ptarray_isccw(line->points), 1);
+ lwline_free(line);
}
static void test_ptarray_contains_point(void)
diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index 036bfd950..ce7632769 100644
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -67,6 +67,15 @@ ptarray_signed_area_add(double *sum, double *compensation, double value)
*sum = t;
}
+static inline void
+ptarray_signed_area_add_abs_or_infinite(double *sum, double value)
+{
+ if (isfinite(value))
+ *sum += fabs(value);
+ else
+ *sum = INFINITY;
+}
+
static inline void
ptarray_signed_area_two_sum(double a, double b, double *sum, double *err)
{
@@ -1370,21 +1379,57 @@ ptarray_signed_area(const POINTARRAY *pa)
*/
if (isfinite(detleft) && isfinite(detright) && isfinite(det))
{
+ ptarray_signed_area_expansion residual = {0};
+ double tail_abs = 0.0;
+ double subtract_roundoff;
+ double left_roundoff;
+ double right_roundoff;
+ double ax_byerr, axerr_by, axerr_byerr;
+ double ay_bxerr, ayerr_bx, ayerr_bxerr;
+
b = det - detleft;
- err = (detleft - (det - b)) - (detright + b);
- err += fma(ax, by, -detleft) - fma(ay, bx, -detright);
+ subtract_roundoff = (detleft - (det - b)) - (detright + b);
+ left_roundoff = fma(ax, by, -detleft);
+ right_roundoff = fma(ay, bx, -detright);
+ ptarray_signed_area_expansion_add(&residual, subtract_roundoff);
+ ptarray_signed_area_expansion_add(&residual, left_roundoff);
+ ptarray_signed_area_expansion_add(&residual, -right_roundoff);
+ ptarray_signed_area_add_abs_or_infinite(&tail_abs, subtract_roundoff);
+ ptarray_signed_area_add_abs_or_infinite(&tail_abs, left_roundoff);
+ ptarray_signed_area_add_abs_or_infinite(&tail_abs, right_roundoff);
/*
* The origin shift itself can round away unit-scale
* offsets beside 1e16-scale coordinates. Include the
* low-order subtraction terms in the determinant so the
- * compensated sum still sees their signed area.
+ * compensated sum still sees their signed area. Use an
+ * expansion so opposite large low-order products do not
+ * erase a smaller residual before cancellation is tested.
*/
- err += (ax * byerr + axerr * by) + axerr * byerr;
- err -= (ay * bxerr + ayerr * bx) + ayerr * bxerr;
+ ax_byerr = ax * byerr;
+ axerr_by = axerr * by;
+ axerr_byerr = axerr * byerr;
+ ay_bxerr = ay * bxerr;
+ ayerr_bx = ayerr * bx;
+ ayerr_bxerr = ayerr * bxerr;
+ ptarray_signed_area_expansion_add_product(&residual, ax, byerr, 1.0);
+ ptarray_signed_area_expansion_add_product(&residual, axerr, by, 1.0);
+ ptarray_signed_area_expansion_add_product(&residual, axerr, byerr, 1.0);
+ ptarray_signed_area_expansion_add_product(&residual, ay, bxerr, -1.0);
+ ptarray_signed_area_expansion_add_product(&residual, ayerr, bx, -1.0);
+ ptarray_signed_area_expansion_add_product(&residual, ayerr, bxerr, -1.0);
+ ptarray_signed_area_add_abs_or_infinite(&tail_abs, ax_byerr);
+ ptarray_signed_area_add_abs_or_infinite(&tail_abs, axerr_by);
+ ptarray_signed_area_add_abs_or_infinite(&tail_abs, axerr_byerr);
+ ptarray_signed_area_add_abs_or_infinite(&tail_abs, ay_bxerr);
+ ptarray_signed_area_add_abs_or_infinite(&tail_abs, ayerr_bx);
+ ptarray_signed_area_add_abs_or_infinite(&tail_abs, ayerr_bxerr);
+ err = ptarray_signed_area_expansion_sum(&residual);
+ if (!isfinite(err))
+ err = 0.0;
if (scaled_term)
- scaled_tail_abs_sum += fabs(err);
+ scaled_tail_abs_sum += tail_abs;
else
- tail_abs_sum += fabs(err);
+ tail_abs_sum += tail_abs;
}
ptarray_signed_area_add(active_sum, active_compensation, det);
diff --git a/regress/core/orientation.sql b/regress/core/orientation.sql
index 54a71e855..f8f3d2af9 100644
--- a/regress/core/orientation.sql
+++ b/regress/core/orientation.sql
@@ -99,3 +99,8 @@ SELECT '439', ST_Area('POLYGON ((1.9029580275885475e-13 2.1370106480605854e-89,
SELECT '440', ST_IsPolygonCW( 'POLYGON ((1.9029580275885475e47 2.1370106480605854e-29, 2.5163957753386382e-33 -3.002966181933756e69, 4.77621563620871e-15 -1e-40, -3.7953949389104665e134 -1e260, -2.0800734986321869 2.051992028072436e5, 3.38586824330807e133 -7.962294473171359e134, 1.9029580275885475e47 2.1370106480605854e-29))');
SELECT '441', ST_IsPolygonCCW('POLYGON ((1.9029580275885475e47 2.1370106480605854e-29, 2.5163957753386382e-33 -3.002966181933756e69, 4.77621563620871e-15 -1e-40, -3.7953949389104665e134 -1e260, -2.0800734986321869 2.051992028072436e5, 3.38586824330807e133 -7.962294473171359e134, 1.9029580275885475e47 2.1370106480605854e-29))');
SELECT '442', ST_Area('POLYGON ((1.9029580275885475e47 2.1370106480605854e-29, 2.5163957753386382e-33 -3.002966181933756e69, 4.77621563620871e-15 -1e-40, -3.7953949389104665e134 -1e260, -2.0800734986321869 2.051992028072436e5, 3.38586824330807e133 -7.962294473171359e134, 1.9029580275885475e47 2.1370106480605854e-29))') > 1e260;
+
+-- Scaled subtraction-tail residual after grouped product-tail cancellation
+SELECT '443', ST_IsPolygonCW( 'POLYGON ((3e120 1.0000000000000002e220, 3e155 -1e-50, 3e155 -1.0000000000000003e140, 3e120 1.0000000000000002e220))');
+SELECT '444', ST_IsPolygonCCW('POLYGON ((3e120 1.0000000000000002e220, 3e155 -1e-50, 3e155 -1.0000000000000003e140, 3e120 1.0000000000000002e220))');
+SELECT '445', ST_Area('POLYGON ((3e120 1.0000000000000002e220, 3e155 -1e-50, 3e155 -1.0000000000000003e140, 3e120 1.0000000000000002e220))') > 1e295;
diff --git a/regress/core/orientation_expected b/regress/core/orientation_expected
index 060a96b4f..1d918cede 100644
--- a/regress/core/orientation_expected
+++ b/regress/core/orientation_expected
@@ -58,3 +58,6 @@
440|t
441|f
442|t
+443|t
+444|f
+445|t
commit ff20f053438d1c082f05bc5834a44eca083aa7bb
Author: Darafei Praliaskouski <me at komzpa.net>
Date: Fri Jun 19 14:22:34 2026 +0400
liblwgeom: preserve scaled area residual tails
(cherry picked from commit 89811473e6cbd8bb6cdfff34ecdb6fec9e2ed928)
diff --git a/liblwgeom/cunit/cu_ptarray.c b/liblwgeom/cunit/cu_ptarray.c
index c8cac7dad..933a5ca4f 100644
--- a/liblwgeom/cunit/cu_ptarray.c
+++ b/liblwgeom/cunit/cu_ptarray.c
@@ -382,6 +382,35 @@ static void test_ptarray_signed_area(void)
CU_ASSERT_EQUAL(ptarray_isccw(line->points), 1);
lwline_free(line);
+ /* Preserve residual tails when scaled large product tails cancel. */
+ line =
+ lwgeom_as_lwline(lwgeom_from_text("LINESTRING(1.9029580275885475e47 2.1370106480605854e-29,"
+ "2.5163957753386382e-33 -3.002966181933756e69,"
+ "4.77621563620871e-15 -1e-40,"
+ "-3.7953949389104665e134 -1e260,"
+ "-2.0800734986321869 2.051992028072436e5,"
+ "3.38586824330807e133 -7.962294473171359e134,"
+ "1.9029580275885475e47 2.1370106480605854e-29)"));
+ area = ptarray_signed_area(line->points);
+ CU_ASSERT(area > 1e260);
+ CU_ASSERT(area < 2e260);
+ CU_ASSERT_EQUAL(ptarray_isccw(line->points), 0);
+ lwline_free(line);
+
+ line =
+ lwgeom_as_lwline(lwgeom_from_text("LINESTRING(1.9029580275885475e47 2.1370106480605854e-29,"
+ "3.38586824330807e133 -7.962294473171359e134,"
+ "-2.0800734986321869 2.051992028072436e5,"
+ "-3.7953949389104665e134 -1e260,"
+ "4.77621563620871e-15 -1e-40,"
+ "2.5163957753386382e-33 -3.002966181933756e69,"
+ "1.9029580275885475e47 2.1370106480605854e-29)"));
+ area = ptarray_signed_area(line->points);
+ CU_ASSERT(area < -1e260);
+ CU_ASSERT(area > -2e260);
+ CU_ASSERT_EQUAL(ptarray_isccw(line->points), 1);
+ lwline_free(line);
+
/* Preserve product tails when determinant products nearly cancel. */
line =
lwgeom_as_lwline(lwgeom_from_text("LINESTRING(10000000000000276 10000000000000280,"
diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index 681f1a79d..036bfd950 100644
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -174,7 +174,7 @@ ptarray_signed_area_product_exponent(double a, double b)
}
static double
-ptarray_signed_area_exact(const POINTARRAY *pa)
+ptarray_signed_area_exact(const POINTARRAY *pa, int coordinate_scale)
{
const POINT2D *P1 = getPoint2d_cp(pa, 0);
const POINT2D *P2 = getPoint2d_cp(pa, 1);
@@ -193,6 +193,17 @@ ptarray_signed_area_exact(const POINTARRAY *pa)
ptarray_signed_area_two_diff(P2->y, y0, &ay, &ayerr);
ptarray_signed_area_two_diff(P3->x, x0, &bx, &bxerr);
ptarray_signed_area_two_diff(P3->y, y0, &by, &byerr);
+ if (coordinate_scale > 0)
+ {
+ ax = ldexp(ax, -coordinate_scale);
+ ay = ldexp(ay, -coordinate_scale);
+ bx = ldexp(bx, -coordinate_scale);
+ by = ldexp(by, -coordinate_scale);
+ axerr = ldexp(axerr, -coordinate_scale);
+ ayerr = ldexp(ayerr, -coordinate_scale);
+ bxerr = ldexp(bxerr, -coordinate_scale);
+ byerr = ldexp(byerr, -coordinate_scale);
+ }
ptarray_signed_area_expansion_add_product(&expansion, ax, by, 1.0);
ptarray_signed_area_expansion_add_product(&expansion, ax, byerr, 1.0);
@@ -206,7 +217,7 @@ ptarray_signed_area_exact(const POINTARRAY *pa)
P2 = P3;
}
- return -ptarray_signed_area_expansion_sum(&expansion) / 2.0;
+ return -ldexp(ptarray_signed_area_expansion_sum(&expansion), 2 * coordinate_scale) / 2.0;
}
int
@@ -1267,8 +1278,12 @@ ptarray_signed_area(const POINTARRAY *pa)
double scaled_tail_compensation = 0.0;
double det_abs_sum = 0.0;
double tail_abs_sum = 0.0;
+ double scaled_det_abs_sum = 0.0;
+ double scaled_tail_abs_sum = 0.0;
double x0, y0;
int area_scale = 0;
+ int unscaled_cancelled;
+ int scaled_cancelled;
uint32_t i;
if (! pa || pa->npoints < 3 )
@@ -1314,6 +1329,7 @@ ptarray_signed_area(const POINTARRAY *pa)
double *active_compensation = &compensation;
double *active_tail_sum = &tail_sum;
double *active_tail_compensation = &tail_compensation;
+ int scaled_term = LW_FALSE;
P3 = getPoint2d_cp(pa, i);
ptarray_signed_area_two_diff(P2->x, x0, &ax, &axerr);
@@ -1338,6 +1354,7 @@ ptarray_signed_area(const POINTARRAY *pa)
active_compensation = &scaled_compensation;
active_tail_sum = &scaled_tail_sum;
active_tail_compensation = &scaled_tail_compensation;
+ scaled_term = LW_TRUE;
detleft = ax * by;
detright = ay * bx;
det = detleft - detright;
@@ -1364,12 +1381,20 @@ ptarray_signed_area(const POINTARRAY *pa)
*/
err += (ax * byerr + axerr * by) + axerr * byerr;
err -= (ay * bxerr + ayerr * bx) + ayerr * bxerr;
- tail_abs_sum += fabs(err);
+ if (scaled_term)
+ scaled_tail_abs_sum += fabs(err);
+ else
+ tail_abs_sum += fabs(err);
}
ptarray_signed_area_add(active_sum, active_compensation, det);
if (isfinite(det))
- det_abs_sum += fabs(det);
+ {
+ if (scaled_term)
+ scaled_det_abs_sum += fabs(det);
+ else
+ det_abs_sum += fabs(det);
+ }
/*
* Keep determinant tails in a separate compensated stream. If
* they are folded into the main determinant stream immediately,
@@ -1391,9 +1416,12 @@ ptarray_signed_area(const POINTARRAY *pa)
}
if (tail_sum != 0.0 && isfinite(tail_sum))
ptarray_signed_area_add(&sum, &compensation, tail_sum);
- if (area_scale == 0 && tail_abs_sum > 0.0 && fabs(sum) < FP_MAX(det_abs_sum, tail_abs_sum) * 1e-12)
+ unscaled_cancelled = tail_abs_sum > 0.0 && fabs(sum) < FP_MAX(det_abs_sum, tail_abs_sum) * 1e-12;
+ scaled_cancelled = scaled_tail_abs_sum > 0.0 &&
+ fabs(scaled_sum + scaled_tail_sum) < FP_MAX(scaled_det_abs_sum, scaled_tail_abs_sum) * 1e-12;
+ if (unscaled_cancelled || scaled_cancelled)
{
- double exact_area = ptarray_signed_area_exact(pa);
+ double exact_area = ptarray_signed_area_exact(pa, scaled_cancelled ? area_scale : 0);
if (isfinite(exact_area))
return exact_area;
}
diff --git a/regress/core/orientation.sql b/regress/core/orientation.sql
index a21fe4972..54a71e855 100644
--- a/regress/core/orientation.sql
+++ b/regress/core/orientation.sql
@@ -94,3 +94,8 @@ SELECT '436', ST_IsPolygonCCW('POLYGON ((0 0, 10000000000000000 1, 1 1, 10000000
SELECT '437', ST_IsPolygonCW( 'POLYGON ((1.9029580275885475e-13 2.1370106480605854e-89, 2.5163957753386382e-93 -3002966181.933756, 4.77621563620871e-75 -1e-100, -3.7953949389104665e74 -1e200, -2.0800734986321869e-60 2.051992028072436e-55, 3.38586824330807e73 -7.962294473171359e74, 1.9029580275885475e-13 2.1370106480605854e-89))');
SELECT '438', ST_IsPolygonCCW('POLYGON ((1.9029580275885475e-13 2.1370106480605854e-89, 2.5163957753386382e-93 -3002966181.933756, 4.77621563620871e-75 -1e-100, -3.7953949389104665e74 -1e200, -2.0800734986321869e-60 2.051992028072436e-55, 3.38586824330807e73 -7.962294473171359e74, 1.9029580275885475e-13 2.1370106480605854e-89))');
SELECT '439', ST_Area('POLYGON ((1.9029580275885475e-13 2.1370106480605854e-89, 2.5163957753386382e-93 -3002966181.933756, 4.77621563620871e-75 -1e-100, -3.7953949389104665e74 -1e200, -2.0800734986321869e-60 2.051992028072436e-55, 3.38586824330807e73 -7.962294473171359e74, 1.9029580275885475e-13 2.1370106480605854e-89))') > 1e140;
+
+-- Scaled large product tail cancellation with a smaller residual tail
+SELECT '440', ST_IsPolygonCW( 'POLYGON ((1.9029580275885475e47 2.1370106480605854e-29, 2.5163957753386382e-33 -3.002966181933756e69, 4.77621563620871e-15 -1e-40, -3.7953949389104665e134 -1e260, -2.0800734986321869 2.051992028072436e5, 3.38586824330807e133 -7.962294473171359e134, 1.9029580275885475e47 2.1370106480605854e-29))');
+SELECT '441', ST_IsPolygonCCW('POLYGON ((1.9029580275885475e47 2.1370106480605854e-29, 2.5163957753386382e-33 -3.002966181933756e69, 4.77621563620871e-15 -1e-40, -3.7953949389104665e134 -1e260, -2.0800734986321869 2.051992028072436e5, 3.38586824330807e133 -7.962294473171359e134, 1.9029580275885475e47 2.1370106480605854e-29))');
+SELECT '442', ST_Area('POLYGON ((1.9029580275885475e47 2.1370106480605854e-29, 2.5163957753386382e-33 -3.002966181933756e69, 4.77621563620871e-15 -1e-40, -3.7953949389104665e134 -1e260, -2.0800734986321869 2.051992028072436e5, 3.38586824330807e133 -7.962294473171359e134, 1.9029580275885475e47 2.1370106480605854e-29))') > 1e260;
diff --git a/regress/core/orientation_expected b/regress/core/orientation_expected
index 8d93fa66d..060a96b4f 100644
--- a/regress/core/orientation_expected
+++ b/regress/core/orientation_expected
@@ -55,3 +55,6 @@
437|t
438|f
439|t
+440|t
+441|f
+442|t
commit c0521bcc74cd4168f495909c4fa31a1cc4cd956e
Author: Darafei Praliaskouski <me at komzpa.net>
Date: Fri Jun 19 13:46:14 2026 +0400
liblwgeom: preserve signed area residual tails
(cherry picked from commit b220922f136353eda41cf2e43a56f234c7961168)
diff --git a/liblwgeom/cunit/cu_ptarray.c b/liblwgeom/cunit/cu_ptarray.c
index d8c82f89a..c8cac7dad 100644
--- a/liblwgeom/cunit/cu_ptarray.c
+++ b/liblwgeom/cunit/cu_ptarray.c
@@ -353,6 +353,35 @@ static void test_ptarray_signed_area(void)
CU_ASSERT_EQUAL(ptarray_isccw(line->points), 0);
lwline_free(line);
+ /* Preserve residual tails when large product tails cancel. */
+ line =
+ lwgeom_as_lwline(lwgeom_from_text("LINESTRING(1.9029580275885475e-13 2.1370106480605854e-89,"
+ "2.5163957753386382e-93 -3002966181.933756,"
+ "4.77621563620871e-75 -1e-100,"
+ "-3.7953949389104665e74 -1e200,"
+ "-2.0800734986321869e-60 2.051992028072436e-55,"
+ "3.38586824330807e73 -7.962294473171359e74,"
+ "1.9029580275885475e-13 2.1370106480605854e-89)"));
+ area = ptarray_signed_area(line->points);
+ CU_ASSERT(area > 1e140);
+ CU_ASSERT(area < 2e140);
+ CU_ASSERT_EQUAL(ptarray_isccw(line->points), 0);
+ lwline_free(line);
+
+ line =
+ lwgeom_as_lwline(lwgeom_from_text("LINESTRING(1.9029580275885475e-13 2.1370106480605854e-89,"
+ "3.38586824330807e73 -7.962294473171359e74,"
+ "-2.0800734986321869e-60 2.051992028072436e-55,"
+ "-3.7953949389104665e74 -1e200,"
+ "4.77621563620871e-75 -1e-100,"
+ "2.5163957753386382e-93 -3002966181.933756,"
+ "1.9029580275885475e-13 2.1370106480605854e-89)"));
+ area = ptarray_signed_area(line->points);
+ CU_ASSERT(area < -1e140);
+ CU_ASSERT(area > -2e140);
+ CU_ASSERT_EQUAL(ptarray_isccw(line->points), 1);
+ lwline_free(line);
+
/* Preserve product tails when determinant products nearly cancel. */
line =
lwgeom_as_lwline(lwgeom_from_text("LINESTRING(10000000000000276 10000000000000280,"
diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index 4942c9dcb..681f1a79d 100644
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -36,6 +36,14 @@
#include "liblwgeom_internal.h"
#include "lwgeom_log.h"
+#define PTARRAY_SIGNED_AREA_EXPANSION_MAX 4096
+
+typedef struct {
+ uint32_t nterms;
+ int overflowed;
+ double terms[PTARRAY_SIGNED_AREA_EXPANSION_MAX];
+} ptarray_signed_area_expansion;
+
static inline void
ptarray_signed_area_add(double *sum, double *compensation, double value)
{
@@ -59,6 +67,84 @@ ptarray_signed_area_add(double *sum, double *compensation, double value)
*sum = t;
}
+static inline void
+ptarray_signed_area_two_sum(double a, double b, double *sum, double *err)
+{
+ double bvirt, avirt, bround, around;
+
+ *sum = a + b;
+ bvirt = *sum - a;
+ avirt = *sum - bvirt;
+ bround = b - bvirt;
+ around = a - avirt;
+ *err = around + bround;
+}
+
+static inline void
+ptarray_signed_area_expansion_add(ptarray_signed_area_expansion *expansion, double value)
+{
+ uint32_t i, n = 0;
+ double sum, err;
+
+ if (value == 0.0)
+ return;
+ if (!isfinite(value) || expansion->overflowed)
+ {
+ expansion->overflowed = LW_TRUE;
+ return;
+ }
+
+ for (i = 0; i < expansion->nterms; i++)
+ {
+ ptarray_signed_area_two_sum(value, expansion->terms[i], &sum, &err);
+ if (err != 0.0)
+ {
+ if (n == PTARRAY_SIGNED_AREA_EXPANSION_MAX)
+ {
+ expansion->overflowed = LW_TRUE;
+ return;
+ }
+ expansion->terms[n++] = err;
+ }
+ value = sum;
+ }
+
+ if (value != 0.0)
+ {
+ if (n == PTARRAY_SIGNED_AREA_EXPANSION_MAX)
+ {
+ expansion->overflowed = LW_TRUE;
+ return;
+ }
+ expansion->terms[n++] = value;
+ }
+ expansion->nterms = n;
+}
+
+static inline void
+ptarray_signed_area_expansion_add_product(ptarray_signed_area_expansion *expansion, double a, double b, double sign)
+{
+ double product = a * b;
+
+ ptarray_signed_area_expansion_add(expansion, sign * product);
+ ptarray_signed_area_expansion_add(expansion, sign * fma(a, b, -product));
+}
+
+static inline double
+ptarray_signed_area_expansion_sum(ptarray_signed_area_expansion *expansion)
+{
+ double sum = 0.0;
+ double compensation = 0.0;
+ uint32_t i;
+
+ if (expansion->overflowed)
+ return NAN;
+
+ for (i = 0; i < expansion->nterms; i++)
+ ptarray_signed_area_add(&sum, &compensation, expansion->terms[i]);
+ return sum;
+}
+
static inline void
ptarray_signed_area_two_diff(double a, double b, double *diff, double *err)
{
@@ -87,6 +173,42 @@ ptarray_signed_area_product_exponent(double a, double b)
return ea + eb;
}
+static double
+ptarray_signed_area_exact(const POINTARRAY *pa)
+{
+ const POINT2D *P1 = getPoint2d_cp(pa, 0);
+ const POINT2D *P2 = getPoint2d_cp(pa, 1);
+ const POINT2D *P3;
+ ptarray_signed_area_expansion expansion = {0};
+ double x0 = P1->x;
+ double y0 = P1->y;
+ uint32_t i;
+
+ for (i = 2; i < pa->npoints; i++)
+ {
+ double ax, ay, bx, by, axerr, ayerr, bxerr, byerr;
+ P3 = getPoint2d_cp(pa, i);
+
+ ptarray_signed_area_two_diff(P2->x, x0, &ax, &axerr);
+ ptarray_signed_area_two_diff(P2->y, y0, &ay, &ayerr);
+ ptarray_signed_area_two_diff(P3->x, x0, &bx, &bxerr);
+ ptarray_signed_area_two_diff(P3->y, y0, &by, &byerr);
+
+ ptarray_signed_area_expansion_add_product(&expansion, ax, by, 1.0);
+ ptarray_signed_area_expansion_add_product(&expansion, ax, byerr, 1.0);
+ ptarray_signed_area_expansion_add_product(&expansion, axerr, by, 1.0);
+ ptarray_signed_area_expansion_add_product(&expansion, axerr, byerr, 1.0);
+ ptarray_signed_area_expansion_add_product(&expansion, ay, bx, -1.0);
+ ptarray_signed_area_expansion_add_product(&expansion, ay, bxerr, -1.0);
+ ptarray_signed_area_expansion_add_product(&expansion, ayerr, bx, -1.0);
+ ptarray_signed_area_expansion_add_product(&expansion, ayerr, bxerr, -1.0);
+
+ P2 = P3;
+ }
+
+ return -ptarray_signed_area_expansion_sum(&expansion) / 2.0;
+}
+
int
ptarray_has_z(const POINTARRAY *pa)
{
@@ -1143,6 +1265,8 @@ ptarray_signed_area(const POINTARRAY *pa)
double scaled_compensation = 0.0;
double scaled_tail_sum = 0.0;
double scaled_tail_compensation = 0.0;
+ double det_abs_sum = 0.0;
+ double tail_abs_sum = 0.0;
double x0, y0;
int area_scale = 0;
uint32_t i;
@@ -1240,9 +1364,12 @@ ptarray_signed_area(const POINTARRAY *pa)
*/
err += (ax * byerr + axerr * by) + axerr * byerr;
err -= (ay * bxerr + ayerr * bx) + ayerr * bxerr;
+ tail_abs_sum += fabs(err);
}
ptarray_signed_area_add(active_sum, active_compensation, det);
+ if (isfinite(det))
+ det_abs_sum += fabs(det);
/*
* Keep determinant tails in a separate compensated stream. If
* they are folded into the main determinant stream immediately,
@@ -1264,6 +1391,12 @@ ptarray_signed_area(const POINTARRAY *pa)
}
if (tail_sum != 0.0 && isfinite(tail_sum))
ptarray_signed_area_add(&sum, &compensation, tail_sum);
+ if (area_scale == 0 && tail_abs_sum > 0.0 && fabs(sum) < FP_MAX(det_abs_sum, tail_abs_sum) * 1e-12)
+ {
+ double exact_area = ptarray_signed_area_exact(pa);
+ if (isfinite(exact_area))
+ return exact_area;
+ }
return -sum / 2.0;
}
diff --git a/regress/core/orientation.sql b/regress/core/orientation.sql
index 725e9e195..a21fe4972 100644
--- a/regress/core/orientation.sql
+++ b/regress/core/orientation.sql
@@ -89,3 +89,8 @@ SELECT '434', ST_Area('POLYGON ((0 0, 1e300 1e300, 1e-200 0, 0 0))') > 4e99;
-- Large determinant cancellation after a small determinant tail
SELECT '435', ST_IsPolygonCW( 'POLYGON ((0 0, 10000000000000000 1, 1 1, 10000000000000000 0, 0 0))');
SELECT '436', ST_IsPolygonCCW('POLYGON ((0 0, 10000000000000000 1, 1 1, 10000000000000000 0, 0 0))');
+
+-- Large product tail cancellation with a smaller residual tail
+SELECT '437', ST_IsPolygonCW( 'POLYGON ((1.9029580275885475e-13 2.1370106480605854e-89, 2.5163957753386382e-93 -3002966181.933756, 4.77621563620871e-75 -1e-100, -3.7953949389104665e74 -1e200, -2.0800734986321869e-60 2.051992028072436e-55, 3.38586824330807e73 -7.962294473171359e74, 1.9029580275885475e-13 2.1370106480605854e-89))');
+SELECT '438', ST_IsPolygonCCW('POLYGON ((1.9029580275885475e-13 2.1370106480605854e-89, 2.5163957753386382e-93 -3002966181.933756, 4.77621563620871e-75 -1e-100, -3.7953949389104665e74 -1e200, -2.0800734986321869e-60 2.051992028072436e-55, 3.38586824330807e73 -7.962294473171359e74, 1.9029580275885475e-13 2.1370106480605854e-89))');
+SELECT '439', ST_Area('POLYGON ((1.9029580275885475e-13 2.1370106480605854e-89, 2.5163957753386382e-93 -3002966181.933756, 4.77621563620871e-75 -1e-100, -3.7953949389104665e74 -1e200, -2.0800734986321869e-60 2.051992028072436e-55, 3.38586824330807e73 -7.962294473171359e74, 1.9029580275885475e-13 2.1370106480605854e-89))') > 1e140;
diff --git a/regress/core/orientation_expected b/regress/core/orientation_expected
index 8502182ff..8d93fa66d 100644
--- a/regress/core/orientation_expected
+++ b/regress/core/orientation_expected
@@ -52,3 +52,6 @@
434|t
435|t
436|f
+437|t
+438|f
+439|t
commit 0e16d8991f1356c0b735d6ab3046e036a7e08740
Author: Darafei Praliaskouski <me at komzpa.net>
Date: Wed Jun 17 00:36:20 2026 +0400
liblwgeom: guard compensated area overflow
(cherry picked from commit 0bae360bbefd0d4ce87898d027fe808201ea0a8d)
diff --git a/NEWS b/NEWS
index 1a7f8db0a..c4a2ebec8 100644
--- a/NEWS
+++ b/NEWS
@@ -171,6 +171,8 @@ To take advantage of all postgis_sfcgal extension features SFCGAL 2.3+ is needed
(Darafei Praliaskouski)
- GH-895, Qualify pg_class lookup during SRID/PROJ validation
(Darafei Praliaskouski)
+ - GH-898, Fix ST_DFullyWithin indexed plans and large-coordinate
+ ring-orientation precision regressions (Darafei Praliaskouski)
- #4828, geometry_columns handles NOT VALID SRID checks without errors (Darafei Praliaskouski)
- #6048, [raster] ST_Clip no longer crashes when clipping sparse band
selections (Darafei Praliaskouski)
diff --git a/liblwgeom/cunit/cu_ptarray.c b/liblwgeom/cunit/cu_ptarray.c
index c906d466e..d8c82f89a 100644
--- a/liblwgeom/cunit/cu_ptarray.c
+++ b/liblwgeom/cunit/cu_ptarray.c
@@ -5,6 +5,7 @@
*
* Copyright (C) 2011 Sandro Santilli <strk at kbt.io>
* Copyright (C) 2008 Paul Ramsey
+ * Copyright (C) 2026 Darafei Praliaskouski <me at komzpa.net>
*
* This is free software; you can redistribute and/or modify it under
* the terms of the GNU General Public Licence. See the COPYING file.
@@ -14,6 +15,7 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
+#include <math.h>
#include "CUnit/Basic.h"
#include "CUnit/CUnit.h"
@@ -343,6 +345,82 @@ static void test_ptarray_signed_area(void)
area = ptarray_signed_area(line->points);
ASSERT_DOUBLE_EQUAL_TOLERANCE(area, 0.5, 0.0000001);
lwline_free(line);
+
+ /* Preserve determinant tails when later large determinants cancel. */
+ line = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(0 0,10000000000000000 1,1 1,10000000000000000 0,0 0)"));
+ area = ptarray_signed_area(line->points);
+ ASSERT_DOUBLE_EQUAL_TOLERANCE(area, 0.5, 0.0000001);
+ CU_ASSERT_EQUAL(ptarray_isccw(line->points), 0);
+ lwline_free(line);
+
+ /* Preserve product tails when determinant products nearly cancel. */
+ line =
+ lwgeom_as_lwline(lwgeom_from_text("LINESTRING(10000000000000276 10000000000000280,"
+ "59999999999999680 59999999999999672,"
+ "89999999999999360 89999999999999360,"
+ "90000000000000112 90000000000000128,"
+ "90000000000000704 90000000000000704,"
+ "10000000000000276 10000000000000280)"));
+ area = ptarray_signed_area(line->points);
+ CU_ASSERT(area < 0.0);
+ CU_ASSERT_EQUAL(ptarray_isccw(line->points), 1);
+ lwline_free(line);
+
+ /* Preserve subtraction tails when the origin shift hides unit-scale offsets. */
+ line =
+ lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-10000000000000000 -10000000000000000,"
+ "1 -1,1 0,1 1,-1 0,"
+ "-10000000000000000 -10000000000000000)"));
+ area = ptarray_signed_area(line->points);
+ CU_ASSERT(area < -1e16);
+ CU_ASSERT(area > -2e16);
+ CU_ASSERT_EQUAL(ptarray_isccw(line->points), 1);
+ lwline_free(line);
+
+ /* Overflowed determinants should remain infinite, not poison the sum with NaN. */
+ line = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(0 0,0 1e200,1e200 1e200,1e200 0,0 0)"));
+ area = ptarray_signed_area(line->points);
+ CU_ASSERT(isinf(area) && area > 0);
+ CU_ASSERT_EQUAL(ptarray_isccw(line->points), 0);
+ lwline_free(line);
+
+ /* Rescale overflowing fan products when the determinant itself is finite. */
+ line =
+ lwgeom_as_lwline(lwgeom_from_text("LINESTRING(0 0,1 1e155,1e155 1e155,"
+ "1e155 1.0000000000000001e155,"
+ "1 1.0000000000000001e155,0 0)"));
+ area = ptarray_signed_area(line->points);
+ CU_ASSERT(isfinite(area));
+ CU_ASSERT(area < -1e294);
+ CU_ASSERT(area > -2e294);
+ CU_ASSERT_EQUAL(ptarray_isccw(line->points), 1);
+ lwline_free(line);
+
+ line =
+ lwgeom_as_lwline(lwgeom_from_text("LINESTRING(0 0,1 1.0000000000000001e155,"
+ "1e155 1.0000000000000001e155,"
+ "1e155 1e155,1 1e155,0 0)"));
+ area = ptarray_signed_area(line->points);
+ CU_ASSERT(isfinite(area));
+ CU_ASSERT(area > 1e294);
+ CU_ASSERT(area < 2e294);
+ CU_ASSERT_EQUAL(ptarray_isccw(line->points), 0);
+ lwline_free(line);
+
+ /* Preserve tiny finite products beside unrelated huge coordinates. */
+ line = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(0 0,1e-200 0,1e300 1e300,0 0)"));
+ area = ptarray_signed_area(line->points);
+ CU_ASSERT(area < -4e99);
+ CU_ASSERT(area > -6e99);
+ CU_ASSERT_EQUAL(ptarray_isccw(line->points), 1);
+ lwline_free(line);
+
+ line = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(0 0,1e300 1e300,1e-200 0,0 0)"));
+ area = ptarray_signed_area(line->points);
+ CU_ASSERT(area > 4e99);
+ CU_ASSERT(area < 6e99);
+ CU_ASSERT_EQUAL(ptarray_isccw(line->points), 0);
+ lwline_free(line);
}
static void test_ptarray_contains_point(void)
diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index 4cae95038..4942c9dcb 100644
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -20,19 +20,73 @@
*
* Copyright (C) 2012-2021 Sandro Santilli <strk at kbt.io>
* Copyright (C) 2001-2006 Refractions Research Inc.
- * Copyright 2026 Darafei Praliaskouski
+ * Copyright (C) 2026 Darafei Praliaskouski <me at komzpa.net>
*
**********************************************************************/
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
+#include <float.h>
+#include <limits.h>
+#include <math.h>
#include "../postgis_config.h"
/*#define POSTGIS_DEBUG_LEVEL 4*/
#include "liblwgeom_internal.h"
#include "lwgeom_log.h"
+static inline void
+ptarray_signed_area_add(double *sum, double *compensation, double value)
+{
+ double y, t;
+
+ /*
+ * Kahan compensation is only valid in finite arithmetic. When a finite
+ * ring is so large that the determinant overflows, preserve the IEEE
+ * infinity instead of manufacturing NaN in the compensation term.
+ */
+ if (!isfinite(value) || !isfinite(*sum))
+ {
+ *sum += value;
+ *compensation = 0.0;
+ return;
+ }
+
+ y = value - *compensation;
+ t = *sum + y;
+ *compensation = (t - *sum) - y;
+ *sum = t;
+}
+
+static inline void
+ptarray_signed_area_two_diff(double a, double b, double *diff, double *err)
+{
+ double bvirt, avirt, bround, around;
+
+ *diff = a - b;
+ bvirt = a - *diff;
+ avirt = *diff + bvirt;
+ bround = bvirt - b;
+ around = a - avirt;
+ *err = around + bround;
+}
+
+static inline int
+ptarray_signed_area_product_exponent(double a, double b)
+{
+ int ea, eb;
+
+ if (!isfinite(a) || !isfinite(b))
+ return DBL_MAX_EXP;
+ if (a == 0.0 || b == 0.0)
+ return INT_MIN;
+
+ frexp(fabs(a), &ea);
+ frexp(fabs(b), &eb);
+ return ea + eb;
+}
+
int
ptarray_has_z(const POINTARRAY *pa)
{
@@ -1083,7 +1137,14 @@ ptarray_signed_area(const POINTARRAY *pa)
const POINT2D *P3;
double sum = 0.0;
double compensation = 0.0;
+ double tail_sum = 0.0;
+ double tail_compensation = 0.0;
+ double scaled_sum = 0.0;
+ double scaled_compensation = 0.0;
+ double scaled_tail_sum = 0.0;
+ double scaled_tail_compensation = 0.0;
double x0, y0;
+ int area_scale = 0;
uint32_t i;
if (! pa || pa->npoints < 3 )
@@ -1093,36 +1154,116 @@ ptarray_signed_area(const POINTARRAY *pa)
P2 = getPoint2d_cp(pa, 1);
x0 = P1->x;
y0 = P1->y;
+ /*
+ * Use one binary area scale for overflowing fan triangles. Some skinny
+ * rings have finite total area even though individual fan triangles
+ * overflow; scaling each overflowing term independently would preserve
+ * those infinities and still lose the final cancellation. Finite raw
+ * triangles are kept in the unscaled sum so unrelated tiny determinants do
+ * not underflow just because another triangle needed scaling.
+ */
+ for (i = 1; i < pa->npoints; i++)
+ {
+ double ax, ay, bx, by, axerr, ayerr, bxerr, byerr;
+ int product_exponent;
+ P3 = getPoint2d_cp(pa, i);
+ ptarray_signed_area_two_diff(P2->x, x0, &ax, &axerr);
+ ptarray_signed_area_two_diff(P2->y, y0, &ay, &ayerr);
+ ptarray_signed_area_two_diff(P3->x, x0, &bx, &bxerr);
+ ptarray_signed_area_two_diff(P3->y, y0, &by, &byerr);
+ product_exponent =
+ FP_MAX(ptarray_signed_area_product_exponent(ax, by), ptarray_signed_area_product_exponent(ay, bx));
+ if (product_exponent > DBL_MAX_EXP - 24)
+ {
+ int coordinate_scale = (product_exponent - (DBL_MAX_EXP - 24) + 1) / 2;
+ area_scale = FP_MAX(area_scale, coordinate_scale);
+ }
+ P2 = P3;
+ }
+
+ P2 = getPoint2d_cp(pa, 1);
for ( i = 2; i < pa->npoints; i++ )
{
- double ax, ay, bx, by;
- double detleft, detright, det, err, b, y, t;
+ double ax, ay, bx, by, axerr, ayerr, bxerr, byerr;
+ double detleft, detright, det, err = 0.0, b;
+ double *active_sum = ∑
+ double *active_compensation = &compensation;
+ double *active_tail_sum = &tail_sum;
+ double *active_tail_compensation = &tail_compensation;
P3 = getPoint2d_cp(pa, i);
- ax = P2->x - x0;
- ay = P2->y - y0;
- bx = P3->x - x0;
- by = P3->y - y0;
+ ptarray_signed_area_two_diff(P2->x, x0, &ax, &axerr);
+ ptarray_signed_area_two_diff(P2->y, y0, &ay, &ayerr);
+ ptarray_signed_area_two_diff(P3->x, x0, &bx, &bxerr);
+ ptarray_signed_area_two_diff(P3->y, y0, &by, &byerr);
detleft = ax * by;
detright = ay * bx;
det = detleft - detright;
- b = det - detleft;
- err = (detleft - (det - b)) - (detright + b);
+ if (area_scale > 0 && (!isfinite(detleft) || !isfinite(detright) || !isfinite(det)))
+ {
+ ax = ldexp(ax, -area_scale);
+ ay = ldexp(ay, -area_scale);
+ bx = ldexp(bx, -area_scale);
+ by = ldexp(by, -area_scale);
+ axerr = ldexp(axerr, -area_scale);
+ ayerr = ldexp(ayerr, -area_scale);
+ bxerr = ldexp(bxerr, -area_scale);
+ byerr = ldexp(byerr, -area_scale);
+ active_sum = &scaled_sum;
+ active_compensation = &scaled_compensation;
+ active_tail_sum = &scaled_tail_sum;
+ active_tail_compensation = &scaled_tail_compensation;
+ detleft = ax * by;
+ detright = ay * bx;
+ det = detleft - detright;
+ }
+ /*
+ * Keep the determinant rounding error, including the product
+ * tails, and sum it with Kahan compensation. Near-collinear
+ * large-coordinate rings can otherwise lose the sign in the
+ * rounded products before the subtraction tail is recovered.
+ * The error tail is meaningful only for finite products; with
+ * overflow it can turn an infinite determinant into NaN through
+ * Inf-Inf cancellation.
+ */
+ if (isfinite(detleft) && isfinite(detright) && isfinite(det))
+ {
+ b = det - detleft;
+ err = (detleft - (det - b)) - (detright + b);
+ err += fma(ax, by, -detleft) - fma(ay, bx, -detright);
+ /*
+ * The origin shift itself can round away unit-scale
+ * offsets beside 1e16-scale coordinates. Include the
+ * low-order subtraction terms in the determinant so the
+ * compensated sum still sees their signed area.
+ */
+ err += (ax * byerr + axerr * by) + axerr * byerr;
+ err -= (ay * bxerr + ayerr * bx) + ayerr * bxerr;
+ }
- y = det - compensation;
- t = sum + y;
- compensation = (t - sum) - y;
- sum = t;
-
- y = err - compensation;
- t = sum + y;
- compensation = (t - sum) - y;
- sum = t;
+ ptarray_signed_area_add(active_sum, active_compensation, det);
+ /*
+ * Keep determinant tails in a separate compensated stream. If
+ * they are folded into the main determinant stream immediately,
+ * a later cancellation between large determinants can erase the
+ * small residual that fixes the exact double-coordinate sign.
+ */
+ ptarray_signed_area_add(active_tail_sum, active_tail_compensation, err);
/* Move forwards! */
P2 = P3;
}
+ if (scaled_sum != 0.0)
+ ptarray_signed_area_add(&sum, &compensation, ldexp(scaled_sum, 2 * area_scale));
+ if (scaled_tail_sum != 0.0)
+ {
+ double scaled_tail = ldexp(scaled_tail_sum, 2 * area_scale);
+ if (isfinite(scaled_tail))
+ ptarray_signed_area_add(&tail_sum, &tail_compensation, scaled_tail);
+ }
+ if (tail_sum != 0.0 && isfinite(tail_sum))
+ ptarray_signed_area_add(&sum, &compensation, tail_sum);
return -sum / 2.0;
}
diff --git a/regress/core/orientation.sql b/regress/core/orientation.sql
index 1ee927fd5..725e9e195 100644
--- a/regress/core/orientation.sql
+++ b/regress/core/orientation.sql
@@ -62,3 +62,30 @@ SELECT '416', ST_IsPolygonCCW('MULTIPOLYGON (((0 0, 1 1, 1 0, 0 0)), ((100 0, 10
SELECT '417', ST_IsPolygonCW( 'MULTIPOLYGON (((0 0, 1 0, 1 1, 0 0)), ((100 0, 101 1, 101 0, 100 0)))');
SELECT '418', ST_IsPolygonCCW('MULTIPOLYGON (((0 0, 1 0, 1 1, 0 0)), ((100 0, 101 1, 101 0, 100 0)))');
+-- Near-collinear exterior ring with large coordinates
+SELECT '419', ST_IsPolygonCW( 'POLYGON ((10000000000000276 10000000000000280, 59999999999999680 59999999999999672, 89999999999999360 89999999999999360, 90000000000000112 90000000000000128, 90000000000000704 90000000000000704, 10000000000000276 10000000000000280))');
+SELECT '420', ST_IsPolygonCCW('POLYGON ((10000000000000276 10000000000000280, 59999999999999680 59999999999999672, 89999999999999360 89999999999999360, 90000000000000112 90000000000000128, 90000000000000704 90000000000000704, 10000000000000276 10000000000000280))');
+
+-- Large origin shift with unit-scale offsets
+SELECT '421', ST_IsPolygonCW( 'POLYGON ((-10000000000000000 -10000000000000000, 1 -1, 1 0, 1 1, -1 0, -10000000000000000 -10000000000000000))');
+SELECT '422', ST_IsPolygonCCW('POLYGON ((-10000000000000000 -10000000000000000, 1 -1, 1 0, 1 1, -1 0, -10000000000000000 -10000000000000000))');
+
+-- Skinny finite-area ring with individually overflowing fan products
+SELECT '423', ST_IsPolygonCW( 'POLYGON ((0 0, 1 1e155, 1e155 1e155, 1e155 1.0000000000000001e155, 1 1.0000000000000001e155, 0 0))');
+SELECT '424', ST_IsPolygonCCW('POLYGON ((0 0, 1 1e155, 1e155 1e155, 1e155 1.0000000000000001e155, 1 1.0000000000000001e155, 0 0))');
+SELECT '425', ST_IsPolygonCW( 'POLYGON ((0 0, 1 1.0000000000000001e155, 1e155 1.0000000000000001e155, 1e155 1e155, 1 1e155, 0 0))');
+SELECT '426', ST_IsPolygonCCW('POLYGON ((0 0, 1 1.0000000000000001e155, 1e155 1.0000000000000001e155, 1e155 1e155, 1 1e155, 0 0))');
+SELECT '427', ST_Area('POLYGON ((0 0, 1 1e155, 1e155 1e155, 1e155 1.0000000000000001e155, 1 1.0000000000000001e155, 0 0))') > 1e294;
+SELECT '428', ST_Area('POLYGON ((0 0, 1 1.0000000000000001e155, 1e155 1.0000000000000001e155, 1e155 1e155, 1 1e155, 0 0))') > 1e294;
+
+-- Tiny finite determinant product beside huge coordinates
+SELECT '429', ST_IsPolygonCW( 'POLYGON ((0 0, 1e-200 0, 1e300 1e300, 0 0))');
+SELECT '430', ST_IsPolygonCCW('POLYGON ((0 0, 1e-200 0, 1e300 1e300, 0 0))');
+SELECT '431', ST_Area('POLYGON ((0 0, 1e-200 0, 1e300 1e300, 0 0))') > 4e99;
+SELECT '432', ST_IsPolygonCW( 'POLYGON ((0 0, 1e300 1e300, 1e-200 0, 0 0))');
+SELECT '433', ST_IsPolygonCCW('POLYGON ((0 0, 1e300 1e300, 1e-200 0, 0 0))');
+SELECT '434', ST_Area('POLYGON ((0 0, 1e300 1e300, 1e-200 0, 0 0))') > 4e99;
+
+-- Large determinant cancellation after a small determinant tail
+SELECT '435', ST_IsPolygonCW( 'POLYGON ((0 0, 10000000000000000 1, 1 1, 10000000000000000 0, 0 0))');
+SELECT '436', ST_IsPolygonCCW('POLYGON ((0 0, 10000000000000000 1, 1 1, 10000000000000000 0, 0 0))');
diff --git a/regress/core/orientation_expected b/regress/core/orientation_expected
index cb8789ff0..8502182ff 100644
--- a/regress/core/orientation_expected
+++ b/regress/core/orientation_expected
@@ -34,3 +34,21 @@
416|f
417|f
418|f
+419|f
+420|t
+421|f
+422|t
+423|f
+424|t
+425|t
+426|f
+427|t
+428|t
+429|f
+430|t
+431|t
+432|t
+433|f
+434|t
+435|t
+436|f
commit 5336ccf1fe35d63672357a11f0082f4b70c0bbf8
Author: Darafei Praliaskouski <me at komzpa.net>
Date: Wed Jun 17 00:12:21 2026 +0400
geometry: avoid extended precision area dependency
(cherry picked from commit 7162829dbea13fb3acc72bbee2884abf26e09994)
diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index 3091efc06..4cae95038 100644
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -1081,9 +1081,9 @@ ptarray_signed_area(const POINTARRAY *pa)
const POINT2D *P1;
const POINT2D *P2;
const POINT2D *P3;
- long double sum = 0.0;
- double x0;
- long double x, y1, y2;
+ double sum = 0.0;
+ double compensation = 0.0;
+ double x0, y0;
uint32_t i;
if (! pa || pa->npoints < 3 )
@@ -1092,19 +1092,38 @@ ptarray_signed_area(const POINTARRAY *pa)
P1 = getPoint2d_cp(pa, 0);
P2 = getPoint2d_cp(pa, 1);
x0 = P1->x;
+ y0 = P1->y;
for ( i = 2; i < pa->npoints; i++ )
{
+ double ax, ay, bx, by;
+ double detleft, detright, det, err, b, y, t;
P3 = getPoint2d_cp(pa, i);
- x = (long double)P2->x - x0;
- y1 = P3->y;
- y2 = P1->y;
- sum += x * (y2-y1);
+
+ ax = P2->x - x0;
+ ay = P2->y - y0;
+ bx = P3->x - x0;
+ by = P3->y - y0;
+
+ detleft = ax * by;
+ detright = ay * bx;
+ det = detleft - detright;
+ b = det - detleft;
+ err = (detleft - (det - b)) - (detright + b);
+
+ y = det - compensation;
+ t = sum + y;
+ compensation = (t - sum) - y;
+ sum = t;
+
+ y = err - compensation;
+ t = sum + y;
+ compensation = (t - sum) - y;
+ sum = t;
/* Move forwards! */
- P1 = P2;
P2 = P3;
}
- return (double)(sum / 2.0);
+ return -sum / 2.0;
}
int
commit 5124ad4a50a4a914e4b188240997ee3f06e258d7
Author: Darafei Praliaskouski <me at komzpa.net>
Date: Tue Jun 16 23:34:07 2026 +0400
geometry: fix planner and precision regressions
(cherry picked from commit 2609ea702dd59fcdfa6331530b3eda5fa5f55f69)
diff --git a/liblwgeom/cunit/cu_ptarray.c b/liblwgeom/cunit/cu_ptarray.c
index 8cf3d8c11..c906d466e 100644
--- a/liblwgeom/cunit/cu_ptarray.c
+++ b/liblwgeom/cunit/cu_ptarray.c
@@ -338,6 +338,11 @@ static void test_ptarray_signed_area(void)
ASSERT_DOUBLE_EQUAL_TOLERANCE(area, -4.0, 0.0000001);
lwline_free(line);
+ /* Preserve small residuals through cancellation in the area sum. */
+ line = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(0 0,10000000000000000 0,1 -1,10000000000000000 -1,0 0)"));
+ area = ptarray_signed_area(line->points);
+ ASSERT_DOUBLE_EQUAL_TOLERANCE(area, 0.5, 0.0000001);
+ lwline_free(line);
}
static void test_ptarray_contains_point(void)
diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index 3baae4da6..3091efc06 100644
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -1081,8 +1081,9 @@ ptarray_signed_area(const POINTARRAY *pa)
const POINT2D *P1;
const POINT2D *P2;
const POINT2D *P3;
- double sum = 0.0;
- double x0, x, y1, y2;
+ long double sum = 0.0;
+ double x0;
+ long double x, y1, y2;
uint32_t i;
if (! pa || pa->npoints < 3 )
@@ -1094,7 +1095,7 @@ ptarray_signed_area(const POINTARRAY *pa)
for ( i = 2; i < pa->npoints; i++ )
{
P3 = getPoint2d_cp(pa, i);
- x = P2->x - x0;
+ x = (long double)P2->x - x0;
y1 = P3->y;
y2 = P1->y;
sum += x * (y2-y1);
@@ -1103,7 +1104,7 @@ ptarray_signed_area(const POINTARRAY *pa)
P1 = P2;
P2 = P3;
}
- return sum / 2.0;
+ return (double)(sum / 2.0);
}
int
diff --git a/postgis/gserialized_supportfn.c b/postgis/gserialized_supportfn.c
index b6fe6dbba..a7c0e1136 100644
--- a/postgis/gserialized_supportfn.c
+++ b/postgis/gserialized_supportfn.c
@@ -433,10 +433,16 @@ Datum postgis_index_supportfn(PG_FUNCTION_ARGS)
rightdatatype = exprType(rightarg);
/*
- * Given the index operator family and the arguments and the
- * desired strategy number we can now lookup the operator
- * we want (usually && or &&&).
- */
+ * Given the index operator family and the arguments and the
+ * desired strategy number we can now lookup the operator
+ * we want (usually && or &&&). ST_DFullyWithin is directional:
+ * for an index on the first argument, overlap with the expanded
+ * second argument is the safe lossy condition. The stronger
+ * contained-by condition is only safe for an index on the second
+ * argument, which must be covered by the expanded first argument.
+ */
+ if (idxfn.index == ST_DFULLYWITHIN_IDX && req->indexarg == 0)
+ idxfn.index = ST_DWITHIN_IDX;
oproid = get_opfamily_member(opfamilyoid,
leftdatatype,
rightdatatype,
diff --git a/regress/core/tickets.sql b/regress/core/tickets.sql
index c8a6e69cb..83f954139 100644
--- a/regress/core/tickets.sql
+++ b/regress/core/tickets.sql
@@ -1612,9 +1612,11 @@ SELECT '#5747', ST_Length('MULTISURFACE (((0 0, 1 0, 1 1, 0 1, 0 0)), CURVEPOLYG
-- #5855
CREATE TEMP TABLE TEST (street text, extent geometry(Polygon,32633));
INSERT INTO test VALUES ('Knosesmauet','0103000020797F0000010000000500000010B2468761BDDFC06390523AA1B0594110B2468761BDDFC030554D9BC3B0594107992AF50799DFC030554D9BC3B0594107992AF50799DFC06390523AA1B0594110B2468761BDDFC06390523AA1B05941');
-SET enable_seqscan=false;
CREATE INDEX test_idx ON test USING GIST ( extent );
+-- These checks protect index-support rewrites, so keep sequential scans
+-- disabled only around the #5855 queries that must exercise GiST paths.
+SET enable_seqscan = off;
SELECT '#5855', street
FROM test
WHERE ST_DFullyWithin(
@@ -1623,6 +1625,21 @@ WHERE ST_DFullyWithin(
1700
);
+CREATE TEMP TABLE dfullywithin_indexed_first_arg (id integer, geom geometry);
+INSERT INTO dfullywithin_indexed_first_arg
+VALUES (1, 'POLYGON((-10 -10,-10 10,10 10,10 -10,-10 -10))');
+CREATE INDEX dfullywithin_indexed_first_arg_gist
+ON dfullywithin_indexed_first_arg USING GIST (geom);
+
+SELECT '#5855.1', id
+FROM dfullywithin_indexed_first_arg
+WHERE ST_DFullyWithin(
+ geom,
+ ST_GeomFromText('POINT(0 0)'),
+ 0
+);
+RESET enable_seqscan;
+
SELECT '#5876', ST_AsText(ST_AddPoint(
'LINESTRING (1 1, 2 2)'::geometry,
'POINT EMPTY'::geometry), 2);
diff --git a/regress/core/tickets_expected b/regress/core/tickets_expected
index 639472965..bed2d3302 100644
--- a/regress/core/tickets_expected
+++ b/regress/core/tickets_expected
@@ -492,6 +492,7 @@ ERROR: Geometry contains invalid coordinates
#5686|0
#5747|0
#5855|Knosesmauet
+#5855.1|1
#5876|LINESTRING(1 1,2 2)
#5082|t
public.test5978.shape SRID:4326 TYPE:POINT DIMS:2
commit 6fff4b58bbf6b671fdb34a0bbaa70d1f8ab11b2c
Author: Darafei Praliaskouski <me at komzpa.net>
Date: Mon Jul 6 11:30:38 2026 +0400
topology: preserve large ids from edge wrappers
(cherry picked from commit 3553756154863a8eebb1c56af7f56e10ca91f7b6)
Closes https://github.com/postgis/postgis/pull/1153
diff --git a/topology/postgis_topology.c b/topology/postgis_topology.c
index 4bbf220ae..41491d080 100644
--- a/topology/postgis_topology.c
+++ b/topology/postgis_topology.c
@@ -4726,7 +4726,7 @@ Datum ST_RemEdgeModFace(PG_FUNCTION_ARGS)
{
text* toponame_text;
char* toponame;
- int ret;
+ LWT_ELEMID ret;
LWT_ELEMID node_id;
LWT_TOPOLOGY *topo;
@@ -4781,7 +4781,7 @@ Datum ST_RemEdgeNewFace(PG_FUNCTION_ARGS)
{
text* toponame_text;
char* toponame;
- int ret;
+ LWT_ELEMID ret;
LWT_ELEMID node_id;
LWT_TOPOLOGY *topo;
@@ -4834,7 +4834,7 @@ Datum ST_ModEdgeHeal(PG_FUNCTION_ARGS)
{
text* toponame_text;
char* toponame;
- int ret;
+ LWT_ELEMID ret;
LWT_ELEMID eid1, eid2;
LWT_TOPOLOGY *topo;
@@ -4888,7 +4888,7 @@ Datum ST_NewEdgeHeal(PG_FUNCTION_ARGS)
{
text* toponame_text;
char* toponame;
- int ret;
+ LWT_ELEMID ret;
LWT_ELEMID eid1, eid2;
LWT_TOPOLOGY *topo;
commit 6ea8421662ed55c8b895e1a6e4d4ef3edefde4c6
Author: Darafei Praliaskouski <me at komzpa.net>
Date: Mon Jul 6 11:30:02 2026 +0400
topology: fix corrupt column fallback variables
(cherry picked from commit 3c6f0e146badc6528c17f7f1a666cf33c71cb694)
Closes https://github.com/postgis/postgis/pull/1149
diff --git a/topology/sql/manage/FixCorruptTopoGeometryColumn.sql.in b/topology/sql/manage/FixCorruptTopoGeometryColumn.sql.in
index 9c1ca04ec..b385d4686 100644
--- a/topology/sql/manage/FixCorruptTopoGeometryColumn.sql.in
+++ b/topology/sql/manage/FixCorruptTopoGeometryColumn.sql.in
@@ -72,10 +72,10 @@ WHERE pg_type.typname = 'topogeometry' AND pga.attname = 'id'
l.feature_type
)::topogeometry
FROM topology.layer AS l
- WHERE l.topology_id = (%3$I).topology_id AND l.layer_id = (%3$I).layer_id AND (%3$I).type <> l.feature_type ', topo_schema, topo_table, topo_column);
+ WHERE l.topology_id = (%3$I).topology_id AND l.layer_id = (%3$I).layer_id AND (%3$I).type <> l.feature_type ', layerSchema, layerTable, layerColumn);
EXECUTE var_sql;
GET DIAGNOSTICS var_row_count = ROW_COUNT;
- result = result || format('%s rows updated for %s.%s.%s column back to integer id type', var_row_count, topo_schema, topo_table, topo_column);
+ result = result || format('%s rows updated for %s.%s.%s column back to integer id type', var_row_count, layerSchema, layerTable, layerColumn);
END IF;
RETURN result;
END
commit 8a03193ebc452f335c0f51c26dbdbb53b1a15904
Author: Darafei Praliaskouski <me at komzpa.net>
Date: Mon Jul 6 11:27:02 2026 +0400
utils: create upgrade-check tempdir safely
(cherry picked from commit e1fc3de47c29025a232738f78da8d330f2281984)
Closes https://github.com/postgis/postgis/pull/1148
diff --git a/utils/check_create_upgrade_replaced_function_errors.sh b/utils/check_create_upgrade_replaced_function_errors.sh
index caa5d375d..df9ebe224 100644
--- a/utils/check_create_upgrade_replaced_function_errors.sh
+++ b/utils/check_create_upgrade_replaced_function_errors.sh
@@ -2,11 +2,9 @@
set -e
-TMPDIR="${TMPDIR:-/tmp}/check_create_upgrade_replaced_function_errors.$$"
+TMPDIR="$(mktemp -d "${TMPDIR:-/tmp}/check_create_upgrade_replaced_function_errors.XXXXXX")"
trap 'rm -rf "${TMPDIR}"' EXIT
-mkdir -p "${TMPDIR}"
-
cat > "${TMPDIR}/postgis.sql" <<'SQL'
-- INSTALL VERSION: 3.7.0dev
-- Availability: 1.0.0
commit d9525f032bba08eefdce60e55f333c101dace23d
Author: Darafei Praliaskouski <me at komzpa.net>
Date: Mon Jul 6 11:26:47 2026 +0400
doc: fix localized EPUB help target
(cherry picked from commit af647934be4c73ed3b4a9a0d1553d0094f7d2c75)
Closes https://github.com/postgis/postgis/pull/1147
diff --git a/doc/Makefile.in b/doc/Makefile.in
index cfeafac25..5288648c6 100644
--- a/doc/Makefile.in
+++ b/doc/Makefile.in
@@ -487,7 +487,7 @@ epub-uninstall:
rm -f $(DESTDIR)$(docdir)/postgis-${POSTGIS_MAJOR_VERSION}.${POSTGIS_MINOR_VERSION}.${POSTGIS_MICRO_VERSION}$(DOCSUFFIX).epub
epub: postgis-${POSTGIS_MAJOR_VERSION}.${POSTGIS_MINOR_VERSION}.${POSTGIS_MICRO_VERSION}$(DOCSUFFIX).epub ## Generate postgis-<version>-en.epub
-pdf-localized: ## Generate a postgis-<version>-<lang>.epub for supported languages
+epub-localized: ## Generate a postgis-<version>-<lang>.epub for supported languages
pdf: postgis-${POSTGIS_MAJOR_VERSION}.${POSTGIS_MINOR_VERSION}.${POSTGIS_MICRO_VERSION}$(DOCSUFFIX).pdf ## Generate postgis-<version>-en.pdf
pdf-localized: ## Generate a postgis-<version>-<lang>.pdf for supported languages
commit 9000ef67ed38616b10e06d8da7bd012f9d866aa0
Author: Darafei Praliaskouski <me at komzpa.net>
Date: Mon Jul 6 11:29:02 2026 +0400
topology: clean up TopoRingIsCCW early returns
(cherry picked from commit 770ace4aeb699a2e55885f2893c60ec30eb5c008)
Closes https://github.com/postgis/postgis/pull/1143
diff --git a/topology/postgis_topology.c b/topology/postgis_topology.c
index 2eb19bb67..4bbf220ae 100644
--- a/topology/postgis_topology.c
+++ b/topology/postgis_topology.c
@@ -5917,7 +5917,9 @@ Datum TopoRingIsCCW(PG_FUNCTION_ARGS)
if ( lwgeom_is_empty(lwgeom) )
{
- PG_RETURN_BOOL(false);
+ lwgeom_free(lwgeom);
+ PG_FREE_IF_COPY(geom, 0);
+ PG_RETURN_BOOL(false);
}
if (lwgeom->type == POLYGONTYPE)
@@ -5930,8 +5932,10 @@ Datum TopoRingIsCCW(PG_FUNCTION_ARGS)
}
else
{
- lwpgerror("Unsupported geometry type passed to TopoRingIsCCW");
- PG_RETURN_NULL();
+ lwgeom_free(lwgeom);
+ PG_FREE_IF_COPY(geom, 0);
+ lwpgerror("Unsupported geometry type passed to TopoRingIsCCW");
+ PG_RETURN_NULL();
}
isCCW = lwt_IsTopoRingCCW(pa);
-----------------------------------------------------------------------
Summary of changes:
NEWS | 5 +
doc/Makefile.in | 2 +-
liblwgeom/cunit/cu_ptarray.c | 164 +++++++++
liblwgeom/ptarray.c | 383 ++++++++++++++++++++-
postgis/gserialized_supportfn.c | 14 +-
regress/core/orientation.sql | 42 +++
regress/core/orientation_expected | 27 ++
regress/core/tickets.sql | 19 +-
regress/core/tickets_expected | 1 +
topology/postgis_topology.c | 18 +-
.../sql/manage/FixCorruptTopoGeometryColumn.sql.in | 4 +-
...heck_create_upgrade_replaced_function_errors.sh | 4 +-
12 files changed, 657 insertions(+), 26 deletions(-)
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list