Faulhaber's formula expresses a sum over some finite number of naturals to the $m^{th}$ power in terms of the Bernoulli numbers $B_{j}$ (using the $B_{1} = 1/2$ convention) or polynomials $\hat{B}_{j}$ as follows

$$\sum_{k=1}^{n} k^{m} = \frac{1}{m+1} \sum_{j=0}^{m} \begin{pmatrix} m+1 \\ j \end{pmatrix} B_{j} n^{m+1-j} = \frac{1}{m+1} \left( \hat{B}_{m+1}(n+1) - \hat{B}_{m+1}(1) \right) \tag 1$$

I recently needed to derive an expression similar to $(1)$ but for a double summation over squares $k_{1}^{2} + k_{2}^{2}$ rather than $k$, summing to arbitrary $n_{1}, n_{2}$ and with $m \in \mathbb{N} \setminus \{1\}$. To this end, using the binomial theorem and then Faulhaber's formula on the resultant product gives

\begin{align} \sum_{k_{1},k_{2}=1}^{n_{1},n_{2}} (k_{1}^{2} + k_{2}^{2})^{m} &= \sum_{k_{1},k_{2}=1}^{n_{1},n_{2}} \sum_{i=0}^{m} \begin{pmatrix} m \\ i \end{pmatrix} k_{1}^{2i} k_{2}^{2m-2i} \\ &= \sum_{i=0}^{m} \begin{pmatrix} m \\ i \end{pmatrix} \sum_{k_{1}=1}^{n_{1}} k_{1}^{2i} \sum_{k_{2}=1}^{n_{2}} k_{2}^{2m-2i} \\ &= \sum_{i=0}^{m} \begin{pmatrix} m \\ i \end{pmatrix} \frac{1}{2i+1} \frac{1}{2m-2i+1} \left[\sum_{j_{1}=0}^{2i} \begin{pmatrix} 2i+1 \\ j_{1} \end{pmatrix} B_{j_{1}} n_{1}^{2i+1-j_{1}} \sum_{j_{2}=0}^{2m-2i} \begin{pmatrix} 2m-2i+1 \\ j_{2} \end{pmatrix} B_{j_{2}} n_{2}^{2m-2i+1-j_{2}} \right] \\ &= \sum_{i=0}^{m} \begin{pmatrix} m \\ i \end{pmatrix} \frac{1}{2i+1} \frac{1}{2m-2i+1} \bigg[\left( \hat{B}_{2i+1}(n_{1}+1) - \hat{B}_{2i+1}(1) \right) \left( \hat{B}_{2m-2i+1}(n_{2}+1) - \hat{B}_{2m-2i+1}(1) \right) \bigg] \end{align}

What I then wanted to do was to use this to find a polynomial $G(l, m, n_{1}, n_{2})$ such that

$$G(l, m, n_{1}, n_{2}) = \frac{\sum_{k_{1},k_{2}=1}^{n_{1},n_{2}} (k_{1}^{2} + k_{2}^{2})^{m+l}}{\sum_{k_{1},k_{2}=1}^{n_{1},n_{2}} (k_{1}^{2} + k_{2}^{2})^{m}} \tag 2$$

for arbitrary values of $(l, m, n_{1}, n_{2})$ with $l \in \mathbb{N}$. I was hoping to be able to find $G$ in terms of the Bernoulli polynomials and numbers themselves, though unfortunately I haven't been able to make any headway.

I was wondering if it is possible to construct a $G$ that satisfies $(2)$ for arbitrary naturals $(l, m, n_{1}, n_{2})$? I don't mind a construction for the special case $n_{1} = n_{2} = n$ either. Or does anyone know any references they could point me to? Any help is appreciated.

  • $\begingroup$ $G(l,m,n_1,n_2)$ as defined by (2) is not a polynomial. $\endgroup$
    – Ira Gessel
    Jan 11, 2022 at 19:17
  • $\begingroup$ @IraGessel I'm not quite following. The RHS of $(2)$ is a rational function if that is what you are alluding to, but I am asking whether we can write the rational function as a polynomial $G$. Sorry if that didn't come through in the question. $\endgroup$ Jan 11, 2022 at 20:13
  • 2
    $\begingroup$ I calculate $$G(1, 2, n, n) = \frac{3(72n^4 + 165n^3 + 70n^2 - 30n + 3)}{7(28n^2 + 33n - 1)}$$ is non-polynomial. $\endgroup$ Jan 11, 2022 at 20:39

1 Answer 1


For the special case $n_1 = n_2 = n$ and $1 \le m \le 9$ it turns out that the double sum yields a polynomial of the form $P(n)(2n+1)(n+1)n^2$ where $P(n)$ is irreducible of degree $2(m-1)$. Obviously the ratio of two irreducible polynomials of different degrees cannot be a polynomial. It may be that for higher values of $m$ the polynomial factors further, but I wouldn't be optimistic about finding polynomial ratios.

My Sage code (online demo):

n = var('n')

def extrapolate_polynomial(ys, x):
    Extrapolates a polynomial using Neville's algorithm.

    ys: the values of the polynomial at 0, 1, 2, ...
    x: the value at which to evaluate the polynomial
    n = len(ys)
    for delta in range(1, n):
        ys = [(x - i) * ys[i+1] - (x - i - delta) * ys[i] for i in range(n - delta)]
    return ys[0] / factorial(n - 1)

def faulhaber(k, n):
    ys = [0] * (k + 2)
    for x in range(1, k + 2):
        ys[x] = ys[x-1] + x**k

    return extrapolate_polynomial(ys, n)

def double_sum(m):
    return expand(sum(binomial(m,i) * faulhaber(2*i, n) * faulhaber(2*m-2*i,n) for i in range(m+1)))

print(factor(double_sum(3) / double_sum(2)))
for m in range(1, 10):
    print(m, "\t", factor(double_sum(m)))
  • $\begingroup$ Thanks Peter, this is exactly the kind of answer I was looking for. $\endgroup$ Jan 12, 2022 at 18:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.