Post

no, i will not do it in ocaml

Posted
By redandgreynotebook
33 min read

but, with the amount of recursive backtracking this thing has, it’d probably perform comparably in ocaml, with a lot less segfaults along the way.

problem

grid

this grid can be partitioned into nine L-shaped “hooks”: the largest 9-by-9 (so, 17 squares total), then 8-by-8, all the way down to 1-by-1 (which is just a square). partition the grid into said hooks, and place nine 9’s in one of the hooks, eight 8’s in another, and so on, such that:

  • the occupied squares must form a connected region, where two squares are connected if they share an edge
  • every 2-by-2 box must contain at least one empty square
  • for each row, the greatest common divisor of the set of numbers formed by concatenating digits (reading left-to-right) is equal to the outside left number
  • for each column, the greatest common divisor of the set of numbers formed by concatenating digits (reading top-to-bottom) is equal to the outside top number

the source, jane street puzzles, june 2023, has an example.

solution

the solution can be split into two parts: narrowing down the number of hook-number configurations, and solving each configuration to see if it’s the one.

reduce search space

the number of possible configurations is way too high to try and solve each one individually. after noticing that 2 has to be assigned to the 2-by-2 hook and 1 has to be assigned to the square, that still leaves around $4^8 \cdot 7!$ (about 330 million) configurations. the $4^8$ comes from the fact that each (non-square) hook can be oriented one of four ways: corner in one of the four of top left, top right, bottom left, bottom right; and the $7!$ comes from permuting the hook assignments of the numbers 3 through 9.

consider a completed configuration config, represented as a 9-by-9 grid of numbers, with each cell represented by the number it’s assigned to by its hook. for example, the example solution’s configuration would look like this:

1
2
3
4
5

4 4 4 4 4
5 5 5 5 4
3 3 3 5 4
1 2 3 5 4
2 2 3 5 4

if a configuration is feasible, then each individual row and column must have a solution that satisfies its corresponding gcd constraint. for example, this configuration fails because there’s no way to fill in row 4 to get a gcd of 123:

1
2
3
4
5

1 2 3 4 5
2 2 3 4 5
3 3 3 4 5
4 4 4 4 5
5 5 5 5 5

thus, if we had functions solve_row and solve_column that tried to find individual row/column solutions (disregarding all others), we could eliminate boards as follows. also, note that for each number i in the hook, there must be at least i occurrences of it. here’s some high-level pseudocode.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

let check_config config = 
    for each i from 1 to 9
        if there are less than i occurrences of i in config 
            return false

    if there exists a row that solve_row is false 
    or there exists a column that solve_column is false
        return false

    return true

let solve_row config row = 
    for each of the 2^9 possible row states
        if it fulfills the row's gcd constraint
            return true

    return false

let solve_column config column = you get the point

but this still means checking all configurations, even if we’re not solving all of them, which is still way too much work. so, let’s incrementally build up our possible search space of configurations: starting from the empty grid, we try and add the 9-by-9 hook-number assignment, then check_config this unfinished configuration, and only continue filling this configuration if it’s valid. we’ll need to change check_config a bit, to account for the fact that config may be unfinished.

as for how to fill in an individual hook, we can keep a bounding_box of the area of config we haven’t filled in yet; at the start bounding_box would just be [(0, 0); (0, 8); (8, 8); (8, 0)]. then, let dir represent which corner the hook covers: 0 for top left, 1 for top right, 2 for bottom right, 3 for bottom left.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

let check_config config = 
    for each i from 1 to 9
        if i is present in config
            if there are less than i occurrences of i in config
                return false

    for each completed row of config
        if solve_row config row is false
            return false
    for each completed column of config
        if solve_column config column is false
            return false
    
    return true

let fill_hook config bounding_box dir num = 
    let k = side length of bounding_box
    if num = 1 and k != 1
        return none
    else if num = 2 and k != 2
        return none
    else if 2 * k - 1 < num
        return none 

    let hook = k-by-k hook covering the corner of bounding_box indicated by dir
    shrink bounding_box appropriately as given by dir
    for each cell in hook
        config[cell] = num

    return config, bounding_box

now to find all valid configurations valid_configs, we can use backtracking. keep a set unused_nums to indicate the numbers we’ve yet to assign; then, when you fill a hook with a number, remove num from unused_nums and recursively call valid_configs as follows.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

