• Đầu tiên, ta cùng tìm hiểu về công thức tính ~f(B)~ từ dãy con ~B~ bất kỳ.

Nhận xét 1: Điều kiện điểm ~(x, y)~ thuộc đường thẳng đi qua ~(B_i, 0)~ và ~(0, B_i)~ tương đương với ~x+y = B_i~

Nhận xét 2: Số lượng lattice point nằm trên đoạn thẳng nối ~(0,0)~ với ~(x, y)~ (~x, y~ nguyên dương) là ~gcd(x, y)~.

Nhận xét 3: Nếu cặp số nguyên dương ~(x, y)~ thỏa mãn ~x+y = B_i~ thì ~gcd(x, y) < B_i~ và ~gcd(x, y)~ phải là một ước của ~B_i~

Nhận xét 4 : Nếu gọi ~p_i~ là ước nguyên tố bé nhất của ~B_i~ thì có đúng ~(p_i-1)~ cách để chọn ra cặp số ~(x, y)~ sao cho ~x+y = B_i~ và ~gcd(x, y)~ là ước lớn thứ 2 của ~B_i~. Mặt khác, từ nhận xét 3 có thể kết luận rằng ~gcd(x, y)~ lớn nhất chính bằng ước lớn thứ 2 của ~B_i~.

Theo nhận xét 1 và 2, dễ thấy ~f(B)~ là số cách chọn ~m~ cặp số nguyên dương sao cho cặp thứ số ~i~ là ~(x_i, y_i)~ thỏa mãn ~x_i+y_i = B_i~ với ~gcd(x_i, y_i)~ lớn nhất có thể. Vậy từ nhận xét 4 ta rút ra được công thức ~f(B[1 \dots m]) = \prod_{i = 1}^{m} (p_i - 1)~.

• Sau đó, hãy tìm công thức tính ~s(A)~ với mảng ~A~ bất kỳ.

Theo định nghĩa thì ~s(A) = \sum_{B \in A, B \neq \emptyset} f(B)~ nhưng ta sẽ đơn giản hóa công thức này bằng toán tổ hợp. Với mỗi phần tử ~A_i~ có hai trường hợp xảy ra, một là ~A_i~ nằm ngoài dãy con ~B~ được chọn, hai là ~A_i~ nằm trong dãy con ~B~ và trường hợp này có ~(p_i-1)~ ~(p_i~ là ước nguyên tố bé nhất của ~A_i)~ cách chọn cặp số ~(x_i, y_i)~ sao cho ~x_i+y_i=A_i~ với ~gcd(x_i, y_i)~ lớn nhất. Vậy mỗi phần tử ~A_i~ có ~(p_i-1+1) = p_i~ cấu hình phù hợp, suy ra ~s(A) = -1 + \prod_{i = 1}^{n} p_i~ tương đương ~s(A)+1 = \prod_{i = 1}^{n} p_i~.

Chú ý : Ta cần cộng ~1~ vào vế trái để loại trừ trường hợp ~B~ là dãy con rỗng, ngoài ra lúc này vế phải của đẳng thức là tích các số nguyên tố nên bài toán đã đơn giản hơn rất nhiều.

• Gọi mảng ~A~ sau khi xóa là ~A'~, cho truy vấn ~[L, R]~, ta cần đếm số giá trị khác nhau của ~s(A')~ với ~A'~ là một dãy con của ~A[L\dots R]~.

Ta có: ~s(A[L \dots R]) + 1 = \prod_{i = L}^{R} p_i~.

