Suppose that \(N\) has \(B\) digits.

We consider the contribution of an \(b\)-digit integer \(m\) to the answer. In other words, we consider the answer where the string appended as \(T\) is replaced by repeated 0 of the same length unless it is \(m\).

First of all, if there are \(i_k\) integers with \(k\) digits after \(m\) in \(P\), then \(m\) contributes to \(f(P)\) by \(m10^{\sum_k ki_k}\). Thus, denoting by \(C_k\) the number of \(k\)-digit integers between \(1\) and \(N\) except for \(m\),

• for each \(k=1,2,\ldots,B\), there are \(\displaystyle \binom{C_k}{i_k}\) ways to choose \(i_k\) integers with \(k\) digits from the \(C_k\) choices;
• there are \(\displaystyle \left(\sum_k i_k\right)!\) ways to arrange the \(\displaystyle \sum_k i_k\) integers after \(m\);
• there are \(\displaystyle \left(N-1-\sum_ki_k\right)!\) ways to arrange the \(\displaystyle N-1-\sum_k i_k\) integers before \(m\);

so \(m\) contributes to the answer by

\[\sum_{0\le i_k\le C_k\ (1\le k\le B)} m10^{\sum_k ki_k}\left(\prod_k \binom{C_k}{i_k}\right)\left(\sum_k i_k\right)!\left(N-1-\sum_ki_k\right)!.\]

Here, \(\displaystyle \sum_{0\le i_k\le C_k\ (1\le k\le B)}\) represents the sum of the function inside the sigma over all the length-\(B\) sequences \(i\) such that \(0\le i_k\le C_k\) for all \(k=1,2,\ldots,B\).

Moreover, denoting by \(\displaystyle F_k(x)=\sum_{i=0}^{C_k}\binom{C_k}{i}10^{ki}x^i=(1+10^k x)^{C_k}\), the contribution of \(m\) to the answer can be written as \(\displaystyle \sum_i f_ix^i=\prod_kF_k(x)\), where \(\displaystyle \sum_i f_ix^i=\prod_kF_k(x)\). Moreover, denoting by \(C_k^0\) the number of \(k\)-digit integers between \(1\) and \(N\), and letting \(\displaystyle F_k^0(x)=\sum_{i=0}^{C_k^0}\binom{C_k^0}{i}10^{ki}x^i=(1+10^k x)^{C_k^0}\), we have \(\displaystyle \sum_i f_ix^i=\frac{F^0(x)}{1+10^b x}\) where \(\displaystyle F^0(x)=\prod_k F_k^0(x)\).

Expanding \(F^0\) requires \((B-1)\) convolutions, but the degree of \(F_k\) is \(C_k^0=O(10^k)\), so it can be done in \(\displaystyle O\left(\sum_{k=2}^B 10^k\log10^k\right)=O(B10^B)\) time by repeatedly picking two polynomials with the smallest degrees and convolving them.

Moreover, \(F^0\) is independent of \(b\). Thus, we can first calculate \(F^0\), and then compute \(\displaystyle \frac{F^0(x)}{1+10^b x}\) for each \(b\), the number of digits in \(m\). This polynomial division can be processed in \(O(N)\) time for each \(b\) in the same manner as the column division.

Since \(m\) has at most \(B\) digits, where \(B=O(\log N)\), the problem can be solved in a total of \(O(N\log N)\) time by implementing the algorithm above appropriately

Sample code (C++)

#include <bits/stdc++.h>
#include <atcoder/all>
using namespace std;
using mint = atcoder::modint998244353;
int main()
{
    int n;
    cin >> n;
    int B = to_string(n).size();
    vector<mint> fact(n + 1), finv(n + 1);
    fact[0] = 1;
    for (int i = 1; i <= n; i++)
        fact[i] = i * fact[i - 1];
    finv[n] = fact[n].inv();
    for (int i = n; i > 0; i--)
        finv[i - 1] = i * finv[i];
    auto comb = [&](int n, int k) -> mint
    { return fact[n] * finv[k] * finv[n - k]; };
    vector<int> ten(B + 1);
    ten[0] = 1;
    for (int i = 1; i <= B; i++)
        ten[i] = 10 * ten[i - 1];
    vector<int> c(B + 1);
    for (int i = 1; i <= n; i++)
        c[to_string(i).size()]++;
    vector<mint> F0 = {1};
    for (int k = 1; k <= B; k++)
    {
        vector<mint> new_f(c[k] + 1);
        mint ten_c = 1;
        for (int i = 0; i < c[k] + 1; i++)
        {
            new_f[i] = comb(c[k], i) * ten_c;
            ten_c *= ten[k];
        }
        F0 = atcoder::convolution(F0, new_f);
    }
    F0.push_back(0);
    mint ans = 0;
    for (int b = 1; b <= B; b++)
    {
        auto f = F0;
        for (int i = 0; i < n; i++)
            f[i + 1] -= ten[b] * f[i];
        mint res = 0;
        for (int i = 0; i < n; i++)
            res += f[i] * fact[i] * fact[n - 1 - i];
        long long le = ten[b - 1], ri = min(ten[b], n + 1);
        ans += res * (ri * (ri - 1) / 2 - le * (le - 1) / 2);
    }
    cout << ans.val() << endl;
}