let valid_configs = valid_configs_h [(0, 0); (0, 8); (8, 8); (8, 0)] [1; 2; 3; 4; 5; 6; 7; 8; 9] (blank 9-by-9 board)

let valid_configs_h bounding_box unused_nums unfinished_config = 
    let k = side length of bounding_box
    if k = 1
        fill the last cell of unfinished_config, which must be bounding_box, with 1
        if check_config unfinished_config is true
            return [unfinished_config]
        else
            return []
            
    let result = []
    for each dir from 0 to 3
        for each num in unused_nums
            if fill_hook unfinished_config bounding_box dir num returned something
                let test_config, test_bounds = result from fill_hook
                if check_config test_config is true
                    remove num from unused_nums
                    append (valid_configs_h test_bounds unused_nums test_config) onto result
                    add num back onto unused_nums          
            
    return result

slight problem. if you run this, you get every one of the $\sim 4^8 \cdot 7!$ possible configurations we mentioned earlier. this is because solve_row, solve_column always return true: simply making the entire row/column blank is a solution. to actually make our gcd constraints restrictive, we have to prove there does not exist an empty row/column.

tougher gcd constraints

this is not entirely rigorous, but it does provide a heuristic as to why no row can be empty.

assume the first row is empty. then, you need to fill the remaining 8-by-9 grid with 45 numbers. if you pack it really tightly, you’ll find you can fit at most 56:

1
2
3
4
5
6
7
8
9

- - - - - - - - - 
X - X - X - X - X
X X X X X X X X X
X - X - X - X - X
X X X X X X X X X
X - X - X - X - X
X X X X X X X X X
X - X - X - X - X
X X X X X X X X X

but now keep in mind a row/column with all filled numbers is its own nine-digit-long gcd, so we need to break those up (it’s best to break along intersections if you can):

1
2
3
4
5
6
7
8
9

- - - - - - - - - 
X - X - X - X - X
X X - X X X X X X
X - X - X - X - X
X X X X - X X X -
X - X - X - X - X
- X X X X X - X X
X - X - X - X - X
X X X X - X X X X

and now keep in mind that the outermost 9-by-9 hook, wherever it might be, can hold at most 9 numbers, so clearly it’s either located in the top left or top right corner (and is assigned the number 7). but then the 8-by-8 hook can hold at most 9 numbers so it has to be located in the same corner as the 9-by-9 hook is; without loss of generality we assume both are located in the top left (and is assigned the number 8). now here’s the problem: how do we fill the 7-by-7 hook? at the very least, the 7-by-7 hook contains at least 11 numbers; without loss of generality assume it’s in the bottom right. we assign it the number 9 and are forced to remove two more:

1
2
3
4
5
6
7
8
9

- - - - - - - - - 
X - X - X - X - X
X X - X X X X X X
X - X - X - X - X
X X X X - X X X -
X - X - X - X - X
- X X X X X - X X
X - X - X - X - -
X X X X - X X - X

you get the point: examining all hooks eventually forces us to have less than 45 filled cells, so a solution is impossible. using a similar argument, we can show that no row/column can have two or less cells filled, implying that each row/column has at least two distinct numbers for the purpose of calculating gcd. this greately helps us narrow down individual solutions to rows/columns:

1
2
3
4
5
6
7
8

let solve_row config row = 
    for each of the 2^9 possible row states
        if it has at least two distinct numbers and fulfills the row's gcd constraint
            return true

    return false

let solve_column config column = you get the point

doing this leaves us with only 1005 valid configurations, out of around 330 million total. it means that around one in every $3 \cdot 10^6$ configurations are valid, a huge cutdown of search space.

solve each configuration

our first approach was simple backtracking: given a configuration, we know that the square had to be filled, so then we could do the following backtracking. note that we must add an adjacency check if we’re trying to fill a cell with a number (and not a space), as to not create additional disconnected components.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

let solve_config config grid = 
    if grid is completely filled
        if the grid passes all constraints
            return true
        return false
    
    for each unfilled cell in grid
        if cell is adjacent to a number
            if check_cell grid cell config[cell] is true
                let grid[cell] = config[cell]
                if solve_config config grid is true
                    return true
                unfill cell
        if check_cell grid cell ' ' is true
            let grid[cell] = ' '
            if solve_config config grid is true
                return true
            unfill cell

    return false

