You’re running into the right idea (use the GF)

but the way you add each $p$-block is $O(M^2)$ overall (and your update of $Q_p$ is also $O(M)$ per $p$, so $O(M^2)$ total). With $M\le 3\cdot10^5$ that will time out.

Below are two optimizations, in order of impact.

1) Do it in $\tilde O(M)$ multiplications via a divide-&-conquer (segment-tree) of polynomials

Key observation: define the (truncated) polynomials

• $H_k(x)= (1-x^k)^{-2} = \sum_{c\ge 0}(c+1)\,x^{ck}$ (cut off at degree $M$),
• $G_p(x)=\dfrac{x^p}{1-x^p} = \sum_{c\ge 1}x^{cp}$ (cut off at degree $M$).

Then

Instead of forming each prefix product and convolving it with $G_p$ one-by-one, build a segment tree over $p\in[1,M]$ and compute everything with $O(M)$ polynomial multiplications (each is NTT/CRT-based, degrees trimmed to $\le M$):

• Build a product tree $\mathcal P$: each node stores $P_{[L,R]}(x)=\prod_{v=L}^R H_v(x)$.
• Build a “weighted sum” tree $\mathcal D$ bottom-up with the identity $$D_{[L,R]}(x) \;=\; D_{[L,mid]}(x)\;+\;P_{[L,mid]}(x)\cdot D_{[mid+1,R]}(x),$$ with leaf $D_{[p,p]}(x)=G_p(x)$.
• Then $D_{[1,M]}(x)=\sum_{p=1}^{M}\Big(\prod_{v=1}^{p-1}H_v(x)\Big)\,G_p(x)=F(x)$.

Every internal node performs one polynomial multiplication for the product tree and one for the D-tree, both truncated to degree $M$. With an NTT this is $O(M\log^2 M)$ total and is easily fast enough.

Implementation tips:

• Leaves:
  • $H_k$ is very sparse: set poly[0]+=1, poly[k]+=2, poly[2k]+=3, … until degree $>M$.
  • $G_p$ is a simple 0/1 periodic vector: set poly[cp]=1 for cp in {p,2p,… ≤ M}.
• Always truncate to degree $M$ after each multiply.
• Use NTT(s). Since the modulus is $10^9+7$ (not NTT-friendly), do MTT/CRT:
  • run NTT under two or three friendly primes (e.g., 998244353, 1004535809, 469762049),
  • Garner / CRT the result back mod $10^9+7$.
• For very small degrees, switch to a naive $O(nm)$ multiply to reduce constants.
• Pre-allocate buffers; reuse vectors to avoid allocations.

This completely eliminates your $O(M^2)$ loops; you never form $Q_p$ explicitly and you never add “multiples of $p$” with nested loops—the tree does the algebra for all $p$ at once.

Sketch (only the structure; drop-in NTT is standard and omitted for brevity):

cpp
// Build leaves H_k and G_p up to degree M.
// Product tree: P[node] = P[left]*P[right]  (truncate to deg M)
// D-tree: D[node] = D[left] + P[left] * D[right]  (truncate to deg M)
// Answer poly = D[root]; print coeffs 1..M.

struct Poly { vector<int> a; }; // coeffs mod 1e9+7 via CRT in mul()

Poly mul(const Poly& A, const Poly& B, int M); // NTT+CRT, truncate to M
void add_inplace(Poly& A, const Poly& B);      // (mod 1e9+7), resize as needed
void trim(Poly& A, int M);

vector<Poly> Hleaf(M+1), Gleaf(M+1);

void build_H_G(int M) {
    for (int k = 1; k <= M; ++k) {
        auto& H = Hleaf[k].a; H.assign(M+1, 0);
        for (int t = 0, deg = 0; deg <= M; ++t, deg = t*k) {
            if (deg > M) break;
            H[deg] = (H[deg] + (t+1)) % MOD;  // (1 - x^k)^{-2}
        }
        auto& G = Gleaf[k].a; G.assign(M+1, 0);
        for (int deg = k; deg <= M; deg += k) G[deg] = 1; // x^k/(1-x^k)
    }
}

struct Node { Poly P, D; };
vector<Node> seg(4*(M+5));

void build(int idx, int L, int R, int M) {
    if (L == R) {
        seg[idx].P.a = Hleaf[L].a;     trim(seg[idx].P, M);
        seg[idx].D.a = Gleaf[L].a;     trim(seg[idx].D, M);
        return;
    }
    int mid = (L+R)>>1;
    build(idx<<1, L, mid, M);
    build(idx<<1|1, mid+1, R, M);
    seg[idx].P = mul(seg[idx<<1].P, seg[idx<<1|1].P, M);
    Poly rightWeighted = mul(seg[idx<<1].P, seg[idx<<1|1].D, M);
    seg[idx].D = seg[idx<<1].D;
    add_inplace(seg[idx].D, rightWeighted);
    trim(seg[idx].D, M);
}

After build(1,1,M,M), seg[1].D.a[n] (for 1 ≤ n ≤ M) are your answers.

2) If you must keep your current shape, at least drop a log-factor (still too slow at 3e5 though)

Right now you do:

cpp
for (int m = p; m <= M; m += p)
  for (int N = m; N <= M; ++N)
    f[N] += Q[N - m];

Rewrite that using residue-class prefix sums (for each r in [0..p-1]):

• let a[i] = Q[r + i*p] for i ≥ 0,
• prefix pref[i] = pref[i-1] + a[i],
• then for j ≥ 1, add pref[j-1] to f[r + j*p].

This makes the B_p update $O(M)$ (down from $O(M\log M)$ per $p$), but the whole algorithm remains $O(M^2)$ and will still TLE at the given limit. Use it only as an incremental improvement if you can’t implement #1 yet.

Complexity & practicality

• Segment-tree D&C + NTT/CRT: ~$O(M \log^2 M)$ polynomial time with small constants; passes comfortably for $M=3\cdot10^5$.
• Memory: store a few hundred MB at most if you aggressively trim to degree $M$ and reuse buffers.

If you want, I can share a full drop-in NTT+CRT template wired to the tree above; it’s a bit long for a single message, but that’s the exact code pattern to make this pass.