Look at the elements that are not in the right position. Which of them must be selected? Can selecting any other element help?

Every element that is not in the right position has to be selected; otherwise, it never moves.

Selecting an element that is already in the right position cannot help. For it to stay correct after the reversal, it has to be paired with itself, because its value occurs nowhere else in the permutation. Thus it can only be the middle selected element, and removing it changes nothing.

So, if the array is already sorted, choose any one element. Otherwise, select exactly the misplaced elements and reverse them. The answer is YES if this sorts the array, and NO otherwise.

If the misplaced positions are listed from left to right, the reversal pairs the first with the last, the second with the second-last, and so on. Check that every pair contains the two values those positions need.

This takes $$$O(n)$$$ time and $$$O(n)$$$ memory.

#include <bits/stdc++.h>

using namespace std;

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int t;
    cin >> t;
    while (t--) {
        int n;
        cin >> n;

        vector<int> p(n);
        for (int &x : p) cin >> x;

        vector<int> bad;
        for (int i = 0; i < n; ++i) {
            if (p[i] != i + 1) {
                bad.push_back(i);
            }
        }

        bool ok = true;
        int k = bad.size();
        for (int i = 0; i < k; ++i) {
            if (p[bad[i]] != bad[k - 1 - i] + 1) {
                ok = false;
                break;
            }
        }

        cout << (ok ? "YES\n" : "NO\n");
    }
}

Write out the score for a few terms and collect the coefficient of each chosen value.

Fix the last chosen position. Which $$$m-1$$$ earlier values are best, and how can we maintain them while moving this position to the right?

Expanding the score gives

Only the last chosen value has a positive coefficient. Once we decide that $$$a_j$$$ is last, we should take the $$$m-1$$$ smallest values before it. Their positions are already in the right order, so they always form a valid subsequence together with $$$a_j$$$.

Scan $$$j$$$ from left to right. Before considering $$$a_j$$$, keep the $$$m-1$$$ smallest values among $$$a_1,\ldots,a_{j-1}$$$ in a max-heap, and let their sum be $$$s$$$. The best score ending at $$$j$$$ is

Then insert $$$a_j$$$ for later positions and remove the maximum if the heap becomes too large. Take the best score over all $$$j$$$. For $$$m=1$$$, the heap is empty and this simply takes the maximum element.

The complexity is $$$O(n\log(m+1))$$$ time and $$$O(m)$$$ extra memory.

#include <bits/stdc++.h>

using namespace std;

using ll = long long;

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int t;
    cin >> t;
    while (t--) {
        int n, m;
        cin >> n >> m;

        vector<ll> a(n);
        for (ll &x : a) cin >> x;

        priority_queue<ll> pq;
        ll sum = 0, ans = LLONG_MIN;

        for (int i = 0; i < n; ++i) {
            if ((int)pq.size() == m - 1) {
                ans = max(ans, m * a[i] - sum);
            }

            pq.push(a[i]);
            sum += a[i];
            if ((int)pq.size() == m) {
                sum -= pq.top();
                pq.pop();
            }
        }

        cout << ans << '\n';
    }
}

After sorting, every parent lies to the right of their child. What does that imply about cycles and the root?

The total cost is a sum of edge costs. Fix one possible edge; how many choices remain for the other children?

Sort the ratings:

The maximum is forced to be the root. Give every other position $$$i$$$ any parent $$$j \gt i$$$. Following parent links always moves right, so there is no cycle and every path ends at the maximum. Thus any such choices form a valid tree.

We cannot enumerate the trees. Instead, consider one possible edge at a time and add its cost once for every tree containing it.

Fix the edge from child $$$i$$$ to parent $$$j \gt i$$$. Its cost is $$$b_j-b_i$$$. Every other child $$$k$$$ still has $$$n-k-1$$$ choices, so this edge appears in

trees. The value $$$W_i$$$ does not depend on $$$j$$$, so all edges leaving $$$i$$$ contribute

A suffix sum gives the sum of ratings to the right of each $$$i$$$. Prefix and suffix products of the numbers $$$n-k-1$$$ give every $$$W_i$$$ without division.