let check_cell grid cell target = 
    let grid[cell] = target

    if two-by-two constraint is broken return false
    for each i from 1 to 9
        if there are more than i occurrences of i in grid return false
    for each completed row and column of grid
        if the gcd constraint is broken return false

    return true

slight problem. for each hook with side length $1 < k \leq 9$, there are at most $\binom{2k - 1}{k} = O\left[\binom{2k}{k}\right] \sim \frac{4^k}{\sqrt{k\pi}}$ possible states of that hook, for an upper bound of $\frac{2^{88}}{\sqrt{9! \cdot \pi^9}} \approx 3 \cdot 10^{21}$ possible board states. yeah no.

let’s go back to finding individual row/column solutions satisfying individual gcd constraints. if $S$ is the set of all row states satisfying that row gcd constraint, then there may be cells common to each element in $S$, thus forcing that cell state in the final solution. of course, we need to check this against the forced cells from column gcds so:

1
2
3
4
5
6
7
8
9
10
11
12
13
14

let forced_grid config = 
    let result = empty 9-by-9 grid
    let forced_rows = the grid gotten from the forced cells from all individual rows
    let forced_columns = the grid gotten from the forced cells from all invidual columns

    for each cell in result
        if forced_rows[cell], forced_columns[cell] disagree
            return none
        if exactly one of forced_rows[cell], forced_columns[cell] is actually forced
            result[cell] = that forced state
        else
            result[cell] stays unfilled

    return result

where forced_rows, forced_columns disagree if one of them says cell says it’s forced to be a number and the other says it’s forced to be a space. on average, around half the cells in each grid are forced, greatly narrowing down the possibilities.

putting it together

so now it’s a bit clearer how it works:

  • call valid_configs to generate all valid configurations
  • for each config of these configurations, use solve_grid config (forced_grid config) to solve it

and running our code yields a solution of

1
2
3
4
5
6
7
8
9

5 5 @ @ @ 5 5 @ @ 
@ 9 9 9 9 @ 9 @ 8 
@ @ 6 @ 4 4 4 7 8 
@ 9 @ 4 3 @ @ @ 8 
@ 9 6 @ 3 1 2 @ @ 
@ 9 @ @ 3 @ 2 7 @ 
@ 9 6 6 6 6 @ 7 8 
@ 7 @ 7 @ 7 @ 7 @ 
5 8 @ 8 8 8 @ @ @

with an answer of $15552$. we are done.

for completeness, here’s our full implementation in c++. it takes around a minute to solve in replit, taking up 100mb of ram. this compares really nicely to a solution that uses a satisfiability solver, so i’m pretty happy with my results. to see just how optimized forced_grid makes it, a solve_config without it takes around 500-1500 seconds, while a solve_config forcing cells first takes at most 1.5 seconds.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639

#include <iostream>
#include <vector>
#include <string>
#include <cctype>
#include <stack>
#include <queue>
#include <algorithm>
#include <tuple>
#include <chrono>
#include <functional>
#include <numeric>

using std::vector;
using std::pair;
using std::tuple;
using std::make_tuple;
using std::get;
using std::string;
using std::stack;
using std::queue;
using std::prev_permutation;
using std::chrono::steady_clock;
using std::chrono::duration;
using std::chrono::milliseconds;
using std::chrono::duration_cast;
using std::function;
using std::accumulate;

typedef vector<vector<char>> board;
typedef pair<int, int> cell;

template <typename T, typename U> pair<T, U> operator+(const pair<T, U> &a, const pair<T, U> &b) {
    return {a.first + b.first, a.second + b.second};
}
template <typename T, typename U> pair<T, U> operator-(const pair<T, U> &a, const pair<T, U> &b) {
    return {a.first - b.first, a.second - b.second};
}
template <typename T, typename U> pair<T, U> operator*(const pair<T, U> &a, const int &b) {
    return {a.first * b, a.second * b};
}

