Solving a Three-Variable Recursion via Generating Functions

Author: vchernoy

404.html

@@ -568,6 +568,8 @@ <h1>Page not found</h1>
   <h2>Latest</h2>
   <ul>
 
+      <li><a href="/post/three-var-recursive/">Solving a Three-Variable Recursion via Generating Functions</a></li>
+    
       <li><a href="/post/">Posts</a></li>
 
       <li><a href="/post/perfect-distribution/">Perfect Distribution: GCD in Disguise</a></li>

content/post/three-var-recursive.md

@@ -0,0 +1,155 @@
++++
+date = "2018-04-10T03:14:00Z"
+math = true
+highlight = true
+tags = ["math", "python", "generating function", "recursion", "combinatorics", "dynamic programming", "multivariate function"]
+title = "Solving a Three-Variable Recursion via Generating Functions"
+draft = false
+
+[header]
+image = ""
+caption = ""
++++
+
+This post extends the generating-function technique from the [two-variable recursion][two-var-recursive-func] to a three-variable case. I originally wrote this as an answer to a [Math Stack Exchange question][math-se]; here it is adapted for the blog with clearer exposition and code.
+
+## The Problem
+
+We want to solve the recurrence
+
+$$a(m,n,k) = 2a(m-1,n-1,k-1) + a(m-1,n-1,k) + a(m-1,n,k-1) + a(m,n-1,k-1)$$
+
+where $m$, $n$, $k$ are nonnegative integers, with boundary conditions:
+
+- $a(0,0,0) = a(1,0,0) = a(0,1,0) = a(0,0,1) = 1$
+- $a(m,0,0) = a(0,m,0) = a(0,0,m) = 0$ for $m > 1$
+- $a(m,n,k)$ is symmetric in $m,n,k$
+
+A subtlety: $a(0,1,1)$ is not defined by the recurrence alone, since it would require values like $a(-1,0,0)$. We take $a(m,n,k) = 0$ whenever any argument is negative.
+
+## The Generating Function
+
+Define
+
+$$\Phi(x,y,z) = \sum_{m,n,k \geq 0} a(m,n,k)\, x^m y^n z^k$$
+
+Substituting the recurrence and collecting terms, we get
+
+$$
+\begin{aligned}
+\Phi(x,y,z)
+&= \sum_{m,n,k} a(m,n,k)\, x^m y^n z^k \\
+&= 2\Phi\, xyz + \Phi\, xy + \Phi\, xz + \Phi\, yz + 1 + x + y + z
+\end{aligned}
+$$
+
+where the boundary terms $1 + x + y + z$ come from $a(0,0,0)$, $a(1,0,0)$, $a(0,1,0)$, $a(0,0,1)$. Solving for $\Phi$:
+
+$$\Phi(x,y,z) = \frac{1 + x + y + z}{1 - 2xyz - xy - xz - yz}$$
+
+## From Generating Function to Closed Form
+
+Using $\frac{1}{1-\rho} = \sum_{i \geq 0} \rho^i$ and the multinomial expansion
+
+$$(x_1+x_2+x_3+x_4)^N = \sum_{k_1+k_2+k_3+k_4=N} \binom{N}{k_1,k_2,k_3,k_4} x_1^{k_1} x_2^{k_2} x_3^{k_3} x_4^{k_4}$$
+
+with $\binom{N}{k_1,k_2,k_3,k_4} = \frac{N!}{k_1!\,k_2!\,k_3!\,k_4!}$, we expand the denominator. Let $\rho = 2xyz + xy + xz + yz$. Then
+
+$$
+\Phi = (1+x+y+z) \sum_{N \geq 0} \rho^N
+$$
+
+Writing $\rho^N$ with terms $(2xyz)^{k_1}(xy)^{k_2}(xz)^{k_3}(yz)^{k_4}$ where $k_1+k_2+k_3+k_4 = N$, and extracting the coefficient of $x^m y^n z^k$, yields the closed form:
+
+$$
+\begin{aligned}
+a(m,n,k)
+&= \sum_{\max(m,n,k) \leq N \leq \frac{m+n+k}{2}}
+   \binom{N}{m+n+k-2N,\, N-m,\, N-n,\, N-k} 2^{m+n+k-2N} \\
+&\quad + \text{three similar sums from the } 1,\, x,\, y,\, z \text{ terms}
+\end{aligned}
+$$
+
+The full expression has four sums corresponding to the four terms in the numerator $1+x+y+z$. The exact formulas are:
+
+$$
+\begin{aligned}
+a(m,n,k)
+&= \sum_{N=\max(m,n,k)}^{\lfloor (m+n+k)/2 \rfloor}
+   \binom{N}{m+n+k-2N,\, N-m,\, N-n,\, N-k} 2^{m+n+k-2N} \\
+&\quad + \sum_{N=\max(m-1,n,k)}^{\lfloor (m+n+k-1)/2 \rfloor}
+   \binom{N}{m+n+k-2N-1,\, N-m+1,\, N-n,\, N-k} 2^{m+n+k-2N-1} \\
+&\quad + \sum_{N=\max(m,n-1,k)}^{\lfloor (m+n+k-1)/2 \rfloor}
+   \binom{N}{m+n+k-2N-1,\, N-m,\, N-n+1,\, N-k} 2^{m+n+k-2N-1} \\
+&\quad + \sum_{N=\max(m,n,k-1)}^{\lfloor (m+n+k-1)/2 \rfloor}
+   \binom{N}{m+n+k-2N-1,\, N-m,\, N-n,\, N-k+1} 2^{m+n+k-2N-1}
+\end{aligned}
+$$
+
+There may be room to simplify this further; the symmetry in $m,n,k$ could help.
+
+## Complexity
+
+- **Recursion with memoization (DP):** $\Theta(m \cdot n \cdot k)$ time and space.
+- **Closed form:** Precompute factorials, then loop over $N$; time and space $\Theta(m+n+k)$, ignoring the cost of arithmetic on large integers.
+
+## Implementation
+
+Both the recursive and closed-form versions in Python:
+
+```python
+import functools
+import math
+
+@functools.lru_cache(maxsize=None)
+def a_rec(m: int, n: int, k: int) -> int:
+    if min(m, n, k) < 0:
+        return 0
+    if m + n + k == 0 or m + n + k == 1:
+        return 1
+    if m + n == 0 or m + k == 0 or n + k == 0:
+        return 0
+    return (
+        2 * a_rec(m - 1, n - 1, k - 1)
+        + a_rec(m - 1, n - 1, k)
+        + a_rec(m - 1, n, k - 1)
+        + a_rec(m, n - 1, k - 1)
+    )
+
+@functools.lru_cache(maxsize=None)
+def _binom4(N: int, a: int, b: int, c: int) -> int:
+    r = N - a - b - c
+    vals = sorted([a, b, c, r])
+    assert vals[0] >= 0
+    return math.factorial(N) // (
+        math.factorial(vals[0]) * math.factorial(vals[1])
+        * math.factorial(vals[2]) * math.factorial(vals[3])
+    )
+
+def a_closed(m: int, n: int, k: int) -> int:
+    if min(m, n, k) < 0:
+        return 0
+    s = 0
+    for N in range(max(m, n, k), (m + n + k) // 2 + 1):
+        s += _binom4(N, N - m, N - n, N - k) * 2 ** (m + n + k - 2 * N)
+    for N in range(max(m - 1, n, k), (m + n + k - 1) // 2 + 1):
+        s += _binom4(N, N - m + 1, N - n, N - k) * 2 ** (m + n + k - 2 * N - 1)
+    for N in range(max(m, n - 1, k), (m + n + k - 1) // 2 + 1):
+        s += _binom4(N, N - m, N - n + 1, N - k) * 2 ** (m + n + k - 2 * N - 1)
+    for N in range(max(m, n, k - 1), (m + n + k - 1) // 2 + 1):
+        s += _binom4(N, N - m, N - n, N - k + 1) * 2 ** (m + n + k - 2 * N - 1)
+    return s
+
+# Sanity check
+r, r1 = a_rec(100, 200, 210), a_closed(100, 200, 210)
+print(f"Recursive: {r}, Closed: {r1}, Match: {r == r1}")
+```
+
+We did not exploit the symmetry $a(m,n,k) = a(\sigma(m,n,k))$ for permutations $\sigma$; it could speed up computation but does not obviously simplify the closed expression.
+
+---
+
+*Originally answered on [Math Stack Exchange][math-se].*
+
+[two-var-recursive-func]: /post/two-var-recursive-func/
+[math-se]: https://math.stackexchange.com/questions/1093271/how-to-solve-this-multivariable-recursion/2730331#2730331

