Here’s a straightforward Java implementation that runs in $O(n^2)$ time and $O(n)$ memory. It uses the combinatorial interpretation

$$W_n \;=\;\sum_{\substack{f:\,S\hookrightarrow\{1,\dots,n\}\\f(i)\le i}} \prod_{i\in S} f(i),$$

where the sum runs over all partial injections $f$ whose domain $S\subseteq\{1,\dots,n\}$ satisfies $f(i)\le i$. We build this up for $i=1,2,\dots,n$ by maintaining:

• g[i]: the total “weight” $W_i$ over partial injections on $\{1,\dots,i\}$;
• w[t]: the total weight of those injections among these that ever use the target value $t$.

Then at step $i\!+\!1$, every old injection either

1. stays an injection of size $i$ (no new map at $i+1$),
2. gains a new map $f(i+1)=t$ for some $t\in[1..i+1]$ not already used

and each such extension multiplies its old weight by $t$. Putting that together you get

java
// g0 = 1, w[*]=0
g[i] = g[i-1]
       +  ∑_{t=1..i}  t * (g[i-1] - w[t])
       (all mod m)

and then

w[t] = w[t]                       // keep all old injections that already used t
       + t * (g[i-1] - w[t])      // plus all the new ones that pick t at step i

The total is then g[n] % m.

That runs in about $\sum_{i=1}^n i=O(n^2)$, which will handle up to a few thousand comfortably. If you need to go higher (e.g.\ $n$ up to $10^5$), you’ll need to speed up the prefix‐sum/point‐updates—for example with a Fenwick tree / BIT or segment tree—so that both the “sum of $t\cdot w[t]$” and the “point‐update on $w[t]$” become $O(\log n)$ instead of $O(n)$.

java
import java.io.*;
import java.util.*;

public class Main {
    public static void main(String[] args) throws IOException {
        // Fast input
        BufferedReader br = new BufferedReader(new InputStreamReader(System.in));
        StringTokenizer st = new StringTokenizer(br.readLine());
        int n = Integer.parseInt(st.nextToken());
        int m = Integer.parseInt(st.nextToken());
        
        // g[i] = total weight for injections on {1..i}
        // w[t] = total weight of those injections that ever use target t
        long[] g = new long[n+1];
        long[] w = new long[n+1];
        
        g[0] = 1;  // empty injection has weight 1
        // w[*] start at 0
        
        // We will need to compute sum_{t=1..i} t * (g[i-1] - w[t]) mod m
        // To do that in O(n^2) we just do the loop directly.
        for (int i = 1; i <= n; i++) {
            long prev = g[i-1];
            
            // 1) compute S = sum_{t=1..i} t * (prev - w[t])
            long S = 0;
            for (int t = 1; t <= i; t++) {
                long diff = (prev - w[t]) % m;
                if (diff < 0) diff += m;
                S = (S + t * diff) % m;
            }
            
            // 2) add the “no‐new‐map” case
            g[i] = (prev + S) % m;
            
            // 3) update w[t] for t=1..i
            //    old w[t] stays (those injections already used t, and pick nothing new)
            //    plus t * (prev - old w[t]) for those that pick t now
            for (int t = 1; t <= i; t++) {
                long old = w[t];
                long diff = (prev - old) % m;
                if (diff < 0) diff += m;
                w[t] = (old + t * diff) % m;
            }
        }
        
        // answer is g[n] mod m
        System.out.println(g[n] % m);
    }
}

---

How it works

1. g[i] is exactly

   $$\sum_{\substack{f:\,S\hookrightarrow\{1,\dots,i\}\\f(k)\le k}} \prod_{k\in S}f(k),$$

   the total “weight” of all partial injections on the first $i$ positions.

2. w[t] accumulates the same sum but only over those injections that use the value $t$ somewhere in their domain.

3. When we go from $i-1$ to $i$:

• No new mapping at position $i$ just carries over all the old weight g[i-1].
• Mapping to $t$ for $1\le t\le i$ can be done exactly from those old injections that haven’t used $t$ yet (that’s why we look at g[i-1] - w[t]); each such extension multiplies the weight by $t$.
  Adding up over $t$ gives the second term

$$\sum_{t=1}^i t\;\bigl(g[i-1]\;-\;w[t]\bigr).$$

4. Finally, we update w[t] by keeping all the old injections that already used $t$ (they still do in the new step) plus those new ones that pick $t$ now (again t * (g[i-1]-w[t])).

Total time is

$$\sum_{i=1}^n\;O(i)\;=\;O(n^2),$$

which is fine up to a few thousand. If you need $n$ up to, say, $10^5$, you can keep a Fenwick Tree (BIT) over the array $w[t]$ to maintain

• the prefix sums $\sum t\,w[t]$, and
• point‐updates on individual $w[t]$,

each in $O(\log n)$ instead of $O(n)$, bringing the total to $O(n\log n)$.

---

Note: Because there are $(n+1)!$ valid sequences, any brute‐force over all of them is completely infeasible for even moderate $n$; this “injection‐DP” is the standard way to collapse the double sum (over sequences $\mathbf a$ and the removal‐choices) down to a single $O(n^2)$ pass.