Ackermann Function

This is a more or less verbatim implementation of the equations.
We memoize to save a bit of computation.

ClearAll[A]
A[0,y_]   :=  y+1;
A[x_,0]   := A[x,0]= A[x-1,1];
A[x_,y_] := A[x,y]= A[x-1,A[x,y-1]];

A table of the first few values at levels 0 through 3.

Table[ A[i,j], {i,0,3}, {j,0,6} ] // TableForm[#,TableHeadings->{Range[0,3],Range[0,6]}]&

Making a Trace

We can extract information about exactly which values of  A(x,y)  were computed so far. Since we are memoizing, this boils down to printing out the hash table associated with the function.
Don't look at the code, this is horrible Mathematica hacking.

all =Cases[First /@( First /@ DownValues[A]/. 
{HoldPattern[A[i_,j_]] -> {i,j}}),{_Integer,_Integer} ];
trace = Flatten[{#,A@@#}]&/@ all;

The size of the hash table, at each level:

Length[trace]
Length[ Cases[trace,{1,__}]]
Length[ Cases[trace,{2,__}]]
Length[ Cases[trace,{3,__}]]

Cases[trace, {3, __}]

Level 4

Things start to become very tedious at level 4.

A[4,0]

$RecursionLimit = 5000;

A[4,1]

Cheating

We can still obtain the value for A(4,1)  by cheating a bit and avoiding calls to the lower levels.
The following implementation returns values directly at levels 3 and lower.

ClearAll[A]

A[0,y_] := y+1;
A[1,y_] := y+2;
A[2,y_] := 2y+3;
A[3,y_] := 2^(y+3)-3;

A[x_,0] := A[x,0]= A[x-1,1];
A[x_,y_] := A[x,y]= A[x-1,A[x,y-1]];

A(4,1) Revealed

Checking small values:

Table[ A[i,j], {i,0,3}, {j,0,6} ] //TableForm[#,TableHeadings->{Range[0,3],Range[0,6]}]&

Not so small:

A[4,1]

In fact, we can even do   A(4,2) = A(3,A(4,1)) = (2^( 65536))-3:

A[4,2]

This is a number with 19729 decimal digits.

% // DigitCount // Total

But   A(4,3) = A(3,A(4,2)) = (2^( A(4,2)+3))-3  would crash the machine.

Closed Forms

The first 3 levels can easily be expressed in terms of well-known arithmetic functions.

ClearAll[a,n]

RSolve[  {a[n+1] ==  a[n] + 1, a[0]==2}, a[n], n ]

RSolve[  {a[n+1] ==  a[n] + 2, a[0]==3}, a[n], n ]

RSolve[  {a[n+1] == 2 a[n] + 3, a[0]==5}, a[n], n ]

But here things fall apart: super-exponentiation is not a standard function.

RSolve[  {a[n+1] == 2^( a[n] + 3)-3, a[0]==13}, a[n], n ]