【min_25筛】筛欧拉函数、莫比乌斯函数、积性函数

传送门:欧拉函数、莫比乌斯函数

花了一下午学了一下min25筛。

一知半解。不过还是学会了套板子。

mu函数要看成质数个数取反,在预处理g的时候不能带负数(我也不知道为什么)

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 1e5 + 233;
ll primes[maxn], cnt;
bool st[maxn];
ll sp1[maxn], sp2[maxn], g1[maxn], g2[maxn], sm[maxn], g[maxn], w[maxn], wcnt;
ll id1[maxn], id2[maxn];
ll inv2, inv6, qn;
ll ksm(ll n, ll k)
{
    ll res = 1;
    while(k)
    {
        if(k & 1) res = res * n;
        n = n * n ;
        k >>= 1;
    }
    return res;
}
void get_primes(int n)
{
    for(int i = 2; i <= n; i++) {
        if(!st[i]) {
            primes[++cnt] = i;
        }
        for(int j = 1; j <= cnt && i * primes[j] <= n; j++) {
            st[i * primes[j]] = 1;
            if(i % primes[j] == 0) 
                break;
        }
    }
    for(int i = 1; i <= cnt; i++) {
        sp2[i] = primes[i] + sp2[i - 1], sp1[i] = sp1[i - 1] + 1;
        sm[i] = i;
    }
}
ll n;
ll dfs(int x, int y)
{
    if(primes[y] >= x) return 0;
    int k = x;
    if(k <= qn) k = id1[x]; else k = id2[n / x];
    ll res = (g2[k] - g1[k]) - (sp2[y] - sp1[y]);
    for(int i = y + 1; i <= cnt && primes[i] * primes[i] <= x; i++)
    {
        ll p = primes[i];
        for(int e = 1; p <= x; e++, p *= primes[i])
        {
            res = (res + (p / primes[i] * (primes[i] - 1)) * (dfs(x / p, i) + (e != 1)));
        }
    }
    return res;
}
ll dfs2(int x, int y)
{
    if(primes[y] >= x) return 0;
    int k = x;
    if(k <= qn) k = id1[x]; else k = id2[n / x];
    ll res = (g[k]) - (-sm[y]);
    for(int i = y + 1; i <= cnt && primes[i] * primes[i] <= x; i++)
    {
        ll p = primes[i];
        res = (res - (dfs2(x / p, i)));
    }
    return res;
}
int main()
{
    get_primes(1e5 + 10);
    int T; cin >> T;
    while(T--)
    {
        cin >> n;
        qn = sqrt(n);
        // memset(g1, 0, sizeof g1);
        // memset(g2, 0, sizeof g2);
        // memset(g, 0, sizeof g);
        // memset(w, 0, sizeof w);
        wcnt = 0;
        for(ll l = 1, r; l <= n; l = r + 1)
        {
            r = n / (n / l);
            w[++wcnt] = n / r;
            ll t = (n / r);
            g1[wcnt] = t - 1; //g(1)排除
            g2[wcnt] = (t + 1) * t / 2 - 1;
            g[wcnt] = t - 1;
            if(n / r <= qn) id1[n / r] = wcnt; else id2[r] = wcnt;
        }

        for(int i = 1; i <= cnt; i++)
        {
            ll pp = primes[i] * primes[i];
            for(int j = 1; j <= wcnt &&  pp <= w[j]; j++)
            {
                ll k = w[j] / primes[i];
                if(k <= qn) k = id1[k]; else k = id2[n / k];    
                g2[j] = g2[j] - primes[i] * (g2[k] - sp2[i - 1]);
                g1[j] = g1[j] - (g1[k] - sp1[i - 1]);
                g[j] = g[j] - (g[k] - sm[i - 1]);
            }
        }
        for(int i = 1; i <= wcnt; i++) g[i] = -g[i];
        // for(int i = 1; i <= wcnt; i++) cerr << g[i] << endl;
        printf("%lld %lld\n", dfs(n, 0) + 1, dfs2(n, 0) + 1);
    }

}

