My handling of the cyclicity was
incorrect above! My sincere apologies. 😢
The definition of a cyclic barycenter among 32 regularly placed weights (be they bits or histogram bins) \$b_n\$, \$n = 0 \dots 31\$, is
by definition$$\left\lbrace \, \begin{aligned}
x &= \sum_{n=0}^{31} b_n \cos\left(\frac{n \pi}{16}\right) \\
y &= \sum_{n=0}^{31} b_n \sin\left(\frac{n \pi}{16}\right) \\
c &= \left\lfloor 32 + \frac{16}{\pi}\operatorname{atan2}(y, x) \right\rceil \mod 32 \\
\end{aligned} \right .$$
where \$c = 0 \dots 31\$ is the barycenter index.
Consider this as an unit circle, with 32 possible "weights" placed regularly on it, at 360°/32 = 11.25° intervals. We use standard 2D barycenter for these 32 points: the weights are either 0 or 1 for binary, or bin counts for histograms. The answer is the 11.25° sector from origin the barycenter ends up in. (The distance from origin to barycenter is not relevant, although very close to origin means the barycenter is not well defined.)
Thus, one very efficient implementation in C would be
cyclic_barycenter(bits),
// SPDX-License-Identifier: CC0-1.0
#include <stdint.h>
// First quadrant arctangent for 32 directions; returns 0 .. 8.
static uint_fast8_t atan32_internal(uint_fast16_t y, uint_fast16_t x)
{
static const struct {
uint16_t x;
uint16_t y;
} points[8] = {
{ .x = 26581, .y = 2618 },
{ .x = 31703, .y = 9617 },
{ .x = 9113, .y = 4871 },
{ .x = 5716, .y = 4691 },
{ .x = 23963, .y = 29199 },
{ .x = 15914, .y = 29773 },
{ .x = 2121, .y = 6992 },
{ .x = 1025, .y = 10407 },
};
uint32_t ppos, pneg;
int_fast8_t imin = 0, imax = 8;
while (1) {
int_fast8_t i = (imin + imax) / 2;
ppos = (uint32_t)y * points[i].x;
pneg = (uint32_t)x * points[i].y;
if (i == imin)
return i + (ppos >= pneg);
if (ppos > pneg)
imin = i;
else
if (ppos < pneg)
imax = i;
else
return i + 1;
}
}
// Verified for x,y in -32768..32767 inclusive:
// atan32(y, x) == (31 & (int)round(atan2((double)y, (double)x) * 16/PI))
uint_fast8_t atan32(int_fast16_t y, int_fast16_t x)
{
if (y > 0) {
if (x > 0) {
return 31 & ( atan32_internal(y, x));
} else
if (x < 0) {
return 31 & (16 - atan32_internal(y, -x));
} else {
return 8;
}
} else
if (y < 0) {
if (x > 0) {
return 31 & (32 - atan32_internal(-y, x));
} else
if (x < 0) {
return 31 & (16 + atan32_internal(-y, -x));
} else {
return 24;
}
} else {
if (x > 0) {
return 0;
} else
if (x < 0) {
return 16;
} else {
return 0; // (0, 0)
}
}
}
int cyclic_barycenter(uint32_t u)
{
// Scale 3209 minimizes error while keeping maximum sum within signed 16-bit numbers
static const int16_t c[40] = {
3147, 2965, 2668, 2269, 1783, 1228, 626, 0, -626, -1228, -1783, -2269, -2668, -2965, -3147, -3209,
-3147, -2965, -2668, -2269, -1783, -1228, -626, 0, 626, 1228, 1783, 2269, 2668, 2965, 3147, 3209,
3147, 2965, 2668, 2269, 1783, 1228, 626, 0,
};
int32_t x = 0, y = 0;
uint_fast8_t i = 32;
while ((i--)|u) {
if (u & 1) {
x += c[i];
y += c[i+8];
}
u >>= 1;
}
return atan32(y, x);
}
Only one loop over the 32 bits are done; then, three conditionals, followed by a binary search over 8 entries. The binary search does exactly 4 iterations, with each iteration doing two 16×16=32-bit (unsigned) multiplications, for a total of 8 16×16=32-bit multiplications. This is why I claim this to be efficient.
Note that the cost is not much different when a histogram is used, although then I would switch to 32-bit signed integers and lookup (and 32×32=64-bit unsigned multiply operations).
Two look-up tables are used. One contains 40 signed 16-bit integers, and the other 8 pairs of signed 16-bit integers, for a total of 56 signed 16-bit integers, or 112 bytes of read-only memory.
I have verified that
atan32(y,x) matches
31 & (int)round(16/3.14159265358979323846*atan2(y,x)) exactly, for all 16-bit signed
x and
y, -32768 .. 32767.
atan32_internal(y,x) takes a point in the positive quadrant, and returns the angle by comparing the direction to eight known points closest to but within the lower boundary of each 11.5° sector, using a binary search and exterior product. Essentially, if \$(x_A, y_A)\$ and \$(x_B, y_B)\$ are two vectors, \$y_A x_B - x_A y_B > 0\$ if the former vector is at a larger (positive) angle than the latter vector. The stored points are at the lower (positive) angle boundary of the 11.5° sector.
cyclic_barycenter(bits) uses a precalculated table of sines and cosines in a funky order: cos(0°)=c[31], sin(0°)=c[39], cos(11.5°)=c[30], sin(11.5°)=c[38], and so on. The scale is 3209, because it minimizes relative errors in the individual integer variables (compared to the exact values of sine/cosine), while keeping the maximum sum combination within signed 16-bit range, -32768..+32767 (in -32581..+32581, actually). Essentially, this minimizes the quantization error in the barycenter coordinates, while keeping the barycenter within 16-bit signed integer range.
I've only lightly tested it, but it seems to produce expected results, within the precision limits.