1 Introduction↩︎

The longest increasing subsequence problem has been a topic of great interest studied from several complementary viewpoints. On static inputs, the classical exact computation problem admits an \(O(n\log n)\) comparison-based algorithm [1]. For permutations of \([n]\), faster word-RAM algorithms are known. In particular, the standard predecessor-based implementation gives \(O(n\log\log n)\) time, and Crochemore and Porat give an \(O(n\log\log \ell)\) bound in terms of the LIS length \(\ell\) [2]. Later work developed streaming and sublinear-time approximation algorithms for \(\operatorname{LIS}\) and distance to monotonicity [3]–[5]. On dynamic inputs, variants of the problem have been studied in sliding-window and fully dynamic models [6]–[9]. On random permutations, the Robinson–Schensted correspondence connects \(\operatorname{LIS}\) to Plancherel-distributed Young diagrams. This viewpoint is rooted in the Ulam–Hammersley problem and Hammersley’s interacting particle process, and leads to the limit-shape results of Logan–Shepp and Vershik–Kerov and the fluctuation theorem of Baik–Deift–Johansson. Subsequent work connects these Plancherel asymptotics with random matrix theory and determinantal processes. See also the surveys of Aldous–Diaconis and Stanley [10]–[19]. For a comprehensive introduction to the mathematics of the longest increasing subsequence, see [20].

Despite this extensive literature, a basic algorithmic question has apparently gone unaddressed: how does one sample a permutation of length \(n\) uniformly at random, subject to the constraint that its \(\operatorname{LIS}\) has a prescribed length \(k\)? For a permutation \(\pi\in S_n\), let \(\operatorname{LIS}(\pi)\) denote the length of its longest increasing subsequence. Given input parameters \((n, k)\), we want to sample exactly from the uniform distribution on \[\Omega_{n,k}:=\{\pi\in S_n:\operatorname{LIS}(\pi)=k\}.\]

Rejection sampling from \(S_n\) does not give a polynomial-time algorithm as the probability that a uniform random permutation has \(\operatorname{LIS}\) equal to a specific \(k\) can be as small as \(1/n!\). One could in principle attempt to sample from \(\Omega_{n,k}\) via Markov chain Monte Carlo, but no Markov chain on \(\Omega_{n,k}\) with a provable polynomial mixing time is known, and our samplers are exact rather than approximate.

The closest works we are aware of concern exact sampling problems in the same Robinson–Schensted and Young-tableau setting, but not uniform sampling of permutations with a prescribed LIS. Thomas and Yong interpret Hecke insertion as an exact polynomial-time sampler for a Plancherel-type measure, and Betea, Boutillier, Bouttier, Chapuy, Corteel, and Vuletić give perfect sampling algorithms for Schur processes [21], [22]. To the best of our knowledge, the exact uniform sampling problem for \(\Omega_{n,k}\) has not been studied.

2 Our Results↩︎

We first give a fast direct sampler for the regime \(k\in\Theta(n)\). The sampler uses rejection sampling on an expanded proposal space of pairs \((\pi,I)\), where \(I\) is a distinguished \(k\)-term increasing subsequence of \(\pi\). A permutation with several LISs may have several representatives in this expanded space, so we use a canonical rule that accepts only the pair in which \(I\) is the leftmost LIS. Thus each permutation with \(\operatorname{LIS}=k\) contributes exactly one accepted proposal. We then prove that, when \(k\in\Theta(n)\), the acceptance probability is bounded below by a positive constant, giving the claimed expected running time.