index.html

@@ -341,7 +341,7 @@
 <meta property="og:description" content="Software Engineer" /><meta property="og:image" content="/media/icon_hu_43cf117bf1a42c34.png" /><meta property="og:locale" content="en-us" />
 
 
-    <meta property="og:updated_time" content="2026-02-08T19:47:43&#43;00:00" />
+    <meta property="og:updated_time" content="2018-04-10T03:14:00&#43;00:00" />
 
 
 
@@ -975,14 +975,23 @@ <h1 class="mb-0">Recent Posts</h1>
   <div class="media-body">
 
     <div class="section-subheading article-title mb-0 mt-0">
-      <a href="/post/perfect-distribution/" >Perfect Distribution: GCD in Disguise</a>
+      <a href="/post/three-var-recursive/" >Solving a Three-Variable Recursion via Generating Functions</a>
     </div>
 
 
-    <a href="/post/perfect-distribution/"  class="summary-link">
+    <a href="/post/three-var-recursive/"  class="summary-link">
       <div class="article-style">
-        <h2 id="perfect-distribution-gcd-in-disguise">Perfect Distribution: GCD in Disguise</h2>
-<p>We discuss an algorithm that distributes $a$ ones among $n$ positions so that the gaps between consecutive ones differ by at most one—a <strong>perfect distribution</strong>. I developed it while working on profiling, stress, and negative testing of a system that needed exactly this kind of uniform spread. I am not aware of prior art; if you know of related work, I would be interested to hear.</p>
+        <p>This post extends the generating-function technique from the <a href="/post/two-var-recursive-func/">two-variable recursion</a> to a three-variable case. I originally wrote this as an answer to a <a href="https://math.stackexchange.com/questions/1093271/how-to-solve-this-multivariable-recursion/2730331#2730331" target="_blank" rel="noopener">Math Stack Exchange question</a>; here it is adapted for the blog with clearer exposition and code.</p>
+<h2 id="the-problem">The Problem</h2>
+<p>We want to solve the recurrence</p>
+<p>$$a(m,n,k) = 2a(m-1,n-1,k-1) + a(m-1,n-1,k) + a(m-1,n,k-1) + a(m,n-1,k-1)$$</p>
+<p>where $m$, $n$, $k$ are nonnegative integers, with boundary conditions:</p>
+<ul>
+<li>$a(0,0,0) = a(1,0,0) = a(0,1,0) = a(0,0,1) = 1$</li>
+<li>$a(m,0,0) = a(0,m,0) = a(0,0,m) = 0$ for $m &gt; 1$</li>
+<li>$a(m,n,k)$ is symmetric in $m,n,k$</li>
+</ul>
+<p>A subtlety: $a(0,1,1)$ is not defined by the recurrence alone, since it would require values like $a(-1,0,0)$. We take $a(m,n,k) = 0$ whenever any argument is negative.</p>
       </div>
     </a>
 