const cell up = {-1, 0}, down = {1, 0}, left = {0, -1}, right = {0, 1};
const vector<cell> corners = { {0, 0}, {0, 1}, {1, 1}, {1, 0} };
const vector<cell> upper_left = {up, left, up + left}, upper_right = {up, right, up + right}, lower_left = {down, left, down + left}, lower_right = {down, right, down + right};
const vector<char> int_to_char = {'1', '2', '3', '4', '5', '6', '7', '8', '9'};
const vector<int> top_constraints = {5, 1, 6, 1, 8, 1, 22, 7, 8};
const vector<int> left_constraints = {55, 1, 6, 1, 24, 3, 6, 7, 2};

board create_board(int);
board copy_board(board &);
bool two_by_two(board &, cell);
bool adjacent(board &, cell);
bool connected(board &);
int answer(board &);
int unfilled(board &);
vector<int> components(board &, function<bool(char)>);
vector<int> get_row_ints(board &, int);
vector<int> get_col_ints(board &, int);
int gcd(int, int);
int gcd_list(vector<int>);
bool check_row_gcd(board &, int, int);
bool check_col_gcd(board &, int, int);
bool fill_cell(board &, cell, char, int, int);
vector<cell> create_hook(cell, cell, cell);
bool fill_hook(board &, vector<cell> &, int, char);
vector<vector<char>> solve_row(board &, board &, int, int, int = 0);
vector<vector<char>> solve_col(board &, board &, int, int, int = 0);
bool check_config(board &, vector<int> &, vector<int> &);
vector<char> forced_row(board &, board &, int, int);
vector<char> forced_col(board &, board &, int, int);
bool forced_board(board &, board &, vector<int> &, vector<int> &);
vector<int> create_hook_reqs(board &);
bool solve_config(board &, board &, vector<int> &, vector<int> &, vector<int> &, int);
bool solve_puzzle(int, vector<int> &, vector<int> &);
void print_board(board &);
void print_solution(board &);

int main() {
    vector<int> top_gcds(top_constraints), left_gcds(left_constraints);

    solve_puzzle(9, top_gcds, left_gcds)

    return 0;
}

// '#' represents an unfilled square, while ' ' represents a space
board create_board(int N) {
    board result(N, vector<char>(N, '#'));
    return result;
}

// makes deep copy of board
board copy_board(board &b) {
    board result;
    for (auto& r : b) result.push_back(vector<char>(r));
    return result;
}

// satisfies the "at least one space in each 2x2 square" constraint
bool two_by_two(board &b, cell pos) {
    int N = b.size();
    if (pos.first != 0) { // can check up
        if (pos.second != 0) { // can check left
            vector<cell> empty;
            for (auto c : upper_left) {
                c = c + pos;
                if (b[c.first][c.second] == ' ' || b[c.first][c.second] == '#') empty.push_back(c);
            }
            // if there's only one empty cell in this block force it to be blank
            if (empty.size() == 0) return false;
        }
        if (pos.second != N - 1) { // can check right
            vector<cell> empty;
            for (auto c : upper_right) {
                c = c + pos;
                if (b[c.first][c.second] == ' ' || b[c.first][c.second] == '#') empty.push_back(c);
            }
            // if there's only one empty cell in this block force it to be blank
            if (empty.size() == 0) return false;
        }
    }
    if (pos.first != N - 1) { // can check down
        if (pos.second != 0) { // can check left
            vector<cell> empty;
            for (auto c : lower_left) {
                c = c + pos;
                if (b[c.first][c.second] == ' ' || b[c.first][c.second] == '#') empty.push_back(c);
            }
            // if there's only one empty cell in this block force it to be blank
            if (empty.size() == 0) return false;
        }
        if (pos.second != N - 1) { // can check right
            vector<cell> empty;
            for (auto c : lower_right) {
                c = c + pos;
                if (b[c.first][c.second] == ' ' || b[c.first][c.second] == '#') empty.push_back(c);
            }
            // if there's only one empty cell in this block force it to be blank
            if (empty.size() == 0) return false;
        }
    }

    return true;
}

// check if cell is adjacent to any others
bool adjacent(board &b, cell pos) {
    int N = b.size();
    if (pos.first != 0) {
        auto c = pos + up;
        if (isdigit(b[c.first][c.second])) return true;
    }
    if (pos.first != N - 1) {
        auto c = pos + down;
        if (isdigit(b[c.first][c.second])) return true;
    }
    if (pos.second != 0) {
        auto c = pos + left;
        if (isdigit(b[c.first][c.second])) return true;
    }
    if (pos.second != N - 1) {
        auto c = pos + right;
        if (isdigit(b[c.first][c.second])) return true;
    }

    return false;
}

