A set of squares painted black corresponds to how the piece advances, so it is sufficient to count the latter.

For now, let us ignore the computational complexity and focus on solving the problem. It is boiled down to DP (Dynamic Programming) represented by the following pseudocode:

dp={1,0,...,0};
for(i=1;i<=N;i++)
  for(j=i+A[i];j<=N;j+=A[i])
    dp[j]+=dp[i];
print(sum(dp));

When \(A_i\) are large, iteration on \(j\) is not repeated many times, so this approach seems reasonable, but if \(A_i=(1,1,\dots,1)\) for example, it costs \(O(N^2)\) time, so this algorithm is hard to pass.

Instead, let us consider how we can update the DP array fast when \(A_i\) are small.

For \(i<j\), a transition from \(dp[i]\) to \(dp[j]\) occurs if and only if \((j-i)\) is a multiple of \(A_i\) (i.e. the remainder when \(i\) is divided by \(A_i\) equal that of \(j\) by \(A_i\)).

Thus, if we define an array like \(dp2[d][r] = \) { remainder divided by \(d\) is \(r\) }, the problem seems solvable (which is indeed true).

So what does it mean by \(A_i\) are “small” or “large”?
Suppose that \(A_i \le b\) is the boundary of whether \(A_i\) is small or large.

• If you perform the operation for small \(A_i\), the time complexity for that part is \(O(Nb)\).
• If you perform the operation for large \(A_i\), the time complexity for that part is \(O(N \times \frac{N}{b})\).

To minimize \(O(Nb + N \times \frac{N}{b})\), it is optimal to let \(b=\sqrt{N}\) (short proof: inequality of arithmetic and geometric means).
When implementing, one can fix \(b\) to an arbitrary constant around \(\sqrt{N}\). Depending on constant factors, the optimal \(b\) might be far from \(\sqrt{N}\), so it might be good idea to try several in some cases.

The route to the correct implementation is left as an exercise for the readers. We show sample code after the editorial.

The technique of choosing a boundary at \(\sqrt{N}\) or splitting into \(\sqrt{N}\) packets to solve subproblem fast is called square root decomposition.

Sample code (C++):

#include<bits/stdc++.h>

using namespace std;

#define mod 998244353
#define BD 500

int dp[200005]={0};
int dp2[505][505]={0};

int main(){
  int n;
  cin >> n;
  dp[0]=1;
  for(int i=0;i<n;i++){
    int a;
    cin >> a;
    for(int d=1;d<=BD;d++){
      dp[i]+=dp2[d][i%d];
      if(dp[i]>=mod){dp[i]-=mod;}
    }

    if(a<=BD){
      dp2[a][i%a]+=dp[i];
      if(dp2[a][i%a]>=mod){dp2[a][i%a]-=mod;}
    }
    else{
      for(int j=i+a;j<n;j+=a){
        dp[j]+=dp[i];
        if(dp[j]>=mod){dp[j]-=mod;}
      }
    }
  }

  int res=0;
  for(int i=0;i<n;i++){
    res+=dp[i];
    if(res>=mod){res-=mod;}
  }
  cout << res << "\n";
  return 0;
}

Note: we avoid division by \(998244353\) using conditional branches.