(* 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=3 *)
semi = 2; d = 3; 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}];
   Do[Do[Do[
        Do[Do[Do[
              part1 = (-1)^(r1 + r2 + r3)*
                semi^(c1 + c2 + 
                   c3)*(n - (r1 + r2 + r3) - (c1 + c2 + c3))!*
                Binomial[n1 - r1, c1]*Binomial[n2 - r2, c2]*
                Binomial[n3 - r3, c3];
              
              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];];
              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]];];
                   z = Flatten[{z1, z2, z3}];
                   
                   Do[m1 = NthSubset[p1, Range[c1 + c2 + c3]]; 
                    j1 = Length[m1]; 
                    If[j1 == 0, s1 = 0, 
                    s1 = j1 + Sum[z[[m1[[i1]]]], {i1, 1, j1}];];
                    rest2 = Complement[Range[c1 + c2 + c3], m1];
                    Do[m2 = NthSubset[p2, rest2]; j2 = Length[m2];
                    
                    If[j2 == 0, s2 = 0, 
                    s2 = j2 + Sum[z[[m2[[i2]]]], {i2, 1, j2}];];
                    s3 = r1 + r2 + r3 - s1 - s2;
                    j3 = c1 + c2 + c3 - j1 - j2;
                    
                    If[n1 - s1 >= 0 && n2 - s2 >= 0 && n3 - s3 >= 0, 
                    suma = 
                    suma + 
                    part1*Binomial[n1 - s1, j1]*j1!*
                    Binomial[n2 - s2, j2]*j2!*Binomial[n3 - s3, j3]*
                    j3!;];
                    , {p2, 0, 2^(c1 + c2 + c3 - j1) - 1}];
                    , {p1, 0, 2^(c1 + c2 + c3) - 1}];
                   , {l3, 1, nl3}];
                 , {l2, 1, nl2}];
               , {l1, 1, nl1}];
              , {c3, 0, Min[r3, n3 - r3]}];
            , {c2, 0, Min[r2, n2 - r2]}];
          , {c1, 0, Min[r1, n1 - r1]}];
        , {r3, 0, n3 - 1}];
      , {r2, 0, n2 - 1}];
    , {r1, 0, n1 - 1}];
   suma, n!]
  , {n, 1, 18}]