// check if the numbers form a connected region
bool connected(board &b) {
    int N = b.size();
    auto test = components(b, [](char c) {return isdigit(c);});
    return test.size() == 1 && test[0] == (N * N + N) / 2;
}

// gets product of empty regions
int answer(board &b) {
    int result = 1;
    auto test = components(b, [](char c) {return c == ' ';});
    for (auto& i : test) result *= i;
    return result;
}

// gets how many unfilled cells there are
int unfilled (board &b) {
    int result = 0;
    for (auto &v : b) {
        for (auto &c : v) result += c == '#';
    }
    return result;
}

// checks if filled regions are all orthogonally connected given a certain condition
vector<int> components(board &b, function<bool(char)> criterion) {
    int N = b.size();
    queue<cell> q;
    vector<int> result;
    vector<vector<bool>> visited(N, vector<bool>(N, false));
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
            if (!visited[i][j] && criterion(b[i][j])) {
                int count = 0;
                q.push({i, j});
                visited[i][j] = true;
                while (!q.empty()) {
                    if (q.front().first != 0) {
                        auto c = q.front() + up;
                        if (!visited[c.first][c.second] && criterion(b[c.first][c.second])) {
                            q.push(c);
                            visited[c.first][c.second] = true;
                        }
                    }
                    if (q.front().first != N - 1) {
                        auto c = q.front() + down;
                        if (!visited[c.first][c.second] && criterion(b[c.first][c.second])) {
                            q.push(c);
                            visited[c.first][c.second] = true;
                        }
                    }
                    if (q.front().second != 0) {
                        auto c = q.front() + left;
                        if (!visited[c.first][c.second] && criterion(b[c.first][c.second])) {
                            q.push(c);
                            visited[c.first][c.second] = true;
                        }
                    }
                    if (q.front().second != N - 1) {
                        auto c = q.front() + right;
                        if (!visited[c.first][c.second] && criterion(b[c.first][c.second])) {
                            q.push(c);
                            visited[c.first][c.second] = true;
                        }
                    }
                    ++count;
                    q.pop();
                }
                result.push_back(count);
            }
        }
    }

    return result;
}

// concatenates digits in the specified row 
// only returns "finalized" integers, with spaces both before and after the integer
// so, the integer isn't surrounded by unfilled cells - filled spaces only
vector<int> get_row_ints(board &b, int row) {
    int N = b.size();
    vector<int> result;
    string acc = "";
    bool valid_int = isdigit(b[row][0]);
    for (int i = 0; i < N; ++i) {
        // start of an integer
        if (b[row][i] == ' ' && !valid_int && i + 1 < N && isdigit(b[row][i + 1])) valid_int = true;
        else if (valid_int && isdigit(b[row][i])) {
            acc += b[row][i];
            // end of an integer
            if (i + 1 < N) {
                if (b[row][i + 1] == ' ' || b[row][i + 1] == '#') {
                    if (b[row][i + 1] == ' ') result.push_back(std::stoi(acc));
                    acc.clear();
                    valid_int = false;
                }
            } else result.push_back(std::stoi(acc));
        }
    }

    return result;
}

// concatenates digits in the specified column
vector<int> get_col_ints(board &b, int col) {
    int N = b.size();
    vector<int> result;
    string acc = "";
    bool valid_int = isdigit(b[0][col]);
    for (int i = 0; i < N; ++i) {
        // start of an integer
        if (b[i][col] == ' ' && !valid_int && i + 1 < N && isdigit(b[i + 1][col])) valid_int = true;
        else if (valid_int && isdigit(b[i][col])) {
            acc += b[i][col];
            // end of an integer
            if (i + 1 < N) {
                if (b[i + 1][col] == ' ' || b[i + 1][col] == '#') {
                    if (b[i + 1][col] == ' ') {result.push_back(std::stoi(acc));}
                    acc.clear();
                    valid_int = false;
                }
            } else result.push_back(std::stoi(acc));
        }
    }

    return result;
}

// greatest common denominator function
int gcd(int a, int b) {
    return a ? gcd(b % a, a) : b;
}

