# Simplifying a rational function in terms of Bernoulli numbers and polynomials

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.

• $G(l,m,n_1,n_2)$ as defined by (2) is not a polynomial. Jan 11, 2022 at 19:17
• @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. Jan 11, 2022 at 20:13
• 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. Jan 11, 2022 at 20:39

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)))
print("")
for m in range(1, 10):
print(m, "\t", factor(double_sum(m)))

• Thanks Peter, this is exactly the kind of answer I was looking for. Jan 12, 2022 at 18:46