We can turn this into a standard digit DP problem. The problem is how we should capture the transition of the No. of the scary substrings. Suppose we are in the $$$i$$$-th (from the beginning) digit of number, we can keep track of "The parity of No. of scary substrings", "The parity of No. of suffixes of $$$s[1,i]$$$ with odd sum", "The parity of No. of suffixes of $$$s[1,i]$$$ with even sum". Then it will not be hard to find the transition.

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
const int B = (int)1e5 + 50, mod = 998244353;

ll dp[B][2][2][2][2];
int bit[B], b;
int pw2[B];
//int A, B;

ll get(int d, int scary, int odd_suf, int even_suf, bool flag, bool st){
    if(d == -1) {
        return scary % 2 == 1;
    }
    if(!flag && ~dp[d][scary][odd_suf][even_suf][st]) return dp[d][scary][odd_suf][even_suf][st];
    int lim = flag ? bit[d] : 9;
    int res = 0;
    for(int i = 0; i <= lim; i++){
        bool nxt_st = st || i != 0;
        int nxt_scary = (scary + (i % 2 == 0 ? even_suf + 1 : odd_suf)) % 2;
        int nxt_odd = (i % 2 == 0 ? odd_suf : even_suf + 1) % 2;
        int nxt_even = (i % 2 == 0 ? even_suf + 1 : odd_suf) % 2;
        if(!nxt_st) {nxt_scary = nxt_odd = nxt_even = 0;}
        res += get(d - 1, nxt_scary, nxt_odd, nxt_even, flag && lim == i, nxt_st);
        res %= mod;
    }
    return flag ? res : dp[d][scary][odd_suf][even_suf][st] = res;
}

ll solve(string str){
    b = 0;
    memset(dp, -1, sizeof(dp));
    for(int i = str.length() - 1; i >= 0; i--) {
        bit[b++] = str[i] - '0';
    }
    return get(b-1, 0, 0, 0, true, false);
}

int main() {
    string L, R;
    cin >> L >> R;
    L.back() -= 1;
    int carry = 0;
    for(int i = L.size() - 1; i >= 0; i--) {
        if(carry) {
            L[i]--;
            carry = 0;
        }
        if(L[i] < '0') {
            L[i] += 10;
            carry++;
        }
    }
    ll res = (solve(R) - solve(L) + mod) % mod;
    cout << res << endl;
}

Try to build the graph for each different query. We can see the sum of the number of nodes in the graph will not exceed $$$O(n \cdot \sqrt{\max f_i})$$$, since each node will appear in at most $$$2 \sqrt{\max f_i}$$$ (the number of its divisors). One can also see the number of edges are bounded. Then for each graph, we run dfs and count the answer.

#include <bits/stdc++.h>
using namespace std;

const int N = (int)1e5 + 50, Q = (int)1e6 + 50;

int n,m,q;
vector<int> G[N];

bool vis[N];
bool gd[N];
int f[N];
vector<int> dvs[Q];
int res[Q];

void dfs(int v) {
    vis[v] = true;
    for(int nxt : G[v]) {
        if(!vis[nxt] && gd[nxt]) dfs(nxt);
    }
}

int solve(int dv) {
    for(int i : dvs[dv]) {
        vis[i] = false;
        gd[i] = true;
    }

    int res = 0;
    for(int i : dvs[dv]) {
        if(!vis[i]) {
            dfs(i);
            res++;
        }
    }
    for(int i : dvs[dv]) {
        gd[i] = false;
    }
    return res;
}

int main() {
    memset(res, -1, sizeof(res));
    cin >> n >> m >> q;
    for(int i = 0; i < n; i++) {
        cin >> f[i];
        int j;
        for(j = 1; j * j < f[i]; j++) {
            if(f[i] % j == 0) {
                dvs[j].push_back(i);
                dvs[f[i]/j].push_back(i);
            }
        }
        if(j * j == f[i]) dvs[j].push_back(i);
    }
    for(int i = 0; i < m; i++) {
        int a, b; cin >> a >> b; a--, b--;
        G[a].push_back(b);
        G[b].push_back(a);
    }
    while(q--) {
        int x;
        cin >> x;
        if(res[x] != -1) {
            cout << res[x] << "\n";
            continue;
        }
        res[x] = solve(x);
        cout << res[x] << "\n";
    }

}

For each pair of neighboring layers, we try to compute the answer by inclusion-exclusion principle and sum them together. Let $$$f(a, b)$$$ be the number ways to connect two layers of size $$$a$$$ and $$$b$$$ respectively. Then we have

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
const int N = (int)5e5 + 50;
const ll mod = 998244353;

int n, num[N];
ll fac[N], facinv[N];
ll inv[N];
ll pw2[N];

ll fp(ll x, ll k){
    if(k == 0) return 1;
    ll hf = fp(x, k/2);
    return k % 2 ? hf * hf % mod * x % mod: hf * hf % mod;
}

ll comb(int n, int k){
    return fac[n] * facinv[k] % mod * facinv[n - k] % mod;
}

ll calc(int a, int b) {
    if(a < b) swap(a, b);

    ll res = 0;
    int sig = 1;
    for(int i = b; i >= 0; i--) {
        res += sig * fp(pw2[i]-1, a) * comb(b, i) % mod;
        sig = -sig;
        res %= mod;
    }
    return res;
}

int main(){
    inv[1] = 1;
    for(int i = 2; i < N; i++) inv[i] = (mod - (mod / i) * inv[mod % i] % mod) % mod;
    fac[0] = 1;
    for(int i = 1; i <= N-1; i++) fac[i] = fac[i-1] * i % mod;
    facinv[N-1] = fp(fac[N-1], mod - 2);
    for(int i = N-1 - 1; i >= 0; i--) facinv[i] = facinv[i+1] * (i+1) % mod;
    pw2[0] = 1;
    for(int i = 1; i < N; i++) pw2[i] = pw2[i-1] * 2 % mod;


    cin >> n;
    for(int i = 0; i < n; i++) cin >> num[i];
    ll res = 1;
    for(int i = 0; i < n - 1; i++) {
        res *= calc(num[i], num[i+1]);
        res %= mod;
    }
    cout << (res + mod) % mod << endl;
}