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

Question

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 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.

Answer 1

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;
            if (optimize_spread) {
                    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
    if(!optimize_spread || max != 2)
    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;
                add_sum(max+1,y,x,b->sum);
                add_placed(1,y,x,b->placed);
                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;
                add_placed(-1,y,x,b->placed);
                add_sum(-(max+1),y,x,b->sum);
                b->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 <= 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;
                    add_sum(max+1,y,x,b->sum);
                    add_sum(1,y+ys[d1],x+xs[d1],b->sum);
                    add_placed(1,y,x,b->placed);
                        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;
                    add_placed(-1,y,x,b->placed);
                    add_sum(-(max+1),y,x,b->sum);
                    add_sum(-1,y+ys[d1],x+xs[d1],b->sum);
                    b->b[y+ys[d1]][x+xs[d1]] = 0;
                    b->b[y][x] = 0;
                }
            }
    }
    // place with 2 new 1-cells
    if(!optimize_spread || max != 2)
    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;
                    add_sum(max+1,y,x,b->sum);
                    add_sum(1,y+ys[d1],x+xs[d1],b->sum);
                    add_sum(1,y+ys[d2],x+xs[d2],b->sum);
                    add_placed(1,y,x,b->placed);
                        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;
                    add_placed(-1,y,x,b->placed);
                    add_sum(-(max+1),y,x,b->sum);
                    add_sum(-1,y+ys[d1],x+xs[d1],b->sum);
                    add_sum(-1,y+ys[d2],x+xs[d2],b->sum);
                    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;
                    add_sum(max+1,y,x,b->sum);
                    add_sum(1,y+ys[d1],x+xs[d1],b->sum);
                    add_sum(1,y+ys[d2],x+xs[d2],b->sum);
                    add_sum(1,y+ys[d3],x+xs[d3],b->sum);
                    add_placed(1,y,x,b->placed);
                        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;
                    add_placed(-1,y,x,b->placed);
                    add_sum(-(max+1),y,x,b->sum);
                    add_sum(-1,y+ys[d1],x+xs[d1],b->sum);
                    add_sum(-1,y+ys[d2],x+xs[d2],b->sum);
                    add_sum(-1,y+ys[d3],x+xs[d3],b->sum);
                    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;
};

Answer 2

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]