Summation Involving Product Of Two Identical Polynomials. Recently I stuck, to a problem. However I rarely think that there is some proper formula for this problem, but here I am in search of algorithm's or theorem that relate to this problem or can solve this problem. 
We have three integers  positive integers $a, b, n$ and we need to compute this summation: 
$$\sum\limits_{k = 1}^n{(k^a - {(k - 1)}^a)(k^b - {(k - 1)}^b)}$$ 
It looks simple, if we try to solve it using computer, but real problem lies in it's constraints: 
$$n < 10^{12}$$
$$a < 10^4$$
$$b < 10^4$$ 
Clearly, we cannot iterate from $k = 1$ to $k = n$, thus we need to think of an algorithm that solve it in constant time
or in complexity in terms of $a$ and $b$ 
One can easily find relation like this $\sum\limits_{k = 1}^n{(k^a - {(k - 1)}^a)} = n^a$ but the problem has product of two such terms. So, Please give me some suggestion to solve this problem.
 A: This is not a full answer because the bounds for $a$ and $b$ given by OP imply that there are situations where $a>>n$ and $b>>n.$ (For example, $n=10, a=1000, b=2000.$)  The following asymptotic formula works well if $n>>a$ and $n>>b.$  Given the square of bounds for the parameters for $a$ and $b$ relative to $n$, the formula actually covers about all but 1 part in 10^16 of the domain.
$$ S:=\sum_{k=1}^n \big(k^a -(k-1)^a)\big(k^b -(k-1)^b) \sim
\frac{n^{a+b-1} a\,b}{a+b-1} \Big( 1 - \frac{(a-1)(b-1)}{12 n^2} \frac{a+b-1}{a+b-3)} \Big) $$
The proof is easily enough constructed by using a few terms of the Faulhaber formula,
$$ (*) \quad \sum_{k=1}^n k^m = \frac{n^{m+1}}{m+1} + \frac{n^m}{2} + m \, \frac{n^{m-1}}{12} + \binom{m+1}{3}\frac{B_3}{m+1}n^{m-2} + ... $$
and the binomial theorem.  Note that the Bernoulli number $B_3$ is 0. To go beyond the two terms given in the asymptotic formula, you'll need more terms in the binomial expansion and the Faulhaber formula.  By a simple series arrangement,
$$ S= - n^{a+b} + \sum_{k=1}^n 2 k^a - k^a(k-1)^b - k^b(k-1)^a $$  Define a symbol $$(a,b)_k := \binom{a}{k} + \binom{b}{k}.$$  Then to a 4th order binomial expansion,
$$S=-n^{a+b} + (a+b)\sum_{k=1}^n k^{a+b-1} - (a,b)_2\sum_{k=1}^n k^{a+b-2} 
+ (a,b)_3\sum_{k=1}^n k^{a+b-3}  - (a,b)_4\sum_{k=1}^n k^{a+b-4} $$ 
Use (*): 
$$ S= - n^{a+b}  + (a+b)\Big(\frac{n^{a+b}}{a+b} + \frac{n^{a+b-1}}{2} +  \frac{a+b-1}{12}n^{a+b-2} + 0 \cdot n^{a+b-3} + ... \Big) $$
$$- (a,2)_n \Big( \frac{n^{a+b-1}}{a+b-1} + \frac{n^{a+b-2}}{2} +  \frac{a+b-2}{12}n^{a+b-3} + ... \Big)$$
$$+(a,3)_n \Big( \frac{n^{a+b-2}}{a+b-2} + \frac{n^{a+b-3}}{2} + ...\Big) -
(a,4)_n \Big( \frac{n^{a+b-3}}{a+b-3} + ...\Big)$$
The ellipses indicate terms that are not shown, because the powers of $n$ become too small to appear in the final answer. The first 2 terms cancel. Organize the rest as a descending series in powers of $n:$
$$ S = n^{a+b-1}\Big( \frac{a+b}{2} - \frac{(a,b)_2}{a+b-1} \Big) + n^{a+b-2}\Big(
\frac{(a+b)(a+b-1)}{12} - \frac{(a,b)_2}{2} + \frac{(a,b)_3}{a+b-2} \Big)$$
$$+n^{a+b-3}\Big(-\frac{(a+b-2)}{12}(a,b)_2 + \frac{(a,b)_3}{2} - \frac{(a,b)_4}{a+b-3} \Big) + ...$$  Simplification of the coefficients gives the answer of the final form.
It can be seen that once $a$ or $b$ become comparable to $n,$ the second term will be as large as the first and so the asymptotic solution begins to fail.  Also, we must have a+b>3.
As an example, for n=5000, a=20, b=60:  The first term of the asymptotic expansion gives about 5 digits precision, and both terms give about 10 digits precision.
A: I think writing as 
$$n + \sum_{k=1}^n{\frac{k-1}{k}^{ab}-\frac{k-1}{k}^{a}-\frac{k-1}{k}^{b}}$$
more pertinent to start with because you have a commun separated structure of $\frac{k-1}{k}^{x}$
A: If we expand the terms we get four terms out of which two terms can be calculated by using Faulhaber's formula(https://en.wikipedia.org/wiki/Faulhaber%27s_formula) and I am facing the same problem right now.I'm not to able to get the summation for the other two terms.
