Problem 916: Restricted Permutations

Project Euler Problem 916 Solution

EulerSolve provides an optimized solution for Project Euler Problem 916, Restricted Permutations, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each $n$, let $P(n)$ be the number of permutations $\pi$ of $\{1,2,\dots,2n\}$ satisfying $$\operatorname{LIS}(\pi)\le n+1,\qquad \operatorname{LDS}(\pi)\le 2.$$ The task is to compute $P(100000000)$ modulo $10^9+7$. A direct count is completely infeasible, because the ambient space has $(2n)!$ permutations. The solution works only because the two subsequence constraints force the associated Young diagram into an extremely small family. Mathematical Approach Define $$P(n)=\#\left\{\pi\in S_{2n}:\operatorname{LIS}(\pi)\le n+1,\ \operatorname{LDS}(\pi)\le 2\right\}.$$ The key move is to stop counting permutations directly and instead count the possible RSK shapes and the number of permutations attached to each such shape. From subsequence bounds to Young diagram shapes The Robinson-Schensted correspondence sends every permutation $\pi\in S_{2n}$ to a partition $\lambda\vdash 2n$ with $$\lambda_1=\operatorname{LIS}(\pi),\qquad \lambda'_1=\operatorname{LDS}(\pi).$$ Here $\lambda_1$ is the first row length, while $\lambda'_1$ is the first column length, equivalently the total number of rows. Therefore $\operatorname{LDS}(\pi)\le 2$ means that $\lambda$ has at most two rows. So every admissible shape has the form $\lambda=(a,b)$ with $a\ge b\ge 0$ and $a+b=2n$. The bound on the longest increasing subsequence gives $a\le n+1$....

Detailed mathematical approach

Problem Summary

For each $n$, let $P(n)$ be the number of permutations $\pi$ of $\{1,2,\dots,2n\}$ satisfying

$$\operatorname{LIS}(\pi)\le n+1,\qquad \operatorname{LDS}(\pi)\le 2.$$

The task is to compute $P(100000000)$ modulo $10^9+7$. A direct count is completely infeasible, because the ambient space has $(2n)!$ permutations. The solution works only because the two subsequence constraints force the associated Young diagram into an extremely small family.

Mathematical Approach

Define

$$P(n)=\#\left\{\pi\in S_{2n}:\operatorname{LIS}(\pi)\le n+1,\ \operatorname{LDS}(\pi)\le 2\right\}.$$

The key move is to stop counting permutations directly and instead count the possible RSK shapes and the number of permutations attached to each such shape.

From subsequence bounds to Young diagram shapes

The Robinson-Schensted correspondence sends every permutation $\pi\in S_{2n}$ to a partition $\lambda\vdash 2n$ with

$$\lambda_1=\operatorname{LIS}(\pi),\qquad \lambda'_1=\operatorname{LDS}(\pi).$$

Here $\lambda_1$ is the first row length, while $\lambda'_1$ is the first column length, equivalently the total number of rows. Therefore $\operatorname{LDS}(\pi)\le 2$ means that $\lambda$ has at most two rows.

So every admissible shape has the form $\lambda=(a,b)$ with $a\ge b\ge 0$ and $a+b=2n$. The bound on the longest increasing subsequence gives $a\le n+1$. But every two-row partition of size $2n$ also satisfies $a\ge (a+b)/2=n$. Hence only two values of $a$ are possible:

$$a\in\{n,n+1\}.$$

This leaves exactly two admissible shapes,

$$\lambda=(n,n),\qquad \lambda=(n+1,n-1).$$

That is the decisive simplification: the original count over $(2n)!$ permutations collapses to a sum over just two Young shapes.

Why a fixed shape contributes the square of a tableau count

RSK is a bijection between permutations and pairs $(P,Q)$ of standard Young tableaux of the same shape. If $f^\lambda$ denotes the number of standard Young tableaux of shape $\lambda$, then the number of permutations with RSK shape $\lambda$ is

$$\#\{\pi\in S_{2n}:\operatorname{shape}(\pi)=\lambda\}=(f^\lambda)^2.$$

So once the two admissible shapes are known, the problem is reduced to evaluating $f^\lambda$ for those two shapes and squaring the resulting tableau counts.

Hook-length evaluation for the two-row family

For a two-row shape $(a,b)$, the hook lengths can be written explicitly. In the first row, the leftmost $b$ cells have hook lengths $a-j+2$ for $1\le j\le b$, and the remaining $a-b$ cells have hook lengths $a-j+1$ for $b<j\le a$. In the second row, the hook lengths are $b-j+1$ for $1\le j\le b$.

Multiplying them gives

$$\prod_{u\in(a,b)} h(u)=\left(\prod_{j=1}^{b}(a-j+2)\right)\left(\prod_{j=b+1}^{a}(a-j+1)\right)\left(\prod_{j=1}^{b}(b-j+1)\right)=\frac{(a+1)!\,b!}{a-b+1}.$$

Now the hook-length formula yields

$$f^{(a,b)}=\frac{(a+b)!}{\prod_{u\in(a,b)} h(u)}=\frac{a-b+1}{a+1}\binom{a+b}{b}.$$