The complexity is $$$O(n\log n)$$$ time and $$$O(n)$$$ memory.

#include <bits/stdc++.h>

using namespace std;

using ll = long long;

const int MOD = 998244353;

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int t;
    cin >> t;
    while (t--) {
        int n;
        cin >> n;
        vector<ll> b(n);
        for (ll &x : b) cin >> x;
        sort(b.begin(), b.end());

        vector<ll> sum(n + 1);
        for (int i = n - 1; i >= 0; --i) {
            sum[i] = sum[i + 1] + b[i];
        }

        vector<ll> pref(n, 1), suff(n, 1);
        for (int i = 0; i < n - 1; ++i) {
            pref[i + 1] = pref[i] * (n - i - 1) % MOD;
        }
        for (int i = n - 2; i >= 0; --i) {
            suff[i] = suff[i + 1] * (n - i - 1) % MOD;
        }

        ll ans = 0;
        for (int i = 0; i < n - 1; ++i) {
            ll diff = sum[i + 1] - (n - i - 1) * b[i];
            ll ways = pref[i] * suff[i + 1] % MOD;
            ans = (ans + diff % MOD * ways) % MOD;
        }
        cout << ans << '\n';
    }
}

Since $$$2\equiv-1\pmod 3$$$, try an alternating prefix sum.

A divisible substring corresponds to two equal prefix remainders. How should the counts of the three remainders be distributed?

With three ones, the prefix remainder is constant on four blocks. Try their boundaries around $$$n/3$$$, $$$2n/3$$$, and the end of the string.

Modulo $$$3$$$, the powers of $$$2$$$ alternate between $$$1$$$ and $$$-1$$$. Define the alternating prefix sum

If $$$v$$$ is the binary value of $$$s_l\ldots s_r$$$, then

Therefore, this substring is divisible by $$$3$$$ exactly when $$$p_{l-1}=p_r$$$.

Suppose the remainders $$$0,1,2$$$ occur $$$c_0,c_1,c_2$$$ times among the $$$n+1$$$ prefix sums. Every pair of equal prefix remainders gives one divisible substring, so

For a fixed sum $$$c_0+c_1+c_2=n+1$$$, this is smallest when the three counts differ by at most one.

Now we only need to build a string with such balanced counts. A zero repeats the current prefix remainder. A one changes it by $$$-1$$$ at an odd position and by $$$+1$$$ at an even position. Thus, three ones divide the prefix sums into four constant blocks. To make the three remainders appear almost equally often, the block boundaries should be near the thirds of the string.

For $$$n\ge3$$$, simply try:

• the first one at $$$\lfloor n/3\rfloor$$$ or $$$\lfloor n/3\rfloor+1$$$;
• the second one at $$$\lfloor2n/3\rfloor$$$ or $$$\lfloor2n/3\rfloor+1$$$;
• either no third one, or a third one at $$$n$$$.

For each of these at most eight strings, count how many prefix sums have each remainder. Output the first one whose three counts differ by at most one.

The searched ranges contain the following construction. Let $$$q=\lfloor n/3\rfloor$$$ and consider a string of length $$$3q+2$$$.

• If $$$q$$$ is odd, put ones at $$$q$$$, $$$2q+1$$$, and $$$3q+2$$$. The four blocks have remainders $$$0,2,1,0$$$ and lengths $$$q,q+1,q+1,1$$$.
• If $$$q$$$ is even, put ones at $$$q+1$$$, $$$2q+1$$$, and $$$3q+2$$$. The blocks have remainders $$$0,2,1,2$$$ and lengths $$$q+1,q,q+1,1$$$.

In both cases, each remainder occurs $$$q+1$$$ times. The real length $$$n$$$ is obtained by cutting off at most two characters. The last two prefix remainders are different, so the remaining counts are still balanced. If position $$$3q+2$$$ was cut off, this is exactly the choice with no third one. Hence, the string is among the eight candidates above.

For $$$n=1$$$ and $$$n=2$$$, output 1 and 10, respectively. For $$$n=2$$$, three distinct prefix remainders are impossible, while 10 has only one divisible substring, so it is optimal.