@@ -1006,7 +1015,7 @@ <h2 id="perfect-distribution-gcd-in-disguise">Perfect Distribution: GCD in Disgu
 
 
 
-    Aug 1, 2017
+    Apr 10, 2018
   </span>
 
 
@@ -1015,7 +1024,7 @@ <h2 id="perfect-distribution-gcd-in-disguise">Perfect Distribution: GCD in Disgu
 
   <span class="middot-divider"></span>
   <span class="article-reading-time">
-    7 min read
+    4 min read
   </span>
 
 
@@ -1063,15 +1072,14 @@ <h2 id="perfect-distribution-gcd-in-disguise">Perfect Distribution: GCD in Disgu
   <div class="media-body">
 
     <div class="section-subheading article-title mb-0 mt-0">
-      <a href="/post/binomial-modulo-prime/" >Binomial Coefficients Modulo a Prime: Fermat&#39;s Theorem and the Non-Adjacent Selection Problem</a>
+      <a href="/post/perfect-distribution/" >Perfect Distribution: GCD in Disguise</a>
     </div>
 
 
-    <a href="/post/binomial-modulo-prime/"  class="summary-link">
+    <a href="/post/perfect-distribution/"  class="summary-link">
       <div class="article-style">
-        <p>In the <a href="/post/efficient-implementation-non-adjacent-selection/">previous post</a>, we implemented the closed form $F_{n,m} = \binom{n-m+1}{m}$ using Python&rsquo;s <code>math.factorial</code>, and with <code>scipy</code> and <code>sympy</code>. Here we cover the common competitive-programming case: computing the answer <strong>modulo a large prime</strong> $M$ (e.g. $M = 10^9+7$).</p>
-<h2 id="why-modulo">Why modulo?</h2>
-<p>In counting problems, the result can be huge even for moderate input. Often the problem asks for the answer modulo a big prime so that it fits in a standard integer type. We could compute the full number and then take the remainder, but that forces expensive long-integer arithmetic. Computing <strong>everything</strong> modulo $M$ from the start is much faster.</p>
+        <h2 id="perfect-distribution-gcd-in-disguise">Perfect Distribution: GCD in Disguise</h2>
+<p>We discuss an algorithm that distributes $a$ ones among $n$ positions so that the gaps between consecutive ones differ by at most one—a <strong>perfect distribution</strong>. I developed it while working on profiling, stress, and negative testing of a system that needed exactly this kind of uniform spread. I am not aware of prior art; if you know of related work, I would be interested to hear.</p>
       </div>
     </a>
 
