That's a very nice challenge!

Here is a Prolog formulation of the task:

:- use_module(library(clpfd)).

n_tour(N, Vs) :-
L #= N*N,
length(Vs0, L),
successors(Vs0, N, 1),
append(Vs0, [_], Vs), % last element is for open tours
circuit(Vs).

successors([], _, _).
successors([V|Vs], N, K0) :-
findall(Num, n_k_next(N, K0, Num), [Next|Nexts]),
foldl(num_to_dom, Nexts, Next, Dom),
Dummy #= N*N + 1, % last element is for open tours
V in Dom \/ Dummy,
K1 #= K0 + 1,
successors(Vs, N, K1).

num_to_dom(N, D0, D0\/N).

n_x_y_k(N, X, Y, K) :- [X,Y] ins 1..N, K #= N*(Y-1) + X.

n_k_next(N, K, Next) :-
n_x_y_k(N, X0, Y0, K),
( [DX,DY] ins -2 \/ 2 % 2 diagonally
; [DX,DY] ins -3 \/ 0 \/ 3, % 3 vertically or horizontally
abs(DX) + abs(DY) #= 3
),
[X,Y] ins 1..N,
X #= X0 + DX,
Y #= Y0 + DY,
n_x_y_k(N, X, Y, Next),
label([DX,DY]).

The key element of this solution is the CPL(FD) constraint circuit/1, which describes a Hamiltonian circuit. It uses a list of finite domain variables to represent solutions: Each element of the list is an integer, denoting the position of the successor in the list. For example, a list with 3 elements admits precisely 2 Hamiltonian circuits, which we can find via search:

?- Vs = [_,_,_], circuit(Vs), label(Vs).
Vs = [2, 3, 1] ;
Vs = [3, 1, 2].

The rest of the code sets up the domains of the involved variables to constrain them to admissible moves. A dummy element is used to allow open tours. The predicate n_x_y_k/4 relates X/Y coordinates to list indices. You can easily adapt this to other variants of the puzzle (e.g., knight's tour) by changing n_k_next/3.

A major attraction of a declarative solution is that it can be used in all directions: We can use the exact same code to test, generate, and to complete solutions. For example, we can use the Prolog program to show that the 5×5 solution that is shown in the article is uniquely defined if we fix just 4 elements in the first row:

?- n_tour(5, Vs),
last(Vs, 1),
Vs = [_,_,18,19,26|_],
label(Vs).
Vs = [4, 5, 18, 19, 26, 21, 22, 20, 24|...] ;
false.

A few additional definitions let us print solutions in a more readable form:

:- set_prolog_flag(double_quotes, chars).

print_tour(Vs0) :-
reverse(Vs0, [First|Vs1]),
reverse(Vs1, Vs),
length(Vs, L),
L #= N*N, N #> 0,
length(Ts, N),
tour_enumeration(Vs, N, First, Es),
phrase(format_string(Ts, 0, 4), Fs),
maplist(format(Fs), Es).

format_(Fs, Args, Xs0, Xs) :- format(chars(Xs0,Xs), Fs, Args).

format_string([], _, _) --> "\n".
format_string([_|Rest], N0, I) -->
{ N #= N0 + I }, % I is textual width of square
"~t~w~", call(format_("~w|", [N])),
format_string(Rest, N, I).

tour_enumeration(Vs, N, First, Es) :-
length(Es, N),
maplist(same_length(Es), Es),
append(Es, Ls),
foldl(vs_enumeration(Vs, Ls), Vs, First-1, _).

vs_enumeration(Vs, Ls, _, V0-E0, V-E) :-
E #= E0 + 1,
nth1(V0, Ls, E0),
nth1(V0, Vs, V).

For example, here is the 5×5 solution from the article again:

?- n_tour(5, Vs),
last(Vs, 1),
Vs = [_,_,18,19,26|_],
label(Vs),
print_tour(Vs).
1 24 14 2 25
16 21 5 8 20
13 10 18 23 11
4 7 15 3 6
17 22 12 9 19

And here is a query that solves the 10×10 instance:

?- n_tour(10, Vs),
time(label(Vs)),
print_tour(Vs).

Yielding:

% 21,852,020 inferences, 3.323 CPU in 3.349 seconds (99% CPU, 6575193 Lips)
63 58 53 64 59 54 65 60 55 66
21 16 11 22 17 12 23 18 13 24
52 49 62 57 50 61 56 67 28 85
10 7 20 15 8 19 14 25 82 72
47 44 51 48 45 68 29 86 69 30
5 91 9 6 92 26 83 73 27 84
42 96 46 43 97 87 70 31 81 71
36 39 93 35 38 74 34 78 75 33
4 90 98 3 89 99 2 88 100 1
41 95 37 40 94 79 76 32 80 77

Since the search for solutions is decoupled from the problem description, we can easily try different search strategies. For example, using the strategy "ff" (first fail) measurably reduces the running time in this case:

?- n_tour(10, Vs),
time(labeling([ff], Vs)),
print_tour(Vs).
% 8,317,298 inferences, 1.344 CPU in 1.355 seconds (99% CPU, 6190382 Lips)
89 84 49 90 85 48 23 86 47 22
81 63 13 82 64 61 68 65 60 69
50 91 88 51 92 87 46 43 24 36
14 83 80 62 12 66 59 70 67 21
27 52 93 26 53 44 25 37 45 42
79 57 15 8 58 71 11 40 72 35
94 31 54 95 32 38 19 33 75 20
28 7 98 29 6 9 73 3 10 41
78 56 16 77 55 17 76 39 18 34
97 30 5 96 99 4 1 100 74 2

Using CLP(FD) constraints typically yields much faster solutions than using brute-force search, since constraints propagate information before and also during the search to prune the search tree. Especially for more complex tasks, the computational cost of constraint propagation is typically negligible in comparison to traversing the entire search tree.

However, if you want, you can also easily obtain a brute-force variation from this program, by posting the constraints after the search. This turns the whole program into a "generate and test" approach, which is typically much slower than constraint-based search.