# RSolve and a Piecewise Recurrence + Runtime

This should really be a comment to this answer but I don't have enough rep. I'd like to know how to do the same thing (efficiently - see below) but with $$\mod 3$$ conditions, i.e. define a[k] to be the $$0\mod 3$$ index terms of the function, b[k] (resp. c[k]) the $$1\mod 3$$ (resp. $$2\mod 3$$) indexes. Copying the method in the linked answer works for some simpler recursive formulas. However, for the more complicated formulas that I'm interested in, the process does not seem to stop, even given an hour + to run. What can I do about this?

Edit: Here's a more specific version of my question. I have something like the following recursive formula for a sequence $$a_n$$: $$a_n=\begin{cases} 3a_{n-1}+3^{n+3}:& n\not\equiv 1\mod 3\\ 3a_{n-1}+5\cdot 3^n:& n\equiv 1\mod 3 \end{cases}$$ and $$a_1=1$$. I want code that gives a closed form for $$a_n$$ (probably will be piecewise defined, that's fine). The answer that I linked does this for $$\mod 2$$ conditions, and anything higher seems to complicate matters so much that it doesn't terminate in time.

• If you just want an answer, the ordinary generating function for the sequence $a_n$ is (x*(1 + 243*x + 729*x^2 + 378*x^3))/((-1 + 3*x)^2*(1 + 3*x + 9*x^2)). First few terms are {1, 246, 1467, 4806, 20979, 82620, 258795}. – Somos Feb 4 at 3:19
• @Somos How did you obtain this? – user234988 Feb 4 at 3:21
• I actually used my very quick PARI/GP code. It would take me much more time to get Mathematica to get similar results. – Somos Feb 4 at 3:22
• Mind sharing? Unless that's against site rules as we're on the mathematica SE. – user234988 Feb 4 at 3:25
• Recent versions of PARI/GP have a function which does the job. You can use t=seralgdep(A,1,4) where A = y + 246*y^2 + ... + 953532*y^8 + O(y^9) and t=(81*y^4 - 27*y^3 -3*y + 1)*x + (-378*y^4 - 729*y^3 - 243*y^2 - y) gives the rational form by solving $t=0$ for $x$ in terms of $y$. – Somos Feb 4 at 3:38

Okay, the problem is easier than it looks. You just have to generate a few terms of the sequence. To define the sequence just use code like this:

a[1] = 1; a[n_] := a[n] = 3 a[n-1] + 3^n If[Mod[n, 3]==1, 5, 3^3];


To find a suitable rational generating function just use code like this:

rule = Solve[0 == CoefficientList[ Normal[Sum[a[n] x^n, {n, 9}]
(1 + q[1] x + q[2] x^2 + q[3] x^3 + q[4] x^4) -
x (1 + p[1] x + p[2] x^2 + p[3] x^3 + p[4] x^4) +
O[x]^10], x],
{p[1], p[2], p[3], p[4], q[1], q[2], q[3], q[4]}][[1]];


which returns the result

{p[1] -> 243, p[2] -> 729, p[3] -> 378, p[4] -> 0,
q[1] -> -3, q[2] -> 0, q[3] -> -27, q[4] -> 81}}


and then your generating function is given by

x (1 + p[1] x + p[2] x^2 + p[3] x^3 + p[4] x^4)/(1 + q[1] x +
q[2] x^2 + q[3] x^3 + q[4] x^4) /. rule