@@ -1095,7 +1103,7 @@ <h2 id="why-modulo">Why modulo?</h2>
 
 
 
-    Jul 8, 2017
+    Aug 1, 2017
   </span>
 
 
@@ -1104,7 +1112,7 @@ <h2 id="why-modulo">Why modulo?</h2>
 
   <span class="middot-divider"></span>
   <span class="article-reading-time">
-    2 min read
+    7 min read
   </span>
 
 
@@ -1152,17 +1160,15 @@ <h2 id="why-modulo">Why modulo?</h2>
   <div class="media-body">
 
     <div class="section-subheading article-title mb-0 mt-0">
-      <a href="/post/efficient-implementation-non-adjacent-selection/" >Efficient Implementation of the Non-Adjacent Selection Formula</a>
+      <a href="/post/binomial-modulo-prime/" >Binomial Coefficients Modulo a Prime: Fermat&#39;s Theorem and the Non-Adjacent Selection Problem</a>
     </div>
 
 
-    <a href="/post/efficient-implementation-non-adjacent-selection/"  class="summary-link">
+    <a href="/post/binomial-modulo-prime/"  class="summary-link">
       <div class="article-style">
-        <p>In the <a href="/post/two-var-recursive-func/">previous post</a>, we derived the closed form for the non-adjacent selection problem:</p>
-<p>$$ F_{n, m} = {n - m + 1 \choose m} $$</p>
-<p>Now we discuss how to implement this efficiently in Python—from a simple factorial-based solution to library implementations. For the common case of computing the answer <strong>modulo a large prime</strong> (e.g. in competitive programming), see the <a href="/post/binomial-modulo-prime/">next post</a>.</p>
-<h2 id="fast-solutions-based-on-binomials">Fast Solutions Based on Binomials</h2>
-<p>We can reflect the closed form in very trivial Python code:</p>
+        <p>In the <a href="/post/efficient-implementation-non-adjacent-selection/">previous post</a>, we implemented the closed form $F_{n,m} = \binom{n-m+1}{m}$ using Python&rsquo;s <code>math.factorial</code>, and with <code>scipy</code> and <code>sympy</code>. Here we cover the common competitive-programming case: computing the answer <strong>modulo a large prime</strong> $M$ (e.g. $M = 10^9+7$).</p>
+<h2 id="why-modulo">Why modulo?</h2>
+<p>In counting problems, the result can be huge even for moderate input. Often the problem asks for the answer modulo a big prime so that it fits in a standard integer type. We could compute the full number and then take the remainder, but that forces expensive long-integer arithmetic. Computing <strong>everything</strong> modulo $M$ from the start is much faster.</p>
       </div>
     </a>
 
@@ -1186,7 +1192,7 @@ <h2 id="fast-solutions-based-on-binomials">Fast Solutions Based on Binomials</h2
 
 
 
-    Jul 7, 2017
+    Jul 8, 2017
   </span>
 
 
@@ -1195,7 +1201,7 @@ <h2 id="fast-solutions-based-on-binomials">Fast Solutions Based on Binomials</h2
 
   <span class="middot-divider"></span>
   <span class="article-reading-time">
-    4 min read
+    2 min read
   </span>
 
 
@@ -1243,19 +1249,17 @@ <h2 id="fast-solutions-based-on-binomials">Fast Solutions Based on Binomials</h2
   <div class="media-body">
 
     <div class="section-subheading article-title mb-0 mt-0">
-      <a href="/post/two-var-recursive-func/" >Cracking Multivariate Recursive Equations Using Generating Functions</a>
+      <a href="/post/efficient-implementation-non-adjacent-selection/" >Efficient Implementation of the Non-Adjacent Selection Formula</a>
     </div>
 
 