// calculates gcd of a list
int gcd_list(vector<int> arr) {
    if (arr.size() == 0) return 0;
    int result = arr[0];
    for (int i = 1; i < arr.size(); ++i) result = gcd(result, arr[i]);
    return result;
}

// checks gcd constraint of a filled row
bool check_row_gcd(board &b, int row, int left_gcd) {
    if (!left_gcd) return true;
    int N = b.size();
    bool filled = true;
    for (int i = 0; i < N; ++i) filled = filled && b[row][i] != '#';
    if (!filled) return true;
    return get_row_ints(b, row).size() > 1 ? gcd_list(get_row_ints(b, row)) == left_gcd : false;
}

// checks gcd constraint of a filled column
bool check_col_gcd(board &b, int col, int top_gcd) {
    if (!top_gcd) return true;
    int N = b.size();
    bool filled = true;
    for (int i = 0; i < N; ++i) filled = filled && b[i][col] != '#';
    if (!filled) return true;
    return get_col_ints(b, col).size() > 1 ? gcd_list(get_col_ints(b, col)) == top_gcd : false;
}

// checks all constraints and fills cell if you can
bool fill_cell(board &b, cell pos, char target, int top_gcd, int left_gcd) {
    // bool result = target == ' ' ? true : adjacent(b, pos);
    bool result = true;
    b[pos.first][pos.second] = target;
    // std::cout << "Target cell: (" << pos.first << ", " << pos.second << "): " << target << std::endl;
    // print_board(b);
    result = result && check_row_gcd(b, pos.first, left_gcd) && check_col_gcd(b, pos.second, top_gcd) && (target == ' ' ? true : two_by_two(b, pos));
    if (!result) b[pos.first][pos.second] = '#';
    // std::cout << result << std::endl;

    return result;
}

// creates a hook given its three endpoints are in clockwise order
vector<cell> create_hook(cell end1, cell corner, cell end2) {
    vector<cell> result = {end1};
    int size = std::max(abs((corner - end1).first), abs((corner - end1).second));
    if (size == 0) return result;
    cell first_leg = corner - end1;
    first_leg.first /= size, first_leg.second /= size;
    while (result.back() != corner) result.push_back(result.back() + first_leg);
    cell second_leg = end2 - corner;
    second_leg.first /= size, second_leg.second /= size;
    while (result.back() != end2) result.push_back(result.back() + second_leg);

    return result;
}

// checks if hook can feasibly be assigned number
bool fill_hook(board &b, vector<cell> &bounding_box, int dir, char target) {
    auto hook = create_hook(bounding_box[(dir + 3) & 3], bounding_box[dir], bounding_box[(dir + 1) & 3]);
    if (hook.size() < target - '0') return false;
    if (target == '1' && hook.size() != 1) return false;
    if (target == '2' && hook.size() != 3) return false;
    switch (dir) {
        case 0:
            bounding_box[0] = bounding_box[0] + down + right;
            bounding_box[3] = bounding_box[3] + right;
            bounding_box[1] = bounding_box[1] + down;
            break;
        case 1:
            bounding_box[0] = bounding_box[0] + down;
            bounding_box[2] = bounding_box[2] + left;
            bounding_box[1] = bounding_box[1] + down + left;
            break;
        case 2:
            bounding_box[3] = bounding_box[3] + up;
            bounding_box[2] = bounding_box[2] + left + up;
            bounding_box[1] = bounding_box[1] + left;
            break;
        case 3:
            bounding_box[3] = bounding_box[3] + up + right;
            bounding_box[2] = bounding_box[2] + up;
            bounding_box[0] = bounding_box[0] + right;
            break;
        default: 
            std::cout << "How did we get here?" << std::endl;
            break;
    }
    for (auto &pos : hook) b[pos.first][pos.second] = target;

    return true;
}

// solve an individual row
vector<vector<char>> solve_row(board &b, board &config, int row, int left_gcd, int step) {
    int N = b.size();
    vector<vector<char>> result;
    if (step == N) {
        if (check_row_gcd(b, row, left_gcd)) result.push_back(vector<char>());
        return result;
    }
    b[row][step] = config[row][step];
    auto partial_solution = solve_row(b, config, row, left_gcd, step + 1);
    for (auto& s : partial_solution) {
        s.insert(s.begin(), config[row][step]);
        result.push_back(s);
    }
    b[row][step] = ' ';
    partial_solution = solve_row(b, config, row, left_gcd, step + 1);
    for (auto& s : partial_solution) {
        s.insert(s.begin(), ' ');
        result.push_back(s);
    }
    b[row][step] = '#';
    
    return result;
}

