(* Vaclav Kotesovec 2011, based on article by R. Tauraso, 2006 *)
<< Combinatorica`
(* semi=2 Rook + Leaper[d,d] *)
(* semi=1 Rook + semi-Leaper[d,d] *)
(* version for d=4 *)
semi = 1; d = 4; Table[If[n >= d, suma = 0;
   n1 = 0; Do[If[Mod[i - 1, d] == 0, n1 = n1 + 1], {i, 1, n}];
   n2 = 0; Do[If[Mod[i - 2, d] == 0, n2 = n2 + 1], {i, 1, n}];
   n3 = 0; Do[If[Mod[i - 3, d] == 0, n3 = n3 + 1], {i, 1, n}];
   n4 = 0; Do[If[Mod[i - 4, d] == 0, n4 = n4 + 1], {i, 1, n}];
   Do[Do[Do[
        Do[r = r1 + r2 + r3 + r4; Do[Do[Do[Do[c = c1 + c2 + c3 + c4;
                 
                 part1 = (-1)^r*semi^c*(n - r - c)!*
                   Binomial[n1 - r1, c1]*Binomial[n2 - r2, c2]*
                   Binomial[n3 - r3, c3]*Binomial[n4 - r4, c4];
                 
                 If[c1 == 0, If[r1 == 0, nl1 = 1, nl1 = 0];, 
                  nl1 = NumberOfCompositions[r1 - c1, c1];];
                 
                 If[c2 == 0, If[r2 == 0, nl2 = 1, nl2 = 0];, 
                  nl2 = NumberOfCompositions[r2 - c2, c2];];
                 
                 If[c3 == 0, If[r3 == 0, nl3 = 1, nl3 = 0];, 
                  nl3 = NumberOfCompositions[r3 - c3, c3];];
                 
                 If[c4 == 0, If[r4 == 0, nl4 = 1, nl4 = 0];, 
                  nl4 = NumberOfCompositions[r4 - c4, c4];];
                 Do[Do[Do[Do[
                    
                    If[c1 == 0, z1 = {};, 
                    z1 = Compositions[r1 - c1, c1][[l1]];];
                    
                    If[c2 == 0, z2 = {};, 
                    z2 = Compositions[r2 - c2, c2][[l2]];];
                    
                    If[c3 == 0, z3 = {};, 
                    z3 = Compositions[r3 - c3, c3][[l3]];];
                    
                    If[c4 == 0, z4 = {};, 
                    z4 = Compositions[r4 - c4, c4][[l4]];];
                    z = Flatten[{z1, z2, z3, z4}];
                    
                    Do[m1 = NthSubset[p1, Range[c]]; j1 = Length[m1]; 
                    If[j1 == 0, s1 = 0, 
                    s1 = j1 + Sum[z[[m1[[i1]]]], {i1, 1, j1}];];
                    rest2 = Complement[Range[c], m1];
                    Do[m2 = NthSubset[p2, rest2]; j2 = Length[m2];
                    
                    If[j2 == 0, s2 = 0, 
                    s2 = j2 + Sum[z[[m2[[i2]]]], {i2, 1, j2}];];
                    rest3 = Complement[rest2, m2];
                    Do[m3 = NthSubset[p3, rest3]; j3 = Length[m3];
                    
                    If[j3 == 0, s3 = 0, 
                    s3 = j3 + Sum[z[[m3[[i3]]]], {i3, 1, j3}];];
                    s4 = r - s1 - s2 - s3;
                    j4 = c - j1 - j2 - j3;
                    
                    If[n1 - s1 >= 0 && n2 - s2 >= 0 && n3 - s3 >= 0 && 
                    n4 - s4 >= 0, 
                    suma = 
                    suma + 
                    part1*Binomial[n1 - s1, j1]*j1!*
                    Binomial[n2 - s2, j2]*j2!*Binomial[n3 - s3, j3]*
                    j3!*Binomial[n4 - s4, j4]*j4!;];
                    , {p3, 0, 2^(c - j1 - j2) - 1}];
                    , {p2, 0, 2^(c - j1) - 1}];
                    , {p1, 0, 2^c - 1}];
                    , {l4, 1, nl4}];
                    , {l3, 1, nl3}];
                    , {l2, 1, nl2}];
                  , {l1, 1, nl1}];
                 , {c4, 0, Min[r4, n4 - r4]}];
               , {c3, 0, Min[r3, n3 - r3]}];
             , {c2, 0, Min[r2, n2 - r2]}];
           , {c1, 0, Min[r1, n1 - r1]}];
          , {r4, 0, n4 - 1}];
        , {r3, 0, n3 - 1}];
      , {r2, 0, n2 - 1}];
    , {r1, 0, n1 - 1}];
   suma, n!]
  , {n, 1, 16}]