SIGFIB - Editorial



Author: Devendra Agarwal
Tester: Praveen Dhinwa and Hiroto Sekido
Editorialist: Lalit Kundu






Find the value for the expression ∑ (6*x*y*z*Fibo[x]*Fibo[y]*Fibo[z]) where x+y+z = N.


Free Advice: Please sit with a pen/pencil and a copy to grasp the solution for the problem.

We will denote Fibo[x] with F[x] , and T[x,y,z] = Fibo[x]*Fibo[y]*Fibo[z] in the whole editorial.

= ∑ 6*x*y*z*T[x,y,z]

= ∑ 6*x*y*(N-x-y)*T[x,y,z] //x+y+z = N

= ∑ 6*x*y*N*T[x,y,z] - 6*x2*y*T[x,y,z] - 6*x*y2*T[x,y,z] //Follows from previous step

From symmetry of x and y, we can say that

∑ 6*x2*y*T[x,y,z] = ∑ 6*x*y2*T[x,y,z] ,where x+y+z = N

Answer = ∑ 6*x*y*N*T[x,y,z] - 2*∑ 6*x*y2*T[x,y,z]

Let us now bifurcate to solve both expressions individually,i.e.

Let E1 = ∑ 6*x*y*T[x,y,z] , E2 = ∑ 6*x*y*N*T[x,y,z]

Answer = E1*N - 2*E2


E1 =

= ∑ (2*x*y + 2*y*z + 2*z*x)*T[x,y,z] //Reason: symmetry

= ∑ ((x+y+z)2-x2-y2-z2)*T[x,y,z] //Formula of 2*x*y + 2*y*z + 2*z*x

= ∑ (N2 - 3*x2)*T[x,y,z] //Reason: x+y+z= N , and from symmetry of x,y,and z.

= N2*∑ T[x,y,z] - ∑ 3*x2*T[x,y,z]


E2 =

= ∑ 6*x2*y*T[x,y,z]

= ∑ (3*x2*y+6*y2*x)*T[x,y,z] //Symmetry

Using the Identity: (x+y+z)3 = x3 + y3 + z3 + 3*x2*y+6*y2*x + 3*(x+y)*z*(x+y+z) , and also using the fact that x+y+z = N , we get :

= ∑ (N3 - (x3 + y3 + z3 + 3*(xy+yz)*N))*T[x,y,z]