The complexity is $$$O(n)$$$ time and $$$O(n)$$$ memory per test case.

#include <bits/stdc++.h>

using namespace std;

bool hasBalancedPrefixes(const string& s) {
    array<int, 3> count{1, 0, 0};
    int remainder = 0;

    for (int i = 0; i < (int)s.size(); ++i) {
        if (s[i] == '1') {
            remainder = (remainder + (i % 2 == 0 ? 2 : 1)) % 3;
        }
        ++count[remainder];
    }

    return *max_element(count.begin(), count.end())
         - *min_element(count.begin(), count.end()) <= 1;
}

string constructString(int n) {
    if (n == 1) {
        return "1";
    }
    if (n == 2) {
        return "10";
    }

    for (int first = n / 3; first <= n / 3 + 1; ++first) {
        for (int second = 2 * n / 3; second <= 2 * n / 3 + 1; ++second) {
            for (int putLastOne = 0; putLastOne <= 1; ++putLastOne) {
                if (first >= second || (putLastOne && second >= n)) {
                    continue;
                }

                string s(n, '0');
                s[first - 1] = '1';
                s[second - 1] = '1';
                if (putLastOne) {
                    s[n - 1] = '1';
                }

                if (hasBalancedPrefixes(s)) {
                    return s;
                }
            }
        }
    }

    assert(false);
    return {};
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int t;
    cin >> t;
    while (t--) {
        int n;
        cin >> n;
        cout << constructString(n) << '\n';
    }
    return 0;
}

An operation never changes the order of two values. If the minimum and maximum become equal, what happens to every value between them?

Let $$$dp[x][y]$$$ be the answer for minimum $$$x$$$ and maximum $$$y$$$. For $$$x \lt y$$$, can we decrease $$$y$$$ while leaving $$$x$$$ unchanged?

After computing every $$$dp[x][y]$$$, count subsequences with minimum $$$x$$$ and maximum $$$y$$$.

First ignore the counting of subsequences. Let $$$dp[x][y]$$$ be the largest value to which the pair $$$[x,y]$$$ can be made equal.

Every operation preserves the order of the values. Therefore, if the minimum and maximum of a subsequence become equal, every value between them becomes equal as well. The contribution of a subsequence $$$b$$$ is simply

We can compute this table directly. Clearly, $$$dp[x][x]=x$$$. Now suppose $$$x \lt y$$$ and try to decrease the larger endpoint.

If some prime divides $$$y$$$ but does not divide $$$x$$$, choosing it changes the pair to $$$[x,y-1]$$$. Also, any sequence that equalizes $$$x$$$ and $$$y$$$ must equalize the value $$$y-1$$$ between them. Thus, this case gives $$$dp[x][y-1]$$$.

Otherwise, every prime divisor of $$$y$$$ also divides $$$x$$$, so choosing any of them changes the pair to $$$[x-1,y-1]$$$.

Why can doing something else first not give a better result? Look at the first operation in any successful sequence that decreases $$$y$$$. Immediately after it, the upper endpoint is $$$y-1$$$, while the lower endpoint is already at most $$$x-1$$$: it either decreased earlier, or the same prime decreases it together with $$$y$$$. The remaining operations equalize these two values, so order preservation means that they also equalize the value $$$x-1$$$ between them. Therefore, this case gives $$$dp[x-1][y-1]$$$.

The transition is

Here the limit is only $$$M=3\,000$$$, so we can build the whole $$$M\times M$$$ table. Let need[y] be the product of the distinct prime divisors of $$$y$$$. Sieve all need[y]; the first case applies exactly when x % need[y] == 0. Process $$$y-x$$$ from small to large and $$$x$$$ from small to large, so both states used by the transition are ready.

Now count subsequences by their minimum and maximum. Let $$$c_v$$$ be the number of occurrences of $$$v$$$.

If both endpoints are $$$v$$$, any nonempty subset of its occurrences works, giving $$$2^{c_v}-1$$$ choices.