传送门:min25模板题

只会套板子。。

// luogu-judger-enable-o2
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 2e5 + 233;
ll mod = 1e9 + 7;
ll primes[maxn], phi[maxn], cnt;
ll w[maxn], wcnt; // 下标
bool st[maxn];//int st[N];

ll s_p1[maxn], s_p2[maxn], g1[maxn], g2[maxn], id1[maxn], id2[maxn];
// f(p) = p*p - p
ll inv2, inv6, qn;
ll ksm(ll n, ll k)
{
    ll res = 1;
    while(k)
    {
        if(k & 1) res = res * n % mod;
        n = n * n % mod;
        k >>= 1;
    }
    return res;
}
void get_primes(int n)
{
    // mu[1] = phi[1] = 1;
    phi[1] = 1;
    for(int i = 2; i <= n; i++) {
        if(!st[i]) {
            primes[++cnt] = i;
            phi[i] = i - 1;
            // mu[i] = -1;
        }
        for(int j = 1; j <= cnt && i * primes[j] <= n; j++) {
            st[i * primes[j]] = 1;
            if(i % primes[j] == 0) {
                phi[i * primes[j]] = phi[i] * primes[j];
                // mu[i * primes[j]] = 0;
                break;
            }else {
                phi[i * primes[j]] = phi[i] * (primes[j] - 1);
                // mu[i * primes[j]] = -mu[i];
            }
        }
    }
    for(int i = 1; i <= cnt; i++) 
        s_p1[i] = (primes[i] + s_p1[i - 1]) % mod, s_p2[i] = (s_p2[i - 1] + primes[i] * primes[i]) % mod;
}
ll n; 
ll dfs(ll x, int y)
{
    if(primes[y] >= x) return 0;
    ll k = x;
    if(k <= qn) k = id1[x]; else k = id2[n / x];
    ll res = ((g2[k] - g1[k]) - (s_p2[y] - s_p1[y])) % mod;
    if(res < 0) res += mod;
    for(int i = y + 1; i <= cnt && primes[i] * primes[i] <= x; i++)
    {
        ll p = primes[i];
        for(int e = 1; p <= x; e++, p *= primes[i])
        {
            ll now = p % mod;
            res = (res + now * (now - 1LL) % mod * (dfs(x / p, i) % mod + (e != 1)) % mod ) % mod;
        }
    }
    return res;
}
int main()
{
    inv2 = ksm(2, mod - 2); inv6 = ksm(6, mod - 2);
    cin >> n;
    qn = sqrt(n);
    get_primes(1e5 + 10);
    for(ll l = 1, r; l <= n; l = r + 1)
    {
        r = n / (n / l);
        w[++wcnt] = n / r;
        ll t = (n / r) % mod;
        g1[wcnt] = (t + 1) * t % mod * inv2 % mod - 1; //g(1)排除
        g2[wcnt] = (2 * t + 1) * (t + 1) % mod * t % mod * inv6 % mod - 1;
        if(n / r <= qn) id1[n / r] = wcnt; else id2[r] = wcnt;
    }
    for(int i = 1; i <= cnt; i++)
    {
        for(int j = 1; j <= wcnt && primes[i] * primes[i] <= w[j]; j++)
        {
            ll k = w[j] / primes[i];
            if(k <= qn) k = id1[k]; else k = id2[n / k];
            g1[j] = (g1[j] - primes[i] * (g1[k] - s_p1[i - 1]) % mod) % mod;
            g2[j] = (g2[j] - primes[i] * primes[i] % mod * (g2[k] - s_p2[i - 1]) % mod) % mod;
            if(g1[j] < 0) g1[j] += mod;
            if(g2[j] < 0) g2[j] += mod;
        }
    }
    printf("%lld\n", dfs(n, 0) + 1LL);
}

发表评论

邮箱地址不会被公开。 必填项已用*标注