Long ago I was amazed to read Zeilberger’s binomial coefficient masterpiece, A = B, and have ever since been surprised that more people don’t know about it. Sometime in the past year I had this feeling again while reading Kauers’ overview of the cylindrical algebraic decomposition algorithm (CAD).

The sales pitch is this: There exists an algorithm that can determine whether any first-order statement about the reals involving polynomial equations and inequalities is true or false. More generally, there exists an algorithm that can turn any such first order statement into an equivalent one without any quantifiers. The precise statement is more technical, and the potential applications are much broader, but that’s the idea. Apparently we’ve known that this algorithm has existed since the 1960’s, and have had an effective implementation since the 1970’s. Well, it was news to me!

Let me demonstrate a recent problem I was trying to solve that is made almost trivial with CAD.

For which integers $a$ and $b$ does the cubic

have a positive minimum on the region $\{a s + bt \geq 1\} \cap [0, 1]^2$?

This is an annoying problem. It is only just possible to do by hand, because finding the critical points of a cubic amounts to finding the roots of a quadratic, which is easy. But the quadratic is very messy, and it isn’t obvious which critical points are extrema without some asymptotic analysis.

To apply CAD, we can phrase our question like this: For which $a, b$ does the formula

hold? (Note that we don’t need another variable for the minimum itself. The region is compact, so $f(s, t) > 0$ is enough.) If we give this statement to a CAD algorithm, then it will turn it into an unquantified statement, meaning that it only depends on $a$ and $b$. Here’s what Mathematica’s built in CAD does:

f = 1 - 2 a s + a^2 s^2 - b s^2 + s^3 + a b s^3 - b t - 3 s t + a b s t + a s^2 t + b^2 s^2 t + a t^2 + 2 b s t^2 + t^3
CylindricalDecomposition[
 ForAll[{s, t},
  Implies[0 <= s <= 1 && 0 <= t <= 1 && a*s + b*t >= 1, f > 0]],
 {b, a}]
