# how can I calculate $a_{n+3}=a_n+a_{n+1}+a_{n+2}$?

recently I need to calculate this: $$a(n+1)=((a(n-2)+a(n-1)+a(n)) \bmod 10000)$$

and get $$a(20000000)$$ (for example).

I know RecurrenceTable,but

RecurrenceTable[{a[n + 1] == Mod[a[n] + a[n - 1] + a[n - 2], 10000],
a[1] == 1, a[2] == 1, a[3] == 1}, a, {n, 1, 20000000}] // Last


require a bit of long time.(and extra space)

using RSolve is also useless.

Nest seems only support $$a_{n+1}=f(a_n)$$

this is also not so useful.

a[1] == 1;
a[2] == 1;
a[3] == 1;
a[n_] := a[n] = Mod[a[n - 3] + a[n - 1] + a[n - 2], 100000];
a[20000000]


So is there a way which can Space Complexity be $$O(1)$$ as well as fast? (I know matrix exponentiating by squaring,but it seems hard to write.)

Nest[] is in fact usable here, if you are ingenious enough and willing to wait a little:

Last[Nest[Append[Rest[#], Mod[{1, 1, 1}.#, 100000]] &, {1, 1, 1}, 20000000 - 3]]
16287


An alternative is to use an undocumented function for modular matrix exponentiation:

Mod[First[AlgebraMatrixPowerMod[{{1, 1, 1}, {1, 0, 0}, {0, 1, 0}},
20000000 - 3, 100000].{1, 1, 1}], 100000]
16287


You need to use Set rather than Equal for the initial conditions.

Clear[a, n]

a[1] = 1;
a[2] = 1;
a[3] = 1;
a[n_] := a[n] = Mod[a[n - 3] + a[n - 1] + a[n - 2], 100000];
Last[a /@ Range[4, 20000000]] // AbsoluteTiming

(* {107.222, 16287} *)


Similar approach to JM...

mat = {{1, 1, 1},
{1, 0, 0},
{0, 1, 0}};

init = {1,1,1};

next[vec_] := Mod[mat.vec, 100000];

First@Nest[next, {1, 1, 1}, 20000000-3] // Timing

(*  {42.0625, 16287}  *)


Inspired by JM's second approach, but using what is documented...still pretty fast. I thought it might overflow, but no problem.

First@Mod[Mod[MatrixPower[mat, 20000000 - 3], 100000].{1, 1, 1}, 100000] // Timing

(* {0.03125, 16287} *)

• You can simplify your second approach by using the action form of MatrixPower[]: First[Mod[MatrixPower[mat, 20000000 - 3, {1, 1, 1}], 100000]]. Commented Jan 24, 2020 at 3:07
• Alternatively there is AlgebraMatrixPowerMod. Unfortunately it does not (so far as I am aware) take that optional vector argument. But at least it limits the intermediate integer swell. Commented Jan 24, 2020 at 17:29
• Yup, I was riffing on JM's answer, which uses it. That one was "too cheap to meter" time-wise. Very nice. Commented Jan 24, 2020 at 18:31