For $$$x \lt y$$$, choose at least one occurrence of both $$$x$$$ and $$$y$$$. Every occurrence with value strictly between them is optional, so the number of choices is

Use prefix sums of the frequencies for the exponent, multiply each count by $$$dp[x][y]$$$, and add everything to the answer.

The one-time DP takes $$$O(M^2)$$$ time and memory. Each test case takes $$$O(n^2)$$$ time and $$$O(n)$$$ additional memory.

#include <bits/stdc++.h>

using namespace std;

const int MAXV = 3000;
const int MOD = 998244353;

int dp[MAXV + 1][MAXV + 1], need[MAXV + 1], p2[MAXV + 1];

void precompute() {
    fill(need, need + MAXV + 1, 1);
    for (int p = 2; p <= MAXV; ++p) {
        if (need[p] != 1) {
            continue;
        }
        for (int value = p; value <= MAXV; value += p) {
            need[value] *= p;
        }
    }

    for (int x = 1; x <= MAXV; ++x) {
        dp[x][x] = x;
    }
    for (int difference = 1; difference < MAXV; ++difference) {
        for (int x = 1; x + difference <= MAXV; ++x) {
            int y = x + difference;
            if (x % need[y] == 0) {
                dp[x][y] = dp[x - 1][y - 1];
            } else {
                dp[x][y] = dp[x][y - 1];
            }
        }
    }

    p2[0] = 1;
    for (int i = 1; i <= MAXV; ++i) {
        p2[i] = 2LL * p2[i - 1] % MOD;
    }
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    precompute();

    int t;
    cin >> t;
    while (t--) {
        int n;
        cin >> n;
        vector<int> freq(n + 1);
        for (int i = 0; i < n; ++i) {
            int value;
            cin >> value;
            ++freq[value];
        }

        vector<int> prefix(n + 1);
        for (int value = 1; value <= n; ++value) {
            prefix[value] = prefix[value - 1] + freq[value];
        }

        long long answer = 0;
        for (int x = 1; x <= n; ++x) {
            if (freq[x] == 0) {
                continue;
            }
            int chooseX = p2[freq[x]] - 1;
            answer = (answer + 1LL * x * chooseX) % MOD;

            for (int y = x + 1; y <= n; ++y) {
                if (freq[y] == 0) {
                    continue;
                }
                int chooseY = p2[freq[y]] - 1;

                int middle = prefix[y - 1] - prefix[x];
                long long ways = 1LL * chooseX * chooseY % MOD;
                ways = ways * p2[middle] % MOD;
                answer = (answer + ways * dp[x][y]) % MOD;
            }
        }
        cout << answer << '\n';
    }
}

Let $$$dp[x][y]$$$ be the answer for minimum $$$x$$$ and maximum $$$y$$$. E1 builds and sums over every pair. Fix $$$x$$$ and scan larger values of $$$y$$$. When can $$$dp[x][y]$$$ differ from $$$dp[x][y-1]$$$?

Let $$$r(y)$$$ be the product of the distinct prime divisors of $$$y$$$. The row for $$$x$$$ can change only at positions where $$$r(y)\mid x$$$. Generate these positions by fixing $$$y$$$ and visiting the smaller multiples of $$$r(y)$$$.

For fixed $$$x$$$, separate the number of subsequences with endpoints $$$(x,y)$$$ into an $$$x$$$-only factor and a $$$y$$$-only factor. Prefix sums of the $$$y$$$-only factors can add a whole constant part of the row at once.

The easy solution builds the full DP table and then loops over every minimum and maximum. We need to compress both quadratic parts.

Let $$$r(y)$$$ be the product of the distinct prime divisors of $$$y$$$. The transition from E1 is

Fix $$$x$$$ and read the row from left to right. Unless $$$r(y)$$$ divides $$$x$$$, the value is copied from the previous position. Thus, the row can change only at positions satisfying $$$r(y)\mid x$$$; call them blockers.

For example, the blockers of row $$$2$$$ are $$$4,8,16,\ldots$$$. At $$$y=4$$$, the row changes from $$$2$$$ to $$$dp[1][3]=1$$$. Every later blocker also reads $$$1$$$ from row $$$1$$$, so the whole row is stored as only $$$(2,2)$$$ and $$$(4,1)$$$.