(b <= -(17/4) && (a < (2 - b)/2 - 1/2 Sqrt[-4 + b^2] || a > Root[27 - 4 b^3 + 18 b #1 - b^2 #1^2 + 4 #1^3 &, 3]))
    || (-(17/4) < b <= -2 && (a < (2 - b)/2 - 1/2 Sqrt[-4 + b^2] || a > (2 - b)/2 + 1/2 Sqrt[-4 + b^2]))
    || -2 < b < 1
    || (1 <= b < 2.08... && a > Root[27 - 4 b^3 + 18 b #1 - b^2 #1^2 + 4 #1^3 &, 1])
    || (b == 2.08... && a > (2 - b)/2 + 1/2 Sqrt[-4 + b^2])
    || (b > 2.08... && a > Root[27 - 4 b^3 + 18 b #1 - b^2 #1^2 + 4 #1^3 &, 1])

The output is a disjoint union of conditions which, together, are equivalent to our input. The full output is a bit messy, so here it is again with only the parts relevant to “integers $b \geq 1$”:

    (1 <= b < 2.08... && a > Root[27 - 4 b^3 + 18 b #1 - b^2 #1^2 + 4 #1^3 &, 1])
      || (b > 2.08... && a > Root[27 - 4 b^3 + 18 b #1 - b^2 #1^2 + 4 #1^3 &, 1])

And now, noticing that the inner condition on $a$ is the same in both clauses, and that we only care about integer values, we see that this whole thing is equivalent to just this:

    a > Root[27 - 4 b^3 + 18 b #1 - b^2 #1^2 + 4 #1^3 &, 1]

This says that our original statement about $f$ is equivalent to “$a$ is larger than the smallest root $x$ of the cubic $27 - 4 b^3 + 18 b x - b^2 x^2 + 4 x^3$.” This is certainly simpler, but we can go further.

Using CAD again, we can determine when this cubic has only one real root. Because it’s a cubic, this is the negation of the statement

where $q$ is the cubic in question. Let’s pass that to Mathematica:

q = 27 - 4  b^3 + 18  b  x - b^2  x^2 + 4  x^3
CylindricalDecomposition[
 Exists[{x, y, z},
  x < y < z &&
   q == 0 && (q /. {x -> y}) == 0 && (x /. {x -> z}) == 0], {b}]

b < 1.89...

If $b > 1$ (we only care about integers), then the cubic only has one root. In that case, because its leading coefficient in $x$ is positive, being larger than the real root is equivalent to saying that the cubic is positive. If $b = 1$, then we can just do CAD again but substitute $b = 1$ first, which reveals that the condition is $a > -1$.

So, all together, we have this: If $b = 1$, then $f(s, t) > 0$ on the region for all positive integers $a$; if $b \geq 2$, then $f(s, t) > 0$ on the region if and only if

This is a very precise, quantitative statement, and it only takes a few minutes of computing time (plus some exploring) to discover.

Puiseux series

In one sense the above inequality tells us everything, but in another sense is tells us very little. I don’t know how to get a feel for what solutions of the inequality look like. Should $a$ be bigger than $b$? Smaller? It would be nice to have an approximate quantitative statement which is more digestible.

If you plot our inequality, then you see something like two quadratics stitched together. And if you compute the critical points of the original cubic $f(s, t)$, you start to see things involving square roots. This all points to some simpler quadratic condition.

Indeed, there is such a condition. It comes from a very classical tool called Puiseux series. It turns out that every algebraic equation in two variables¹ has a power series solution where the exponents might have fractions in them. Our equation

is “solved”² by a power series that looks like

So $a$ needs to be smaller than roughly $b^2 / 4$ for our condition to work. You can compute this expansion using the following Maple snippet:

p := 27 + 4 * a^3 + 18 * a * b - a^2 * b^2 - 4 * b^3:
# We want coefficients that shrink as b grows,
# so make it a power series in 1 / b.
q := numer(subs(b=1/b, p)):

with(gfun):
expr := algeqtoseries(q, b, a, 5)[1]:
# Get back to a power series in b.
expr := expand(subs(b=1/b, expr));

In fact, we can compute as many terms as we want:

A nice thing about Puiseux series is that their coefficients satisfy recurrences. If we write our series solution as

then the generating function for $s(n)$ (the series itself!) is algebraic, and every sequence with an algebraic generating function is D-finite. We can compute a recurrence for $s(n)$ by turning our algebraic equation into a differential equation, then reading the recurrence from that. Maple can do all of this for us:

p := 27 + 4 * a^3 + 18 * a * b - a^2 * b^2 - 4 * b^3:
q := numer(subs(b=1/b, p)):

with(gfun):
diffeqn := algeqntodiffeq(subs(a=a(b), q), a(b)):
rec := diffeqntorec(diffeqn, a(b), s(n));

It turns out that the coefficients $s(n)$ eventually satisfy the recurrence

with initial conditions

This implies $s(3) = 0$, so there shouldn’t be any $b^{-3}$ in the expansion, and also

so the next nonzero term should be $4 b^{-4}$. That’s right! In fact, looking back, we see that all of the nonzero terms have coefficients that are $1$ mod $3$ (except for the $b^2$ term), which checks out with this recurrence.

Our recurrence is order three, but it only involves two terms, so there’s a good chance we can get a closed form. The sequence $S(n) = s(3n + 1)$ satisfies the recurrence

Unrolling this recurrence yields

which gives us something to work with. That product turns out to be

This last factorial expression can be written lots of ways. Let’s settle for

Putting everything together and doing some algebra, we get

In summary, we have

Integrality

There is something unusual about our Puiseux series:

Why are the coefficients all integers? It isn’t obvious from the formula

or the recurrence

So there is some explaining to do.

There’s a nice fact about binomial coefficient divisibility: If $\gcd(n + 1, k + 1) = g$, then $k + 1$ divides $g{n \choose k}$. This gives a non-combinatorial proof that the Catalan numbers are integers, for example, because $2n + 1$ and $n + 1$ are relatively prime.

This condition applies directly to the term $2n + 1$ and the binomial coefficient ${3n + 2 \choose n + 2} = {3n + 2 \choose 2n}$. Since $\gcd(2n + 1, 3n + 3) = \gcd(n - 1, 3)$, the condition tells us that $2n + 1$ divides our binomial coefficient if $n$ is not $1$ mod $3$, and it divides $3$ times our binomial coefficient otherwise. This explains one division from the formula, and also why the $3$ in the numerator is necessary.

For the other division, note that $\gcd(2n + 1, 2n + 3) = 1$, so we just need to know when $2n + 3$ divides ${3n + 2 \choose 2n}$. We can write

and apply our condition on the binomial coefficient on the right. It says that $2n + 3$ divides that binomial coefficient if $n$ is not divisible by $3$, and it divides $3$ times the coefficient otherwise. Because $\gcd(n + 2, 2n + 3) = 1$, this translates back to the binomial coefficient on the left.

In summary, it turns out that $2n + 1$ and $2n + 3$ are relatively prime, they both divide

and they both need that $3$ in disjoint cases based on the value of $n \bmod 3$. That explains why $S(n)$ is always an integer.

By the way, I should point out that $S(n) / (-4)^n$ is essentially A000139 in the OEIS. It’s very surprising that our sequence, which arose from studying a difficult calculus problem, is even tangentially related to combinatorics. It would be interesting to see if there is any deeper connection. (There probably isn’t, but who knows!)

The original context for the calculus problem was thinking about some complicated thing involving third order recurrences. It would also be interesting to see if any of what we did here makes sense for higher order recurrences.