

$IterationLimit = $RecursionLimit = 10^7;

(**
 * define f such that f[i,j,n] evaluates to the number of walks in the
 * quarterplane starting at the origin, doing n admissible steps, and
 * ending up at point i,j.
 **)

ClearAll[c];
(* gessel walk carefully *)
c[i_Integer, j_Integer, k_Integer, n_Integer] :=
 Which[
   Min[i,j,k,n] < 0, 0, (* less than no steps *)
   i==j==k==n==0, 1, (* initial value *)
   n>k || (j+n)>(i+k), 0, (* out of quarter plane *)
   True, c[i,j,k,n] = c[i-1,j,k,n] + c[i,j-1,k,n] + c[i,j,k-1,n] + c[i,j,k,n-1]
 ];

(* Kreweras 3d *)
ClearAll[c];
c[i_Integer, j_Integer, k_Integer, n_Integer] :=
 Which[
   !(i<=j<=k), c@@Append[Sort[{i,j,k}], n],
   i < 0 || k > n || 2k > n+i || 3k > n+i+j || OddQ[n+i+j+k], 0,
   n == 0, KroneckerDelta[i, j, k, n],
   True, c[i,j,k,n] = 
     c[i+1,j,k,n-1] + c[i,j+1,k,n-1] + c[i,j,k+1,n-1] + c[i-1,j-1,k-1,n-1]
 ];

(* W, S, NE, in n steps to point (i,j); Kreweras *)
ClearAll[c];
c[i_Integer, j_Integer, n_Integer] :=
 Which[
   i > j, c[j, i, n],
   Max[i, j] > n || i < 0 || j < 0 || n < 0 || 2j > n+i, 0,
   n <= 1, KroneckerDelta[i, j, n],
   True, c[i,j,n] = Which[
     i == 0, c[1,j,n-1]+c[0,j+1,n-1],
     j == 0, c[i,1,n-1]+c[i+1,0,n-1],
     True, c[i,j+1,n-1]+c[i+1,j,n-1]+c[i-1,j-1,n-1]
   ]
 ];

(* NE, E, W, SW, in n steps to point (i,j); Gessel *)
ClearAll[c]; (* checked *)
c[i_Integer, j_Integer, n_Integer] :=
  Which[
    Max[i, j] > n || OddQ[ n + i ] || i < 0 || j < 0, 0,
    i == n, Binomial[n, j],
    j == n, KroneckerDelta[n, i], 
    i == j == n == 0, 1,
    (* 
    i == j == 0, 
     c[0, 0, n] = c[1, 1, n - 1] + c[1, 0, n - 1], *)
    i == 0, (* contains j = 0 *)
     c[0, j, n] = c[1, j, n - 1] + c[1, j+1, n - 1],
    j == 0, 
     c[i, 0, n] =
        c[i+1, 1, n - 1] + c[i-1, 0, n - 1]
       +c[i+1, 0, n - 1],
    True, 
     c[i, j, n] = (*
       -c[-2+i, -2+j, n] - c[-2 + i, -1 + j, n] 
       +c[-1 + i, -1 + j, 1 + n] - c[i, -1 + j, n] *)
        c[i-1, j, n - 1] + c[i-1, j-1, n - 1]
       +c[i+1, j+1, n - 1] + c[i+1, j, n - 1]
  ];
c[n_Integer] := c[n]=Sum[c[i,j,n]x^i y^j,{i,0,n},{j,0,n}];

ig[x_,y_,t_,n_Integer] := 
  Nest[
   Cancel[Together[1+t*(x+x*y)*#+t/x*(#-(#/.x->0))
                  +t/x/y*(#-(#/.x->0)-(#/.y->0)+(#/.x->0/.y->0))]]&
  , 1, n];


(* NE, SW, N, in n steps to point (i,j) *)
(*
ClearAll[c];
c[i_Integer, j_Integer, n_Integer] :=
  Which[
    i + j > 2*n, 0,
    i == j == n == 0, 1,
    j == 0, (* includes i=j=0 *)
      c[i, 0, n] = c[i + 1, 1, n - 1],
    i == 0,
      c[0, j, n] = c[0, j - 1, n - 1] + c[1, j + 1, n - 1],
    True,
      c[i, j, n] = 
         c[i + 1, j + 1, n - 1] + c[i - 1, j - 1, n - 1]
       + c[i, j - 1, n - 1]
  ];
*)



