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.

  • $\begingroup$ 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}. $\endgroup$
    – Somos
    Commented Feb 4, 2019 at 3:19
  • $\begingroup$ @Somos How did you obtain this? $\endgroup$
    – user234988
    Commented Feb 4, 2019 at 3:21
  • $\begingroup$ I actually used my very quick PARI/GP code. It would take me much more time to get Mathematica to get similar results. $\endgroup$
    – Somos
    Commented Feb 4, 2019 at 3:22
  • $\begingroup$ Mind sharing? Unless that's against site rules as we're on the mathematica SE. $\endgroup$
    – user234988
    Commented Feb 4, 2019 at 3:25
  • $\begingroup$ 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$. $\endgroup$
    – Somos
    Commented Feb 4, 2019 at 3:38

1 Answer 1


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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.