// solve an individual column
vector<vector<char>> solve_col(board &b, board &config, int col, int top_gcd, int step) {
    int N = b.size();
    vector<vector<char>> result;
    if (step == N) {
        if (check_col_gcd(b, col, top_gcd)) result.push_back(vector<char>());
        return result;
    }
    b[step][col] = config[step][col];
    auto partial_solution = solve_col(b, config, col, top_gcd, step + 1);
    for (auto& s : partial_solution) {
        s.insert(s.begin(), config[step][col]);
        result.push_back(s);
    }
    b[step][col] = ' ';
    partial_solution = solve_col(b, config, col, top_gcd, step + 1);
    for (auto& s : partial_solution) {
        s.insert(s.begin(), ' ');
        result.push_back(s);
    }
    b[step][col] = '#';
    
    return result;
}

// check a (potentially in-progress) configuration to see if it violates gcd constraints
bool check_config(board &config, vector<int> &top_constraints, vector<int> &left_constraints) {
    int N = config.size();
    for (int i = 0; i < N; ++i) {
        bool row_filled = true, col_filled = true;
        for (int j = 0; j < N; ++j) {
            row_filled = row_filled && config[i][j] != '#';
            col_filled = col_filled && config[j][i] != '#';
        }
        auto b = create_board(N);
        if (row_filled && !solve_row(b, config, i, left_constraints[i]).size()) return false;
        if (col_filled && !solve_col(b, config, i, top_constraints[i]).size()) return false;
    }

    return true;
}

// find common intersections of all solutions to a row
vector<char> forced_row(board &b, board &config, int row, int left_gcd) {
    int N = b.size();
    vector<char> result(N, '#');
    auto possibilities = solve_row(b, config, row, left_gcd);
    if (!possibilities.size()) return vector<char>();
    for (int i = 0; i < N; ++i) {
        bool all_filled = true, all_space = true;
        for(auto &p : possibilities) {
            all_filled = all_filled && isdigit(p[i]);
            all_space = all_space && p[i] == ' ';
        }
        if (all_filled) result[i] = config[row][i];
        else if (all_space) result[i] = ' ';
    }

    return result;
}

// find common intersections of all solutions to a column
vector<char> forced_col(board &b, board &config, int col, int top_gcd) {
    int N = b.size();
    vector<char> result(N, '#');
    auto possibilities = solve_col(b, config, col, top_gcd);
    if (!possibilities.size()) return vector<char>();
    for (int i = 0; i < N; ++i) {
        bool all_filled = true, all_space = true;
        for(auto &p : possibilities) {
            all_filled = all_filled && isdigit(p[i]);
            all_space = all_space && p[i] == ' ';
        }
        if (all_filled) result[i] = config[i][col];
        else if (all_space) result[i] = ' ';
    }

    return result;
}

// find common intersections to row / column solutions
// returns false if the intersections break the other feasibility rules
bool forced_board(board &b, board &config, vector<int> &top_constraints, vector<int> &left_constraints) {
    int N = b.size();
    auto forced_rows = create_board(N), forced_cols = create_board(N);
    for (int i = 0; i < N; ++i) {
        auto partial_row = forced_row(b, config, i, left_constraints[i]), partial_col = forced_col(b, config, i, top_constraints[i]);
        if (partial_row.size() == 0 || partial_col.size() == 0) return false;
        for (int j = 0; j < N; ++j) {
            forced_rows[i][j] = partial_row[j];
            forced_cols[j][i] = partial_col[j];
        }
    }
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
            bool filled = isdigit(forced_rows[i][j]) || isdigit(forced_cols[i][j]);
            bool space = forced_rows[i][j] == ' ' || forced_cols[i][j] == ' ';
            if (filled && space) return false;
            else if (filled) b[i][j] = config[i][j];
            else if (space) b[i][j] = ' ';
        }
    }
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
            if (isdigit(b[i][j]) && !two_by_two(b, {i, j})) return false;
        }
    }
    auto hook_reqs = create_hook_reqs(b);
    for (auto &i : hook_reqs) if (i < 0) return false;

    return true;
}