We can generate all blockers without checking every pair. For each $$$y$$$, visit

and add $$$y$$$ to the blocker list of each visited $$$x$$$. Processing $$$y$$$ in increasing order keeps every list sorted.

Now build the rows in increasing order of $$$x$$$. Store a row as breakpoints $$$(s,z)$$$: its value becomes $$$z$$$ starting at position $$$s$$$. Begin row $$$x$$$ with $$$(x,x)$$$. At a blocker $$$y$$$, advance a pointer in row $$$x-1$$$ while the next breakpoint starts at or before $$$y-1$$$. Its current value is $$$dp[x-1][y-1]$$$. Add a new breakpoint only when this value differs from the previous one.

For the fixed limit $$$M=300\,000$$$, generating the blocker lists visits $$$5\,970\,168$$$ pairs. Unfortunately, a bound on the number of generated pairs is not easy to calculate for general $$$n$$$, so the best way to solve the problem in-contest is to experimentally verify that the number of pairs is small enough for the given $$$n$$$. After equal neighboring parts are merged, all rows together contain $$$640\,704$$$ breakpoints, and no row contains more than eight.

We still need to avoid iterating over every pair of endpoints. Let $$$c_v$$$ be the frequency of $$$v$$$, and write

For $$$x \lt y$$$, the number of subsequences with minimum $$$x$$$ and maximum $$$y$$$ is

Call the two factors $$$L_x$$$ and $$$R_y$$$. Build prefix sums of $$$R_y$$$. Intersect each compressed part with $$$x \lt y\le n$$$. If the remaining interval is $$$u\le y\le v$$$ and the row value there is $$$z$$$, its whole contribution is

which is one prefix-sum query. Skip values with $$$c_x=0$$$. The case $$$x=y$$$ contributes $$$xC_x$$$ separately.

The radical sieve takes $$$O(M\log\log M)$$$ time. The remaining preprocessing is linear in the $$$5\,970\,168$$$ generated pairs and $$$640\,704$$$ stored breakpoints. Afterward, each present minimum scans at most eight parts, so the additional work is virtually linear in the total input size. The blocker lists and compressed rows have the sizes given above; together with the $$$O(M)$$$ arrays, they fit comfortably in memory.

#include <bits/stdc++.h>

using namespace std;

const int MAXN = 300000;
const int MOD = 998244353;

struct Breakpoint {
    int pos;
    int value;
};

array<int, MAXN + 1> freq, prefixCount, prefix, p2, ip2;
vector<vector<Breakpoint>> rows;

void buildRows() {
    vector<int> rad(MAXN + 1, 1);
    for (int p = 2; p <= MAXN; ++p) {
        if (rad[p] != 1) {
            continue;
        }
        for (int x = p; x <= MAXN; x += p) {
            rad[x] *= p;
        }
    }

    vector<vector<int>> blockers(MAXN + 1);
    for (int y = 2; y <= MAXN; ++y) {
        for (int x = rad[y]; x < y; x += rad[y]) {
            blockers[x].push_back(y);
        }
    }

    rows.resize(MAXN + 1);
    rows[0].push_back({0, 0});
    for (int x = 1; x <= MAXN; ++x) {
        auto &row = rows[x];
        auto &previous = rows[x - 1];
        row.push_back({x, x});

        int ptr = 0;
        for (int y : blockers[x]) {
            while (ptr + 1 < (int)previous.size() &&
                   previous[ptr + 1].pos < y) {
                ++ptr;
            }
            int value = previous[ptr].value;
            if (value != row.back().value) {
                row.push_back({y, value});
            }
        }
    }
}

void buildPowers() {
    const int inv2 = (MOD + 1) / 2;
    p2[0] = ip2[0] = 1;
    for (int i = 1; i <= MAXN; ++i) {
        p2[i] = 2LL * p2[i - 1] % MOD;
        ip2[i] = 1LL * inv2 * ip2[i - 1] % MOD;
    }
}

