Description：

定義 f(n) 表示 1 到 n 這 n 個數，同 n 的最大公約數的和。
例如 f(1)=1,f(2)=3,f(3)=5,f(6)=15 ，比如 n=6 時 1,2,3,4,5,6 同 6 的最大公約數分別為 1,2,3,2,1,6 ，它們的和是 15 。
小Q遇到了一個問題，Nod老師要求他回家計算 g(n)=∑d|nf(d) 這個函式模 1000000007 的值，例如 g(6)=f(1)+f(2)+f(3)+f(6)=24 。
小Q想了一下，如果有 n≤109 的條件，那麼他是會做的，但是現在這個 n 非常的大，他覺得有點不太會算了，你能幫他計算一下 g(n) 模 1000000007 的值嗎？
由於 n 非常大，所以 n 被表示成了 m 個正整數 x1,x2,⋯,xm 的乘積，即 n=∏mi=1xi 。
由於 m 也非常大，所以對於 x

基準時間限制：1.3 秒 空間限制：131072 KB

題解：

先隨便反演一下。
g(n)=∑d|nf(d)
=∑d|n∑d′|dd′∗φ(d/d′)
=∑d′|nd′∗∑d|nd′d∗φ(ndd′)
=∑d′|nd′∗(n/d′)
=∑d′|nn
=n∗σ0(n)

m雖然很大，但是c很小，所以可以找x的迴圈節。

現在需要快速將10000000個數分解質因數，可以線篩每個數的最小質因子和指數及對應的次冪，然後暴力統計，注意這題時限卡的很緊，一開始我這麼做，然後T飛了。

於是去百度上膜了下標，發現只需要記錄最小質因子，然後倒著掃一遍即可。

注意空間卡得也緊，所以需要重複利用陣列。

Code：

#include<cstdio>
#include<cstring>
#define ll long long
#define fo(i, x, y) for(int i = x; i <= y; i ++)
#define fd(i, x, y) for(int i = x; i >= y; i --)
using namespace std;

const int N = 1e7, mo = 1e9 + 7;

int
 
 p[700000], u[N + 5], x[N + 5], bz[N + 5];
int a, b, c;
ll m, ans;

void Build() {
    fo(i, 2, N) {
        if(!bz[i]) p[++ p[0]] = i, u[i] = i;
        fo(j, 1, p[0]) {
            int k = i * p[j];
            if(k > N) break;
            bz[k] = 1; u[k] = p[j];
            if(i % p[j] == 0) break;
        }
    }
}

ll ksm(ll x, ll y) {
    ll s = 1;
    for(; y; y >>= 1, x = x * x % mo)
        if(y & 1)
            s = s * x % mo;
    return s;
}

void add(int &x, int y) {
    x += y, x = x >= mo ? x - mo : x;
}

int main() {
    scanf("%lld %d %d %d %d", &m, &x[1], &a, &b, &c);
    Build(); memset(bz, 0, sizeof(bz));
    int k = 1; bz[x[1]] = 1;
    while(1) {
        x[k + 1] = ((ll)x[k] * a + b) % c + 1; k ++;
        if(bz[x[k]]) break;
        bz[x[k]] = k;
    }
    int last = bz[x[k]] - 1; k --;
    int round = k - last;
    memset(bz, 0, sizeof(bz));
    ans = 1; int tem = (m - last) / round % mo;
    fo(i, last + 1, k) {
        ans = ans * x[i] % mo;
        add(bz[x[i]], tem);
    }
    ans = ksm(ans, (m - last) / round % (mo - 1));
    fo(i, 1, last)
        ans = ans * x[i] % mo, add(bz[x[i]], 1);
    tem = (m - last) % round;
    fo(i, last + 1, last + tem)
        ans = ans * x[i] % mo, add(bz[x[i]], 1);
    fd(i, N, 2)
        if(bz[i] && u[i] != i)
            add(bz[u[i]], bz[i]), add(bz[i / u[i]], bz[i]);
    fo(i, 1, p[0]) ans = ans * (bz[p[i]] + 1) % mo;
    printf("%lld", ans);
}