Strategies for simplifying recurrences with sums

The problem:

Count the number of permutations of n distinct objects that leave none of them fixed.

This is actually a well-studied problem with pretty well-established terminology (derangements). Mathematica even has a built-in function for it, Subfactorial[n], which it is perfectly capable of producing if you feed it a simple enough recurrence:

RSolve[{a[n] == (n - 1) (a[n - 1] + a[n - 2]),
a[0] == 1, a[1] == 0}, a[n], n] // FullSimplify
(* output: {{a[n] -> Subfactorial[n]}} *)


But when I attacked the problem, I was unable to come up with this simple recurrence, and only produced recurrences that RSolve couldn't handle.

Demonstrate how one could successfully use the analytical features (simplification, solving, reduction, etc.) of Mathematica to obtain the closed-form solution Subfactorial[n], beginning from one of my recurrences below.

(note: the intent of the last restriction above is to weed out responses that merely derive the simpler formula directly from the problem. It is okay for a solution to not meet this restriction, so long as it doesn't require me to have a Jimmy Neutron-style "brain blast!" to see it)

I do not mind if a bit of guessing, trial and error, or manipulation-by-hand is required. I am merely seeking techniques applicable to this problem which may also be applicable to others.

My first attempt: Count ways to partition into cycles of size $> 1$.

Suppose that you have $N$ objects whose fates have yet to be decided. From these, pick elements to form $m$ cycles ($m\ge1$) of size $b>1$, dividing by $m!$ to account for the order in which the cycles can be chosen. To ensure the same $b$ is never considered twice, a second parameter $b_\text{max}$ is used to ensure that $b$ strictly decreases between recursive evalutaions.

$$\text{RSolve}\left[\left\{a[n,~b_\text{max}]=\sum _{b=2}^{b_\text{max}} \sum _{m=1}^{\left\lfloor \frac{n}{b}\right\rfloor } \frac{n! a[n-b m,~b-1]}{m! b^m (n-b m)!},~a[0,~b]=1\right\},~a[j,~k],~\{j,~k\}\right]$$

(* % // InputForm *)
RSolve[{a[n, bmax] == Sum[(a[-(b*m) + n, -1 + b]*n!)/(b^m*m!*(-(b*m) + n)!),
{b, 2, bmax}, {m, 1, Floor[n/b]}], a[0, b] == 1}, a[j, k], {j, k}]


Unfortunately, this is not even a valid way to invoke RSolve (it will merely spit back out what you typed); RSolve doesn't like multivariable recurrences.

My second attempt: A single valued recurrence.

I tried another approach using the inclusion-exclusion principle, which eventually turned into this:

To count the ways to arrange $N$ objects with no fixed points, subtract from $N!$ the number of arrangements with at least one fixed point. For the latter, simply pick $m\ge1$ fixed points, and then the remaining $N-m$ have no fixed points (giving a recurrence):

$$\text{RSolve}\left[\left\{f[n]=n!-\sum _{m=1}^n \binom{n}{m} f[n-m],~f[0]=1\right\},~f[n],~n\right]$$

RSolve[{f[n] == n! - Sum[Binomial[n, m]*f[-m + n], {m, 1, n}], f[0] == 1}, f[n], n]


But again, Mathematica just spits the same form back out.

• Its strange on my computer I have no problem to obtain a solution in terms of BesselK functions a[n] -> (BesselI[n, -2] BesselK[1, 2] + BesselI[1, 2] BesselK[n, 2])/( BesselI[1, 2] BesselK[0, 2] + BesselI[0, 2] BesselK[1, 2]). What is your Mathematica version ? Commented Sep 19, 2016 at 4:06
• @cyrille.piatecki my version is 10.4 Commented Sep 19, 2016 at 14:44

Not precisely solving the recurrences from the question, but a different approach to getting the formula from the very definition.

Make a condition such that elem is on a kth position; will be used to exclude such a situation:

cond[k_, elem_] := Table[_, k - 1]~Join~{elem}~Join~{___}
a[n_] := Alternatives @@ Table[cond[i, i], {i, 1, n}] /; n > 1
a[1] := {_}


Let's take for example lists of length up to max = 10:

seq = Length /@ Table[DeleteCases[Permutations[Range[n]], a[n]], {n, 1, max}]


{0, 1, 2, 9, 44, 265, 1854, 14833, 133496, 1334961}

This is the same as Subfactorial:

Subfactorial /@ Range[max]


{0, 1, 2, 9, 44, 265, 1854, 14833, 133496, 1334961}

One can use FindSequenceFunction to guess what's the origin of sequence seq:

FindSequenceFunction[seq]


Subfactorial

In fact, max = 4 is sufficient for FindSequenceFunction to return the correct solution.

EDIT: I'd like to emphasize that FindSequenceFunction is a guess on what the underlying formula might be. As an illustrative example, consider a polynomial

f[x_] := -7 + 12 x - 6 x^2 + x^3


Its first 4 values are

f /@ Range[4]


{0, 1, 2, 9}

which happen to coincide with Subfactorial. On the other hand, f /@ Range[5] is {0, 1, 2, 9, 28}, and applying FindSequenceFunction yields (correctly) f[x]. Hence, FindSequenceFunction, and a somewhat related FindFormula, should be used with caution and treated more as hints rather than a definite solution.

• Huh, never would've expected there to be a feature that can guess analytical solutions from brute force values! I have clarified the "from my recurrences" criterion, because this does meet the rule in spirit. Commented Sep 19, 2016 at 15:08