int rangeSum(int left, int right) {
    if (left > right) {
        return 0;
    }
    return (prefix[right] - prefix[left - 1] + MOD) % MOD;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    buildRows();
    buildPowers();

    int t;
    cin >> t;
    while (t--) {
        int n;
        cin >> n;
        vector<int> values;
        for (int i = 0; i < n; ++i) {
            int x;
            cin >> x;
            if (freq[x] == 0) {
                values.push_back(x);
            }
            ++freq[x];
        }

        int before = 0, ans = 0;
        prefixCount[0] = prefix[0] = 0;
        for (int x = 1; x <= n; ++x) {
            int ways = p2[freq[x]] - 1;
            int rightWeight = 1LL * ways * p2[before] % MOD;
            prefix[x] = prefix[x - 1] + rightWeight;
            if (prefix[x] >= MOD) {
                prefix[x] -= MOD;
            }
            ans = (ans + 1LL * x * ways) % MOD;
            before += freq[x];
            prefixCount[x] = before;
        }

        for (int x : values) {
            int ways = p2[freq[x]] - 1;
            int leftWeight = 1LL * ways * ip2[prefixCount[x]] % MOD;

            int rowSum = 0;
            auto &row = rows[x];
            for (int i = 0; i < (int)row.size(); ++i) {
                int left = max(x + 1, row[i].pos);
                if (left > n) {
                    break;
                }
                int right = n;
                if (i + 1 < (int)row.size()) {
                    right = min(n, row[i + 1].pos - 1);
                }
                rowSum = (rowSum +
                          1LL * row[i].value * rangeSum(left, right)) % MOD;
            }
            ans = (ans + 1LL * leftWeight * rowSum) % MOD;
        }

        cout << ans << '\n';
        for (int x : values) {
            freq[x] = 0;
        }
    }
}

Try to make an expression for $$$n!$$$ as a warm-up

Try to use round and division to make an expression which is $$$0$$$ or $$$1$$$ at specific $$$n$$$.

Look for a formula relating $$$n!$$$ and the number of derangements.

Note that $$$10\,000$$$ looks like a lot of characters for an expression, but for $$$n$$$ up to $$$50$$$, due to the particular restrictions of the Deranged Calculator language, this limit is quite tight. So for everything we come up with, we have to be sure we are not wasting too many characters. The model solution uses $$$5969$$$ characters, and there are even better methods that achieve (way) fewer characters.

Let's first try to solve an easier problem: what if, instead of the number of derangements of $$$n$$$ elements, we try to calculate the number of permutations, $$$n!$$$, instead?

We can immediately notice that if we want to use integer constants in our program, we have to make them using workarounds such as n/n or (n+n+n)/n. To save some length, we can precompute a dynamic programming table, which tells us a small construction for the constant $$$i$$$ for $$$1\leq i\leq50$$$. Transitions can be coded to use + and * operators, considering all ways to make $$$i=x+y$$$ or $$$i=x\cdot y$$$, with some appropriate use of parentheses and /n. It turns out that if the the other parts of the expression are well optimized, even writing constants naively in the form $$$(n+n+\dots+n)/n$$$ can work. In less optimized expressions, optimizing expressions for constants is maybe required

Next up, we have to face the issue that we cannot use recursion or loops to calculate the factorial. The one saving grace is that the round function can help us to make some kind of step functions so we can filter out specific $$$n$$$ for our needs.

It would be great if we could make the following program:

(1*[n>=2]+1)*(2*[n>=3]+1)*...*(49*[n>=50]+1)

In such a program, we basically compute $$$50!$$$, but we modulate which terms we want to multiply and which we replace with $$$1$$$. Here [n>=2] is an indicator expression that should output $$$1$$$ if the given inequality is true and $$$0$$$ otherwise.

So now our goal has shifted to building an indicator expression [n>=x] for a given constant $$$x$$$. Here lots of tricks can be used involving the round function. One particularly nice one is the following: round(n/(n+x)). For a given constant $$$x$$$ (for which we have a short expression precomputed), this function gives exactly $$$1$$$ once $$$n\geq x$$$, and $$$0$$$ for $$$2\leq n \lt x$$$. So this solves calculating factorials in the DC language in a reasonable number of characters (this depends on how optimized your round-based "if statement" is, and how well optimized the DP is).

