All expressions containing \(x\) below is formal power series.

Pentagonal Number Theorem

The pentagonal number theorem claims:

\[\displaystyle \prod_{i=1}^{\infty} (1-x^i) = \sum_{n=-\infty}^{\infty} (-1)^nx^{\frac{n(3n-1)}{2}}.\]

We do not introduce the proof here. If you are interested, see Wikipedia ( https://en.wikipedia.org/wiki/Pentagonal_number_theorem ).

Solution of this problem

The sought value is expressed by:

\[[x^M]\ \prod_{i=1}^{N} (x+2x^2+3x^3+\dots+ix^i+ix^{i+1}+ix^{i+2}+\dots).\]

We can transform the factor to obtain

\[[x^M]\ \prod_{i=1}^{N} \frac{x(1-x^i)}{(1-x)^2},\]

whose derivation is explained later.

Furthermore, we have

\[ \begin{aligned} [x^M]\ \prod_{i=1}^{N} \frac{x(1-x^i)}{(1-x)^2} &= [x^{M-N}]\ \prod_{i=1}^{N} \frac{(1-x^i)}{(1-x)^2}\\ &= [x^{M-N}]\ \frac{\prod_{i=1}^{\infty} (1-x^i)}{(1-x)^{2N}}, \end{aligned} \]

where we used \(M-N \leq N\) in the last equality.

Now the expression on the left hand side of the pentagonal number theorem appears as is in the numerator on the right hand side of the expression above. Thus, there are only \(O(\sqrt{N})\) terms with degrees less than or equal to \({M-N}\) in the expansion of the infinite product on the left hand side in the numerator. Denoting them by \(a_1 x^{e_1},a_2 x^{e_2},\dots,a_k x^{e_k}\), we have

\[ \begin{aligned} [x^{M-N}]\ \frac{\prod_{i=1}^{\infty} (1-x^i)}{(1-x)^{2N}} &=[x^{M-N}]\ \frac{\sum_{i=1}^{k} a_kx^{e_k}}{(1-x)^{2N}}\\ &=\sum_{i=1}^{k} a_k \binom{M+N-{e_k}-1}{2N-1}. \end{aligned} \]

Therefore, we can solve this problem in a total of \(O(\sqrt{N}\log_{\mod}{N} +\mod)\) time using Lucas’ theorem to evaluate the binomial coefficients.

The Derivation that is “explained later”

A formal power series of this form can be often transformed into a simple fractional form by considering the progression of differences. Finding the differences is equivalent to multiplying \((1-x)\) to the formal power series. Indeed, if we let

\[f_i(x) = (x+2x^2+3x^3+\dots+ix^i+ix^{i+1}+ix^{i+2}+\dots),\]

then we have

\[f_i(x)(1-x) = (x+x^2+x^3+\dots+x^i),\]

\[f_i(x)(1-x)^2 = (x-x^{i+1}),\]

so

\[f_i(x) = \frac{x(1-x^i)}{(1-x)^2}.\]

Sample Code (C++):

#include<bits/stdc++.h>
#include<atcoder/all>

using namespace std;
using namespace atcoder;

using ll = long long;
using mint = modint;

const int mod = 200003;

mint fact[mod];
mint ifact[mod];

// binom(n, k) mod 200003
mint binom(ll n, ll k) {
    if (k < 0 or k > n) return 0;
    mint res = 1;
    while (n) {
        ll n0 = n % mod;
        ll k0 = k % mod;
        if (n0 < k0) return 0;
        res *= fact[n0] * ifact[k0] * ifact[n0 - k0];
        n /= mod;
        k /= mod;
    }
    return res;
}

int main() {
    mint::set_mod(mod);
    fact[0] = 1;
    for (int i = 1; i < mod; i++) fact[i] = fact[i - 1] * i;
    ifact[mod - 1] = fact[mod - 1].inv();
    for (int i = mod - 1; i >= 1; i--) ifact[i - 1] = ifact[i] * i;
    
    ll n, m;
    cin >> n >> m;
    
    mint ans = 0;
    for (ll i = 0; i * (3 * i + 1) / 2 <= m - n; i++) {
        ll j = i * (3 * i + 1) / 2;
        ans += binom(m + n - j - 1, 2 * n - 1) * (i % 2 ? -1 : 1);
    }
    for (ll i = 1; i * (3 * i - 1) / 2 <= m - n; i++) {
        ll j = i * (3 * i - 1) / 2;
        ans += binom(m + n - j - 1, 2 * n - 1) * (i % 2 ? -1 : 1);
    }
    cout << ans.val() << endl;
}