# Extend the most recent "nice" OEIS sequence: stepping stone puzzle on a grid

Today Neil Sloane of the OEIS sent out an email asking for a confirmation of the current terms, and computation of some larger terms of the latest OEIS sequence A337663 with the keyword "nice".

Here's how this sequence works: You label $$\n\$$ cells on the infinite square grid with $$\1\$$s, and then

place the numbers $$\2,3,4,\dots,m\$$ in order, subject to the rule that when you place $$\k\$$, the sum of its [horizontal, vertical, and diagonal] neighbors must equal $$\k\$$. Then $$\A337663(n)\$$ is the maximum $$\m\$$ that can be achieved over all initial placings of the $$\1\$$-cells.

Here's an illustration of $$\A337663(2) = 16\$$:

+----+----+----+----+----+----+
|  9 |  5 | 10 | 11 |    |    |
+----+----+----+----+----+----+
|    |  4 |  1 |    |    |    |
+----+----+----+----+----+----+
| 12 |  8 |  3 |  2 |    | 16 |
+----+----+----+----+----+----+
|    |    |    |  6 |  1 | 15 |
+----+----+----+----+----+----+
|    |    | 13 |  7 | 14 |    |
+----+----+----+----+----+----+


Notice that the $$\2\$$-cell has two $$\1\$$-cells as neighbors (and $$\1 + 1 = 2\$$); the $$\3\$$-cell has a $$\1\$$-cell and a $$\2\$$-cell as neighbors (and $$\1 + 2 = 3\$$); etc.

# The challenge.

Like this previous challenge, the goal of this is to compute as many terms as possible in this sequence, which begins 1, 16, 28, 38 and where the $$\n\$$-th term is $$\A337663(n)\$$.

Run your code for as long as you'd like. The winner of this challenge will be the user who posts the most terms of the sequence, along with their code to generate it. If two users post the same number of terms, then whoever posts their last term earliest wins.

• If there is no time limit on this challenge, how is the winning condition objective? Oct 7, 2020 at 21:06
• @LiefdeWen Why would it not be objective? You either compute a new term of the sequence or you do not. Oct 7, 2020 at 23:26
• @Sisyphus You're right, soz Oct 8, 2020 at 8:03
• Here is an interesting question. Let D(i) be the smallest Manhattan distance between the i-th 1 and any other 1. What is the maximum value of D(i) possible? So far the largest I have found is 7 in the optimal N=5 solution: oeis.org/A337663/a337663.txt Can we get distances larger than 7? Oct 24, 2020 at 1:36
• This sequence is the subject of the latest Numberphile video: youtube.com/watch?v=m4Uth-EaTZ8 Jan 12 at 22:29

# C + OpenMP, N=5

New version with a proper handling of the 112-1113 case.

a(5) = 49
0 46 26  0  0  0  0  0  0  0  0 35  0

0 20  0  6 28 48  0  0  0  0 34  1 36

39 19  1  2  3 17  0 30  0  0 33  0 37

0  0 18  7  1  4  9  0 21 32  0  0  0

0 40  0  8 38  5 43 10 11  0 44  0  0

0  0 22  0 13  0 15  0  1 12  0  0  0

47 23  0 14 27  0 31 16 29  0  0  0  0

0 24  1  0 41  0  0  0 45  0  0  0  0

49 25  0 42  0  0  0  0  0  0  0  0  0


### How it works

The program will only work for N=5, for higher numbers you'd need some adjustments. First lets take a look how an easier approach for N=4 would look like. We need at least 112 next to each other in some arrangement. Because there are only two 1s left, every other number cannot be made only by new 1s.

So starting from the six possible starting positions for 112:

1 1   1 2 1   1 _ 1   1 _   1 _ _   1 _ _
2 _           _ 2 _   2 1   _ 2 1   _ 2 _
_ _ 1