For general \(k\), the acceptance probability of the rejection sampler can be too small and so we must take an entirely different approach. The Robinson–Schensted correspondence gives a bijection between permutations of \([n]\) and pairs of Young tableaux, often written as a pair \((P, Q)\). It also identifies \(\operatorname{LIS}(\pi)\) with the first row length of the associated partition shape \(\lambda\). Figure 1 shows an example of the correspondence. The number of permutations of shape \(\lambda\) equals \((f^\lambda)^2\), where \(f^\lambda=|\operatorname{SYT}(\lambda)|\) is the number of standard Young tableaux of shape \(\lambda\) [23]. Thus uniform sampling from \(\Omega_{n,k}\) is equivalent to sampling a partition \(\lambda\vdash n\) from \[\mathbf{P}_{n,k}(\lambda)\propto (f^\lambda)^2\,\mathbf{1}_{\{\lambda_1=k\}},\] and then sampling a uniform pair of tableaux of shape \(\lambda\). Given a pair of tableaux, we can then recover the sampled permutation by the inverse Robinson–Schensted correspondence, implemented by reverse insertion. Unfortunately a naive dynamic program over Young diagrams may have superpolynomially many states, since there are \(\exp(\Theta(\sqrt n))\) partitions of \(n\). The Robinson–Schensted reduction and the step of randomly sampling tableaux of a fixed shape are standard, leaving open the question of how to sample the partition \(\lambda\) efficiently.

Our main contribution is an exact sampler for the conditioned shape distribution. The sampler can be viewed as an exact dynamic program over shapes, but with the large completion sums compressed into determinant coefficients. First, we rewrite the shape problem as a weighted fixed-sum problem over strict partitions, that is partitions with distinct parts. Then we expose the sequence one coordinate at a time, choosing each candidate with probability proportional to the total remaining weight of all completions below that candidate. These completion weights are not computed by enumerating completions. Instead, the Cauchy–Binet formula turns each completion sum, which contains a squared Vandermonde factor, into a coefficient of the determinant of a matrix whose entries are polynomials in an auxiliary variable. This determinant representation is what makes the general sampler polynomial-time. The simplest implementation evaluates this determinant oracle separately for each candidate at a stage, which already gives a polynomial-time exact sampler. The direct implementation result is stated in Theorem 2 and proved in Section 5.

The direct implementation above is self-contained apart from standard determinant computation and already gives a polynomial-time exact sampler for all \(1\le k\le n\). We also give a faster implementation of the same sequential sampler. This speedup does not change the distributional argument. It only changes how the completion scores at one stage are evaluated. The faster implementation uses a superfast Hankel determinant primitive as a structured-linear-algebra black box.

For the word-RAM analysis, the integer completion scores are computed modulo word-size primes and reconstructed by Chinese remaindering. This keeps all field operations at unit word cost, and the number of primes contributes the \(\tilde{O}(kn)\) overhead accounted for in the running-time bounds below. Theorem 3 gives the final result.

The key point is that the completion sums driving the sequential sampler, which can be exponentially large in \(n\), are never enumerated. Each is evaluated exactly as a single determinant coefficient via the Cauchy–Binet formula, and at each stage all candidate determinants are computed together by exploiting their shared Hankel structure.

3 A direct sampler when \(k \in \Theta(n)\)↩︎

Our direct sampler works directly with permutations and uses rejection sampling. The proposal distribution is supported on an expanded space of pairs \((\pi,I)\), where \(I\) is a \(k\)-term increasing subsequence of the permutation \(\pi\). The acceptance rule is designed so that, for each permutation with \(\operatorname{LIS}=k\), exactly one associated pair is accepted: the pair for which \(I\) is the leftmost LIS. The role of the leftmost LIS is simply as a canonical choice of LIS that can be computed rapidly.

Throughout this section let \(m:=n-k\). The construction and its exactness are valid for every \(1\le k\le n\). The constant acceptance-probability bound, and hence the final expected-time bound, will be proved only for sequences \(k=k(n)\) with \(k\in\Theta(n)\).

Define the expanded proposal space \[\mathcal{E}_{n,k}:= \{(\pi,I): \pi\in S_n,\; I=(i_1<\cdots<i_k),\; \pi(i_1)<\cdots<\pi(i_k)\}.\] Thus \(I\) is a list of positions whose entries form an increasing subsequence of \(\pi\) of length \(k\).

