# RecursionLimit error appearing if solve the recursion formula

f[0, 0] = 1;
f[a_, b_] := (Solve[b f[a, b] - f[a, b - 1] == 0, f[a, b]] //
Flatten)[[1, 2]] /; b != 0
f[a_, b_] := (Solve[a f[a, b] - f[a - 1, b] == 0, f[a, b]] //
Flatten)[[1, 2]] /; a != 0


if we calculate f[0,2], there is a RecursionLimit error. Of course, in this case, it can be solved by writing the recursion formula explicitly, e.g.

f[0, b_] := f[0, b - 1]/b /; b != 0
f[a_, b_] := f[a - 1, b]/a /; a != 0


But if the recursive formula requires solving some equations and the results are lengthy, is there an alternative way besides explicitly writing out the recursive formula?

You can not use the symbol f[a,b] in defining f[a,b]. The following code works:

f[0, 0] = 1;
f[a_, b_] := (Solve[b  x - f[a, b - 1] == 0, x] // Flatten)[[1, 2]] /;
b != 0
f[a_, b_] := (Solve[a  x - f[a - 1, b] == 0, x] // Flatten)[[1, 2]] /;
a != 0

• Solve inside a delayed assignment is not a good idea. Commented May 27 at 6:23
• An immediate assignment is much more efficient: Clear[f]; f[0, 0] = 1; f[a_, b_ /; b > 0] = SolveValues[b x - f[a, b - 1] == 0, x][[1]]; f[a_ /; a > 0, b_] = SolveValues[a x - f[a - 1, b] == 0, x][[1]]; and adding memoization (see @BobHanlon's answer) speeds it up further. Commented May 27 at 6:37
• @Roman, I completely agree with you. I just pointed out what caused the infinite loop. Commented May 28 at 7:41

The obvious solution with RSolve or RecurrenceTable do not work (bug?):

RSolve[{f[0, 0] == 1, b f[a, b] == f[a, b - 1], a f[a, b] == f[a - 1, b]}, f[a, b], {a, b}]
(*    RSolve: There are fewer dependent variables than equations, so the system is overdetermined.    *)

RecurrenceTable[{f[0, 0] == 1, b f[a, b] == f[a, b - 1], a f[a, b] == f[a - 1, b]}, f[a, b], {a, 0, 3}, {b, 0, 3}]
(*    RecurrenceTable: There are fewer dependent variables than equations, so the system is overdetermined.    *)


But we can solve it one dimension at a time: first, for fixed $$a$$,

RSolve[{b fa[b] == fa[b - 1]}, fa[b], b]
(*    {{fa[b] -> C[1]/Pochhammer[2, -1 + b]}}    *)


So let's set a formula for $$f_a(b)=f(a,b)$$:

f[a_, b_] = c[a]/Pochhammer[2, -1 + b];


and solve for $$c(a)$$ to get the full formula:

RSolve[{f[0, 0] == 1, a f[a, b] == f[a - 1, b]}, c[a], a]
(*    {{c[a] -> 1/Pochhammer[2, -1 + a]}}    *)


The full solution is

f[a, b] /. First[%] // FullSimplify
(*    1/(Gamma[1 + a] Gamma[1 + b])    *)


which is simply $$f(a,b)=\frac{1}{a!b!}$$.

\$Version

(* "14.0.0 for Mac OS X ARM (64-bit) (December 13, 2023)" *)

Clear["Global*"]


The approach offered by A. Kato:

f[0, 0] = 1;
f[a_, b_] := (Solve[
b   x - f[a, b - 1] == 0, x] // Flatten)[[1, 2]] /; b != 0
f[a_, b_] := (Solve[
a   x - f[a - 1, b] == 0, x] // Flatten)[[1, 2]] /; a != 0

calc1 = AbsoluteTiming[(
tab1 = Table[f[a, b], {a, 0, 4}, {b, 0, 4}]) // Grid]


Recursion with memoization

f2[0, 0] = 1;
f2[0, b_Integer?Positive] := f2[0, b] = f2[0, b - 1]/b;
f2[a_Integer?Positive, b_Integer?NonNegative] := f2[a, b] = f2[a - 1, b]/a;

calc2 = AbsoluteTiming[(
tab2 = Table[f2[a, b], {a, 0, 6}, {b, 0, 6}]) // Grid]


Both methods produce the same result.

tab1 === tab2

(* True *)


However, the second approach is far more efficient

calc1[[1]]/calc2[[1]]

(* 67.3106 *)
`