= ∑ (N3 - (3*x3 + 3*(2*x*y)*N)*T[x,y,z] //Symmetry Arguments

= N3* ∑T[x,y,z] - 3*∑ x3*T[x,y,z] - N* ∑ 6*x*y*T[x,y,z]

= N3* ∑T[x,y,z] - 3*∑ x3*T[x,y,z] - N*E1

For calculating E2 , you need E1 , ∑T[x,y,z] and ∑ x3*T[x,y,z]. For calculating E1 , you need ∑ T[x,y,z] and ∑ x2*T[x,y,z].

Overall for calculating answer , all we need is to calculate these three terms : ∑T[x,y,z] , ∑ x2*T[x,y,z] and ∑ x3*T[x,y,z].

Let us try to solve ∑T[x,y,z] :

First off all , you should write a small code to find the pisano period of all numbers from 1 to 100000, you will note that the maximum pisano period is 375000. Maximum value of Period / M <= 6.
Now , you must have realised that the solution is aiming at Order(period of fibo mod M )
Let f denotes the fibonacci series and let P denotes it’s period.
Consider the polynomial G(x):
(f[0] + f1∗x + f2∗x2 + … + f[P-1]∗xP-1 )2
Coefficient of xi [ i< P ] is f[0]*f[i]+f1*f[i-1]+…+f[i]*f[0]
Coefficient of xi [ i ≤ P and i < 2*P ] : f[P-1]*f[i-(P+1)] +…+f[i-(P+1)]*f[P-1]
You can easily find the coefficients in O(N) time instead of using FFT algorithm.
Coefficient[i]=(Coefficient[i+2]-Coefficient[i+1]-f[i+1])mod M for i ≤ P − 3
Try to find out for i>=P.

Now,try to find out the solution to ∑T[x,y,z] using the polynomial.

let al = N mod P , I = N/P , NP = N mod P
for(int i=0;i<P;i++)
	s = (N-i)/P
	ns1 = ( ns1+  f[i]\*z1\*normal(s) )%M;		//normal(s) is the series : Sigma(s-t) , t goes from 0 to s-1
		ns3 = (ns3+f[i]\*z2)%M;
	al=(al-1+P) % P;

s1 = (s1+s2\*I+s3) % M;

s1 contains the value of ∑T[x,y,z]. You need to start from small N and M , try to find the pattern and code it.You may find the pseudo code to be beneficial.

For solving the other two version , it is quite straight forward from here. You will get a different series to solve in both cases and the overall structure remains almost the same.

Below is the pseudo code for ∑ x2*T[x,y,z]

let al = N mod P , I = N/P , NP = N mod P
for(int i=0;i<P;i++)
	s = (N-i)/P
	ss1 = ( ss1+  f[i]\*z1\*square( i , s) )%M;        //square(i,s) is the series : Sigma( ((i+t\*p)^2)\*(s-t)) , t goes from 0 to s-1
	ss2 = ( ss2+ fz2\*squaring(i,s) )				//squaring(x,s) is x^2 + (x+P)^2 + ....s terms
		ss3 = (ss3 + x<sup>2</sup>\*f[i]\*z2)
	al=(al-1+P) % P;

ss1 = (ss1+ss2+ss3) % M;

I hope you can find for ∑ x3*T[x,y,z] yourself.

I saw some awesome solutions , I would like to invite every participant to share their awesome ideas with us :slight_smile: .


O(period) per test case = O(M) per test case.


Author’s solution
Tester’s solution


Awesome problem!!! I guess there can be a lot of ways to solve this problem. In my case I derived a formula for the nth term by solving set of simultaneous equations of 12 variables taking a clue from this paper:

F(N) = ∑ (6xyzFibo[x]*Fibo[y]*Fibo[z]) = ((5N^5 + 35.N^3 - 4.N).Fibo(N) - 36.N^2.Lucas(N))/500

The problem with the above approach was that for M which have factors which are powers of 2 and 5, mod inverse of 500 does not exist. For these cases I used a different approach, the idea again from the same paper. I used the recurrence relations for the rth Fibonacci convolution sequences. I pre-computed 5th Fibonacci convolution sequence for the highest powers of 2 and 5 possible. The generating function for the 5th Fibonacci convolution sequence is (x/(1-x-x^2))^6, whereas our sequence has a generating function of 6.(x.(1+x^2)/(1-x-x^2)^2)^3 or 6.(x^-3 + 3.x^-1 + 3.x^1 + x^3).(x/(1-x-x^2))^6.
So we can obtain F(N) from the 5th Fibonacci convolution sequence by using the formula F(N) = 6. (F5(N-3) + 3.F5(n-1) + 3.F5(N+1) +F5(N+3)) , where F5(N) is the 5th Fibonacci convolution sequence as pre-computed above.

Finally, to address how we arrived at the generating function of our sequence, we notice that our function is the result of convolution of the sequence N.Fibo(N) with itself 2 times. Generating function of the sequence N.Fibo(N) is the first derivative of the Fibonacci sequence shifted by 1 position whose generating function is x/(1-x-x^2) and hence has a generating function x.(1+x^2)/(1-x-x^2)^2
Hence the generating function for our sequence is 6 times the cube of this as we are convolving 2 times.

this is how i solve it(chinese)


I went the same route but got stuck on the case where M 2 ==0 || M 5 == 0 Could you please provide the link of your solution?

omg awesome explanation!!! I really get blown away by how truly talented some people are…and even though I’ll never reach this level, just having the chance to learn a few things is already awesome for me :smiley: Thanks for this


Here’s my approach (first look at my solution here: ). It seems I’m pulling a lot of numbers out of nowhere, so I’ll try to why all these numbers appear.

The first step I did was try to find a generating function for this sequence. This is because it is very likely to correlate to a recurrence, and if we have a recurrence, we can get the Nth term in log N time using matrix exponentiation (if you know how to compute fib(n) in log n time, then this will kind of explain how to do this). So, let’s try to find the generating function.

As @n2n_ has already pointed out, the generating function for this sequence is (x(x^2+1) / (1-x-x^2)^2)^3 (read his explanation on how to derive it, I pretty much did the same thing). Now knowing how to convert generating functions into a recurrence is very useful. In general, if a sequence a(n) has generating function 1/p(x) where p(x)=1 - c1 * x - c2 * x^2 - … - ck * x^k, then the corresponding recurrence is a(n) = c1 * a(n-1) + c2 * a(n - 2) + … + ck * a_(n-k). Now, what about sequences with generating functions of the form q(x)/p(x). Well, let’s say q(x) = b0 + b1x + b2x^2 + … + bm x^m. Then, if we have already derived a(n) above from just 1/p(x), our actual value of the nth term a’(n) = b0 a(n) + b1 a(n - 1) + … + bm a(n-m).

Thus, we have a general form for the recurrence given that the numerator and denominator are polynomials, which they luckily are. We fully expand these out, to get the recurrence relationship as described above. This recurrence will have each term rely on the 12 terms before this. As for the values of a(0) to a(11), we can brute force the taylor expansion of 1/p(x) to get these. Thus given the recurrence, we can solve each case in log N time. This is the basic idea of my solution, but I had some trouble getting this to just pass because of a constant factor, so I did some optimizations (pre computing some periods) to get it to eventually pass.


If you use Cayley-Hamilton theorem, you could replace matrix multiplication by polynomial multiplication, which will give a 12X speedup in this case.


Thanks for the tip! I’ll be sure to read up on that when I have more time since I don’t quite understand that now.

Access denied in testers soln.

I got the generating function, and I found a 12th order recursion that described the sequence. Sadly, I couldn’t do necessary matrix multiplications within the time limit. Alas, it seemed like a decent approach at the time.

1 Like

Great work :slight_smile:
I was going for this problem, did a bit of mathematics but other problems took a lot of time. A good way to learn, Thanks :smiley:

That’s really a fancy solution!

@djdolls, I’ve been thinking on how to use Cayley-Hamilton Thm to obtain the speedup, but couldn’t find out. Could you explain how you would do it? Thanks!

I’ll try to explain. Let’s try to motivate this first. Suppose we have a matrix A and we want to find A^1000. Suppose we also know that p(A) = 0 for some polynomial p. To demonstrate, let’s suppose p(A) = A^2 - 2A + 1 = 0. Now, we know A^1000 - 2A^999 + A^998 = 0, so A^1000 = 2A^999 - A^998. Then, we can write A^999 and A^998 in terms of smaller powers and so on until we just get things in terms of A^1, and 1. Now, this is much easier to evaluate, and is very fast (a 12x speedup in this case). Now, how do we find this characteristic polynomial? We just use the Cayley-Hamilton theorem.

Do you have references to learn how to convert generating functions to recurrence relations? I asked in SE.math but got no answer.

I have a simple, but manual, way of doing it. I’ll make an example with the fibonacci generating function.

First put the GF in that form: p(x)/(1-f(x)), so: 1/(1-(z+z^2)).

The recurrence relation is pretty simple now: a[n]=a[n-1]+a[n-2], had it been 3x/(1-(z+2z^2+3z^3)) it would have been a[n]=a[n-1]+2a[n-2]+3a[n-3].

Put the GF in wolframalpha to do the Taylor expansion and obtain the first terms, in the fibonacci series you would obtain: 1+z+2z^2+… so first terms would be a[0]=1,a[1]=1.

In case you don’t have access to internet you will have to play a little bit: p(x)/(1-f(x)) = GF(x) try to guess the first coefficients of GF(x) so that when you multiply if by (1-f(x)) you get p(x).

@lg5293 wouldn’t your algorithm be O(n)? Sure the constant would be smaller but in this case n=10^18.

please see

I have explained how to use Cayley-Hamilton theorem to solve linear recurrence faster than matrix exponentiation method.

nice explanation, thanks!