But what about derangements? It turns out there are some nice formulas relating factorials and derangements. One particular fact is that

where $$$e\approx2.71828$$$ is Euler's number, and $$$D_n$$$ is the number of derangements for this $$$n$$$. This means we can basically reuse all our work for the factorial case, but now we have to find a way to get a sufficiently precise approximation of $$$e$$$ in our program. There are again many possibilities, but as we already know how to create all integer constants up to $$$50$$$ with a small number of characters, one way is to use the Taylor expansion of $$$e^x$$$, evaluated at $$$x=1$$$:

We can cut off the series at some point, and as we are dividing $$$50!$$$ by $$$e$$$, it makes sense that we need roughly the first $$$50$$$ terms of the expansion, as we need a relative error of $$$\sim \frac{1}{50!}$$$. Summing all the terms from $$$k=0$$$ to $$$k=50$$$ can be confirmed to work. Naively writing a sum of the inverse factorials will use $$$\sim50^2/2$$$ constants, and this will probably not fit. Instead, we can rewrite the sum as

making use of the common terms in the factorials. This only uses $$$\sim50$$$ constants, and some $$$1$$$'s, which are very cheap.

Taking this all together, we can write our final solution as round(n!/e), with our expressions for $$$n!$$$ and $$$e$$$ substituted.

The model solution tries to do some more tricks to optimize the precomputation DP for the constants, but they are not needed:

• Use a bigger DP table $$$DP_{i,j}$$$, for $$$0\leq j \lt 4$$$, which calculates a small-length expression for getting an expression with value $$$i\cdot n^j$$$, and also consider dividing using /n as transitions.
• Think about when exactly you can skip parentheses.
• Try to also use the - operator in your DP. Due to both - and / transitions, the DP is no longer acyclic, but as the table is very small, we can just keep iterating through all possible transitions multiple times, and eventually the DP table stays unchanged.

As a bonus, what is the lowest you can push the number of characters?

N = 50
L = 4
dp = [[None]*L for i in range(N+1)]
def upd(l,n,s):
    if not dp[n][l] or len(s)<len(dp[n][l]): dp[n][l]=s
for i in range(1,L):
    upd(i,1,'n' + (i-1)*'*n')

def parens(s):
    if len(s)>1: return '(' + s + ')'
    return s

for iter in range(1,N):
    for l1 in range(0,L):
        for l2 in range(0,L):
            for i in range(1,N+1):
                for j in range(1,N+1):
                    if dp[i][l1] and dp[j][l2]:
                        a = dp[i][l1]
                        b = dp[j][l2]
                        if i<j and l1==l2: upd(l1,j-i,b + '-' + parens(a))
                        if i+j<=N and l1==l2: upd(l1,i+j,f"{a}+{b}")
                        if i*j<=N and l1+l2<L: upd(l1+l2,i*j,f"{parens(a)}*{parens(b)}")
                        if i%j==0 and l1-l2>=0: upd(l1-l2, i//j, f"{parens(a)}/{parens(b)}")

e = "n/n"
for i in range(N,0,-1): e = f"n/n+n/({dp[i][1]})*({e})"

def geqn(i): return f"round(n/({dp[i][0]}+n))"

def factorial():
    res = "n/n"
    for j in range(2,N+1):
        res += f"*({geqn(j)}*({dp[j-1][0]})+n/n)"
    return res
fact = factorial()
ans = f"round({fact}/({e}))"
print(ans)

def n(x):
    if x==1: return "n/n"
    else: return '(' + "+".join("n"*x) + ')/n'
N, e, fact = 50, n(1),n(1)
def geqn(i): return f"round(n/({n(i)}+n))"
for j in range(2,N+1): fact += f"*({geqn(j)}*({n(j-1)})+n/n)"
for i in range(N,0,-1): e = f"n/n+n/({'+'.join('n'*i)})*({e})"
print(f"round({fact}/({e}))")