Polynomial Integer Sums as Functional Equations

The sum of natural numbers

Can be described as a recursive function

Or

We seek $f(x)$. We will start with a guess-and-check method starting with $f(x) = x^2$.

This guess for $f(x)$ does not satisfy the constraint and must be updated. A more clear method for doing this will be explained later, for now, halve current guess:

Now, the guess for $f(x)$ is a lot closer to satisfying the equation. Update guess again

A Bit more formally

This method, as previously shown, works by incrementally adding and removing polynomials to $f(x)$ until

Where $h(x)$ is some polynomial function. A table has been made to better understand the effect which adding a power of $x$ to the guess for $f(x)$ has on the overall equation.

The above table indicates that, if $x^2$ is added to the guess for $f(x)$, then $f(x)-f(x-1)$ will have $2x-1$ added to it as a result.

Using this table, the method will again be applied to the sum of squares, where $h(x) = x^2$. It is not necessary to use the table, but the algebra becomes pretty murky without it.

The starting guess for $f(x)$ is informed by the highest power of $x$ in $h(x)$. Because $h(x)=x^2$, the first guess will be $f(x) = \frac{1}{3}x^3$. The specific $\frac{1}{3}x^3$ is chosen because it will make $f(x) - f(x-1)$ contain $x^2$ which is of course desirable.

Now that the expression contains $x^2$, the unnecessary components $-x + \frac{1}{3}$ must be removed. Remove $-x$ by adding $\frac{1}{2}x^2$ to the guess for $f(x)$

Next, $-\frac{1}{6}$ is removed by adding $\frac{1}{6}x$ to the guess for $f(x)$

Thus, the sum of squares formula is

Solver Algorithm Implementation

A sample algorithm implementation in Python using NumPy is provided:

import fractions as F
from numpy.polynomial import polynomial as P

# Possibly unexpected behavior:
assert P.polytrim([0]) == [0]

def polyset(p, n, v):
    while len(p) <= n:
        p.append(0)
    p[n] = v

# Given f(x) - f(x-1) = g(x) Where `g` is given value
given = [F.Fraction(0), F.Fraction(1)] # x^2

# Populate table
table = []
base = [F.Fraction(-1), F.Fraction(1)]
tmp = base.copy()
for i in range(len(given) + 1):
    x = tmp[:-1]
    x = P.polymul(x, [-1])
    leading_coef = x[-1]
    table.append((leading_coef, x))
    tmp = P.polymul(tmp, base)
tmp = None
base = None

# move given to left side of equation: f(x)-f(x-1) - g(x) = 0
working = P.polymul(given, -1)
out = []

# Work until a satisfying function $f$ is found
# such that the equation is zero on both sides
while len(working) > 1 or working[-1] != 0:
    high_pwr = len(working) - 1
    hp_c = working[high_pwr]
    assert hp_c != 0
    lc, tp = table[high_pwr]
    factor = hp_c / lc * -1
    polyset(out, high_pwr + 1, factor)
    tp = P.polymul(tp, [factor])
    working = P.polyadd(working, tp)
    working = P.polytrim(working)

print("The sum of: ", P.Polynomial(given))
print("Has the closed form: ", P.Polynomial(out))

Additionally, I’ve ported it to JavaScript so you can use it interactively. All polynomials must be ine expanded form.