Dễ thấy nếu ~A'~ là một dãy không rỗng của ~A[L\dots R]~ thì ~s(A')+1 ~ sẽ là một ước của ~s(A[L\dots R]) + 1~. Mặt khác, hiển nhiên tồn tại ít nhất một dãy con ~A'~ ứng với mỗi ước của ~s(A[L\dots R]) + 1~ hay ước của ~\prod_{I = L}^{R} p_i~. Vì vậy ta chỉ cần đếm số ước của ~\prod_{I = L}^{R} p_i~ là ra số giá trị phân biệt của ~s(A')~. Công việc này có thể thực hiện bằng cách đếm số lần xuất hiện của từng số nguyên tố trong ~\prod_{I = L}^{R} p_i~. Đáp án cuối cùng sẽ là ~-1 + \prod (cnt_i + 1)~ với ~cnt_i~ là số lần xuất hiện của số nguyên tố thứ ~i~ ở trong đoạn ~p[L\dots R]~ ~(-1~ là để loại đi trường ~A'~ là dãy con rỗng~)~.

Ta sẽ sử dụng thuật toán Mo để duy trì mảng đếm ~cnt_i~ (số lần xuất hiện của số nguyên tố thứ ~i~ trong đoạn đang xét). Do các số nguyên tố có thể rất lớn nên cần nén số các số nguyên tố lại.

• Tính ~p_i~ bằng cách sử dụng sàng nguyên tố trên đoạn (do ~max(A) - min(A) \le 10^7~). Ta có thể tăng tốc tính toán bằng cách dùng nhận xét ~A_i \le 10^{12}~ suy ra nếu ~p_i > \sqrt{10^{12}} = 10^6~ thì ~A_i = p_i~ (Vì ~p_i~ là ước nguyên tố bé nhất, nếu có ước nguyên tố khác ~ > p_i~ thì rõ ràng khiến ~A_i > 10^{12}~). Nén các ~p_i~ tìm được theo thứ tự từ ~1~ trở đi để thuận tiện sử dụng mảng đếm ~cnt~ thông thường.

Chú ý : sử dụng phép chia dư cho ~10^9 + 7~ ở bước tính kết quả một cách cẩn thận.

Độ phức tạp: ~O((N + Q) \sqrt{N})~.

Code mẫu

#include <bits/stdc++.h>

using namespace std;

#define int long long 
#define fo(i, a, b) for(int i = a; i <= b; i++)
#define _fo(i, a, b) for(int i = a; i >= b; i--)
#define foa(i, a) for (auto &i : a)
#define sz(a) ((int) a.size())
#define all(a) begin(a), end(a)
#define fi first
#define se second
#define pb(x) push_back(x)
#define mk(x, y) make_pair(x, y)

typedef long long ll;
typedef pair<ll, ll> pll;
typedef vector<ll> vl;
typedef pair<int, int> pii;
typedef vector<int> vi;

const int LOG = 22; 
const int MAX = 1e5+5;
const int MOD = 1e9+7;
const int SQRT = 316;
const ll INF = 1e12;
const ll lon = 1e18;

ll n, q;
ll a[100005];

vector<ll> smp;
ll sieve[10000005];
map<ll, ll> compress;

ll ways = 1;
ll cnt[100005], inv[100005];
vector<pair<pll, ll>> offline[505];
ll ans[100005];

ll add(ll a, ll b) { return (a+b) % MOD; }
ll mul(ll a, ll b) { return (a*b) % MOD; }
ll sub(ll a, ll b) { return (a+MOD-b) % MOD; }

ll pw(ll val, ll temp) {
    ll curr = 1;

    while(temp > 0) {
        if(temp & 1) curr = mul(curr, val);
        val = mul(val, val);
        temp >>= 1;
    }

    return curr;
}

void precompute() {
    ll low = a[1];
    fo(i, 1, n) low = min(low, a[i]);
    low--;

    ll lim = 1e6;
    for(ll p = 2; p <= lim; p++) {
        if(sieve[p] == 0) {
            smp.pb(p);
            for(ll i = p*p; i <= lim; i += p) sieve[i] = 1;
        }
    }
    fill(sieve+1, sieve+lim+1, 0);

    lim = 1e7;
    foa(p, smp) {
        ll temp = p*(low/p+1);
        for(ll i = temp; i <= low+lim; i += p) {
            if(sieve[i-low] == 0) sieve[i-low] = p; 
        }
    }

    fo(i, 1, n) {
        ll temp = sieve[a[i]-low];
        if(temp == 0) temp = a[i];
        a[i] = temp;
        compress[temp] = 0;
    }
    ll cnt = 0;
    foa(elem, compress) {
        cnt++;
        elem.se = cnt;
    }
    fo(i, 1, n) a[i] = compress[a[i]];  

    fo(i, 1, 100001) inv[i] = pw(i, MOD-2);
}

void update(ll p, ll delta) {
    ways = mul(ways, inv[cnt[p]+1]);
    cnt[p] += delta;
    ways = mul(ways, cnt[p]+1);
}

ll get(ll l, ll r) {
    ll res = 1;
    fo(i, l, r) update(a[i], 1);
    res = sub(ways, 1);
    fo(i, l, r) update(a[i], -1);
    return res;
}

void solve() {
    fo(i, 1, 500) sort(all(offline[i]));
    for(ll i = 1; i <= n; i += SQRT) {
        ll block = (i-1) / SQRT, ptr = i+1;
        while(!offline[block].empty()) {
            if(offline[block].back().fi.fi != ptr) {
                ptr--;
                update(a[ptr], 1);
                continue;
            }
            ll l = offline[block].back().fi.fi, r = offline[block].back().fi.se, id = offline[block].back().se;
            fo(j, i+1, r) update(a[j], 1);
            ans[id] = sub(ways, 1);
            fo(j, i+1, r) update(a[j], -1);
            offline[block].pop_back();
        }
        fo(j, ptr, i) update(a[j], -1);
    }
}

bool check() {
    ll mx = a[1];
    fo(i, 1, n) mx = max(mx, a[i]);
    ll mn = a[1];
    fo(i, 1, n) mn = min(mn, a[i]);
    ll lim = 1e7;
    if(mx-mn >= lim) return false;
    return true;
}

signed main() {
    ios_base::sync_with_stdio(0); cin.tie(0); cout.tie(0);

    cin >> n >> q;
    fo(i, 1, n) cin >> a[i];
    if(check()) {
    precompute();

    fo(i, 1, q) {
        ll l, r;
        cin >> l >> r;
        if(r-l+1 > SQRT) {
            ll block = (r-1) / SQRT;
            offline[block].push_back(mk(mk(l, r), i));
        }
        else ans[i] = get(l, r);
    }
    solve();

    fo(i, 1, q) cout << ans[i] << "\n";     
    }
    else cout << "cac\n";
}