[Faulhaber's formula](https://en.wikipedia.org/wiki/Faulhaber%27s_formula) expresses a sum over some finite number of naturals to the $m^{th}$ power in terms of the [Bernoulli numbers](https://en.wikipedia.org/wiki/Faulhaber%27s_formula#Summae_Potestatum) $B\_{j}$ (using the $B\_{1} = 1/2$ convention) or [polynomials](https://en.wikipedia.org/wiki/Faulhaber%27s_formula#Alternate_expressions) $\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.



Simplifying a rational function in terms of Bernoulli numbers and polynomials
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](https://sagecell.sagemath.org/?z=eJx9UU2PgjAUPGvif5hwsSC6wt5M9rjX_QMuMVWKNLbFlOLCv99X8APdZCFQ-t7MdN5g8IELt2xu5uFs6u9cFBCts_xcKe7EjpbOVFpyxbo6RhtuZtNJEAT0_nzAanA8kGhqaY74EheplJhTUx0rK12pV_6ISVdv4EpBJ6uGqFXR70Z87rCOkcRIY6xWRJq0I4Zv_5TyUMJVEL5EDl4k7iYNTaiEIfM04KSoLHKhHIc0sNwcBaNTTD8U2SLslrVYQoaI0NVbuUgy2g41enrurZfBy8mHlLkjMtKzwjXWeOQ6wxsKfnAUAuXoYcko7oI3quR7Ydnp5mWwQrwI7IQF0pv59sn40Lqa37YZkfy6JNMLtFF06uO-Gvnnr5qRm7xq9krs6kYz7aVB113hzE3OfGsvr3Qd91k9Zkgj6QVfanrp61R-iUwvkjC8nn620jg25MRGNt5Dim-0T3vGgA4C-vSa-imYZH2zPsB0jODbBTH-qmtS-wUybuVN&lang=sage&interacts=eJyLjgUAARUAuQ==)):



```
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)))

```

