I've spent some time reading the source code of silk and figuring out what's going on in advance.py and advance.hpp. I thought I'd share in case others are interested. My starting point is the following paraphrased comments by APG, from Discord:
No cell in a still life can have over 6 nbrs, so we can represent the known info about stable state with 8 bitplanes: d0, d1, d2, l2, l3, d4, d5, and d6. d/l = dead/live, n = on nbrs. These bitplanes are stored inverted: 1 = impossible and 0 = possible. We refer to a cell as low/stable/high if current - stable is -1/0/1 respectively. i.e. low = on in stable, off in current; high = opposite. From the low/stable/high info about a cell's 8 neighbors, we compute 18 booleans, 2 per vertical line below
-6 -5 | -4 | -3 | -2 | -1 | 0 | +1 | +2 | +3 | +4 +5 +6 +7 +8
specifying whether or not we know actual difference is left/right of the line. We can then deduce info about out the next generation based on these 18 booleans. (We don't need to distinguish between -5 and -6: at most 6 nbrs in stable, and 1, 0 < 2 hence underpopulated in the current generation. Likewise, if the difference is at least +4, the cell is certainly overpopulated.)
That's the high-level overview: low/high/stable info for current gen => 18 booleans per cell => low/high/stable info for next gen. What follows is the details of how one actually goes from one to the other.
In the code and in what follows, there's certain naming conventions: f = forced, d/l = dead/live, us = unstable. not_{low/high/stable} indicates forced to be not {low/high/stable} [as opposed to possibly-not-{low/high/stable}]. Also, silk stores the 64x64 torus in 4 quadrants: the kth warp stores the kth row in each square, as the components of a uint4, a vector of 4 uint32_t's. So keep in mind that this isn't just one uint4, but one per warp.
First, going from low/high/stable for current gen to the 18 booleans. The difference of between the current population and the stable population is the number of high neighbors, minus the number of low neighbors. For brevity, I'll omit the "neighbors" in what follows.
curr_pop - stable_pop = high - low
min(high - low) <= curr_pop - stable_pop <= max(high - low)
min(high) - max(low) <= curr_pop - stable_pop <= max(high) - min(low)
min(high) + min(8 - low) <= 8 + curr_pop - stable_pop <= max(high) + max(8-low).
Here, I added 8 to both sides, and used the fact that maximizing c - x is the same as minimizing x, and vice versa. This is nice because 8-low is the size of the complement of low. In the code, we see
Code: Select all
uint4 lb_bin = kc::sum16<false>(not_low, not_low & not_stable);
uint4 ub_bin = kc::sum16<true>(not_high, not_high & not_stable);
kc::sum16 here "returns the sum of the 16 booleans obtained by taking the values of a and b on each of the 8 neighbours." The "true" here means "subtract one." lb_bin is the sum over the 8 neighbors of (must_be_high + must_be_not_low), the same lower bound I derived above. The upper bound I derived above is might_be_high + might_be_not_low = 16 - (must_be_not_high + must_be_low). Therefore ub_bin = -1 + must_be_not_high + must_be_low = 15 - (upper bound above).
Next, we need to compute the 18 booleans from lb_bin and ub_bin. Relevant comments from the advance.hpp:
Code: Select all
// lb[i] := "we know that (curr_pop - stable_pop >= i - 4)"
uint32_t lb0 = lb_bin.w | lb_bin.z;
...compute lb0-lb8...
// ub[i] := "we know that (curr_pop - stable_pop < i - 4)"
uint32_t ub8 = ub_bin.w | ub_bin.z;
...compute ub0-ub8...
For the lb{i}s, lb_bin <= 8 + curr_pop - stable_pop, hence our lower bound for curr_pop - stable_pop is lb_bin - 8. lb{i} is true iff curr_pop-stable_pop >= i - 4 iff lb_bin - 8 >= i - 4 iff lb_bin >= i + 4. And if you reference that against the omitted binary arithmetic, that does make sense: the digits are w, z, y, x, most significant to least. E.g. lb0 means lb_bin >= 4 which means either w or z is true.
For the ub{i}'s, 8 + curr_pop - stable_pop <= 15 - ub_bin < 16 - ub_bin, hence our strict upper bound for curr_pop - stable_pop is 8 - ub_bin. ub{i} is true iff curr_pop - stable_pop < i - 4 iff 8 - ub_bin <= i - 4 iff ub_bin >= 12 - i. Again, this seems to match the binary arithmetic.
That finishes going from the high/stable/low to the 16 booleans.
Second, going from the 18 booleans to which cells are forced to be on/off/stable in the next generation. Consider some cell that might be of a given state, among d0, d1, d2, l2, l3, d4, d5, and d6. (Potential confusion: in what follows, "state" means "one of d0, d1, etc" not off/on.) Let n be its number of live neighbors. If that cell is indeed of that state, then it is:
Off in next gen <=> (n < 2) | (n < 3 & dead) | (n >= 4)
On in next gen <=> (n >= 2) & (n >= 3 | live) & (n < 4)
unstable in next <=> {on/off} in next where we're considering a {d/l} state
To get the overall forced off/on/unstable, intersect over all possible states for the cell. In advance.hpp, this looks like
Code: Select all
{ // d0
uint32_t fd_next = (ub6 | (ub7 & not_high) | lb8) | ad0;
uint32_t fl_next = (lb6 & (lb7 | not_stable) & ub8) | ad0;
unstable_d0 = fl_next;
gnot_high &= fd_next;
gnot_stable &= fl_next;
}
...
{ // l2
uint32_t fd_next = (ub4 | (ub5 & not_stable) | lb6) | al2;
uint32_t fl_next = (lb4 & (lb5 | not_low) & ub6) | al2;
unstable_l2 = fd_next;
gnot_stable &= fd_next;
gnot_low &= fl_next;
}
...
Here, the "| a{state}" corresponds to "possible state for the cell." (Keep in mind that the a{state} bitplanes are inverted: 1 = impossible.) A 1 in fd_next means not quite "forced dead next," but rather "forced dead next or forced not-of-{state}." gnot_{low/high/stable} start out as all 1's: the &='s remove the cells that are possibly-of-{state}-and-{low/high/stable}, by setting that bit of gnot_{low/high/stable} to 0. The relationship between gnot_{low/high/stable} and f{l/d}_next is different for l states vs d states. E.g. d = off in stable, so on = unstable = high, off = stable = not-high.
The rest of advance.hpp is fairly straightforward: compute which cells are forced stable based off search parameters, and mark those as forced-not-high and forced-not-low. Cells that are forced to be stable cannot be of a state that makes them unstable (this is why we calculated unstable_d0). Check for improvements to knowledge about stable catalyst by comparing new state bit plane to old one. Finally, update the 8 state bit planes, as well as not_low, not_high, and not_stable. Return T/F for improved info/not.
Lengthy post. Please like if you made it to the end and found it useful: writing this all up took some time. Also, APG/mvr/others more familiar with the code, if I'm wrong or missing something, feel free to chime in.