leijurv left a comment (osm2pgsql-dev/osm2pgsql#2503)
I apologize I have to correct my above comment, the paragraph beginning `If the
Lua gives a double which is stored in PG as a real, currently text mode
implicitly does fp64 -> text -> fp32.`, where I claimed this didn't practically
happen, that was incorrect.
The exact scenario, where this PR causes a different result than master, is
when a float64 is sent to Postgres and its mantissa takes the following form:
`(1)???????????????????????10000000000000000000000000000`
That's the implicit leading 1 bit, then 23 irrelevant bits, then a one, then 28
zeros. The odds of this happening to a uniformly random bits would be 1 in
2^29, which is why I (approximately) claimed as such above.
I also "verified" it by ensuring that the fp64->text and text->fp32 was
consistent across all polygons in the planet:
```
gis=# select count(*) from planet_osm_polygon;
count
-----------
854755728
(1 row)
gis=# select count(*) from planet_osm_polygon where st_area(way)::text::real =
st_area(way)::real;
count
-----------
854755728
(1 row)
```
Note that I can't check it on `way_area` because this rounding has already
happened so these bits have been lost. So that's why I used `st_area(way)`. So
this is why I thought we had no practical differences from master to this PR.
Unfortunately, that was misleading, because PostGIS computes areas differently
to Lua (which uses boost).
```
gis=# select st_area(way), st_area(way)::real, way_area from planet_osm_polygon
where osm_id=258622451;
st_area | st_area | way_area
--------------------+-----------+-----------
236.50345608387994 | 236.50345 | 236.50346
(1 row)
```
That st_area has the following mantissa bits:
```
(1)1101100100000001110001001111111011101110011110010010
```
But, for this particular way, Lua gives us `236.50345611572266`, which has a
mantissa in binary of:
```
(1)1101100100000001110001010000000000000000000000000000
```
Note exactly 28 trailing zeros indeed.
```
gis=# select double precision '236.50345611572266'::real, double precision
'236.50345611572266'::text::real;
float4 | float4
-----------+-----------
236.50345 | 236.50346
(1 row)
```
And, oh no, it does indeed round differently, we ended up precisely on the
boundary between two float32s, we get `236.50345` in binary mode versus
`236.50346` in text mode. πΏπΏπΏ
The explanation here is that we have way_area as
`object:as_polygon():transform(3857):area()` which internally calls
`boost::geometry::area()`. That uses
[shoelace](https://en.wikipedia.org/wiki/Shoelace_formula) to compute the area,
and a lot of precision is lost to cancellation as the coordinates here are
already rather large (x β β8.4e6 and y β 4.7e6 metres). We usually don't care
about this because we're going to downcast to a fp32 anyway. The relevant code
can be found
[here](https://github.com/boostorg/geometry/blob/boost-1.83.0/include/boost/geometry/strategy/cartesian/area.hpp#L92-L112).
They just compute `sum += (x1 + x2) * (y1 - y2)`. However, PostGIS uses a
simple trick: they subtract out the first coordinate, to recenter all the math
around zero, like: `x = P2->x - x0; sum += x * (y2 - y1)`. Relevant code:
[here](https://github.com/postgis/postgis/blob/3.6.4/liblwgeom/lwpoly.c#L433-L459)
and
[here](https://github.com/postgis/postgis/blob/3.6.4/liblwgeom/ptarray.c#L1151-L1180).
The exact line where they center at zero is subtracting out `x0`
[here](https://github.com/postgis/postgis/blob/3.6.4/liblwgeom/ptarray.c#L1170).
Due to this loss of precision, this rare case of 28 trailing zeros happens
appears to be about 1 in 2,000 rather than 1 in 500,000,000.
Here's a histogram of how many trailing bits of the mantissa are zero, across
all polygons [in
Delaware](https://download.geofabrik.de/north-america/us/delaware.html).
<img width="1600" height="880" alt="way_area_mantissa_bits"
src="https://github.com/user-attachments/assets/15231f4d-bedc-4490-8665-a8b266aaa94d";
/>
The highlighted column alone is what has the difference in rounding (any more
or fewer trailing zeros rounds equally).
**To fix this** and exactly match master, we could check for this case on every
fp64 we send to postgres, and fall back to `strtof(fmt::to_string(d))` if it
matches this bit pattern. Or, we could accept that way_areas will be rounded a
bit differently now. I'm happy either way!
--
Reply to this email directly or view it on GitHub:
https://github.com/osm2pgsql-dev/osm2pgsql/pull/2503#issuecomment-5849223880
You are receiving this because you are subscribed to this thread.
Message ID: <osm2pgsql-dev/osm2pgsql/pull/2503/[email protected]>_______________________________________________
Tile-serving mailing list
[email protected]
https://lists.openstreetmap.org/listinfo/tile-serving
