How to make DifferenceRoot and RecurrenceTable useful for non-numeric difference equations?

Question

In answering a physics forum question this morning, I ran into really bad performance of DifferenceRoot and RecurrenceTable compared to calculating the expressions by naively taking derivatives of an exponential generating functional. A very small amount of digging showed that DifferenceRoot and RecurrenceTable do not simplify expressions as they go.

For example, look at the following output of RecurrenceTable and how it simplifies by just Expanding the result:

In[1]:= RecurrenceTable[f[n] == a f[n - 1] + (a - 1) f[n - 2] && 
                        f[0] == 0 && f[1] == 1, 
                        f, {n, 6}]
% // Expand

Out[1]= {0, 1, a, -1+a+a^2, -a+a^2+a (-1+a+a^2), 1-a-a^2+a (-1+a+a^2)+a (-a+a^2+a (-1+a+a^2))}
Out[2]= {0, 1, a, -1+a+a^2, -2 a+2 a^2+a^3, 1-2 a-2 a^2+3 a^3+a^4}

This quickly gets out of hand, as the leaf count of the 20th iteration (calculated using DifferenceRoot) shows:

dr[k_] := DifferenceRoot[Function[{f, n}, 
          {f[n] == a f[n - 1] + (a - 1) f[n - 2], f[0] == 0, f[1] == 1}]][k]

In[2]:= dr20 = dr[20]; // Timing
        dr20Exp = Expand[dr20]; // Timing
Out[2]= {0.26, Null}
Out[3]= {2.39, Null}

In[4]:= {LeafCount[dr20], LeafCount[dr20Exp]}
Out[4]= {1188383, 92}

Which can be compared to the memoized implementation

In[1]:= mem[n_] := a mem[n-1] + (a-1) mem[n-2] // Expand
        mem[0] = 0; mem[1] = 1;

In[3]:= mem20 = mem[20];//Timing
        LeafCount[mem20]
Out[3]= {0.48, Null}
Out[4]= 92

So my question is: Are there any options/tricks to get DifferenceRoot and RecurrenceTable to apply a (simplifying) function as they go and thus make them useful for non-numeric work?

Edit: A Sjoerd pointed out below, I foolishly chose an example with a RSolveable closed form solution. In this question I'm primarily concerned with the behaviour of DifferenceRoot and RecurrenceTable. If it helps, imagine the the f[n-2] term is multiplied by n, so that there is no simple closed form solution.

Answer 1

I can't really help with your question as I haven't used those functions until now, and the docs don't give a clue. But why don't you just use RSolve here? It gives a closed form solution for each of the table's elements:

sol = f /. RSolve[f[n] == a f[n - 1] + (a - 1) f[n - 2] && 
                  f[0] == 0 && f[1] == 1, f, n
           ][[1, 1]]

sol@Range[6] // Simplify