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 &quot;nice&quot;.
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 code-challenge 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.
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 -&gt; 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 &gt; 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&gt;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&gt;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 &lt;stdint.h&gt;
#include &lt;stddef.h&gt;
#include &lt;stdio.h&gt;
#include &lt;string.h&gt;
#include &lt;stdlib.h&gt;

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

#define MAX(x,y) (x&gt;y?x:y)
#define MIN(x,y) (x&lt;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&lt;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&lt;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(&quot;%d cells, %d max\n&quot;, c, max);
    for(int y=b-&gt;y0;y&lt;=b-&gt;y1;y++){
        for(int x=b-&gt;x0;x&lt;=b-&gt;x1;x++) printf(&quot;%*d&quot;, 3, b-&gt;b[y][x]);
        puts(&quot;\n&quot;);
    }
}

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&lt;bl;i++) {
            if (bs[i].y0 &lt; 2 || bs[i].y1 &gt;= S - 2 || bs[i].x0 &lt; 2 || bs[i].x1 &gt;= S - 2) {
                    st_print(c, max, &amp;bs[i]);
                printf(&quot;board too small %d %d %d %d&quot;, bs[i].y0, bs[i].y1, bs[i].x0, bs[i].x1);
                exit(1);
            }
    }
    // new best
    if (c == MAX_CELL) {
        int score = 0;
            if (optimize_spread) {
                    for (int i=0;i&lt;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 &lt; score) {
                for (int i=0;i&lt;bl;i++)
                        memcpy(&amp;best_b[i], &amp;bs[i], sizeof(board_t));
                *best_max = score;
            }
    }
    // place with 0 new 1-cells
    if(!optimize_spread || max != 2)
    for(int i=0;i&lt;bl;i++) {
            board_t *b=bs+i;
            for(int y=b-&gt;y0;y&lt;=b-&gt;y1;y++)
            for(int x=b-&gt;x0;x&lt;=b-&gt;x1;x++) 
            if(b-&gt;sum[y][x] == max + 1 &amp;&amp; !b-&gt;b[y][x]) {
                b-&gt;b[y][x] = max + 1;
                add_sum(max+1,y,x,b-&gt;sum);
                add_placed(1,y,x,b-&gt;placed);
                int y0o = b-&gt;y0, y1o = b-&gt;y1, x0o = b-&gt;x0, x1o = b-&gt;x1;
                b-&gt;y0 = MIN(b-&gt;y0, y-2);
                b-&gt;y1 = MAX(b-&gt;y1, y+2);
                b-&gt;x0 = MIN(b-&gt;x0, x-2);
                b-&gt;x1 = MAX(b-&gt;x1, x+2);
                step(c, max + 1, bs, bl, best_max, best_b, optimize_spread);
                b-&gt;y0 = y0o, b-&gt;y1 = y1o, b-&gt;x0 = x0o, b-&gt;x1 = x1o;
                add_placed(-1,y,x,b-&gt;placed);
                add_sum(-(max+1),y,x,b-&gt;sum);
                b-&gt;b[y][x] = 0;
            }
    }
    // sorry for the repetition, couldn't get clang to optimize it otherwise
    // place with 1 new 1-cells
    if(!optimize_spread || max != 2)
    if(c + 1 &lt;= MAX_CELL) 
    for(int i=0;i&lt;bl;i++) {
            board_t *b=bs+i;
            for(int y=b-&gt;y0;y&lt;=b-&gt;y1;y++)
            for(int x=b-&gt;x0;x&lt;=b-&gt;x1;x++) 
            if(b-&gt;sum[y][x] == (max + 1) - 1 &amp;&amp; !b-&gt;b[y][x]) {
                for(int d1=0;d1&lt;8;d1++) {
                    if (b-&gt;placed[y+ys[d1]][x+xs[d1]]) continue;
                    
                    b-&gt;b[y+ys[d1]][x+xs[d1]] = 1;
                    b-&gt;b[y][x] = max + 1;
                    add_sum(max+1,y,x,b-&gt;sum);
                    add_sum(1,y+ys[d1],x+xs[d1],b-&gt;sum);
                    add_placed(1,y,x,b-&gt;placed);
                        int y0o = b-&gt;y0, y1o = b-&gt;y1, x0o = b-&gt;x0, x1o = b-&gt;x1;
                        b-&gt;y0 = MIN(b-&gt;y0, y-2);
                        b-&gt;y1 = MAX(b-&gt;y1, y+2);
                        b-&gt;x0 = MIN(b-&gt;x0, x-2);
                        b-&gt;x1 = MAX(b-&gt;x1, x+2);
                        step(c + 1, max + 1, bs, bl, best_max, best_b, optimize_spread);
                        b-&gt;y0 = y0o, b-&gt;y1 = y1o, b-&gt;x0 = x0o, b-&gt;x1 = x1o;
                    add_placed(-1,y,x,b-&gt;placed);
                    add_sum(-(max+1),y,x,b-&gt;sum);
                    add_sum(-1,y+ys[d1],x+xs[d1],b-&gt;sum);
                    b-&gt;b[y+ys[d1]][x+xs[d1]] = 0;
                    b-&gt;b[y][x] = 0;
                }
            }
    }
    // place with 2 new 1-cells
    if(!optimize_spread || max != 2)
    if(c + 2 &lt;= MAX_CELL) 
    for(int i=0;i&lt;bl;i++) {
            board_t *b=bs+i;
            for(int y=b-&gt;y0;y&lt;=b-&gt;y1;y++)
            for(int x=b-&gt;x0;x&lt;=b-&gt;x1;x++) 
            if(b-&gt;sum[y][x] == (max + 1) - 2 &amp;&amp; !b-&gt;b[y][x]) {
                for(int d1=0;d1&lt;8-1;d1++) {
                    if (b-&gt;placed[y+ys[d1]][x+xs[d1]]) continue;
                for(int d2=d1+1;d2&lt;8;d2++) {
                    if (b-&gt;placed[y+ys[d2]][x+xs[d2]]) continue;
                    
                    b-&gt;b[y+ys[d1]][x+xs[d1]] = 1;
                    b-&gt;b[y+ys[d2]][x+xs[d2]] = 1;
                    b-&gt;b[y][x] = max + 1;
                    add_sum(max+1,y,x,b-&gt;sum);
                    add_sum(1,y+ys[d1],x+xs[d1],b-&gt;sum);
                    add_sum(1,y+ys[d2],x+xs[d2],b-&gt;sum);
                    add_placed(1,y,x,b-&gt;placed);
                        int y0o = b-&gt;y0, y1o = b-&gt;y1, x0o = b-&gt;x0, x1o = b-&gt;x1;
                        b-&gt;y0 = MIN(b-&gt;y0, y-2);
                        b-&gt;y1 = MAX(b-&gt;y1, y+2);
                        b-&gt;x0 = MIN(b-&gt;x0, x-2);
                        b-&gt;x1 = MAX(b-&gt;x1, x+2);
                        step(c + 2, max + 1, bs, bl, best_max, best_b, optimize_spread);
                        b-&gt;y0 = y0o, b-&gt;y1 = y1o, b-&gt;x0 = x0o, b-&gt;x1 = x1o;
                    add_placed(-1,y,x,b-&gt;placed);
                    add_sum(-(max+1),y,x,b-&gt;sum);
                    add_sum(-1,y+ys[d1],x+xs[d1],b-&gt;sum);
                    add_sum(-1,y+ys[d2],x+xs[d2],b-&gt;sum);
                    b-&gt;b[y+ys[d1]][x+xs[d1]] = 0;
                    b-&gt;b[y+ys[d2]][x+xs[d2]] = 0;
                    b-&gt;b[y][x] = 0;
                }
                }
            }
    }
    // place with 3 new 1-cells
    if(c + 3 &lt;= MAX_CELL) 
    for(int i=(optimize_spread &amp;&amp; max == 2);i&lt;bl;i++) {
            board_t *b=bs+i;
            for(int y=b-&gt;y0;y&lt;=b-&gt;y1;y++)
            for(int x=b-&gt;x0;x&lt;=b-&gt;x1;x++)
            if(b-&gt;sum[y][x] == (max + 1) - 3 &amp;&amp; !b-&gt;b[y][x]) {
                for(int d1=0;d1&lt;8-2;d1++) {
                    if (b-&gt;placed[y+ys[d1]][x+xs[d1]]) continue;
                for(int d2=d1+1;d2&lt;8-1;d2++) {
                    if (b-&gt;placed[y+ys[d2]][x+xs[d2]]) continue;
                for(int d3=d2+1;d3&lt;8;d3++) {
                    if (b-&gt;placed[y+ys[d3]][x+xs[d3]]) continue;
                    
                    b-&gt;b[y+ys[d1]][x+xs[d1]] = 1;
                    b-&gt;b[y+ys[d2]][x+xs[d2]] = 1;
                    b-&gt;b[y+ys[d3]][x+xs[d3]] = 1;
                    b-&gt;b[y][x] = max + 1;
                    add_sum(max+1,y,x,b-&gt;sum);
                    add_sum(1,y+ys[d1],x+xs[d1],b-&gt;sum);
                    add_sum(1,y+ys[d2],x+xs[d2],b-&gt;sum);
                    add_sum(1,y+ys[d3],x+xs[d3],b-&gt;sum);
                    add_placed(1,y,x,b-&gt;placed);
                        int y0o = b-&gt;y0, y1o = b-&gt;y1, x0o = b-&gt;x0, x1o = b-&gt;x1;
                        b-&gt;y0 = MIN(b-&gt;y0, y-2);
                        b-&gt;y1 = MAX(b-&gt;y1, y+2);
                        b-&gt;x0 = MIN(b-&gt;x0, x-2);
                        b-&gt;x1 = MAX(b-&gt;x1, x+2);
                        step(c + 3, max + 1, bs, bl, best_max, best_b, optimize_spread);
                        b-&gt;y0 = y0o, b-&gt;y1 = y1o, b-&gt;x0 = x0o, b-&gt;x1 = x1o;
                    add_placed(-1,y,x,b-&gt;placed);
                    add_sum(-(max+1),y,x,b-&gt;sum);
                    add_sum(-1,y+ys[d1],x+xs[d1],b-&gt;sum);
                    add_sum(-1,y+ys[d2],x+xs[d2],b-&gt;sum);
                    add_sum(-1,y+ys[d3],x+xs[d3],b-&gt;sum);
                    b-&gt;b[y+ys[d1]][x+xs[d1]] = 0;
                    b-&gt;b[y+ys[d2]][x+xs[d2]] = 0;
                    b-&gt;b[y+ys[d3]][x+xs[d3]] = 0;
                    b-&gt;b[y][x] = 0;
                }
                }
                }
            }
    }
}

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

    for(int y=1;y+1&lt;S;y++)
    for(int x=1;x+1&lt;S;x++)
    for(int yd=-1;yd&lt;=1;yd++)
    for(int xd=-1;xd&lt;=1;xd++)
    if(yd!=0||xd!=0)
        b-&gt;sum[y][x] += b-&gt;b[y+yd][x+xd];
    for(int y=1;y+1&lt;S;y++)
    for(int x=1;x+1&lt;S;x++)
    for(int yd=-1;yd&lt;=1;yd++)
    for(int xd=-1;xd&lt;=1;xd++)
        b-&gt;placed[y][x] += b-&gt;b[y+yd][x+xd] &gt; 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&lt;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(&amp;bs[0], i);
        step(2, 2, bs, 2, &amp;best_max[i], best_b[i], 1);
    }

    int best_i=0, mm = 0;
    for(int i=0;i&lt;6;i++)
        if (best_max[i] &gt; mm)
            mm = best_max[i],
            best_i = i;
    printf(&quot;most spread of distant 112 and 1113: %d\n&quot;, best_max[best_i]);
    st_print(MAX_CELL, best_max[best_i], &amp;best_b[best_i][0]);
    st_print(MAX_CELL, best_max[best_i], &amp;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&lt;6;i++) {
        board_t bs[] = {(board_t){y0 - bb, y0 + bb, x0 - bb, x0 + bb, {0}, {0}, {0}},};
            set_starting_board(&amp;bs[0], i);
        step(2, 2, bs, 1, &amp;best_max[i], best_b[i], 0);
    }

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

    return 0;
};