we can take a look at every space placed two spots away and check their sum (Note: with some proper case handling, you should be fine to check the direct neighbors, though I took the safe route).

      0 0 0 0 0 0
0 1 2 2 1 0


1 1 -> 0 3 . . 1 0 2 _ 0 3 . 4 1 0 0 2 2 2 0 0 0 0 0 0 0 0

For every spot: check if the sum is the next needed number (in this case 3) or if we can still place some 1s: is the sum plus some newly added 1s the next needed number. In the latter case, we need to make sure that the new 1s don't interfere with existing numbers > 1, e.g.

3 1
1 1 1
2


wouldn't be valid as the 2-placement would have been illegal, but

1 1
2   3 1
1


would be fine. Note that I only increase the bounding box for two spots around non-1 numbers. So for the lower right corner, the next spots to try are as following:

1 _ _ _
_ 3 1 _
_ 1 _ _
_ _ _ _
x


The x spot wouldn't get checked, as its number would only neighbor new 1s – and for N=4 this is not possible as mentioned before.

For N>4 this get a little more complicated: it is not guaranteed that every number will be connected to the first 112. Another cluster might start independently: 1113. But after that every number cannot be made only of new 1s, thus will be connected to either 1113 or 112. Note that we don't have to handle anything else in the N=5 case (but would need for N>5): having two clusters with 1 and 11114 will already be handled, as 2 and 3 must also be placed in 11114; so every 11114 will already be checked by 112 or 1113.

So we need to get a bounding box to find out how close 112 and 1113 can be placed. For this we run two boards that cannot touch, scoring them by the sum of the distances they managed to get away from the starting position. This is the best they manage:

  0  0  0  0  0  0  0  0  0

0  0  0  0  0  0  0  0  0

0  0  0  0 15  0  0  0  0

0  0 11 10  5  0  0  0  0

0  0  0  1  4 12  0  0  0

0  0  0  0  2  1 13  0  0

0  0  0  0  0  0 14  0  0

0  0  0  0  0  0  0  0  0

0  0  0  0  0  0  0  0  0

…

0  0  0  0  0  0  0  0

0  0  0  0  0  0  0  0

0  0 17  9  0  3  1  0

0  0  0  8  1  6  1  0

0  0  0 16  7  0  0  0

0  0  0  0  0  0  0  0

0  0  0  0  0  0  0  0


They cannot leave 5 tiles. So if we place the initial 3 within a 20x20 (+ a padding of 4 for off-by-one-errors :-)) field centered around the 2, we either get two disconnected clusters that have a score independent on where they are exactly, or two clusters that eventually will join up. So anything up to

1 1 _ _ _ _ _ _ _ _ _ _ _ 1 1
_ 2 a b c d e _ e d c b a 3 1


will be checked with 11 spaces in between; enough that they cannot meet up.

With all this, then just recursively try out all the possibilities in a depth-first-search. Always modifying only one board, we only need memory for a(N) recursive steps.

OMP is only used to check the initial boards in parallel. This is far from a balanced workload; the final position needs about twice as long as the others. However, it is the easiest to implement. :-)

### The program

Compiled with clang -O3 -o main main.c -fopenmp and ran with time OMP_NUM_THREADS=4 ./main.

#include <stdint.h>
#include <stddef.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>

typedef uint8_t mem_t;
typedef uint16_t sum_t;
#define S 64
const int startx = S/2, starty = S/2;
// for N > 5, adjust the unrolled loops in step
#define MAX_CELL 5
#define MAX_BOARDS 2

#define MAX(x,y) (x>y?x:y)
#define MIN(x,y) (x<y?x:y)

const int ys[8] = {0,1,1,1,0,-1,-1,-1};
const int xs[8] = {1,1,0,-1,-1,-1,0,1};

inline
void add_sum(sum_t v, int y, int x, sum_t sum[S][S]) {
for(int d=0;d<8;d++)
sum[y+ys[d]][x+xs[d]] += v;
}
inline
void add_placed(sum_t v, int y, int x, mem_t placed[S][S]) {
for(int d=0;d<8;d++)
placed[y+ys[d]][x+xs[d]] += v;
}