We sample a pair \((\pi,I)\in\mathcal{E}_{n,k}\) as follows.

1. Choose a \(k\)-subset \(V\subseteq[n]\) uniformly.

2. Choose a \(k\)-subset \(I\subseteq[n]\) uniformly.

3. Place the values of \(V\) in increasing order into the positions \(I\).

4. Fill the remaining \(m\) positions with the remaining \(m\) values in a uniformly random order.

We identify \(I\) with its increasing list of positions.

For a permutation \(\pi\) with \(\operatorname{LIS}(\pi)=k\), let \(L^\star(\pi)\) denote the leftmost LIS of \(\pi\): among all increasing subsequences of length \(k\), choose the one whose list of positions is lexicographically smallest.

3.1 Testing one proposal in \(O(n\log \log n)\) time↩︎

We first show that an individual proposal can be accepted or rejected in \(O(n\log \log n)\) time. For a permutation \(\pi\in S_n\) and an index \(i\in[n]\), let \(r_i\) be the length of the longest increasing subsequence of \(\pi\) that starts at position \(i\).

This is the standard patience-sorting dynamic program with the ordered set of current tail values implemented by an integer predecessor structure, such as a van Emde Boas structure. See, e.g., Crochemore and Porat [2].

Assume now that \(\operatorname{LIS}(\pi)=k\). We construct the leftmost LIS greedily. Starting with the dummy values \(i_0:=0\) and \(v_0:=0\), suppose that \(i_1,\dots,i_{t-1}\) have already been chosen and write \(v_{t-1}:=\pi(i_{t-1})\). We choose \(i_t\) to be the smallest index \(i>i_{t-1}\) such that \(\pi(i)>v_{t-1}\) and such that an increasing subsequence of length \(k-t+1\) can still be completed starting from \(i\), that is, \[r_i\ge k-t+1.\] We then set \(v_t:=\pi(i_t)\).

3.2 Acceptance probability↩︎

