spot_finding: accumulate spot centroids in integers

The photon-weighted position sums were floats, so the centroid's last bit
depended on the build rather than on the data: gcc contracts the multiply-add
in AddPixel into an FMA under -march=x86-64-v3 and cannot at the baseline,
and MSVC does not contract at all under /fp:precise. The GPU extractor had to
match with __fmaf_rn, and the parity test still needed a two-ulp slack for
hosts that do not fuse.

Column, line and the per-pixel count are all integral, so the sums are exact
in int64 and both implementations reach the same bits with nothing to match.
The parity test now demands exact equality unconditionally and gets it,
including on a baseline build.

ConvertToImageCoordinates keeps the sums integral too: the raw -> image map
is a signed axis swap plus an integer translation, so it is applied to the
sums instead of to the centroid.

Drops the SpotToSave constructor, which had no callers and could not have
been converted without quantising the stored centroid.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
2026-08-02 19:03:35 +02:00
co-authored by Claude Opus 5
parent cd16053c2b
commit 8f1b0b2281
5 changed files with 56 additions and 67 deletions
+5 -25
View File
@@ -195,26 +195,11 @@ public:
}
};
// How many representable floats apart two values are.
int64_t UlpDistance(float a, float b) {
int32_t ia, ib;
memcpy(&ia, &a, sizeof(ia));
memcpy(&ib, &b, sizeof(ib));
if (ia < 0) ia = INT32_MIN - ia; // map to a monotone ordering across the sign
if (ib < 0) ib = INT32_MIN - ib;
return std::abs(static_cast<int64_t>(ia) - static_cast<int64_t>(ib));
}
// Everything that decides which spots exist and what they weigh is compared EXACTLY: the number of
// spots, their order, and each one's pixel count, photon sum and maximum. Those are integers, and a
// difference in any of them is a difference in the partition.
//
// The centroid is a float sum, and its last bit is a property of the BUILD rather than of either
// implementation: gcc contracts DiffractionSpot::AddPixel into an FMA whenever the flags allow it
// (the CI -march=x86-64-v3 does), while a baseline -march, or MSVC with its default /fp:precise,
// cannot. The extractor uses __fmaf_rn, so it is bit-exact against a host that fuses and one ulp
// off one that does not. An ulp bound catches a real divergence - which moves a centroid by pixels,
// not by 1e-4 of one - while staying true whichever way the host was built.
// Everything is compared EXACTLY: the number of spots, their order, and each one's pixel count,
// photon sum, maximum and centroid. All of them are built out of integer sums on both sides, so
// there is nothing here that a build flag can move - which is the point of accumulating in integers
// rather than in float, where gcc contracted the multiply-add under -march=x86-64-v3 and not at the
// baseline and left the last bit of the centroid a property of how the host was built.
void RequireIdentical(const std::string &what,
const std::vector<DiffractionSpot> &cpu,
const std::vector<DiffractionSpot> &gpu) {
@@ -228,13 +213,8 @@ void RequireIdentical(const std::string &what,
// RawCoord divides the sums by the photon count, so comparing it compares the sums; a spot
// whose photons sum to zero reports (0,0) on both sides by the same branch.
const Coord c = cpu[i].RawCoord(), g = gpu[i].RawCoord();
#ifdef __FMA__
REQUIRE(memcmp(&g.x, &c.x, sizeof(float)) == 0);
REQUIRE(memcmp(&g.y, &c.y, sizeof(float)) == 0);
#else
REQUIRE(UlpDistance(g.x, c.x) <= 2);
REQUIRE(UlpDistance(g.y, c.y) <= 2);
#endif
}
}