typedef struct board {
int y0, y1, x0, x1;
mem_t b[S][S], placed[S][S];
sum_t sum[S][S];
} board_t;

void st_print(int c, int max, board_t *b) {
printf("%d cells, %d max\n", c, max);
for(int y=b->y0;y<=b->y1;y++){
for(int x=b->x0;x<=b->x1;x++) printf("%*d", 3, b->b[y][x]);
puts("\n");
}
}

void step(int c, mem_t max, board_t *bs, int bl, mem_t *best_max, board_t best_b[MAX_BOARDS], int optimize_spread) {
// check board size
for(int i=0;i<bl;i++) {
if (bs[i].y0 < 2 || bs[i].y1 >= S - 2 || bs[i].x0 < 2 || bs[i].x1 >= S - 2) {
st_print(c, max, &bs[i]);
printf("board too small %d %d %d %d", bs[i].y0, bs[i].y1, bs[i].x0, bs[i].x1);
exit(1);
}
}
// new best
if (c == MAX_CELL) {
int score = 0;
for (int i=0;i<bl;i++)
score += MAX(starty - bs[i].y0,
MAX(bs[i].y1 - starty,
MAX(startx - bs[i].x0,
bs[i].x1 - startx)));
} else {
score = max;
}
if (*best_max < score) {
for (int i=0;i<bl;i++)
memcpy(&best_b[i], &bs[i], sizeof(board_t));
*best_max = score;
}
}
// place with 0 new 1-cells
for(int i=0;i<bl;i++) {
board_t *b=bs+i;
for(int y=b->y0;y<=b->y1;y++)
for(int x=b->x0;x<=b->x1;x++)
if(b->sum[y][x] == max + 1 && !b->b[y][x]) {
b->b[y][x] = max + 1;
int y0o = b->y0, y1o = b->y1, x0o = b->x0, x1o = b->x1;
b->y0 = MIN(b->y0, y-2);
b->y1 = MAX(b->y1, y+2);
b->x0 = MIN(b->x0, x-2);
b->x1 = MAX(b->x1, x+2);
step(c, max + 1, bs, bl, best_max, best_b, optimize_spread);
b->y0 = y0o, b->y1 = y1o, b->x0 = x0o, b->x1 = x1o;
b->b[y][x] = 0;
}
}
// sorry for the repetition, couldn't get clang to optimize it otherwise
// place with 1 new 1-cells
if(c + 1 <= MAX_CELL)
for(int i=0;i<bl;i++) {
board_t *b=bs+i;
for(int y=b->y0;y<=b->y1;y++)
for(int x=b->x0;x<=b->x1;x++)
if(b->sum[y][x] == (max + 1) - 1 && !b->b[y][x]) {
for(int d1=0;d1<8;d1++) {
if (b->placed[y+ys[d1]][x+xs[d1]]) continue;

b->b[y+ys[d1]][x+xs[d1]] = 1;
b->b[y][x] = max + 1;
int y0o = b->y0, y1o = b->y1, x0o = b->x0, x1o = b->x1;
b->y0 = MIN(b->y0, y-2);
b->y1 = MAX(b->y1, y+2);
b->x0 = MIN(b->x0, x-2);
b->x1 = MAX(b->x1, x+2);
step(c + 1, max + 1, bs, bl, best_max, best_b, optimize_spread);
b->y0 = y0o, b->y1 = y1o, b->x0 = x0o, b->x1 = x1o;
b->b[y+ys[d1]][x+xs[d1]] = 0;
b->b[y][x] = 0;
}
}
}
// place with 2 new 1-cells
if(c + 2 <= MAX_CELL)
for(int i=0;i<bl;i++) {
board_t *b=bs+i;
for(int y=b->y0;y<=b->y1;y++)
for(int x=b->x0;x<=b->x1;x++)
if(b->sum[y][x] == (max + 1) - 2 && !b->b[y][x]) {
for(int d1=0;d1<8-1;d1++) {
if (b->placed[y+ys[d1]][x+xs[d1]]) continue;
for(int d2=d1+1;d2<8;d2++) {
if (b->placed[y+ys[d2]][x+xs[d2]]) continue;

b->b[y+ys[d1]][x+xs[d1]] = 1;
b->b[y+ys[d2]][x+xs[d2]] = 1;
b->b[y][x] = max + 1;
int y0o = b->y0, y1o = b->y1, x0o = b->x0, x1o = b->x1;
b->y0 = MIN(b->y0, y-2);
b->y1 = MAX(b->y1, y+2);
b->x0 = MIN(b->x0, x-2);
b->x1 = MAX(b->x1, x+2);
step(c + 2, max + 1, bs, bl, best_max, best_b, optimize_spread);
b->y0 = y0o, b->y1 = y1o, b->x0 = x0o, b->x1 = x1o;
b->b[y+ys[d1]][x+xs[d1]] = 0;
b->b[y+ys[d2]][x+xs[d2]] = 0;
b->b[y][x] = 0;
}
}
}
}
// place with 3 new 1-cells
if(c + 3 <= MAX_CELL)
for(int i=(optimize_spread && max == 2);i<bl;i++) {
board_t *b=bs+i;
for(int y=b->y0;y<=b->y1;y++)
for(int x=b->x0;x<=b->x1;x++)
if(b->sum[y][x] == (max + 1) - 3 && !b->b[y][x]) {
for(int d1=0;d1<8-2;d1++) {
if (b->placed[y+ys[d1]][x+xs[d1]]) continue;
for(int d2=d1+1;d2<8-1;d2++) {
if (b->placed[y+ys[d2]][x+xs[d2]]) continue;
for(int d3=d2+1;d3<8;d3++) {
if (b->placed[y+ys[d3]][x+xs[d3]]) continue;

b->b[y+ys[d1]][x+xs[d1]] = 1;
b->b[y+ys[d2]][x+xs[d2]] = 1;
b->b[y+ys[d3]][x+xs[d3]] = 1;
b->b[y][x] = max + 1;
int y0o = b->y0, y1o = b->y1, x0o = b->x0, x1o = b->x1;
b->y0 = MIN(b->y0, y-2);
b->y1 = MAX(b->y1, y+2);
b->x0 = MIN(b->x0, x-2);
b->x1 = MAX(b->x1, x+2);
step(c + 3, max + 1, bs, bl, best_max, best_b, optimize_spread);
b->y0 = y0o, b->y1 = y1o, b->x0 = x0o, b->x1 = x1o;
b->b[y+ys[d1]][x+xs[d1]] = 0;
b->b[y+ys[d2]][x+xs[d2]] = 0;
b->b[y+ys[d3]][x+xs[d3]] = 0;
b->b[y][x] = 0;
}
}
}
}
}
}