-    <a href="/post/two-var-recursive-func/"  class="summary-link">
+    <a href="/post/efficient-implementation-non-adjacent-selection/"  class="summary-link">
       <div class="article-style">
-        <p>In this post, we return back to the combinatorial problem discussed in <a href="/post/intro-to-dp/">Introduction to Dynamic Programming and Memoization</a> post.
-We will show that generating functions may work great not only for single variable case (see <a href="/post/gen-func-art/">The Art of Generating Functions</a>),
-but also could be very useful for hacking two-variable relations (and of course, in general for multivariate case too).</p>
-<p>For making the post self-contained, we repeat the problem definition here.</p>
-<h2 id="the-problem">The Problem</h2>
-<blockquote>
-<p>Compute the number of ways to choose $m$ elements from $n$ elements such that selected elements in one combination are not adjacent.</p>
+        <p>In the <a href="/post/two-var-recursive-func/">previous post</a>, we derived the closed form for the non-adjacent selection problem:</p>
+<p>$$ F_{n, m} = {n - m + 1 \choose m} $$</p>
+<p>Now we discuss how to implement this efficiently in Python—from a simple factorial-based solution to library implementations. For the common case of computing the answer <strong>modulo a large prime</strong> (e.g. in competitive programming), see the <a href="/post/binomial-modulo-prime/">next post</a>.</p>
+<h2 id="fast-solutions-based-on-binomials">Fast Solutions Based on Binomials</h2>
+<p>We can reflect the closed form in very trivial Python code:</p>
       </div>
     </a>
 
@@ -1279,7 +1283,7 @@ <h2 id="the-problem">The Problem</h2>
 
 
 
-    Jul 6, 2017
+    Jul 7, 2017
   </span>
 
 
@@ -1288,7 +1292,7 @@ <h2 id="the-problem">The Problem</h2>
 
   <span class="middot-divider"></span>
   <span class="article-reading-time">
-    3 min read
+    4 min read
   </span>
 
 
@@ -1336,20 +1340,19 @@ <h2 id="the-problem">The Problem</h2>
   <div class="media-body">
 
     <div class="section-subheading article-title mb-0 mt-0">
-      <a href="/post/gen-func-art/" >The Art of Generating Functions</a>
+      <a href="/post/two-var-recursive-func/" >Cracking Multivariate Recursive Equations Using Generating Functions</a>
     </div>
 
 
-    <a href="/post/gen-func-art/"  class="summary-link">
+    <a href="/post/two-var-recursive-func/"  class="summary-link">
       <div class="article-style">
-        <p>The notion of generating functions and its application to solving recursive equations are very well-known.
-For reader who did not have a chance to get familiar with the topics,
-I recommend to take a look at very good book:
-<a href="https://en.wikipedia.org/wiki/Concrete_Mathematics" target="_blank" rel="noopener">Concrete Mathematics: A Foundation for Computer Science, by Ronald L. Graham, Donald E. Knuth, Oren Patashnik</a>.</p>
-<p>Generating functions are usually applied to single variable recursive equations.
-But actually, the technique may be extended to multivariate recursive equations, or even to a system of recursive equations.
-Readers who are familiar with one-variable case, may jump directly to the next post:
-<a href="/post/two-var-recursive-func/">Cracking Multivariate Recursive Equations Using Generating Functions</a>.</p>
+        <p>In this post, we return back to the combinatorial problem discussed in <a href="/post/intro-to-dp/">Introduction to Dynamic Programming and Memoization</a> post.
+We will show that generating functions may work great not only for single variable case (see <a href="/post/gen-func-art/">The Art of Generating Functions</a>),
+but also could be very useful for hacking two-variable relations (and of course, in general for multivariate case too).</p>
+<p>For making the post self-contained, we repeat the problem definition here.</p>
+<h2 id="the-problem">The Problem</h2>
+<blockquote>
+<p>Compute the number of ways to choose $m$ elements from $n$ elements such that selected elements in one combination are not adjacent.</p>
       </div>
     </a>
 
@@ -1373,7 +1376,7 @@ <h2 id="the-problem">The Problem</h2>
 
 
 
-    Jul 5, 2017
+    Jul 6, 2017
   </span>
 
 
@@ -1382,7 +1385,7 @@ <h2 id="the-problem">The Problem</h2>
 
   <span class="middot-divider"></span>
   <span class="article-reading-time">
-    4 min read
+    3 min read
   </span>