Finding a recurrence for a sum

Question

I am trying to implement the following sum using a programming language: $$\sum_{i=1}^N a^i i^r$$ where $N$, $a$ and $r$ are integers.

The problem is, I cannot find a suitable way to do this. Considering $S_r$ as the above sum, I found the following recurrence: $$S_r=\frac{a^{N+1}(N+1)^r-a\left(1+\displaystyle \sum_{j=0}^{r-1} {r\choose j}S_j\right)}{a-1}$$

But if I use the above recurrence, I cannot evaluate the sum under the given time constraints. So I need to find another recurrence that could help me solve the task but I don't know how to find it. I think that the term: $$\sum_{j=0}^{r-1} {r\choose j}S_j$$ can be simplified.

To clarify a bit, I am looking for a way to evaluate $S_r$ in $O(r)$ time complexity, the current time complexity is $O(r^2)$.

Any help is appreciated. Thanks!

Answer 1

I sense here the Eulerian numbers and Eulerian Polynomials at work.

Consider the following small GNU Maxima script

display2d : false;

load("simplify_sum");

for r : 1 thru 10 do block(
    h : ratsimp(factor(simplify_sum(sum(a^i*i^r,i,1,n))*(a-1)^(r+1))),
    can : ratcoef(h,a^n),
    print(r,"|",ratexpand(factor((h-can*a^n)/a*(-1)^(r+1))))
);

which writes $S_r$ from above as $$ S_r = \left[p_r(a,N)\:a^{N+1} + a \: (-1)^{r+1} A_{r-1}(a)\right]/(a-1)^{r+1} $$

with p_r(a,N) polynomials of degree $2r$ and $A_{r-1}(a)$ Eulerian polynomials.

See e.g. https://de.wikipedia.org/wiki/Euler-Zahlen as a reference for the Eulerian polynomials.

Another reference would be http://oeis.org/A008292.

I would first try to evaluate these $A_{r-1}(a)$ in linear time.

Answer 2

First, note that even if you obtained an algorithm of complexity $O(r)$, this would still not be an efficient one, because "fast" means $O(\log ^p r)$ for some integer $p$. $O(r)$ still means exponential time, because $r = 2^{\log r}$ which is exponential in $\log r$. As customary in computer science, when I write $\log$ I mean $\log _2$.

Second, your problem is complicated by the fact that neither of $r$, $a$, $N$ is fixed, i.e. they are all free to grow, and you probably want your algorithm to be fast with respect to all the three of them (not just $r$). Not an easy task.

In order to devise an algorithm for computing your sum, let us first try to compute it modulo some prime number $p$ (that we shall later choose to be "small").

Performing Euclidean division, there exist $A, B$ such that $N = Ap + B$ with $0 \le B < p$. If $\hat a$ is the remainder of $a$ modulo $p$ and $\hat r$ the remainder of $r$ modulo $p-1$ (yes, $p-1$, not $p$), then your sum becomes

$$\sum \limits _{i=1} ^p a^i i^r + \sum \limits _{i = p+1} ^{2p} a^i i^r + \sum \limits _{i = 2p+1} ^{3p} a^i i^r \dots + \sum \limits _{i = (A-1)p+1} ^{Ap} a^i i^r + \sum \limits _{i = Ap+1} ^{Ap+B} a^i i^r = \\ \sum \limits _{i=1} ^p a^i i^r + \sum \limits _{i=1} ^p a^{i+p} (i+p)^r + \sum \limits _{i=1} ^p a^{i+2p} (i+2p)^r + \dots + \sum \limits _{i=1} ^p a^{i+(A-1)p} (i+(A-1)p)^r + \sum \limits _{i=1} ^B a^{i+Ap} (i+Ap)^r \equiv \\ \sum \limits _{i=1} ^{p-1} \hat a^i i^ {\hat r} (1 + \hat a ^p + \hat a ^{2p} + \dots + \hat a ^{(A-1)p}) + \hat a ^{Ap} \sum \limits _{i=1} ^B \hat a ^i i^{\hat r} = \\ \frac {\hat a ^{Ap} -1} {\hat a^p -1} \sum \limits _{i=1} ^{p-1} \hat a ^i i^{\hat r} + \hat a ^{Ap} \sum \limits _{i=1} ^B \hat a ^i i^{\hat r} .$$

Let's explain a number of things:

• at the end of the formula spanning lines 2 and 3 there is a reduction modulo $p$

• the upper bound in the sum has changed from $p$ to $p-1$ because terms with $i=p$ are $0$ modulo $p$

• the above is valid if $\hat a ^p \ne 1 \mod p$ (which in fact means $\hat a \ne 1 (\mod p)$). If $\hat a = 1 \mod p$ then replace the fraction by $A$

• the jump from $i^r$ to $i^{\hat r}$ is justified by Fermat's little theorem: $x^{p-1} \equiv 1 (\mod p)$ whenever $x \ne 0 (\mod p)$ (and this last condition is fulfilled since $i$ stops at $p-1$). Since by Euclidean division there exist $q$ such that $r = q(p-1) + \hat r$, then $i^r = i^{q(p-1) + \hat r} = i^{\hat r}$

• for exponentials use binary exponentiation performing a reduction modulo $p$ at each step (or any other fast modular exponentiation that you feel comfortable with).

The two sums on the last line should be fast to compute because all the numbers involved $p, \hat a, \hat r, B$ are "small" (i.e. less than $p$ which, itself, is "small").

The question is now how to construct the value of your sum knowing it modulo various "small" primes $p$. We shall use the Chinese remainder theorem, that I shall assume you to be familiar with (if not, it is widely described in various places on the internet). The crucial fact here is that the result is supposed to be computed modulo $10^9 + 7$.

Find $d$ prime numbers $p_1 \le \dots \le p_d$ such that $p_1 p_2 \dots p_d \ge 10^9 + 7$. How large must $d$ be? Well, impose the condition $p_1 ^d \ge 10^9 + 7$, and if you take $p_1 = 2$ you may choose $d = \lfloor \log (10^9 + 7) \rfloor +1 = 30$. Since the $30$-th prime number is $113$, this amounts to choosing $p_1 = 2, p_2 = 3, p_3 = 5, \dots, p_{30} = 113$. Of course, you may choose a different set of prime numbers, making them fewer but larger. The time taken by a concrete program will probably vary, but not its complexity class.

For each $1 \le j \le 30$, do your computation modulo $p_j$ as described above. Finally, apply the Chinese remainder theorem to the $d=30$ results thus obtained. This will produce a (unique) result between $0$ and $p_1 p_2 p_3 \dots p_{30}$. Finally, take the remainder of this result modulo $10^9 + 7$. This is the answer to the problem.

Some final comments: as you see, not requiring the solution modulo some number ($10^9 + 7$, which happens to be prime but this doesn't help us) would have made the problem intractable for large numbers. The idea above was to decompose the problem into many similar problems of small size, solve each of these and combine their results. This is a general idea in programming, and the approach above is often used when dealing with large integer computations. Keep it in mind for future problems.

PS: Are you cheating and looking here for answers of problems from Project Euler?