Substituting the two admissible shapes gives

$$f^{(n,n)}=\frac{1}{n+1}\binom{2n}{n},$$

$$f^{(n+1,n-1)}=\frac{3}{n+2}\binom{2n}{n-1}=\frac{3n}{(n+1)(n+2)}\binom{2n}{n}.$$

Therefore

$$P(n)=\left(\frac{1}{n+1}\binom{2n}{n}\right)^2+\left(\frac{3n}{(n+1)(n+2)}\binom{2n}{n}\right)^2 \pmod{10^9+7}.$$

The first term is the square of the Catalan number $C_n$. The second term comes from the neighboring shape $(n+1,n-1)$, and there are no other contributions because no other shapes survive the subsequence bounds.

Reducing everything to one central binomial coefficient

It is convenient to set

$$B_n=\binom{2n}{n}.$$

Then the answer becomes

$$P(n)=\left(B_n (n+1)^{-1}\right)^2+\left(3n\,B_n (n+1)^{-1}(n+2)^{-1}\right)^2 \pmod{10^9+7}.$$

The identity

$$\binom{2n}{n-1}=B_n\frac{n}{n+1}$$

is what allows the second shape to be expressed using the same central binomial coefficient. This is exactly why the implementations need only one factorial ratio rather than two unrelated binomial computations.

For the target value $n=100000000$, we have $2n=200000000<10^9+7$. So $n!$, $n+1$, and $n+2$ are all nonzero modulo the prime $10^9+7$, and every denominator in the closed form is invertible.

Worked example: $n=2$

When $n=2$, the admissible shapes are $(2,2)$ and $(3,1)$. The hook-length formula gives

$$f^{(2,2)}=\frac{1}{3}\binom{4}{2}=2,\qquad f^{(3,1)}=\frac{3}{4}\binom{4}{1}=3.$$

Hence

$$P(2)=2^2+3^2=13.$$

So among the $4!=24$ permutations of $\{1,2,3,4\}$, exactly 13 satisfy $\operatorname{LIS}(\pi)\le 3$ and $\operatorname{LDS}(\pi)\le 2$. This is the smallest nontrivial checkpoint used to confirm the closed form.

How the Code Works

The C++, Python, and Java implementations all evaluate the closed form above. None of them constructs Young tableaux or searches through permutations for the target case.

Single-pass factorial accumulation

The implementation multiplies the integers $1,2,\dots,2n$ modulo $10^9+7$. When the loop reaches $n$, it records the current product, so the same pass produces both $n!$ and $(2n)!$. From these two quantities it forms

$$\binom{2n}{n}\equiv (2n)!\,(n!)^{-2}\pmod{10^9+7}.$$

Recovering the two tableau contributions

Next the implementation computes the modular inverses of $n!$, $n+1$, and $n+2$ by fast exponentiation and Fermat's little theorem,

$$a^{-1}\equiv a^{10^9+5}\pmod{10^9+7}.$$

With those inverses available, it evaluates the contribution of the shape $(n,n)$, the contribution of the shape $(n+1,n-1)$, squares both values, adds them, and reduces modulo $10^9+7$. The whole computational core is just a linear sweep plus a constant number of modular exponentiations.

Tiny exhaustive checks

One compiled implementation also performs brute-force verification for very small $n$, checking that direct enumeration agrees with the closed form and that known small answers are reproduced. Those checks validate the mathematics, but they are separate from the large-$n$ computation.

Complexity Analysis

The dominant cost is the single loop from $1$ to $2n$, so the runtime is $O(n)$. The modular inverses require only a constant number of exponentiations, each costing $O(\log(10^9+7))$, which does not change the overall asymptotic bound. Memory usage is $O(1)$, because only a few scalar accumulators are stored. The tiny brute-force checks run only for fixed small inputs and therefore do not affect the asymptotic complexity.

Footnotes and References

Problem page: Project Euler 916
Robinson-Schensted correspondence: Wikipedia - Robinson-Schensted correspondence
Hook-length formula: Wikipedia - Hook-length formula
Standard Young tableaux: Wikipedia - Young tableau
Longest increasing subsequence: Wikipedia - Longest increasing subsequence
Catalan number: Wikipedia - Catalan number

Problem 916 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <vector>