void set_starting_board(board_t* b, int i) {
int x0 = startx, y0 = starty;
b->b[y0][x0] = 2;
if (i == 0) b->b[y0-1][x0-1] = 1,
b->b[y0+1][x0+1] = 1;
if (i == 1) b->b[y0-1][x0-1] = 1,
b->b[y0][x0+1] = 1;
if (i == 2) b->b[y0][x0-1] = 1,
b->b[y0][x0+1] = 1;
if (i == 3) b->b[y0-1][x0] = 1,
b->b[y0][x0+1] = 1;
if (i == 4) b->b[y0-1][x0-1] = 1,
b->b[y0-1][x0+1] = 1;
if (i == 5) b->b[y0-1][x0] = 1,
b->b[y0-1][x0+1] = 1;

for(int y=1;y+1<S;y++)
for(int x=1;x+1<S;x++)
for(int yd=-1;yd<=1;yd++)
for(int xd=-1;xd<=1;xd++)
if(yd!=0||xd!=0)
b->sum[y][x] += b->b[y+yd][x+xd];
for(int y=1;y+1<S;y++)
for(int x=1;x+1<S;x++)
for(int yd=-1;yd<=1;yd++)
for(int xd=-1;xd<=1;xd++)
b->placed[y][x] += b->b[y+yd][x+xd] > 1;

}

