Hi everyone!

Some time ago the following "simple math problem" was shared in a Discord chat:

As a lover of simple math problems, I couldn't just pass by this problem. It turned out much harder than any other genfunc problem that I solved before, as the structure of the answer depends on the parity of $$$n$$$ and $$$m$$$, and it's not very natural to track it through genfuncs. It took me few months, I even called for help from higher powers (i.e. sent a pm to Elegia) but I finally have a solution that I somewhat like.

Finding a closed-form genfunc

First of all, it's really inconvenient to sum up from $$$0$$$ to $$$\lfloor m/2 \rfloor$$$ or $$$\lfloor n/2 \rfloor$$$. What we can do is to introduce variables $$$a=n-2j$$$ and $$$b=m-2i$$$, after which the expression under the sum changes to

As we want to sum it up over all $$$i,j,a,b$$$ such that $$$a + 2j = n$$$ and $$$b + 2i = m$$$, we can keep track of these equations with two formal variables $$$x$$$ and $$$y$$$. Then, the sum that we're interested in expresses as

This is the product of two generating functions. The first one is standard and collapses with $$$s=a+b$$$ substitution:

The second generating function is trickier, and collapses into

Unfortunately, I don't have a simple explanation to it, as it is obtained via reduction to a well-known bivariate generating function for Legendre polynomials (see this Mathematics Stack Exchange question for details). So, the problem reduces to finding

Adding series multisection

As I mentioned earlier, since the second function in the product only depends on $$$x^2$$$ and $$$y^2$$$, it would make sense to do a series multisection (aka roots of unity filter) to the first function to go back to $$$x$$$ and $$$y$$$ instead of $$$x^2$$$ and $$$y^2$$$. The multisection is generally done as

where the first summand only has even powers of $$$x$$$, and the second summand only has odd powers of $$$x$$$. In our case, we need to apply it first to the variable $$$x$$$, and the to $$$y$$$. However, we can do a little shortcut, as we know that the common denominator of the resulting expression should be

meaning that the numerator must be

which expands to

Solving for parities

The $$$4$$$ distinct summands in the numerator correspond to all possible combinations of parities of powers of $$$x$$$ and $$$y$$$. But what is really striking here is that the denominator expands to

meaning that it is exactly the expression under the square root of the second function multiplier. Therefore, the full problem reduces to finding (and computing an appropriate linear combination of) $$$[x^n y^m]$$$ of the function

This function is very similar to

about which we already know that $$$[x^n y^m] G(x, y) = \binom{n+m}{n}^2$$$. Indeed, we can notice that

therefore

from which it follows that

This allows to solve the problem in $$$O(1)$$$ per $$$(n, m)$$$ pair with $$$O(n+m)$$$ pre-computation:

def F(n, m):
    if n < 0 or m < 0:
        return 0
    return ((n+1)*sympy.binomial(n+m+1, n+1)**2 + (m+1)*sympy.binomial(n+m+1, m+1)**2) // 2

def T(n, m):
    if n % 2 and m % 2:
        return 2*F(n//2, m//2)
    if n % 2 == 0 and m % 2:
        n, m = m, n
    if n % 2 and m % 2 == 0:
        return F(n//2, m//2) - F(n//2-1, m//2) + F(n//2, m//2-1)
    return F(n//2, m//2) - F(n//2-1, m//2) - F(n//2, m//2-1)

Alternative approaches

The answer to the whole problem can be further simplified as

Thanks to Endagorion for telling me about this form! One can prove it by inspecting how $$$f(n, m)$$$ relates to $$$f(n-1,m)$$$ and $$$f(n,m-1)$$$, but apparently there is also a bijective proof. You may also refer to this publication for further details and alternative proofs.

Follow-up questions to interested audience

1. Is there a sufficiently simple way to solve this problem without just "observing" some weird facts?
2. Is there a more direct way to find the closed form on the second function? (see this question).

UPD: I found out how to compute the genfunc for $$$\binom{i+j}{i}^2$$$ directly, without using Legendre polynomials.

Consider a bivariate polynomial $$$Q_n(x, y)$$$ defined as

We need to sum it up over all $$$n$$$, and we know from Legendre polynomials that $$$Q_n(x) = (y-x)^n P_n(\frac{y+x}{y-x})$$$.

But let's analyze $$$Q_n(x, y)$$$ on its own:

This expands into

Note that

hence we want to compute the sum

Let's sum up over $$$n$$$ for each fixed $$$k$$$:

Thus, we want to compute

On the other hand we know that

thus the sum above compacts into

which then expands into the same expression in the denominator.