While solving Project Euler problem #402, I (among many others who soled or tried) have stumbled upon the problem of calculating the sum of the powers of first Fibonacci numbers.

Wikipedia article on Fibonacci numbers gives two formulas from this family:

The sum of the first powers (i.e. of the numbers themselves):
$$\sum_{i=0}^n f_i = f_{n+2}-1$$

And the sum of the squares:
$$\sum_{i=0}^n f_i^2 = f_n f_{n+1}$$

The simplicity of these expressions suggests that similar closed form solutions exist for other powers too, and it is indeed so. Solving the above-mentioned problem, I wrote a Python3 + Sympy program, which can derive such expressions for arbitrary large powers (well, there are limits of course: computer memory and your time, but powers around 30 can be resolved in reasonable time). It seems that these formulas are not widely represented in the Internet, so I decided to share them. Who knows, maybe I'll save somebody's 5 minutes?

Here are the expressions for the sums of the powers from 1 to 9. (as image)

$$\sum_{i=0}^n f_i = f_{n+2}-1$$

$$\sum_{i=0}^n f_i^2 = \frac{f_{2n} + f_{2n+2} - \left(-1\right)^n}{5}$$

$$\sum_{i=0}^n f_i^3 = \frac{f_{3n+2}-6(-1)^n f_{n-1}+5}{10}$$

$$\sum_{i=0}^n f_i^4 = \frac{1}{25} \left(f_{4n+2}-4(-1)^n f_{2n+1}+6n+3\right)$$

$$\sum_{i=0}^n f_i^5 = \frac{1}{550}\left(4f_{5n+2}+2f_{5n+4}-55(-1)^n f_{3n+1}+220f_{n+2}-175\right)$$

$$\sum_{i=0}^n f_i^6 = \frac{1}{500}\left(f_{6n+2}+f_{6n+4}-8(-1)^n\left(f_{4n+1}+f_{4n+3}+5\right)+60f_{2n}+60f_{2n+2}\right)$$

$$\sum_{i=0}^n f_i^7 = \frac{1}{79750}\left(22f_{7n+2}+88f_{7n+4}-406(-1)^n\left(55f_{n-1}+f_{5n+1}+2f_{5n+3}\right) + \right. $$

$$\left.+6699f_{3n+2}+17375\right)$$

$$\sum_{i=0}^n f_i^8 = \frac{1}{1875}\left(f_{8n+4}-12(-1)^n\left(14f_{2n+1}+f_{6n+3}\right)+84f_{4n+2}+210n+105\right)$$

$$\sum_{i=0}^n f_i^9 = \frac{1}{7576250}\left(1276f_{9n+4}+319f_{9n+5}-18810(-1)^n f_{7n+3}-3762(-1)^n f_{7n+4}+\right.$$

$$\left. + 119016f_{5n+2}+39672f_{5n+3}-509124(-1)^n f_{3n+1}+1527372f_{n+2}-1173125\right)$$

Because there are algorithms that can calculate n'th Fibonacci number in O(log(n)) steps, the above formulas can be used for effective calculation of such sums for very large n, where direct summation is not feasible.

## Derivation

To derive these closed-form representations of \(sum_{i=0}^n{f_i^p}\), I used the following universal approach:

- Apply Binet's formula t: \(f_n=\left(\phi^n-\psi^n\right)/\sqrt{5}\), where \(\phi=-1/\psi=\left(1+\sqrt{5}\right)/2\).
- Expand the power using the binomial formula. \(f_i^p\) will become sum of several exponents.
- Apply the formula for the sum of exponents : \(\sum_{i=0}^n a^i=(a^{i+1}-1)/(a-1)\)
- Go back from the powers of \(\phi,\psi\) to Fibonacci numbers, using the efollowing qualities: \(\phi^n=f_{n+1}-f_n\psi\), \(\psi^n=f_{n+1}-f_n\phi\).

## Source code

All scripts are in the Github repository: dmishin/fibsums.

The repository contains several programs, main of them is `fibsums_general.py`

, which is the script for deriving a formula. For any positive power, it gives two closed-form representations: one expresses the sum via the products and powers of \(f_n, f_{n+1}\), and another - via \(f_{kn}\). Usually, the second one is simpler.

Example of derivation of \(\sum_{i=0}^n f_i^4\):

```
$python fibsums_general.py 4
==================================
Sum F(i)^4 =
= (-4*(-1)**n*F(n)**2 - 4*(-1)**n*F(n + 1)**2 + 6*n - F(n)**4 +
4*F(n)**3*F(n + 1) + 4*F(n)*F(n + 1)**3 + F(n + 1)**4 + 3)/25 =
= (-4*(-1)**n*F(2*n + 1) + 6*n - F(4*n + 3) + F(4*n + 4) + 3)/25
```

This "raw" result can be simplified even more, using known equalities, giving \((f_{4n+2}-4(-1)^nf_{2n+1}+6n+3)/{25}\).
For small powers, computation takes a few seconds, and for powers near 25, calculation requires a minute or two.

In conclusion, here are raw results for powers 1 - 25, in a text file: fibsums-25.txt.