int get_bounding_box() {
int x0 = startx, y0 = starty;
board_t best_b[6][3] = {0};
mem_t best_max[6] = {0};

#pragma omp parallel for
for(int i=0;i<6;i++) {
board_t bs[] = {(board_t){y0 - 3, y0 + 3, x0 - 3, x0 + 3, {0}, {0}, {0}},
(board_t){y0, y0, x0, x0, {0}, {0}, {0}}};
set_starting_board(&bs[0], i);
step(2, 2, bs, 2, &best_max[i], best_b[i], 1);
}

int best_i=0, mm = 0;
for(int i=0;i<6;i++)
if (best_max[i] > mm)
mm = best_max[i],
best_i = i;
printf("most spread of distant 112 and 1113: %d\n", best_max[best_i]);
st_print(MAX_CELL, best_max[best_i], &best_b[best_i][0]);
st_print(MAX_CELL, best_max[best_i], &best_b[best_i][1]);

return best_max[best_i] + 4;
}

int main(int argc, char **argv) {
int bb = get_bounding_box();

int x0 = startx, y0 = starty;
board_t best_b[6][3] = {0};
mem_t best_max[6] = {0};

#pragma omp parallel for
for(int i=0;i<6;i++) {
board_t bs[] = {(board_t){y0 - bb, y0 + bb, x0 - bb, x0 + bb, {0}, {0}, {0}},};
set_starting_board(&bs[0], i);
step(2, 2, bs, 1, &best_max[i], best_b[i], 0);
}

int best_i=0, mm = 0;
for(int i=0;i<6;i++)
if (best_max[i] > mm)
mm = best_max[i],
best_i = i;
st_print(MAX_CELL, best_max[best_i], &best_b[best_i][0]);

return 0;
};

• Will this find a solution that starts with: 1 131 .... much space .... 121 ? Oct 8, 2020 at 20:00

# C (Perl) n=6

My first^Wsecond pass at this is available on github ; I think this should in principle be able to calculate up to a(8), but that'll take a while even now it has been recoded in C.

On my machine it takes 42s for a(4) and 14ks for a(5), traversing 63,200,517 and 18,371,175,865 board positions respectively; rewriting in C gave about a 250x speedup from the initial Perl prototype.

Solution found for a(5) = 49:

  .  . 39  .  .  . 47  . 49
46 20 19  . 40  . 23 24 25
26  .  1 18  . 22  .  1  .
.  6  2  7  8  . 14  . 42
. 28  3  1 38 13 27 41  .
. 48 17  4  5  .  .  .  .
.  .  .  9 43 15 31  .  .
.  . 30  . 10  . 16  .  .
.  .  . 21 11  1 29 45  .
.  .  . 32  . 12  .  .  .
. 34 33  . 44  .  .  .  .
35  1  .  .  .  .  .  .  .
. 36 37  .  .  .  .  .  .


(Oh, that's a symmetry of xash's solution, I somehow expected it to be different.)

Confirming a(6) = 60 took around 10 CPU-weeks (manually sharded) and traversed 4.57e12 positions. Solution found:

  . 56 42  . 60  .  .  .  .  .  .  .  .  .
