返回顶部

模板 - min25筛

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

const int MAXN = 100000 + 5; //MAXN = sqrt(MAXINPUTN)
const ll MOD = 1000000007, inv2 = 500000004, inv3 = 333333336;

int num;
const int MAXP = 10000; //MAXN以内质数的数量,至多与MAXN一样大
ll prime[MAXP], sp1[MAXP], sp2[MAXP];   //在调用init1()之后,可以精确计算出这里需要的大小

int tot;
const int DMAXN = (MAXN<<1) + 5; //整除分块需要的空间
ll n, SQRT, g1[DMAXN], g2[DMAXN], w[DMAXN], ind1[DMAXN], ind2[DMAXN];

void init1() {  //线性筛预处理
    SQRT = sqrt(n), num = 0 ;
    bitset < MAXN+5 > notprime;
    notprime[1] = 1;
    for(int i = 2; i <= SQRT; ++i) {
        if(!notprime[i]) {
            prime[++num] = i;
            sp1[num] = (sp1[num - 1] + i) % MOD;
            sp2[num] = (sp2[num - 1] + 1ll * i * i) % MOD;
            //spx[num] = 前num个质数的 x 次方和
        }
        for(int j = 1; j <= num && prime[j]*i <= SQRT; ++j) {
            notprime[i * prime[j]] = 1;
            if(i % prime[j] == 0)
                break;
        }
    }
    printf("prime num=%d\n", num);    //传入最大的n,算出num之后给MAXP赋值
}

void init2() {
    tot = 0;
    for(ll l = 1, r, t; l <= n; l = r + 1) {
        t = n / l, r = n / t;
        w[++tot] = t;
        g1[tot] = w[tot] % MOD;
        g2[tot] = (((g1[tot] * (g1[tot] + 1)) >> 1) % MOD * (2 * g1[tot] + 1) % MOD * inv3 % MOD) - 1;
        g1[tot] = ((g1[tot] * (g1[tot] + 1)) >> 1) % MOD - 1;
        if(t <= SQRT)
            ind1[t] = tot;
        else
            ind2[r] = tot;
    }
    //gx(n,j) 表示 [1,n]中,i是质数或者i的最小质因子>pj的数的 x 次方和,n滚动省略
    //ind1和ind2用来记录这个数在数组中的位置
    for(int i = 1; i <= num; i++) {
        for(int j = 1; j <= tot && prime[i]*prime[i] <= w[j]; j++) {
            ll k = w[j] / prime[i] <= SQRT ? ind1[w[j] / prime[i]] : ind2[n / (w[j] / prime[i])];
            g1[j] -= prime[i] * (g1[k] - sp1[i - 1] + MOD) % MOD;
            g2[j] -= prime[i] * prime[i] % MOD * (g2[k] - sp2[i - 1] + MOD) % MOD;
            g1[j] %= MOD, g2[j] %= MOD;
            if(g1[j] < 0)
                g1[j] += MOD;
            if(g2[j] < 0)
                g2[j] += MOD;
        }
    }
}

ll S(ll x, int y) {  //第二部分,不要记忆化,记忆化又慢又卡
    if(prime[y] >= x)
        return 0;
    int k = (x <= SQRT) ? ind1[x] : ind2[n / x];
    ll ans = (g2[k] - g1[k] + MOD - (sp2[y] - sp1[y]) + MOD) % MOD;
    for(int i = y + 1; i <= num && prime[i]*prime[i] <= x; i++) {
        ll pe = prime[i];
        for(int e = 1; pe <= x; ++e, pe = pe * prime[i]) {
            ll xx = pe % MOD;
            ans = (ans + xx * (xx - 1) % MOD * (S(x / pe, i) + (e != 1))) % MOD;
        }
    }
    return ans;
}

int main() {
#ifdef Yinku
    freopen("Yinku.in", "r", stdin);
#endif // Yinku
    scanf("%lld", &n);
    init1();
    init2();
    printf("%lld\n", (S(n, 0) + 1) % MOD);  //加上f(1)
    return 0;
}
posted @ 2019-10-27 16:14  Inko  阅读(...)  评论(...编辑  收藏