// specify how much of each integer needs to be filled in each hook
vector<int> create_hook_reqs(board &b) {
    int N = b.size();
    vector<int> result;
    for (int i = 1; i <= N; ++i) result.push_back(i);
    for (auto &v : b) {
        for (auto &c : v) if (isdigit(c)) result[c - '1']--;
    }
    return result;
}

// solve a particular permutation by backtracking
bool solve_config(board &b, board &config, vector<int> &top_constraints, vector<int> &left_constraints, vector<int> &hook_reqs, int steps_left) {
    if (steps_left == 0) {
        if (!connected(b)) return false;
        for (auto &i : hook_reqs) if (i > 0) return false;
        std::cout << "product of areas is " << answer(b) << "\n" << std::endl;
        print_solution(b);
        return true;
    }
    
    int N = b.size();
    if (accumulate(hook_reqs.begin(), hook_reqs.end(), 0) > steps_left) return false;

    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
            if (b[i][j] == '#') {
                if (hook_reqs[config[i][j] - '1'] && fill_cell(b, {i, j}, config[i][j], top_constraints[j], left_constraints[i])) {
                    hook_reqs[config[i][j] - '1']--;
                    if (solve_config(b, config, top_constraints, left_constraints, hook_reqs, steps_left - 1)) return true;
                    hook_reqs[config[i][j] - '1']++;
                    b[i][j] = '#';
                }
                if (fill_cell(b, {i, j}, ' ', top_constraints[j], left_constraints[i])) {
                    if (solve_config(b, config, top_constraints, left_constraints, hook_reqs, steps_left - 1)) return true;
                    b[i][j] = '#';
                }
                b[i][j] = '#';
                return false;
            }
        }
    }

    return false;
}

// solves the entire puzzle
bool solve_puzzle(int N, vector<int> &top_constraints, vector<int> &left_constraints) {
    auto start = steady_clock::now();

    stack<tuple<board, vector<cell>, int, vector<bool>>> s;
    s.push(make_tuple(create_board(N), vector<cell> { {0, 0}, {0, N - 1}, {N - 1, N - 1}, {N - 1, 0} }, N, vector<bool>(N, true)));
    int valid_boards = 0;

    while (!s.empty()) {
        auto curr = s.top();
        s.pop();
        for (int dir = 0; dir < 4; ++dir) {
            for (int i = 0; i < N; ++i) {
                if (get<3>(curr)[i]) {
                    auto new_board = copy_board(get<0>(curr));
                    auto new_bounding_box = vector<cell>(get<1>(curr));
                    if (fill_hook(new_board, new_bounding_box, dir, int_to_char[i])) {
                        if (check_config(new_board, top_constraints, left_constraints)) {
                            if (get<2>(curr) == 1) {
                                ++valid_boards;
                                auto b = create_board(N);
                                if (forced_board(b, new_board, top_constraints, left_constraints)) {
                                    auto hook_reqs = create_hook_reqs(b);
                                    auto cells_left = unfilled(b);
                                    if (solve_config(b, new_board, top_constraints, left_constraints, hook_reqs, cells_left)) {
                                        std::cout << "solution found after " << duration_cast<milliseconds>(steady_clock::now() - start).count() / 1000.0 << " seconds" << std::endl;
                                        return true;
                                    }
                                }
                            } else {
                                auto new_perm = vector<bool>(get<3>(curr));
                                new_perm[i] = false;
                                s.push(make_tuple(new_board, new_bounding_box, get<2>(curr) - 1, new_perm));
                            }
                        }
                    }
                }
            }
        }
    }

    return false;
}

// debug print
void print_board(board &b) {
    for (auto& r : b) {
        for (auto &c : r) std::cout << c << " ";
        std::cout << std::endl;
    }
    std::cout << std::endl;
}

// it's easier to see spaces in the solution
void print_solution(board &b) {
    for (auto& r : b) {
        for (auto &c : r) std::cout << (c == ' ' ? '@' : c) << " ";
        std::cout << std::endl;
    }
    std::cout << std::endl;
}
puzzle
This post is licensed under CC BY 4.0 by the author.