.  . 14 28 32  .  .  .  .  .  .  .  .  .
. 29 10  4  . 35  .  .  .  .  .  .  .  .
. 44  5  1  3 46  .  .  .  .  .  .  .  .
.  .  . 31  2  6  . 37  .  .  .  .  .  .
55  .  . 11  9  1  7 30  .  .  .  .  .  .
54  1 12 45  . 25  8 15  .  .  .  .  .  .
27 26 13  .  . 33  . 40 16 34 51  .  .  .
53  . 39 52  .  .  .  .  1 17  .  .  .  .
.  .  .  .  .  .  . 57 18  . 36  .  .  .
.  .  .  .  .  .  .  . 38 19  .  .  .  .
.  .  .  .  .  .  .  . 58  1 20 41  .  .
.  .  .  .  .  .  .  . 59  . 21  .  . 47
.  .  .  .  .  .  .  .  .  . 43 22 23 24
.  .  .  .  .  .  .  .  .  .  .  .  1 48
.  .  .  .  .  .  .  .  .  .  . 50 49  .


Finding a(7) would, by extrapolation, take 200-250 times as long as a(6). I don't plan to attempt this.

The approach is a) to insert the 1s lazily as needed, and b) to store unconnected groups separately, coalescing them as needed.

Extending beyond a(8) would require allowing for the possibility that we need to simultaneously coalesce 3 or more groups. I won't bother trying to solve that unless I get the speed of a(8) down to under a day or so.

The core work is done by the Board->try function (C: try_board), which tries each possible way to place the next number in the current board, then recurses.

The Group->coalesce (C: coalesce_group) function was the last and trickiest part to write: given two groups, the location within each that will form the common point at which the new value will be inserted, and the number of additional 1s that must be placed around it, this algorithm:

• fixes the orientation of the first group, and tries out each of the 8 possible orientations of the second group;
• first checks the immediate neighbourhood of the common location, looking for orientations that allow the two groups to coexist and leave room for enough additional 1s;
• then tries to overlay one group on the other, checking for further clashes;
• finally generates the k of n combinations of the n available cells around the common location into which the k additional 1s requested can be placed.

The hardest bit is going to be finding bugs, since there are so few data points to check against. I've added more tests, but I don't have confidence that I've found all bugs.

Hugo

[2020-10-10: added precise timings and position counts]

[2020-10-13: progress in C, a(5) found]

[2020-11-05: a(6) = 60 confirmed]

• Can I ask if you are working on n=6? My (fast but may miss some unlikely possibilities) solver has been running for a couple of days and I have found m=60. Oct 15, 2020 at 11:36
• I'm still working on trying to make the code fast enough to try n=6; my computer won't stay up long enough for a multi-day run (though if someone else wants to run it they are welcome). I may try to shard the work, in which case I could also get some partial results early. Oct 16, 2020 at 12:38
• @user1502040 m=60 is impressive for n=6. My best so far is only m=57. Are you able to share this solution? Also what are you getting for n=7? Oct 24, 2020 at 1:47
• Here's a 60: . 56 42 . 60 . . . . . . . . .; . . 14 28 32 . . . . . . . . .; . 29 10 4 . 35 . . . . . . . .; . 44 5 1 3 46 . . . . . . . .; . . . 31 2 6 . 37 . . . . . .; 55 . . 11 9 1 7 30 . . . . . .; 54 1 12 45 . 25 8 15 . . . . . .; 27 26 13 . . 33 . 40 16 34 51 . . .; 53 . 39 52 . . . . 1 17 . . . .; . . . . . . . 57 18 . 36 . . .; . . . . . . . . 38 19 . . . .; . . . . . . . . 58 1 20 41 . .; . . . . . . . . 59 . 21 . . 47; . . . . . . . . . . 43 22 23 24; . . . . . . . . . . . . 1 48; . . . . . . . . . . . 50 49 . Oct 27, 2020 at 1:29
• I've traversed a bit over 3e12 positions for n=6 now, of an estimated 3.6e12 total positions; if that estimate is correct, I'd expect about 8.1e14 for n=7 (which would take me 1-2 years at this rate). Oct 31, 2020 at 1:58