Firstly, we should analyse the sample input. As we can see for two holes on each column, we have 56 ways to connect the given 4 holes with the infinite long shoelace. By watching the 9 possible ways to connect the holes in the case where the order does not matter, we can easily observe that by permuting the given array, we will obtain new arrays which have the same properties as the initial one: it connects at least one hole on each side and also each hole is connected once only.

The 9 possible solutions when the order does not count are these:

1-3 1-4 2-3 2-4
1-3-4 2-3-4 1-2-3 1-2-4
1-2-3-4

Thus, if we permute each of them we obtain: 4*2! + 4*3! + 1*4! = 8 + 24 + 24 = 56

So our first guess is that a valid formula is:
sum of (N choose X) * (N choose Y) * (X+Y)! where X = 1,N and Y = 1,N.

However, as the number of holes on a column is N < 4*10⁶ a O(N²) or even worse O(N²logMod) algorithm will exceed the time limit. Thus we have to think for a different solution.

We can use the inclusion and exclusion principle to divide all the possible types of solutions in two categories: all the possible solutions (which do not repeat a hole) and the incorrect solutions (which do not repeat a hole but do not respect the other criteria). Thus, the number of correct solutions equals the difference between these two.

Lets define Aranj(N,K) the number (N choose K) * K!. Thus the number of possible solutions (correct and incorrect - but not by repeating a hole) is the sum of Aranj(2*N,X) where X=1,2*N. If we force now our picking system to break the criteria that "the shoelace connects at least one hole from the left side column and one hole from the right side column" we have sum of 2*Aranj(N,Y) where Y=1,N ways to tie the shoelace incorrectly (it is symmetrical for both left column and right column). Thus, if we define F(K)= sum of Aranj(K,X) with X=1,K; the O(N) formula is:

F(2*N) - 2*F(N)

C++ solution (Gabriel Robert Inelus)

    1 #include <iostream>
    2 #define M 666013
    3 int F(int x){
    4     int a = 0,v = 1;
    5     for(int i = x; i >= 1; --i){
    6         v = (1LL * v * i) % M;
    7         a = (a + v) % M;
    8     }
    9     return a%M;
   10 }
   11 int main(){
   12     int N;
   13     std::cin>>N;
   14     std::cout<<((F(2*N)%M - 2*F(N))%M+M)%M;
   15     return 0;
   16 }

C solution (Marius Gavrilescu)

    1 #include <stdio.h>
    2 #define MOD 666013
    3 
    4 int main(void){
    5     int n, i;
    6     long long cur = 1, total = 0;
    7     scanf("%d", &n);
    8     for(i = n ; i >= 1 ; i--)
    9         cur = cur * i % MOD, total = (total - 2 * cur + 2 * MOD) % MOD;
   10     cur = 1;
   11     for(i = 2*n ; i >= 1 ; i--)
   12         cur = cur * i % MOD, total = (total + cur) % MOD;
   13     printf("%lld", total);
   14     return 0;
   15 }

Perl solution (Marius Gavrilescu)

    1 #!/usr/bin/perl
    2 use v5.14;
    3 use warnings;
    4 
    5 use constant MOD => 666013;
    6 
    7 my $nr = <>;
    8 
    9 my ($total, $mcur, $pcur) = (0, 1, 1);
   10 for (reverse 1 .. $nr) {
   11     $pcur = $pcur * $_ % MOD;
   12     $total -= 2 * $pcur;
   13 }
   14 for (reverse 1 .. 2*$nr) {
   15     $mcur = $mcur * $_ % MOD;
   16     $total += $mcur;
   17 }
   18 
   19 say (($total + MOD) % MOD);