For a partition \(\nu\) and a box \(u=(i,j)\) in its Young diagram, let \(h_\nu(u)\) denote the hook length of \(u\), namely the number of boxes consisting of \(u\) itself, the boxes to the right of \(u\) in the same row, and the boxes below \(u\) in the same column. Equivalently, \[h_\nu(i,j)=\nu_i-j+\nu'_j-i+1.\] Write \(H_\nu:=\prod_{u\in\nu}h_\nu(u)\). We use the hook-length formula of Frame–Robinson–Thrall, \[f^\nu=\frac{|\nu|!}{H_\nu},\] where \(f^\nu=|\operatorname{SYT}(\nu)|\) is the number of standard Young tableaux of shape \(\nu\) [24].

For a partition \(\lambda=(\lambda_1,\lambda_2,\dots)\), let \(\lambda'\) denote the conjugate partition, whose parts are the column lengths of the Young diagram of \(\lambda\). For the empty partition we use the convention \(\mu_1=0\) and interpret empty products as \(1\).

Any shape with first row of length \(k\) can be written as \(\lambda=(k,\mu)\), where \(\mu\vdash m\) is the sub-diagram below the first row. The factor \(a_k(\mu)\) measures the multiplicative loss in \(f^{(k, \mu)}\) caused by the elongation of the first-row hooks relative to a bare row of length \(k\). In column \(j\), the hook length in an isolated row would be \(k-j+1\). Placing \(\mu\) underneath elongates it by the column height \(\mu'_j\) to \(k-j+1+\mu'_j\). Columns \(j>\mu_1\) are unaffected. For \(\mu\vdash m\), define \[a_k(\mu):= \begin{cases} \displaystyle \prod_{j=1}^{\mu_1}\frac{k-j+1}{k-j+1+\mu'_j}, & \mu_1\le k,\\[1.2em] 0,& \mu_1>k. \end{cases}\] The hook-length formula then yields \(f^{(k,\mu)}=\binom{n}{m}f^\mu\,a_k(\mu)\), where the binomial coefficient accounts for the choice of which \(m\) labels descend below the first row. Squaring and summing over \(\mu\) will express \(|\Omega_{n,k}|\) as a Plancherel expectation, from which the acceptance probability of the rejection sampler follows directly. To that end, define \[A_{n,k}:=\mathbf{E}_{\operatorname{Pl}_m}[a_k(\mu)^2], \qquad \operatorname{Pl}_m(\mu)=\frac{(f^\mu)^2}{m!}.\]

To bound \(A_{n,k}\) below by a positive constant, we need that \(\mu_1\) is not too large under Plancherel measure.

Lemma 4 gives the tail bound on \(\mu_1\) needed to control \(a_k(\mu)\) on a set of positive Plancherel probability, from which the constant lower bound follows.

4 A polynomial-time algorithm for general \(k\)↩︎

We now turn to the general sampler. By the Robinson–Schensted correspondence, sampling a uniform permutation from \(\Omega_{n,k}\) reduces to two steps. First, sample a partition \(\lambda\vdash n\) from \[\mathbf{P}_{n,k}(\lambda)\propto (f^\lambda)^2\,\mathbf{1}_{\{\lambda_1=k\}},\] and second, sample two independent uniform tableaux of shape \(\lambda\) and apply the inverse correspondence. The second step is standard and is deferred to Section 4.4. The first step is the main algorithmic challenge.

4.0.0.1 Overview of the shape sampler.

The central challenge is to sample the Robinson–Schensted shape conditioned on first row length \(k\), where the target weight \((f^\lambda)^2\) is given by the hook-length formula. Direct enumeration over the up to \(\binom{n}{k}\) feasible shapes is not polynomial. Our algorithm overcomes this through three steps, each removing a distinct obstacle.

Step 1: Reparameterisation separates the factorial weights. We conjugate the partition and apply the shift \(x_i:=\rho_i+k-i\), rewriting the hook-length formula as \[(f^\lambda)^2\propto \Delta(x)^2\prod_{i=1}^k x_i!^{-2}.\] Thus the factorial part becomes a product of per-coordinate weights, while the remaining coordinate coupling is isolated in the squared Vandermonde factor \(\Delta(x)^2\). This separation is what makes the later Cauchy–Binet determinant oracle applicable.

Step 2: Sequential sampling reduces the problem to completion scores. The separated weight admits a one-coordinate-at-a-time sampler. Given a prefix \(p=(p_1,\dots,p_r)\), sample the next coordinate \(t\) with probability proportional to the total weight of all completions of \((p_1,\dots,p_r,t)\). These completion scores give the exact conditional distribution by construction. The obstacle is computing them. Each score sums over all valid suffixes, and the number of possible suffixes can be exponential in \(n\).

Step 3: Cauchy–Binet collapses the suffix sum. Hankel structure makes the oracle fast. Every suffix contribution contains \(\Delta(A)^2\), the squared Vandermonde product over the suffix elements. By the Cauchy–Binet formula this entire sum equals a single coefficient of the determinant of a polynomial matrix, replacing an exponential enumeration by a polynomial-time determinant. A further observation drives the faster implementation. When that matrix is evaluated at a fixed scalar \(\xi\), its \((i,j)\) entry depends only on \(i+j\), making it Hankel. A Hankel determinant of size \(L\) can be computed in \(\tilde{O}(L)\) arithmetic operations rather than \(O(L^3)\). Combined with an incremental sweep over candidates \(t\), this reduces the per-stage cost from \(\tilde{O}(n^4 k L^3)\) to \(\tilde{O}(n^3 k L^2)\) word operations, and the total sampling cost from \(\tilde{O}(n^4 k^5)\) to \(\tilde{O}(n^3 k^4)\).

Algorithm 3 gives the overall structure. The rest of this section formalises the weight formula, defines the completion scores, gives the determinant oracle, and proves correctness of the sequential shape sampler. Subsection 4.4 completes the reduction to permutations. Sections 5 and 6 analyse the direct and fast implementations respectively.

4.1 Reparameterisation↩︎

We work with the conjugate partition \(\rho=\lambda'\) rather than \(\lambda\) directly, its parts \(\rho_i\) are the column heights of \(\lambda\), and it is these that appear naturally in the hook-length formula. Since \(\lambda_1=k\), the partition \(\rho\) has exactly \(k\) positive parts, where \(\rho_i\) counts the rows of \(\lambda\) of length at least \(i\), satisfying \[\rho_1\ge\rho_2\ge\cdots\ge\rho_k\ge 1,\qquad \rho_1+\cdots+\rho_k=n.\] Set \(x_i:=\rho_i+k-i\) and \(N:=n+\binom{k}{2}\). The shift turns the weakly decreasing sequence \(\rho\) into the strictly decreasing sequence \[n\ge x_1>x_2>\cdots>x_k\ge 1,\qquad x_1+\cdots+x_k=N,\] and the map \(x\mapsto\rho\) given by \(\rho_i=x_i-k+i\) is its inverse. Thus the shapes \(\lambda\vdash n\) with \(\lambda_1=k\) are in bijection with the strict partitions \(x\) of \(N\) into exactly \(k\) parts, each at most \(n\). This is the same shift Johansson used to express Plancherel measure as a squared Vandermonde factor times independent weights [25]. Here we repurpose it for purely algorithmic purposes and apply it to the conjugate partition.

For a strict partition \(x=(x_1,\dots,x_k)\), write \(\Delta(x):=\prod_{1\le i<j\le k}(x_i-x_j)\).

The following lemma makes the factored form of the target weight explicit. It is the reason the \(x\) coordinates are useful. The hook-length formula, when rewritten in terms of \(x\), separates into a Vandermonde factor and independent factorial weights, yielding a target distribution that admits an efficient sequential sampler.

We sample the strict sequence \(x_1>\cdots>x_k\) one coordinate at a time. At each step, having fixed a prefix \(p=(p_1,\dots,p_r)\) with \(p_1>\cdots>p_r\ge 1\), we choose the next coordinate \(t\) with probability proportional to the total \(\operatorname{wt}\)-mass of all valid extensions of \((p_1,\dots,p_r,t)\) to a full strict sequence. These masses are the completion scores \(S_p(t)\), which we define and compute below.

At a prefix \(p=(p_1,\ldots,p_r)\) of length \(r<k\), define the upper bound \[B(p):= \begin{cases} n, & r=0,\\ p_r-1, & r>0. \end{cases}\] The admissible candidates for the next coordinate are then \[t\in\{1,\ldots,B(p)\}.\]

For such an admissible \(t\), write \[L:=k-r-1,\quad s:=N-\textstyle\sum_{i=1}^r p_i-t,\quad E_t:=\{1,\dots,t-1\}.\] A valid suffix is an \(L\)-element subset \(A\subseteq E_t\) with \(\sum A=s\). For a finite set \(A\), write \[\Delta(A):=\prod_{\substack{a,b\in A\\a<b}}(b-a),\] with the convention \(\Delta(\varnothing)=1\). For a weight function \(\psi\) on \(E\), the completion counter is \[Z_{L,s}(E,\psi) :=\sum_{\substack{A\subseteq E\\|A|=L,\;\sum A=s}} \Delta(A)^2\prod_{a\in A}\psi(a),\] with \(Z_{0,s}(E,\psi)=\mathbf{1}_{s=0}\) since the only suffix of length zero is \(A=\varnothing\).

Since the sampler uses only relative probabilities at a fixed prefix, we may scale all completion masses by any common positive factor. We clear the factorial denominators of the unrevealed coordinates by multiplying through by \((n!)^{2(L+1)}\). For a possible later value \(a\), define \[G_p(a):=\Bigl(\frac{n!}{a!}\Bigr)^2\prod_{i=1}^r(p_i-a)^2,\] recording the scaled contribution of \(a\) and its pairwise interactions with the prefix. Fixing the next coordinate to \(t\) introduces a further pairwise factor \((t-a)^2\) for each suffix element \(a\), giving \[\psi_{p,t}(a):=(t-a)^2\,G_p(a).\] The completion score is then \[S_p(t):=G_p(t)\,Z_{L,s}(E_t,\psi_{p,t}), \qquad 1\le t\le B(p).\]

4.2 A determinant oracle for completion counts↩︎

Recall that the completion score \(S_p(t)\) involves the counter \(Z_{L,s}(E_t,\psi_{p,t})\), a sum of \(\Delta(A)^2\)-weighted terms over all \(L\)-element subsets \(A\subseteq E_t\) with \(\sum A=s\). Enumerating these subsets directly is infeasible as there are \(\binom{t-1}{L}\) such subsets, which can be exponential in \(n\).The following lemma shows that \(Z_{L,s}(E,\psi)\) is instead recoverable as a single coefficient of the determinant of an \(L\times L\) polynomial matrix, and a determinant can be computed in polynomial-time.

The key mechanism is the factorisation \(M_{L,E,\psi}(q)=VD(q)V^\top\), where \(V\) is a Vandermonde matrix and \(D(q)\) is diagonal. The Cauchy–Binet formula expands the determinant as a sum over \(L\)-element index sets \(J\), each contributing \(\det(V_J)^2=\Delta(A)^2\) and a monomial \(q^{\sum A}\). The coefficient of \(q^s\) then collects exactly the terms in \(Z_{L,s}(E,\psi)\).

Let \(E=\{e_1<\cdots<e_m\}\) be a finite set of positive integers, \(R\) a commutative ring, and \(\psi:E\to R\) a weight function. For \(L\ge 0\), define the following \(L\times L\) matrix over \(R[q]\): \[M_{L,E,\psi}(q):=\sum_{y\in E}\psi(y)q^y\,v(y)v(y)^\top, \quad v(y):=(1,y,\dots,y^{L-1})^\top\in R^L.\]

By Lemma 8, computing \(Z_{L,s}(E_t,\psi_{p,t})\) requires only extracting the coefficient \[[q^s]\det M_{L,E_t,\psi_{p,t}}(q),\] and the completion score is then \[S_p(t)=G_p(t)\,[q^s]\det M_{L,E_t,\psi_{p,t}}(q).\] Thus, for every prefix \(p\) and every candidate next coordinate \(t\), the total weight of all completions below \(t\) can be computed exactly by a determinant coefficient. We now use these scores to assemble the one-coordinate-at-a-time shape sampler and prove that its output has the required distribution. The running-time analysis of the direct and faster implementations is deferred to Sections 5 and 6.

4.3 The sequential shape sampler↩︎

We now assemble the one-coordinate-at-a-time sampler. Start with the empty prefix \(p=\varnothing\). At each stage, compute the completion scores \(S_p(t)\) for all admissible \(t\), then draw the next coordinate with probability \[\mathbf{P}(x_{r+1}=t\mid p) = \frac{S_p(t)}{\sum_{u=1}^{B(p)}S_p(u)}.\] For an admissible candidate \(t\), write \[s(t):=N-\sum_{i=1}^r p_i-t.\] A suffix of length \(L\) chosen from \(E_t=\{1,\ldots,t-1\}\) can have total sum only between the sum of the \(L\) smallest positive integers and the sum of the \(L\) largest elements below \(t\). Hence \(S_p(t)=0\) unless \[\binom{L+1}{2}\le s(t)\le Lt-\binom{L+1}{2}.\] Candidates failing this test may therefore be skipped before invoking the determinant oracle. After \(k\) coordinates have been chosen, recover \(\rho_i=x_i-k+i\) and set \(\lambda=\rho'\).

For a prefix \(p\), let \(W(p)\) denote the total \(\operatorname{wt}\)-mass of all valid full strict sequences extending \(p\).

4.4 From sampled shapes to sampled permutations↩︎

The remaining steps are standard. We use the Robinson–Schensted bijection between permutations of \([n]\) and pairs \((P,Q)\) of standard Young tableaux of the same shape, together with its inverse reverse-insertion algorithm [27]. We also use the Greene–Nijenhuis–Wilf hook walk, which samples a uniform standard Young tableau of any fixed shape [28].

5 The running time of a direct implementation of the general sampler↩︎

We now derive the running time of a direct implementation of the general sampler. We assume the word-RAM model with word size \(\Theta(\log n)\). All exact random choices are made from unbiased random bits¹.

Fix an already sampled prefix \(p=(p_1,\dots,p_r)\) of \(x\), and write \[L:=k-r-1.\] Recall the reusable stage weights \[G_p(a)=\Bigl(\frac{n!}{a!}\Bigr)^2\prod_{i=1}^r(p_i-a)^2, \qquad a=1,\dots,n,\] and \[\psi_{p,t}(a)=(t-a)^2G_p(a).\] If the next chosen coordinate is \(u\), then \[G_{(p,u)}(a)=G_p(a)(u-a)^2.\] Thus the array \(\{G_p(a):1\le a\le n\}\) may be updated from one stage to the next by a single pass over \(a=1,\dots,n\). This update cost is lower order in the bounds below.

We now explain how these integer completion scores are computed exactly in the word-RAM model. By Lemma 10, every completion score arising at a stage with suffix length \(L\) has \(O((L+1)n\log n)\) bits. Hence, for a sufficiently large absolute constant \(C\), every completion score arising anywhere in the algorithm has at most \[B:=Ckn\log n\] bits.

Let \[d_{\max}:=n(k-1)+1\le n^2+1 .\] All interpolation steps below use at most \(d_{\max}\) evaluation points. We therefore fix a set \(\mathcal{P}\) of distinct \(\Theta(\log n)\)-bit primes, each larger than \(d_{\max}\), such that \(\prod_{\mathfrak p\in\mathcal{P}}\mathfrak p>2^B\).

Since each prime contributes \(\Theta(\log n)\) bits to the product, it is enough to take \(|\mathcal{P}|=\tilde{O}(kn)\) such primes. For each \(\mathfrak p\in\mathcal{P}\), we run the same computation over \(\mathbf{F}_{\mathfrak p}\). Afterwards, each integer completion score is recovered exactly by Chinese remaindering.

In the proofs below we use the following CRT accounting convention. We first count field operations over a single word-size prime field \(\mathbf{F}_{\mathfrak p}\). Since the computation is repeated for \(\tilde{O}(kn)\) primes, this field operation count is then multiplied by \(\tilde{O}(kn)\). The final integer completion scores are recovered by Chinese remaindering.

6 Speeding up the completion-score oracle↩︎

The direct implementation computes each candidate completion score separately. We now give a faster implementation of the same oracle. The sequential sampler, the completion scores \(S_p(t)\), and the correctness proof (Theorem 8) are unchanged. Only the cost of evaluating the completion scores decreases.

Fix a prefix \(p\), a candidate next coordinate \(t\), and an evaluation point \(\xi\). The polynomial matrix from Lemma 8 becomes \[M_{L,E_t,\psi_{p,t}}(\xi) = \sum_{a<t}(t-a)^2G_p(a)\xi^a\,v(a)v(a)^\top .\] Using zero-based indices \(0\le i,j<L\), its \((i,j)\) entry is \[\sum_{a<t}(t-a)^2G_p(a)\xi^a a^{i+j}.\] Thus the entry depends only on \(i+j\), and so the evaluated matrix is Hankel. Equivalently, it is determined by the \(2L-1\) anti-diagonal quantities corresponding to \(i+j=0,1,\ldots,2L-2\), rather than by all \(L^2\) entries. These anti-diagonal quantities are weighted power sums, and we refer to them as moments.

We use the following structured linear-algebra primitive.

We now describe how the primitive is used at one sampling stage. Fix a prefix \(p\) and suffix length \(L\ge 1\). By the CRT accounting convention of Section 5, the computation is performed over the fixed set \(\mathcal{P}\) of word-size prime fields chosen there. At this stage, the polynomial \(D_t(q)\) has degree at most \(Ln\), so we use \[d:=Ln+1\] distinct nonzero evaluation points for interpolation. Since every \(\mathfrak p\in\mathcal{P}\) is larger than \[d_{\max}=(k-1)n+1\ge d,\] the field \(\mathbf{F}_{\mathfrak p}\) contains enough nonzero evaluation points at every stage.

We first count field operations over a single word-size prime field \(\mathbf{F}_{\mathfrak p}\), for \(\mathfrak p\in\mathcal{P}\). The same computation is then repeated for the \(\tilde{O}(kn)\) primes in \(\mathcal{P}\), and the exact integer completion scores are recovered by Chinese remaindering. Thus an arithmetic bound of \(T\) field operations for one stage gives a \(\tilde{O}(Tkn)\) word-RAM bound.

For a fixed evaluation point \(\xi\in\mathbf{F}_{\mathfrak p}\), define \[m^{(\gamma)}_\ell(t,\xi) := \sum_{a<t} a^{\gamma+\ell}G_p(a)\xi^a, \qquad \gamma=0,1,2, \qquad 0\le \ell\le 2L-2.\] These are the anti-diagonal moments of the cumulative Hankel matrices \[A_\gamma^{(<t)}(\xi) := \sum_{a<t}a^\gamma G_p(a)\xi^a v(a)v(a)^\top, \qquad \gamma=0,1,2.\] Indeed, since \(v(a)v(a)^\top\) has zero-based \((\alpha,\beta)\) entry \(a^{\alpha+\beta}\), the \((\alpha,\beta)\) entry of \(A_\gamma^{(<t)}(\xi)\) is \[\sum_{a<t}a^{\gamma+\alpha+\beta}G_p(a)\xi^a = m^{(\gamma)}_{\alpha+\beta}(t,\xi),\] which depends only on \(\alpha+\beta\).

Using \[(t-a)^2=t^2-2ta+a^2,\] we obtain \[M_{L,E_t,\psi_{p,t}}(\xi) = t^2A_0^{(<t)}(\xi) - 2tA_1^{(<t)}(\xi) + A_2^{(<t)}(\xi).\] Therefore \(M_{L,E_t,\psi_{p,t}}(\xi)\) is Hankel with moment sequence \[h_\ell(t,\xi) = t^2m^{(0)}_\ell(t,\xi) - 2t\,m^{(1)}_\ell(t,\xi) + m^{(2)}_\ell(t,\xi), \qquad 0\le \ell\le 2L-2.\] Writing \[D_t(q):=\det M_{L,E_t,\psi_{p,t}}(q),\] we therefore have \[D_t(\xi) = \det\bigl(h_{\alpha+\beta}(t,\xi)\bigr)_{0\le \alpha,\beta<L}.\]

References↩︎

---

1. Throughout, exact weighted choices are made as follows. Given nonnegative integer weights \(w_1,\dots,w_M\) with sum \(W\), we draw a uniform integer from \(\{1,\dots,W\}\) by rejection sampling from a power-of-two range and select the corresponding interval in the cumulative sums. This costs expected \(O(\log W)\) random bits and arithmetic quasilinear in the bit length of \(W\).↩︎