namespace {

using i64 = std::int64_t;

constexpr i64 kMod = 1'000'000'007LL;

i64 mod_pow(i64 a, i64 e) {
    i64 r = 1;
    i64 x = a % kMod;
    i64 p = e;
    while (p > 0) {
        if (p & 1LL) {
            r = static_cast<i64>((__int128)r * x % kMod);
        }
        x = static_cast<i64>((__int128)x * x % kMod);
        p >>= 1LL;
    }
    return r;
}

i64 mod_inv(i64 a) {
    return mod_pow((a % kMod + kMod) % kMod, kMod - 2);
}

i64 P_formula(std::int64_t n) {
    const std::int64_t two_n = 2 * n;

    i64 fact_n = 1;
    i64 fact_2n = 1;
    for (std::int64_t i = 1; i <= two_n; ++i) {
        fact_2n = static_cast<i64>((__int128)fact_2n * (i % kMod) % kMod);
        if (i == n) {
            fact_n = fact_2n;
        }
    }

    const i64 inv_fact_n = mod_inv(fact_n);
    const i64 c2n_n = static_cast<i64>((__int128)fact_2n * inv_fact_n % kMod * inv_fact_n % kMod);

    const i64 inv_n1 = mod_inv(n + 1);
    const i64 inv_n2 = mod_inv(n + 2);

    const i64 f1 = static_cast<i64>((__int128)c2n_n * inv_n1 % kMod);
    const i64 c2n_n1 = static_cast<i64>((__int128)c2n_n * (n % kMod) % kMod * inv_n1 % kMod);
    const i64 f2 = static_cast<i64>((__int128)3 * c2n_n1 % kMod * inv_n2 % kMod);

    i64 ans = static_cast<i64>((__int128)f1 * f1 % kMod);
    ans += static_cast<i64>((__int128)f2 * f2 % kMod);
    if (ans >= kMod) {
        ans -= kMod;
    }
    return ans;
}

int lis_len(const std::vector<int>& p) {
    const int n = static_cast<int>(p.size());
    std::vector<int> dp(n, 1);
    int best = 0;
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < i; ++j) {
            if (p[j] < p[i]) {
                dp[i] = std::max(dp[i], dp[j] + 1);
            }
        }
        best = std::max(best, dp[i]);
    }
    return best;
}

int lds_len(const std::vector<int>& p) {
    const int n = static_cast<int>(p.size());
    std::vector<int> dp(n, 1);
    int best = 0;
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < i; ++j) {
            if (p[j] > p[i]) {
                dp[i] = std::max(dp[i], dp[j] + 1);
            }
        }
        best = std::max(best, dp[i]);
    }
    return best;
}

std::int64_t brute_count(int n) {
    const int m = 2 * n;
    std::vector<int> p(m);
    for (int i = 0; i < m; ++i) {
        p[i] = i + 1;
    }

    std::int64_t cnt = 0;
    do {
        if (lis_len(p) <= n + 1 && lds_len(p) <= 2) {
            ++cnt;
        }
    } while (std::next_permutation(p.begin(), p.end()));

    return cnt;
}

void validate() {
    for (int n = 1; n <= 4; ++n) {
        assert(P_formula(n) == brute_count(n));
    }

    assert(P_formula(2) == 13);
    assert(P_formula(10) == 45'265'702);
}

}  // namespace

int main() {
    validate();
    std::cout << P_formula(100'000'000) << '\n';
    return 0;
}

Python

MOD = 1000000007

def mod_pow(a, e):
    return pow(a, e, MOD)

def mod_inv(a):
    return pow(a % MOD, MOD - 2, MOD)

def P_formula(n):
    two_n = 2 * n

    fact_n = 1
    fact_2n = 1
    for i in range(1, two_n + 1):
        fact_2n = (fact_2n * i) % MOD
        if i == n:
            fact_n = fact_2n

    inv_fact_n = mod_inv(fact_n)
    c2n_n = (fact_2n * inv_fact_n % MOD * inv_fact_n) % MOD

    inv_n1 = mod_inv(n + 1)
    inv_n2 = mod_inv(n + 2)

    f1 = (c2n_n * inv_n1) % MOD
    c2n_n1 = (c2n_n * (n % MOD) % MOD * inv_n1) % MOD
    f2 = (3 * c2n_n1 % MOD * inv_n2) % MOD

    ans = (f1 * f1) % MOD
    ans = (ans + f2 * f2) % MOD
    return ans

def solve():
    return str(P_formula(100000000))

if __name__ == "__main__":
    print(solve())

Java

import java.math.BigInteger;

public class Euler916 {

    static final long MOD = 1000000007L;

    static long modPow(long a, long e) {
        long r = 1;
        long x = a % MOD;
        long p = e;
        while (p > 0) {
            if ((p & 1L) != 0) {
                r = (r * x) % MOD;
            }
            x = (x * x) % MOD;
            p >>= 1L;
        }
        return r;
    }

    static long modInv(long a) {
        return modPow((a % MOD + MOD) % MOD, MOD - 2);
    }

    static long PFormula(long n) {
        long twoN = 2 * n;

        long factN = 1;
        long fact2N = 1;
        for (long i = 1; i <= twoN; ++i) {
            fact2N = (fact2N * (i % MOD)) % MOD;
            if (i == n) {
                factN = fact2N;
            }
        }

        long invFactN = modInv(factN);
        long c2nN = (fact2N * invFactN % MOD * invFactN) % MOD;

        long invN1 = modInv(n + 1);
        long invN2 = modInv(n + 2);

        long f1 = (c2nN * invN1) % MOD;
        long c2nN1 = (c2nN * (n % MOD) % MOD * invN1) % MOD;
        long f2 = (3 * c2nN1 % MOD * invN2) % MOD;

        long ans = (f1 * f1) % MOD;
        ans = (ans + f2 * f2) % MOD;
        return ans;
    }

    public static String solve() {
        return Long.toString(PFormula(100000000));
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}